# Coarse grained intermolecular interactions on quantum processors

Lewis W. Anderson<EMAIL_ADDRESS>Clarendon Laboratory,
University of Oxford, Parks Road, Oxford OX1 3PU, UK  Martin Kiffner  Centre
for Quantum Technologies, National University of Singapore, 3 Science Drive 2,
Singapore 117543  Clarendon Laboratory, University of Oxford, Parks Road,
Oxford OX1 3PU, UK  Panagiotis Kl. Barkoutsos  IBM Quantum, IBM Research
Zurich, Säumerstrasse 4, 8803 Rüschlikon, Switzerland  Ivano Tavernelli  IBM
Quantum, IBM Research Zurich, Säumerstrasse 4, 8803 Rüschlikon, Switzerland
Jason Crain  IBM Research Europe, The Hartree Centre STFC Laboratory, Sci-Tech
Daresbury, Warrington WA4 4AD, UK  Department of Biochemistry, University of
Oxford, Oxford OX1 3QU, UK  Dieter Jaksch  Institut für Laserphysik,
Universität Hamburg, 22761 Hamburg, Germany  Clarendon Laboratory, University
of Oxford, Parks Road, Oxford OX1 3PU, UK  Centre for Quantum Technologies,
National University of Singapore, 3 Science Drive 2, Singapore 117543

###### Abstract

Variational quantum algorithms (VQAs) are increasingly being applied in
simulations of strongly-bound (covalently bonded) systems using full molecular
orbital basis representations. The application of quantum computers to the
weakly-bound intermolecular and non-covalently bonded regime however has
remained largely unexplored. In this work, we develop a coarse-grained
representation of the electronic response that is ideally suited for
determining the ground state of weakly interacting molecules using a VQA. We
require qubit numbers that grow linearly with the number of molecules and
derive scaling behaviour for the number of circuits and measurements required,
which compare favourably to traditional variational quantum eigensolver
methods. We demonstrate our method on IBM superconducting quantum processors
and show its capability to resolve the dispersion energy as a function of
separation for a pair of non-polar molecules - thereby establishing a means by
which quantum computers can model Van der Waals interactions directly from
zero-point quantum fluctuations. Within this coarse-grained approximation, we
conclude that current-generation quantum hardware is capable of probing
energies in this weakly bound but nevertheless chemically ubiquitous and
biologically important regime. Finally, we perform simulations on an extended
version of the model which includes multipole interactions and anharmonicity -
the consequences of the latter are unexamined in large systems using classical
computational methods but can be incorporated here with low computational
overhead.

##  I Introduction

Computer simulation of matter at the atomic and molecular scale has advanced
to the point where it is now playing an increasingly significant role in the
design and discovery of new functional materials. However, these simulations
all face the common challenge of how to strike inevitable compromises between
the physical realism of an underlying model and the practicality of its
implementation across the deep hierarchy of intermolecular forces responsible
for cohesion in condensed phases.

This hierarchy spans a continuum of energy scales ranging from greater than 1
eV for intra-molecular bonds, extending down to $10^{-1}-10^{-2}$ eV for the
inter-molecular regime Saito (2020). The weakest of these are the dispersion
interactions arising from instantaneous quantum fluctuations orf the electron
distribution. These give rise to correlated induced multipole moments with the
leading-order contribution coming from dipolar coupling – attenuated by the
familiar $1/R^{6}$ asymptotic decay with distance. Dispersion interactions
therefore form part of the non-local correlation energy of the system and are
not pair-wise additive. Although generally weaker than electrostatic
interactions or hydrogen bonds, their long-ranged character and ubiquity means
that they exert influence to some extent in all non-covalent associations.

There is now clear evidence that dispersion effects extend beyond cohesive
energy and have impact on structural and mechanical properties Reilly and
Tkatchenko (2015); Shtogun and Woods (2010) as well as barrier heights and
phase transitions Pastorczak  _et al._ (2017). They can therefore play
governing roles in a wide range of condensed matter phenomena including liquid
state physics, molecular crystal stability, low-dimensional and layered system
stacking, hybrid organic–inorganic interface properties and biophysical
interactions such as drug-ligand binding. Moreover, as these interactions
scale with system size, their importance is expected to grow for larger
molecules and supramolecular assemblies Hermann  _et al._ (2017).

Principally, there are two distinct strategies employed to capture this
molecular interaction hierarchy in computer models: Firstly, there are
empirical potential or force-field methods. Here the system energy is
typically decomposed into bonded (stretches, bends and torsions) and non-
bonded (electrostatic and dispersive) terms. These are parameterised to fit
certain experimental reference quantities. The motion of the system can then
be evolved by classical mechanics (molecular dynamics), and trajectories
averaged to obtain static properties. These methods have the advantage of
extreme computational efficiency meaning that they can be applied to
relatively large systems. However they generally lack explicit electronic
responses and, consequently, transferability to thermodynamic or environmental
conditions outside those of the parameterisation range is limited.

Secondly, there are so-called first principles approaches in which the
Schrodinger equation is solved to various levels of approximation. Here the
prevailing methods for molecular dynamics applications are based on Kohn-Sham
density functional theory. In these models, there are no empirical force laws:
Instead, forces are computed “on-the-fly” directly from the electronic
structure according to the Hellmann-Feynman Hellman (1937); Feynman (1939)
theorem with trajectories evolved on the Born-Oppenheimer surface Born and
Oppenheimer (1927). These techniques are computationally resource-intensive
and tend to be limited to relatively small systems. Also, common density
functional approximations do not properly capture non-local correlations
involved in dispersion forces and there are active efforts to devise
dispersion-corrected density functionals French  _et al._ (2010); von
Lilienfeld  _et al._ (2004); Tkatchenko and Scheffler (2009); Grimme (2011).
High-level quantum chemistry methods (such as configuration interaction (CI)
Slater (1929); Boys (1950); Szalay  _et al._ (2012) or canonical coupled
cluster (CC) Čížek (1966); Bartlett and Musiał (2007) with perturbative
(single, double, triple) excitations) can account for dispersion interactions
explicitly within various levels of approximation albeit with even more severe
limitations imposed by system size scaling Pieniazek  _et al._ (2007);
Oliveira and Kraka (2017); Bistoni (2020).

Quantum computing offers uniquely promising routes for efficient simulations
of correlated electronic systems at the first-principles level. Materials
modelling at atomic and molecular scales is therefore one of its most
anticipated application areas. One such route is to use variational quantum
algorithms (VQAs) Preskill (2018); Cerezo  _et al._ (2021) for solving the
Schrödinger equation for chemical systems McArdle  _et al._ (2020). VQAs have
the benefits of requiring relatively short coherence times and modest numbers
of qubits which is achieved by using a shallow parameterised quantum circuit
embedded within a classical optimisation routine to find the ground state of a
given problem Hamiltonian Peruzzo  _et al._ (2014). Since the first proposal
of using variational methods - namely the variational quantum eigensolver
(VQE) for chemical systems Peruzzo  _et al._ (2014), VQAs have seen
application to a significant number of areas including combinatorial
optimisation Farhi  _et al._ (2014); Farhi and Harrow (2016), strongly
correlated materials Bauer  _et al._ (2016); Kreula  _et al._ (2016a, b) and
non-linear partial differential equations Lubasch  _et al._ (2020); Lloyd  _et
al._ (2020). Current capabilities to model molecular properties from first
principles in large systems using VQAs are however limited by the prohibitive
computational cost of the high-level, wavefunction methods typically deployed
in quantum chemistry. For example, solution of the Schrödinger equation at the
full configuration interaction (so-called “full CI”) level for the ground-
state wave function for molecular systems exhibits factorial scaling
Ollitrault  _et al._ (2020a).

Chemical applications of VQAs have consequently been confined to the strongly-
bound intramolecular regime of the interaction hierarchy and to small systems
with examples now demonstrated for dissociation curves of small molecules
Grimsley  _et al._ (2019); Kandala  _et al._ (2017), ground state energies for
intramolecular bonds Nam  _et al._ (2020), molecular excitation energies
Ollitrault  _et al._ (2020a), covalent binding energies of hydrogen chains,
isomerisation barriers Arute  _et al._ (2020), molecular vibrations McArdle
_et al._ (2019); Ollitrault  _et al._ (2020b) and, two-site DMFT calculations
Rungger  _et al._ (2019); Jaderberg  _et al._ (2020). All except the last two
examples consider Hamiltonians and states within full electronic basis sets
such as Slater-type or Gaussian-type orbitals Helgaker  _et al._ (2014). For
the sake of comparison with our method, we will refer to these as orbital-
based VQE methods.

$-\delta$
---
$+\delta$
---
$-\delta$
---
$+\delta$
---
$\gamma$
---
IBM QPU
---
(VQE: $\bm{\theta}\rightarrow\bm{\theta^{\prime}}$)
---
$\langle\hat{H}(\bm{\theta})\rangle$
---
CPU
---
$\bm{\theta}$
---
$\bm{\theta}$
---
$\bm{\theta}$
---
Figure 1: Schematic of the VQE optimisation procedure and output states
showing instantaneous assymetry of the molecular charge distribution. These
correlated charge density fluctuations lead to dispersion interaction which,
to the lowest order, is the familiar $1/R^{6}$ London attraction.

In contrast, the intermolecular (non-covalent) interaction regime has been
largely unexplored using quantum algorithms and is a gap which will need to be
filled in order to address realistic material science problems using near-
term, noisy, quantum computers. Due to the size of the molecular complexes of
interest and the necessary inclusion of core electrons for dispersive forces,
the required number of electrons and orbitals implies a large number of basis
functions and qubits Sherrill (1996, 2017). This makes the direct scale up of
currently used quantum algorithms for electronic structure calculations very
impractical and therefore strategies are needed to reduce the number of
electrons and orbitals treated explicitly at the highest levels of accuracy.

To address this difficulty, we develop a coarse grained model of electronic
interactions that is ideally suited for VQA and can easily be extended to
represent a wide range of possible interactions. Our approach is inspired by
and extends a maximally coarse grained model used in classical computations
Cipcigan  _et al._ (2019). In this model the responsive (polarisable) part of
the molecular charge distribution is modelled as a quantum harmonic oscillator
embedded in the molecular frame and its properties can be tuned to reproduce
reference values of polarisability or dispersion coefficients for real
molecular species. Zero-point multipolar fluctuations are present by
construction and so dispersion arises naturally. While the model parameters
are defined empirically to fit molecular properties there are no interaction
potentials or force laws between atoms or molecules defined a priori. These
are instead determined from the coarse grained electronic structure as in ab-
initio approaches.

This model blends therefore features from both empirical potential and first
principles methodologies and has already been solved, with model limitations,
using a path integral approach on classical processors. In particular, it has
been successfully applied to noble gas liquids and solids (pure dispersion)
and to water (where both polarisation and dispersion are present) across its
phase diagram where it has correctly predicted vaporisation enthalpy,
temperature of maximum density, temperature dependences of the dielectric
constant and surface tension Cipcigan  _et al._ (2019, 2018); Sokhan  _et al._
(2015a) and revealed molecular signatures of Widom line crossing in the
supercritical fluid Sokhan  _et al._ (2015b). These classical methods have
however been developed with the specific aim of deriving force laws and are
restricted to harmonic confining potentials where certain order interactions
are missing and others are related by simple algebraic relationships that only
hold approximately for real systems Cipcigan  _et al._ (2016).

By using quantum processors to model a coarse grained representation, our
method gives a more general framework than classical approaches. This allows
for the extension to anharmonic oscillators, which would capture these missing
interactions as well as allowing for the introduction of more complex physics
through non-linear interaction terms. Fig. 1 shows a schematic representation
of the algorithm and correlated charge fluctuations captured by the model.

###  I.1 Summary of results

Here we summarize our main findings: Firstly, the number of qubits required to
simulate the model on a quantum computer grows linearly with the number of
molecules. Secondly, the number of unique circuits required to measure our
Hamiltonian with all-to-all dipolar interactions scales at most linearly with
the number of oscillators. This compares favourably to the $O(N^{3})$ unique
circuits that must be measured for a full orbital-based VQE Hamiltonian with
$N$ orbitals Gokhale  _et al._ (2020). Furthermore, we present simulated
results showing how this model can be extended to include anharmonic electron-
nuclear potentials. For our algorithm, relaxation of the harmonic
approximation incurs negligible experimental overhead. Classical methods, such
as path integral and Monte Carlo techniques, typically rely on the efficient
sampling Westbroek  _et al._ (2018) of two-point correlators for Gaussian
states generated by harmonic Hamiltonians. However, the non-Gaussian ground
states of anharmonic Hamiltonians are much more costly to sample. Thus,
anharmonic potentials have not been studied much using classical computational
approaches.

Figure 2: London dispersion energy ($\Delta E$) calculated using variational
quantum algorithm for one-dimensional quantum Drude oscillators computing
interacting I2 dimers. The black line gives the dispersion energy obtained
analytically. Red squares correspond to noise-free evaluation (via state
vector simulations) of the experimentally optimized coupled ground state. Blue
circles correspond to energy directly obtained from experimentally optimized
ground states (error bars represent one-standard deviation on the mean). The
horizontal line marks zero dispersion interaction corresponding to uncoupled
oscillators at $R=\infty$ and vertical dashed line corresponds to the I2 van
der Waals diameter.

We also present proof of principle experiments using IBM superconducting
quantum processors demonstrating the solubility of our model on current
generation hardware. We focus on two linear non-polar molecules where the
electronic responses are described by embedded quantum oscillators. Here the
interactions arise from London dispersion forces for which the leading
contribution comes from dipolar quantum fluctuations. Looking ahead to the
main experimental results of this paper, Fig. 2 shows the result of our VQA;
combined with a simple method to extract the binding energy; for a pair of
these interacting oscillators with parameters chosen to represent an Iodine
(I2) molecular dimer. This demonstration successfully recreates the leading
order London dispersion interaction, thereby giving the first demonstration of
the use of quantum processors to model dispersion forces directly arising from
quantum fluctuations.

Our coarse grained model allows for the use of quantum computers to
efficiently model (with favourable scaling) the inter-molecular interactions
between large molecules, which are important in characterising a wide range of
chemical and biological systems. This fills an essential gap in the chemical
interaction hierarchy and allows us to extend current coarse grained classical
methods to capture a richer set of physical behaviours with the necessary
accuracy. Given the reduced resource requirements (i.e., number of qubits,
gate operations and circuit measurements) when compared to a full orbital-
based electronic structure calculations, the methodology presented here may be
more readily scalable to demonstrating an application of quantum advantage.

The layout of this article is as follows. In Sec. II we describe the coarse
grained model and how it can be solved using a variational quantum algorithm.
This includes: Sec. II.1 describing the coarse grained model of interacting
oscillators; Sec. II.2 giving a description of how the this system is
represented on a quantum computer; Sec. II.3 which makes use of results from
graph theory and measure theory to give a scheme for grouping operators needed
to measure the Hamiltonian as well as derive upper and lower bounds for the
number of shots required to measure the Hamiltonian to a given accuracy; and
Sec. II.4 which describes the variational ansatz used in our experiments. In
Sec. III, we present proof of principle results modelling the London
dispersion interaction of an I2 dimer and a novel technique to subtract the
effect of device noise to accurately measure the binding energy. In Sec. IV we
describe how the one-dimensional model can be extended to three dimensions and
to include anharmonic and non-linear terms to model richer physics than those
accessible to efficient classical methods. In Sec. V we discuss our results
and describe further research directions.

##  II Coarse grained modelling on digital quantum computers

###  II.1 Coarse grained model

In the coarse grained model collective electronic responses are described by a
charged quantum Drude oscillators (QDOs) tethered to an atomic nucleus or
embedded in a molecular frame. For neutral species the Drude (quasi-)particle
and the nucleus to which it is bound have charges $-q$ and $+q$, respectively.
The correlated displacements of these QDOs, which correspond to instantaneous
charge density fluctuations in molecules, allows us to capture dispersive
interactions between non-polar molecules. It is an electronically coarse-
grained model which nonetheless exhibits the full hierarchy of electronic
responses at long range: Many-body polarisation and dispersion to all orders
arise naturally and, in principle, intermolecular forces can then be computed
“on the fly” from the coarse-grained electron distribution of coulomb-coupled
oscillators without the need for a priori force laws or empirical interaction
potentials. One further feature of this model, is that it obeys
distinguishable particle statistics – Fermi statistics are not required. The
model must therefore be augmented by a simple short ranged repulsion term. The
basic method has been demonstrated on classical processors in the case of
noble gas solids and single-component liquids Jones  _et al._ (2013a, b).
Using only the properties of the isolated molecule in the model parameters,
emergent properties of these systems have been predicted across a range of
their respective phase diagrams. For a recent review of these methods and
results see Cipcigan  _et al._ (2019).

For simplicity, for most of this work we will focus on the case of one-
dimensional QDOs with dipole order coupling, however much of the methodology
developed here can be applied straightforwardly to systems of three-
dimensional oscillators and we give a brief description of this extension in
Sec. IV.1. Even the simple one-dimensional model recreates the London
dispersion interaction proportional to $R^{-6}$ and $\alpha^{2}$, where
$\alpha$ is the dipole polarisability and $R$ the physical separation, up-to a
constant factor for spherical non-polar molecules and exactly for a pair of
linear non-polar molecules in the limit of infinitely anisotropic
polarisability (See Appendix A for further details.)

The non-dimensional Hamiltonian for $N$ one-dimensional QDOs with quadratic
coupling is

$\hat{H}=\sum_{i=1}^{N}(\hat{x}_{i}^{2}+\hat{p}_{i}^{2})+\sum_{j>i}\gamma_{i,j}\hat{x}_{i}\hat{x}_{j},$
(1)

where $\hat{x}_{i}$ and $\hat{p}_{i}$ are the bosonic position and momentum
operators for oscillator $i$ and $\gamma_{i,j}\in\mathbb{R}$ is the coupling
between oscillators $i$ and $j$. Throughout this work we will implicitly work
in energy units of $\hbar\omega/2$, for simplicity of notation assuming that
all oscillators share the same trapping frequency $\omega$. For each pair of
identical one-dimensional oscillators aligned parallel to each other and along
the inter-oscillator axis, separated by distance $R_{i,j}$, the coupling
constant is given by
$\gamma_{i,j}=-2(q^{2}/\mu\omega^{2})R_{i,j}^{-3}=-4\alpha R_{i,j}^{-3}$,
where the effective charge of the Drude particle $q$, the effective mass $\mu$
and the oscillator frequency $\omega$ are model parameters. The dipole
polarisability is given by $\alpha=q^{2}/\mu\omega^{2}$. If instead the
oscillators are aligned parallel to each other and perpendicular to the inter-
oscillator axis, the coupling is given by $\gamma_{i,j}=2\alpha R_{i,j}^{-3}$.

One advantage of this harmonically bound coarse grained model is that
efficient classical methods exist which can be used for benchmarking quantum
implementations. However, the price paid for this harmonic approximation is
that the relative strengths of interactions within the model are governed by
Gaussian statistics meaning that there are simple algebraic relations relating
interactions at different orders. This leads to the absence of terms, such as
the second hyperpolarisability (dipole-dipole-dipole-dipole) Jones  _et al._
(2013b), due to certain symmetries that emerge in the harmonic model.

In what follows we present a variational quantum algorithm for calculating
intermolecular interactions between molecules using a one-dimensional QDO
model and outline how to extend this to a full three dimensional
representation. Such a quantum algorithm is not limited to the inclusion of
harmonic Hamiltonian terms and can be extended to include anharmonic trapping
potentials as well as non-linear interaction terms that may allow for a richer
set of physics beyond Gaussian statistics and the study of more realistic
chemical systems Jones  _et al._ (2013b).

###  II.2 Algorithm description

In order to find the ground state of the Hamiltonians of interest we use the
common procedure for a Variational Quantum Eigensolver (VQE) algorithm Peruzzo
_et al._ (2014). The main goal of VQE is to produce quantum states
$\ket{\psi(\bm{\theta})}$ using a quantum circuit parameterised by real
parameters $\bm{\theta}=\set{\theta_{0},\theta_{1},\dots}$ which can be
varied. Throughout the algorithm, the parameters are updated using a classical
feedback loop in order to minimise the expected values of the Hamiltonian on
the particular quantum state. Guided from the variational principle the
optimisation procedure can provide a quantum state that approximates the exact
ground state of the system by solving the minimisation problem

$\min_{\bm{\theta}}\braket{\psi(\bm{\theta})}{\hat{H}}{\psi(\bm{\theta})},$
(2)

provided the variational ansatz is able to represent a state close to the
ground state within the required precision. In the following we outline how we
encode the states of harmonic QDOs into qubit states and how the Hamiltonian
is measured efficiently on a quantum computer before defining our variational
ansatz later.

We represent the state of the QDOs in the Fock basis which for Bosonic
oscillators form an infinite basis. We are limited to a finite number of
qubits and so must restrict ourselves to the space of a finite subset of Fock
states. This will naturally lead to solutions that approximately match the
true states in full Fock space. Since dispersive interaction energies are
expected to be much smaller than the bare energy scale $\hbar\omega$, the
groundstates of realistic Hamiltonians will be dominated by the lowest energy
Fock states and thus the truncation of the Fock space will not severely impact
the accuracy of solutions. We use a compact encoding Sawaya  _et al._ (2020),
which was proposed in earlier work on vibrational modes on a quantum computer
Dumitrescu  _et al._ (2018); McArdle  _et al._ (2019). For this encoding, the
first $d$ Fock states are mapped onto $m=\log_{2}d$ qubits as

$\mathcal{H}_{d-\text{Fock}}\rightarrow\mathcal{H}_{2}^{\otimes
m}:\;\ket{\underline{n}}\rightarrow\ket{\text{bin}_{m}(n)}$ (3)

where $\text{bin}_{m}(n)$ denotes the $m$-length binary representation of
integer $n\in\set{0,1,\dots,d-1}$. Throughout this work, underlined integers
within a $\ket{\cdot}$ will denote a Fock state and all other states will be
written in the computational basis of the qubits. The basis for multiple QDOs
is then given by the tensor product of the single QDO basis such that for a
total of $N$ $2^{m}$-dimensional QDO we will use $M=mN$ qubits.

Although we are restricted to the lowest $d$ Fock states, we expect that the
ground states required will be sufficiently close to the uncoupled ground
state $\ket{\underline{0}}^{\otimes N}$ such that a finite number of basis
states will give a good approximation to the ground state energy, particularly
in the weakly-coupled regime. Since we require logarithmically many qubits,
this qubit encoding will be efficient provided the number of Fock states
required for a chosen accuracy grows less than double-exponentially.

To measure the cost function for a given variational state we must decompose
each term in the Hamiltonian in (1), restricted to a finite Fock basis, into
tensor product Pauli operators $\hat{P}\in\set{\mathds{1},X,Y,Z}^{\otimes M}$
that can be measured on a quantum computer. For the non-interacting and
coupling terms we show this for a single QDO and a single pair of QDOs
respectively and note that the same decomposition can be used for all terms
within the Hamiltonian just acting on different subsets of $m$ and $2m$
qubits. For the non-interacting contribution of a single QDO, there are $m$
non-identity $Z$ operators that must be measured:

$\hat{x}^{2}+\hat{p}^{2}=2^{m}\mathds{1}^{\otimes
m}-\sum_{i=0}^{m-1}2^{i}Z_{i},$ (4)

where $Z_{i}$ is the Pauli Z operator acting on qubit $i$. For a pair of $m$
qubit QDOs the coupled part of the Hamiltonian leads to
$\left(\frac{1}{2}d\log_{2}d\right)^{2}$ Pauli terms. The form of the coupling
terms for arbitrary $d$, as well as a proof for the number of terms, is given
in Appendix C. As an example, the decompositions for $m=1$ and $m=2$ are given
by

$\hat{x}_{1}\hat{x}_{2}\rvert_{m=1}=X_{1}X_{2},$ (5)
$\begin{split}\hat{x}_{1}\hat{x}_{2}\lvert_{m=2}=&\frac{\sqrt{3}+2}{4}X_{1}X_{3}+\frac{\sqrt{3}+1}{4\sqrt{2}}X_{1}X_{2}X_{3}+\frac{\sqrt{3}+1}{4\sqrt{2}}X_{1}X_{3}X_{4}+\frac{1}{4}X_{1}X_{2}X_{3}X_{4}+\frac{\sqrt{3}+1}{4\sqrt{2}}Y_{1}Y_{2}X_{3}+\frac{1}{4}Y_{1}Y_{2}X_{3}X_{4}\\\
&+\frac{\sqrt{3}+1}{4\sqrt{2}}X_{1}Y_{3}Y_{4}+\frac{1}{4}X_{1}X_{2}Y_{3}Y_{4}-\frac{\sqrt{3}-1}{4\sqrt{2}}Y_{1}Y_{2}X_{3}Z_{4}-\frac{\sqrt{3}-1}{4\sqrt{2}}X_{1}Z_{2}Y_{3}Y_{4}-\frac{1}{4}X_{1}Z_{2}X_{3}\\\
&-\frac{\sqrt{3}-1}{4\sqrt{2}}X_{1}Z_{2}X_{3}X_{4}-\frac{1}{4}X_{1}X_{3}Z_{4}-\frac{\sqrt{3}-1}{4\sqrt{2}}X_{1}X_{2}X_{3}Z_{4}+\frac{2-\sqrt{3}}{4}X_{1}Z_{2}X_{3}Z_{4}+\frac{1}{4}Y_{1}Y_{2}Y_{3}Y_{4}.\end{split}$
(6)

The interaction Hamiltonian for $N$ QDOs can then be constructed by combining
these pair-wise interactions for all $N(N-1)/2$ pairs of QDOs with the correct
coupling constant.

###  II.3 Measurement cost

Measure over states | Ungrouped | Grouped
---|---|---
| Circuits | Shots | Circuits | Shots
Uniform spherical | $O\left(N^{2}(\log_{2}d)^{2}d^{2}\right)$ a |  $O\left(\frac{\gamma^{2}N^{4}d^{8}}{\epsilon^{2}}\right)$ b | $O(Nd^{2})$ c |  $O\left(\frac{\gamma^{2}N^{3}(\log_{2}d)^{4}d^{6}}{\epsilon^{2}}\right)$ d
Pure uncoupled |  $O\left(\frac{\gamma^{2}N^{4}d^{4}}{\epsilon^{2}}\right)$ e |  $O\left(\frac{\gamma^{2}N^{3}d^{4}}{\epsilon^{2}}\right)$ f

For full expressions and proofs see Appendix D: aTheorem 1, bTheorem 7,
cTheorem 4, dTheorem 6, eTheorem 9 and, fTheorem 8.

Table 1: Number of unique circuits required using our exact grouping as well
as expected number of shots required to achieve absolute measurement error
$\epsilon$ for a system of $N$ $d$-level one-dimensional QDOs. Here $\gamma$
is a typical scale for the inter-oscillator couplings. The expectations for
number of shots are calculated with respect to distribution over the uniform
spherical measure of states (Uniform spherical) as well as the pure state of
the uncoupled QDOs (Pure uncoupled). For ungrouped results all Pauli operators
in the Hamiltonian are measured separately, for grouped results Pauli
operators are measured in groups according to the method described in Appendix
C Theorem 4. Figure 3: Estimated number of shots required to measure
$\langle\hat{H}(R_{\text{vdW}})\rangle$ with error
$\epsilon=10\%\times\Delta{E}(R_{\text{vdW}})$ for a range of different
molecules. Here $R_{\text{vdW}}$ is the van der Waals diameter of the
respective molecules. The expected number of shots is calculated with respect
to the spherically uniform measure [(a) and (b)] as well as with the uncoupled
state [(c) and (d)] with and without grouping commuting operators. The number
of shots is shown both as a function $N$ for fixed $d=4$ and as a function of
$d$ for fixed $N=2$. Note that in (d), the dashed and solid lines are the
indistinguishable from one another. Values for $R_{\text{vdW}}$ and $\alpha$
are taken from Refs. Maroulis  _et al._ (1997); Yamamoto  _et al._ (2006);
Schropp and Tavan (2008); Huang  _et al._ (2013).

When combining the single oscillator non-interacting and two-oscillator
coupling terms together to measure the total Hamiltonian, if all
$\gamma_{i,j}$ are non-zero there will be a total of
$(\log_{2}d)N+\frac{1}{8}d^{2}(\log_{2}d)^{2}N(N-1)$ non-identity Pauli
operators. Naively $\hat{H}$ can be measured by measuring each of these
operators using its own circuit, however the number of measurements required
can be reduced significantly by measuring mutually-commuting operators
concurrently. Many heuristic methods exist for partitioning the terms into
commuting subsets, see for example Refs. Jena  _et al._ (2019); Crawford  _et
al._ (2021); Verteletskyi  _et al._ (2020); Izmaylov  _et al._ (2019); Zhao
_et al._ (2020); Yen  _et al._ (2020); Gokhale  _et al._ (2020). From the
structure of our Hamiltonian we can derive an exact method for grouping terms
into subsets that qubit-wise commute and therefore obtain an upper bound on
the number of unique measurements that need to be performed. We find an
overall upper bound of $(N-1)\left(\frac{d}{2}\log_{2}d\right)^{2}+1$ or
$N\left(\frac{d}{2}\log_{2}d\right)^{2}+1$ circuits, for even and odd $N$
respectively, needed to measure $\hat{H}$. In general, for a Hamiltonian with
all-to-all $k$th order coupling, there is an upper bound of
$O\left((d\log_{2}d)^{k}N^{k-1}\right)$ scaling with the number of QDOs. Using
heuristic sorting methods the scaling with respect to $d$ can be improved
beyond that given by these upper bounds, however an exact scaling is not
known. For a description of the grouping method, proof of these scalings as
well as examples of the result of using heuristic grouping algorithms, see
Appendix C.

Since we expect that, for physically relevant systems, the ground states will
be dominated by low-lying Fock states, we do not expect the need to increase
$d$ significantly to model non-trivial systems. Therefore, we are primarily
concerned with the scaling with respect to the number of QDOs $N$. Appendix B
shows that this is the case for the pair of interacting QDOs and that $d=4$ is
a good approximation to the true ground state until very large and unphysical
couplings.

Furthermore, using this known grouping, we can estimate the expected number of
shots needed to measure $\langle\hat{H}\rangle$ with absolute error
$\epsilon$. Without knowing the ground state of an arbitrary Hamiltonian, we
calculate the expected number of shots with respect to two different
distributions of states in Hilbert space. The first distribution is the
uniform spherical measure on the full Hilbert space which assumes no knowledge
of the distribution of possible states (Watrous, 2018, Chapter 7). The second
is the pure state of the uncoupled ground state
$\ket{\psi_{0}}=\ket{0}^{\otimes mN}$ which we expect to be close to the
ground state of the system with small non-zero coupling. The scaling of the
expected number of shots required for these two distributions is given in
Table 1 for measuring all Pauli operators individually or by grouping them as
described. We see that in both distributions, grouping the measurements
improves the scaling from $O(N^{4})$ to $O(N^{3})$ with respect to $N$ and
significantly reduces the polynomial scaling with respect to $d$ in the case
of the uniform spherical measure. In all cases, the exact values as well as
proofs for the number of shots required can be found in Appendix D (the number
of shots is an upper bound in the case of the uniform spherical measure and
tight for the pure uncoupled state).

To demonstrate the number of shots required for modelling real molecules of
interest, Fig. 3 highlights the number of shots required to achieve a
measurement error on $\langle\hat{H}\rangle$ of $\epsilon=10\%\times\Delta
E(R_{\text{vdW}})$. Here we consider $\Delta E(R_{\text{vdW}})$ to be the
dispersion interaction between a pair of one-dimensional QDOs representing the
molecule evaluated at $R_{\text{vdW}}$ the Van der Waals diameter of the
molecules. Assuming distributions of quantum states of the uniform spherical
measure as well as a pure uncoupled state, we see that the grouping method
derived in this work can reduce the required number of shots by multiple
orders of magnitude. Further reductions may be possible if heuristic methods
are used to find more optimal groups.

For the prototypical orbital-based VQE methods there are in general $O(N^{4})$
Pauli operators that must be measured where $N$ is the size of the molecular
orbital basis set used. The best known scaling for the number of commuting
operator collections is $O(N^{3})$ Gokhale  _et al._ (2020) which, when
compared to the order $O(N)$ scaling in terms of the number of QDOs for the
dipole coupled Hamiltonian here, suggests that the system we consider in this
work may be significantly more scalable than those considered in orbital-based
VQE methods. To our knowledge, no attempt has been made to find a scaling for
the number of shots required for these.

###  II.4 Variational ansatz

For the results presented in this paper, we use a variational ansatz
consisting of layers of two-qubit operators acting on nearest-neighbour
qubits. Since the problem Hamiltonian is positive, the ground state will be
real-valued (up to a global phase) and we restrict ourselves to using real-
valued, two-qubit operators based on the optimal construction of general
$\text{SO}(4)$ operators found by Vatan and Williams Vatan and Williams
(2004). We removed two parameterised $R_{z}$ gates from their construction to
reduce the number of variational parameters. This naturally restricts the
available operations implemented by each block to a subset of $\text{SO}(4)$
but was found to still give reasonable results in simulation and experiment
while reducing the number of circuits required for gradient measurements. The
circuit used for two QDOs is shown in Fig. 4a which has a total of twelve
variational parameters, the two-qubit circuit blocks are given in Fig. 4b.

(a)

@C=1em @R=1.5em &—0⟩ 1
—0⟩ 1
—0⟩ 1
—0⟩ 121.5em1.2emQDO 1 341.5em1.2emQDO 2

(b)

@C=1.0em @R=1.5em & S 1 R_y(θ_1) R_z(θ_2) 1 S^†
S R R_y(θ_3) R_z(θ_4) R^† S^†

Figure 4: Circuits used in experiments. (a) Four qubit variational ansatz
representing the state of two $d=4$ QDOs. The upper (lower) qubit within the
two-qubit register for each QDO corresponds to the most (least) significant
bit in the binary encoding of the corresponding QDO Fock states. Each two-
qubit block consists of a four-parameter real-valued operator. (b) Composition
of the four-parameter real-valued operators. $R_{z}$ and $R_{y}$ gates in the
centre are Pauli Z and Y rotation gates, each one contains a single
variational parameter. The remaining gates are fixed with
$S=R_{z}\left(\frac{\pi}{2}\right)$ and $R=R_{y}\left(\frac{\pi}{2}\right)$.
This block is based on the general SO(4) operator from Ref. Vatan and Williams
(2004) with two parameterised $R_{z}$ gates removed.

##  III Experimental results using IBM quantum computers

Figure 5: Experimental measurements used for error subtraction. (a)
Measurements of $\langle\hat{H}(R)\rangle_{\rho(\bm{\theta}_{R})}$ where
$\bm{\theta}_{R}$ are the optimal parameters found for each separation $R$
evaluated on ibm_lagos quantum computer. As described in the main text, each
of the 200 experiments is performed for each value of $R$ before moving onto
the next one. (b) Mean values of the measurements for each $R$. (c) Exact
calculation of $\hat{H}$ with no noise for the optimised parameters. Solid
curve shows the true ground state. (d) Absolute difference between exact value
of $\langle\hat{H}\rangle$ and mean of noisy measurements. This shows that,
although the effect of noise varies with each experiment, by measuring each
$R$ in the order described, the effect of noise is approximately equal for
each data point. Data points at $R=\infty$ correspond to the uncoupled
Hamiltonian and parameters found using the optimisation procedure.

To demonstrate the feasibility of calculating weak intermolecular interactions
from the coarse grained model we represent two linear, non-polar molecules at
distance $R$ apart using two one-dimensional QDOs. We choose model parameters
matching those of I2 , which is known to produce a dimer molecular liquid
above its melting point of 114°C and a purely van der Waals liquid above its
boiling point of 180°C - in which dispersion forces dominate the cohesion.

For each QDO, we use a $d=4$ level Fock basis, leading to a total of four
qubits required for the two-molecule system. We model the pair of I2
molecules, both aligned with their molecular axis along the intermolecular
axis such that the one-dimensional QDO models their induced polarisation along
the molecular axis. We use a value for the polarisabilty of
$\alpha=\alpha_{zz}=14.5\text{\AA}^{3}$ the polarisability coefficient along
the molecular axis Maroulis  _et al._ (1997) and a harmonic trapping frequency
of $\hbar\omega=9.61\text{eV}$. For the coarse grained electronic model, it is
usual to choose the $\omega$ parameter (or alternatively the point charge $q$
and mass $\mu$ if working with raw Drude particle parameters). This is done to
capture a certain chosen interaction, such as the $C_{6}$ London dispersion
coefficient, and then further terms (such as $C_{8}$ and $C_{9}$ terms) arise
naturally from the model. In this case, we choose $\hbar\omega$ such that the
full three-dimensional QDO in this setup would give $C_{6}$ coefficient
matching experimental measurements of I2 Kendrick  _et al._ (2012). For the
pair of interacting molecules, the full three-dimensional model and the total
London dispersion energy for the pair of interacting molecules can be found
from repetitions of the one-dimensional model using the correct possibilities
for each remaining axis a discussion of which is deferred until Sec. IV.1. As
proof of principle, we only consider polarisability along the intermolecular
axis. To generate the dispersion energy as a function of inter-molecular
separation $R$, we perform the VQE procedure to calculate the ground state
energy of the coupled QDOs for a range of possible coupling constants.

For the experimental optimisation procedure we used, for $R=\infty$, the
quantum computer ibmq_santiago and, for all other values of $R$, ibmq_montreal
(For further details about the hardware used in this work, see Appendix E).
Each evaluation of $\langle\hat{H}\rangle$ is performed by running ten
circuits with 8192 shots each where the 21 Pauli operators needed to estimate
$\langle\hat{H}\rangle$ are collected into ten groups of operators that qubit-
wise commute as described in Sec. II.2. We evaluated gradients using the
parameter shift rule Mitarai  _et al._ (2018) and performed the optimisation
using the ADAM optimiser (learning rate 0.25, decay rates 0.9, 0.99) Kingma
and Ba (2014). For each value of $R$, the VQE procedure was run for 200
optimisation steps with the final set of parameters taken as the output.

Once the optimal set of parameters was found, we extracted the binding energy
from the noisy quantum computer using a simple scheme to subtract the effect
of device noise. The binding energy due to the London dispersion interaction
at a certain $R$ is then given by

$\Delta E(R)=E(R)-E(R=\infty),$ (7)

where $E(R)$ is the ground state energy at $R$ and with $E(R=\infty)$
corresponding to the ground state energy of the completely uncoupled state. By
modelling the output of the noisy quantum computers as a unitary channel
depending on variational parameters $\bm{\theta}$ along with depolarising
noise with error rate $\lambda$, the output state for a given set of
variational parameters is

$\rho(\bm{\theta})=(1-\lambda){|\psi(\bm{\theta})\rangle\\!\langle\psi(\bm{\theta})|}+\frac{\lambda}{2^{mN}}\mathds{1},$
(8)

where $\ket{\psi{(\bm{\theta})}}$ is the output state from a noiseless ansatz
with parameters $\bm{\theta}$. Assuming $\ket{\psi{(\bm{\theta})}}$ is the
required ground state, this corresponds to a noisy expectation value for the
Hamiltonian

$\langle\hat{H}(R)\rangle_{\rho(\bm{\theta})}=(1-\lambda)E(R)+\frac{\lambda}{2^{mN}}\operatorname*{Tr}[\hat{H}(R)].$
(9)

Therefore, assuming $\lambda$ is independent of the variational parameters and
does not change between experiments, the dispersion energy is given by

$\Delta
E(R)=\frac{\langle\hat{H}(R)\rangle_{\rho(\bm{\theta})}-\langle\hat{H}(R=\infty)\rangle_{\rho(\bm{\theta_{\infty}})}}{(1-\lambda)},$
(10)

where $\bm{\theta}_{\infty}$ are variational parameters found for the ground
state of $\hat{H}(R=\infty)$, here we have used that
$\operatorname*{Tr}[\hat{H}(R)]=\operatorname*{Tr}[\hat{H}_{\text{non-
int}}]=dN\;\forall R$. We use (10) to calculate the dispersion energy from
noisy measurements, estimating $\epsilon$ by measuring the fidelity of the
uncoupled ground state on the noisy computer
$\mathcal{F}=\operatorname*{Tr}\left[\rho(\bm{\theta}_{\infty}){|0\rangle\\!\langle
0|}^{\otimes mN}\right]=(1-\lambda)+\frac{\lambda}{2^{mN}}$. This can be done
without additional cost by reusing the counts measured when calculating
$\langle\hat{H}\rangle_{\rho(\bm{\theta}_{\infty})}$. We find that, for the
devices and circuits we use in this work, this method to subtract the device
noise gives accurate values for the dispersion energy.

To measure the dispersion energy for the optimised states we use the ibm_lagos
quantum computer. One practical consideration when performing this
subtraction, is that due to the significant time dependence of the noise
levels of the device, it is important to evaluate the expectations for
different values of $R$ in (10) with the same qubits on the same device as
close to simultaneously as possible. For the results presented here we took
the mean of 200 repetitions of each measurement of $\langle\hat{H}\rangle$
(with 8192 shots per circuit), performing the measurements for each $R$ one-
after-another before moving onto the next of the 200 repetitions. As many
measurements were performed as possible in a single job and all jobs were
submitted simultaneously to minimise the time between successive jobs. Fig. 5
shows how the level of noise varies significantly throughout the 200
evaluations but the drift seen in the expectation values is approximately
matched for each of the different values of $R$ when measured in this way.

The spectrum achieved by the optimisation procedure in this proof of principle
experiment is shown in Fig. 2. From this, we see that by using the variational
procedure combined with our ansatz, optimal parameters can be found using
noisy gradients calculated on the real device that correspond to energy values
close to the true ground state when evaluated without device noise.
Furthermore, although the raw experimental energies are much larger than even
the uncoupled energy of the system, using the simple subtraction method
described above, we are able to successfully measure dispersion energies on
the device closely matching the exact London dispersion interaction. From Fig.
2 we see that for all values of $R$ tested, the binding energy is within two
standard deviations of the exact value.

##  IV Extending the model

###  IV.1 Three-dimensional QDOs

As described in Sec. II.1, for a complete description of the dipole-dipole
interaction and higher order terms, full three-dimensional QDOs are required.
Here we demonstrate how our representation and variational algorithm can be
extended to three-dimensional QDO systems.

The Hamiltonian for such a system is given by

$\begin{split}\hat{H}=&\sum_{i=1}^{N}\left(\frac{1}{2\mu}\mathbf{\hat{p}}_{i}^{2}+\frac{1}{2}\mu\omega^{2}\mathbf{\hat{x}}_{i}^{2}\right)\\\
&+q^{2}\sum_{j>i}\frac{3(\mathbf{\hat{x}}_{i}\cdot\mathbf{R}_{i,j})(\mathbf{\hat{x}}_{j}\cdot\mathbf{R}_{i,j})-\mathbf{\hat{x}}_{i}\cdot\mathbf{\hat{x}}_{j}R_{i,j}^{2}}{R_{i,j}^{5}},\end{split}$
(11)

where $\mathbf{\hat{x}}_{i}=(\hat{x}_{i},\hat{y}_{i},\hat{z}_{i})$ and
$\mathbf{\hat{p}}_{i}=(\hat{p}_{x_{i}},\hat{p}_{y_{i}},\hat{p}_{z_{i}})$ are
the position and momentum operators for QDO $i$ in each of the three Cartesian
coordinates. $\mathbf{R}_{i,j}$ is the vector separating the centres of the
pairs of QDOs $i$ and $j$ and $\mu$, $\omega$ and $q$ represent the effective
mass, natural frequency and charge for the QDO particles which, for
simplicity, are assumed to be isotropic and equal for all oscillators.

To represent the state of these three-dimensional oscillators on a quantum
computer, we treat each of the three spatial dimensions as a separate harmonic
oscillator with position and momentum operators
$\set{(\hat{x},\hat{p}_{x}),(\hat{y},\hat{p}_{y}),(\hat{z},\hat{p}_{z})}$. As
in Sec. II.2, each of these one-dimensional oscillators can be represented
using a binary encoding in $m$ qubits. Thus, the state of each three
dimensional oscillator state is encoded in a register of $3m$ qubits:

$\begin{split}&\mathcal{H}_{d-\text{Fock}}^{\otimes
3}\rightarrow\mathcal{H}_{2}^{\otimes m}\otimes\mathcal{H}_{2}^{\otimes
m}\otimes\mathcal{H}_{2}^{\otimes m}:\\\
&\ket{\underline{n}_{x},\underline{n}_{y},\underline{n}_{z}}\rightarrow\ket{\text{bin}_{m}(n_{x}),\text{bin}_{m}(n_{y}),\text{bin}_{m}(n_{z})}.\end{split}$
(12)

Following Sec. II.2, an upper bound for the number of circuits that are
required to measure $H$ is $(N-1)\left(\frac{d}{2}\log_{2}d\right)^{2}+1$ or
$N\left(\frac{d}{2}\log_{2}d\right)^{2}+1$, for even and odd $N$ respectively.
See Appendix C for details.

Extending to three dimensions and two interacting species, the Hamiltonian
simplifies to

$\hat{H}=\sum_{u\in\set{x,y,z}}(\hat{u}_{1}^{2}+\hat{p}_{u_{1}}^{2})+(\hat{u}_{2}^{2}+\hat{p}_{u_{2}}^{2})+\gamma_{uu}\hat{u}_{1}\hat{u}_{2},$
(13)

where $\gamma_{xx}=\nicefrac{{2\alpha_{xx}}}{{R^{3}}}$,
$\gamma_{yy}=\nicefrac{{2\alpha_{yy}}}{{R^{3}}}$ and
$\gamma_{zz}=\nicefrac{{-4\alpha_{zz}}}{{R^{3}}}$, $\alpha_{xx}$,
$\alpha_{yy}$, $\alpha_{zz}$ represent the dipole polarisability along each of
the three axes (which we now allow to be anisotropic, corresponding to $q$
being an anisotropic tensor) and the intermolecular axis is in the $z$
direction; along which the physical separation is $R$. We have ignored
dimensions to match the structure of (1). We see from this that the
Hamiltonian is equivalent to three independent one-dimensional oscillators
with the Hamiltonian matching (1) for $N=2$ with the appropriate choices of
coupling constants. Therefore, in the previous section where we modelled the
interaction of two linear molecules aligned along an axis, an extension to a
full three dimensional system follows trivially.

###  IV.2 Induction, anharmonic and non-linear terms

The quadratic Hamiltonians that we have presented so far can easily be
extended to include higher-order and anharmonic terms. Here we briefly
describe three extensions of the model giving rise to additional interactions
and phenomena. We also discuss how the decomposition and grouping strategies
described in Sec. II.3 can be applied to these richer systems.

The cheapest and simplest modification is the inclusion of induction effects
by applying a uniform electric field $E$ to the oscillators. This simple
modification can, for two or more oscillators, lead to non-additive effects
Cipcigan  _et al._ (2019). For one-dimensional oscillators this leads to an
interaction of the form
$\hat{H}_{\text{E-field}}=\sum_{i=1}^{N}q_{i}E\hat{x}_{i}$ for $N$ Drude
oscillators with charges $q_{i}$. These can be measured using the same
decompositions for $\hat{x}$ as were used previously since these interactions
are guaranteed to qubit-wise commute with terms of the form
$\hat{x}\otimes\hat{x}$ incurring no additional circuits to measure the total
Hamiltonian.

For more realistic induction models, we can instead consider oscillators in a
Coulomb potential caused by one or more charges. In this case, the Hamiltonian
will have the form
$\hat{H}_{\text{coulomb}}=\sum_{i=0}^{N}\sum_{j=0}^{M}\phi_{j}(\hat{x}_{i})$
for $N$ oscillators in Coulomb potentials represented by $\phi_{j}$ caused by
$M$ fixed charges. To measure a Hamiltonian of this form, each of the Coulomb
potentials must be Taylor expanded up to some finite order $w$:

$\phi_{j}(\hat{x}_{i})\sim\frac{1}{|r_{j}-\hat{x}_{i}|}\approx\frac{1}{r_{j}}+\frac{\hat{x}_{i}}{r_{j}^{2}}+\frac{\hat{x}_{i}^{2}}{r_{j}^{3}}+\dots+\frac{\hat{x}_{i}^{(w-1)}}{r_{j}^{w}},$
(14)

where $r_{j}\gg|\hat{x}_{i}|$ is the position of the point charge relative to
the centre of the QDO. Again, since each term of the form $\hat{x}_{i}^{w}$
acts on distinct oscillators they can all be measured simultaneously. The
expected contribution to the number of shots will be $O(\sqrt{N})$ with
respect to $N$ thus leaving the asymptotic scaling unchanged when added to the
linear dipolar Hamiltonian of (1).

Using a similar addition to the Hamiltonian we are able to relax the
assumption of a harmonic trapping potential and consider anharmonic trapping
potentials that have not been accessible for large systems using classical
computation methods Jones  _et al._ (2013b). Additional terms of the form
$\hat{H}_{\text{anharm}}=\omega_{3}\hat{x}_{i}^{3}+\dots+\omega_{w}\hat{x}_{i}^{w}$
with appropriate model parameters defining coefficients
$\set{\omega_{3},\omega_{4},\dots,\omega_{w}}$ can be used to choose a
desired, potentially asymmetric, trapping potential. The measurement cost will
scale in the same way as terms in (14). Relaxing the harmonic assumption is
required for inclusion of some higher order terms (such as four-dipole
hyperpolarisabilty Jones  _et al._ (2013b); Jones (2010)) and leads to non-
linear optical effects such as the Kerr effect, second- and third- harmonic
generation and, forms of optical mixing Boyd (2020).

Finally, we introduce non-linear interaction terms. As an example we consider
interactions of the form $\hat{H}_{\text{non-
lin}}=\sum_{i,j,k=0}^{N}\beta_{\langle
i,j,k\rangle}\hat{x}_{i}\hat{x}_{j}\hat{x}_{k}$. Terms similar to this have
been considered for use in quantum algorithms for vibrational spectra, noting
that these are difficult to deal with classical computational methods Sawaya
and Huh (2019); McArdle  _et al._ (2019). If all $O(N^{3})$ possible terms in
$\hat{H}_{\text{non-lin}}$ are non-zero, using the exact grouping that we have
described will contribute an $O(\frac{d^{9}N^{5}}{\epsilon^{2}})$ cost to the
number of shots required. In practice the number of shots may be considerably
lower as heuristic sorting methods will make use of additional qubit-wise
commutation not considered in the grouping method we have used in this work.
For gigher order interactions, these addition commuting pairs will likely be
more common. Additionally, many physical systems may not include all-to-all-
to-all interactions or will have interactions that decay strongly with
distance such that many of the interaction terms can be ignored. Terms of this
form, as well as higher order non-linear interaction terms will allow for a
large amount of flexibility in the model and will allow for the investigation
of a wide range of physical and chemical phenomena.

Figure 6: Ground state energy found using variational quantum algorithm for
two oscillator Hamiltonian extended beyond dipolar model. For VQA results,
optimisation was performed using a shot based simulation, without any
simulated device noise. The values of $E$ plotted were evaluated using state-
vector simulation, i.e., containing no shot statistics. True spectrum for this
Hamiltonian is shown in black dashed line, found by directly diagonalising the
Hamiltonian matrix. The energy for the simple dipole model is shown as the
black solid line, also found by direct diagonalisation.

As an example, we consider (in units of $\frac{1}{2}\hbar\omega$) the two
oscillator Hamiltonian

$\hat{H}=\hat{H}_{\text{non-
int}}+\hat{H}_{\text{int}}+\hat{H}_{\text{extension}},$ (15)

where $\hat{H}_{\text{non-
int}}=\hat{x}_{1}^{2}+\hat{p}_{1}^{2}+x_{2}^{2}+\hat{p}_{2}^{2}$,
$\hat{H}_{\text{int}}=\gamma\hat{x}_{1}\hat{x}_{2}$ are the same as the one-
dimension dipole system of (1) (with $N=2$) and

$\begin{split}\hat{H}_{\text{extension}}=&\frac{\gamma}{2}\hat{x}_{1}^{3}+\frac{\gamma}{2}\hat{x}_{2}^{3}+\frac{\gamma}{4}\hat{x}_{1}^{2}\hat{x}_{2}+\frac{\gamma}{4}\hat{x}_{1}\hat{x}_{2}^{2}\\\
&+\frac{1}{10}\hat{x}_{1}^{4}+\frac{1}{10}\hat{x}_{2}^{4},\end{split}$ (16)

which extends the model to include anharmonic onsite potentials (terms 1, 2, 5
and 6) and non-linear interactions (terms 3 and 4). Fig. 6 shows a simulated
result of using our variational algorithm with $d=4$ to find the ground state
energy. This simulation was performed using the same optimisation method,
settings and variational ansatz as in Sec. III.

##  V Discussion and conclusions

In this work we have developed a variational quantum algorithm for modelling
intermolecular interactions using a coarse grained framework. We have
introduced an encoding scheme, variational ansatz and variational algorithm
for digital quantum computers assessing measurement cost and scaling. We
presented an exact method for grouping mutually commuting operators that can
be measured simultaneously which gave bounds on the number of different
circuits and number of shots required to measure a system with dipole-dipole
interactions. We have shown evidence that this scaling may be more favourable
than for existing VQE methods. We have also presented a pathway for adding
anharmonic and higher order terms that are not easily accessible to classical
computational methods. This, combined with the favourable scaling, suggests
that our algorithm may be a more promising candidate for showing quantum
advantage than other variational methods.

We have further experimentally demonstrated that van der Waals dispersion
energies are resolvable on current generation quantum computing hardware and
that they lie within physically realistic ranges found in molecular materials.
This represents the first use of existing quantum computing capability to
access the weak but influential intermolecular energy scales arising from non-
covalent interactions (including correlated quantum fluctuations) relevant in
condensed states of molecular assemblies. While the specific model used in
this experiment can also be solved on classical processors, the impetus to
develop and demonstrate an efficient and well characterised quantum
implementation includes the following:

Firstly, it expands the palette of problems that can can be accessed on
quantum computing systems by providing a further level of electronic coarse
graining while preserving a complete set of responses, fluctuations and
correlation phenomena at long range. This methodology could therefore be
combined with other electronic structure strategies allowing only selected
electronic subsystems to be treated at the highest level of accuracy. Such
combined hybrid approaches are likely to be required to treat complex or
multiscale problems even if full CI computations can be realised in polynomial
time complexity. Having established that dispersion interactions are
detectable, extensions to the current implementation can be included to move
up the interaction hierarchy to capture induction and many-body polarisation
thereby opening opportunities for versatile molecular simulation problems on
quantum architecture.

Secondly, this framework may also be used to enrich existing simple models
such as bead based descriptions of protein and branched polymer folding which
can approach classically intractable combinatorial complexity and are already
being actively explored for implementation in quantum architecture Fingerhuth
_et al._ (2018); Robert  _et al._ (2021). Including coarse-grained
representations of the long-ranged interactions will result in more realistic
treatments of model folding problems.

###### Acknowledgements.

LWA, MK, DJ acknowledge support from the EPSRC National Quantum Technology Hub
in Networked Quantum Information Technology (EP/M013243/1) and the EPSRC Hub
in Quantum Computing and Simulation (EP/T001062/1). MK and DJ acknowledge
financial support from the National Research Foundation, Prime Ministers
Office, Singapore, and the Ministry of Education, Singapore, under the
Research Centres of Excellence program. The authors would like to acknowledge
the use of the University of Oxford Advanced Research Computing (ARC) facility
in carrying out this work ric . This work was also supported by the Hartree
National Centre for Digital Innovation programme, funded by the Department for
Business, Energy and Industrial Strategy, United Kingdom.

##  Appendix A London dispersion in linear molecules and one-dimensional
quantum Drude oscillators

This section gives shows how the interaction between two axially symmetric
linear molecules can be captured by one-dimensional QDOs. The results and
derivations given here can be found in the literature but are combined and
included here for the reader.

In the framework developed by London London (1942), the interaction between
two anisotropic linear molecules that are axially symmetric about their
respective principle intramolecular bond can be described by ellipsoids of
charge. For a pair of identical linear molecules that share components of the
dipole polarisability $\alpha_{\parallel}$ $(\alpha_{\perp})$ and
characteristic angular frequencies $\omega_{\parallel}$ $(\omega_{\perp})$
parallel (perpendicular) to their intramolecular axis, the London Dispersion
energy at large separation $R$ is

$\begin{split}\Delta
E=-\frac{1}{R^{6}}\bigg{[}&(C_{\parallel,\parallel}+C_{\perp,\perp}-C_{\parallel,\perp})\left\\{\sin^{2}(\theta_{a})\sin^{2}(\theta_{b})\cos(\phi_{a}-\phi_{b})-2\cos(\theta_{a})\cos(\theta_{b})\right\\}^{2}\\\
&+3(C_{\parallel,\perp}-C_{\perp,\perp})\left\\{\cos^{2}(\theta_{a})+\cos^{2}(\theta_{b})\right\\}+2(C_{\parallel,\perp}+4C_{\perp,\perp})\bigg{]},\end{split}$
(17)

where $\theta_{a},\phi_{a},\theta_{b},\phi_{b}$ describe the orientation of
molecules $a$ and $b$ relative to the intramolecular axis. $C_{\cdot,\cdot}$
coefficients are given by

$\begin{split}C_{\parallel,\parallel}=\frac{1}{8}\alpha_{\parallel}^{2}\hbar\omega_{\parallel},\quad&C_{\perp,\perp}=\frac{1}{8}\alpha_{\perp}^{2}\hbar\omega_{\perp},\\\
C_{\parallel,\perp}=\frac{1}{4}\alpha_{\parallel}&\alpha_{\perp}\hbar\frac{\omega_{\parallel}\omega_{\perp}}{\omega_{\parallel}+\omega_{\perp}}.\end{split}$
(18)

For further details about the geometry and derivation see Ref Hirschfelder
_et al._ (1964). When the two molecules are aligned parallel to each other and
the intermolecular axis $\theta_{a}=\theta_{b}=0$ and $\phi_{a}=\phi_{b}$
(which can be set to $0$ without loss of generality) in the completely
anisotropic limit $\alpha_{\parallel}\gg\alpha_{\perp}$ the London dispersion
interaction is given by

$\Delta E=-\frac{\alpha_{\parallel}^{2}\hbar\omega_{\parallel}}{2R^{6}}$ (19)

We now show how two interacting one-dimensional QDOs can recreate this
dispersion interaction. The Hamiltonian for pair of one-dimensional QDO is

$\hat{H}=\frac{\hat{p}_{1}^{2}}{2\mu}+\frac{\hat{p}_{2}^{2}}{2\mu}+\frac{1}{2}\mu\omega^{2}(\hat{x}_{1}^{2}+\hat{x}_{2}^{2}+\gamma\hat{x}_{1}\hat{x}_{2}).$
(20)

where $\gamma=-4\frac{q^{2}}{\mu\omega^{2}}\frac{1}{R^{3}}$. Alternatively
$\gamma=-4\frac{\alpha}{R^{3}}$ when written in terms of the dipole
polarisability $\alpha=\frac{q^{2}}{\mu\omega}$. Defining the transformed
variables

$\hat{p}_{\pm}=\frac{1}{\sqrt{2}}(\hat{p}_{1}\pm\hat{p}_{2}),\qquad\hat{x}_{\pm}=\frac{1}{\sqrt{2}}(\hat{x}_{1}\pm\hat{x}_{2}),$
(21)

the Hamiltonian can be rewritten in this new uncoupled basis as

$\begin{split}\hat{H}=&\left(\frac{\hat{p}_{+}^{2}}{2\mu}+\frac{1}{2}\mu\omega^{2}\left(1+\frac{\gamma}{2}\right)\hat{x}_{+}^{2}\right)\\\
&+\left(\frac{\hat{p}_{-}^{2}}{2\mu}+\frac{1}{2}\mu\omega^{2}\left(1-\frac{\gamma}{2}\right)\hat{x}_{-}^{2}\right).\end{split}$
(22)

Therefore the ground state energy is given by

$\begin{split}E&=\frac{1}{2}\hbar\omega\left(\sqrt{1+\frac{\gamma}{2}}+\sqrt{1-\frac{\gamma}{2}}\right)\\\
&=\frac{1}{2}\hbar\omega\left(2-\frac{\gamma^{2}}{16}+O(\gamma^{4})\right).\end{split}$
(23)

Using the definition $\gamma=-\frac{4\alpha}{R^{3}}$ gives the $R^{6}$ term
corresponding for the oscillator pairs

$\Delta
E=-\frac{1}{32}\gamma^{2}\hbar\omega=-\frac{1}{2}\frac{\alpha^{2}}{R^{6}}\hbar\omega.$
(24)

##  Appendix B Truncated Fock space

For our one-dimensional QDO system, we want to be sure that using a finite
dimensional Fock space for each QDO still gives a reasonable approximation to
the ground state energy of the true infinite-dimensional system. Fig. 7 shows
how the ground state energy of the system deviates for different values of the
local QDO dimension $d$. We can see that even for $d=4$, the ground state
energy is approximated well until very large couplings $\gamma$. For
reference, $\gamma=2$ corresponds to a Hamiltonian that is no longer positive-
semi definite, at this point the charged Drude particles dissociate from the
nuclei and the system breaks down. For physically relevant systems, we expect
to require value of $\gamma$ much lower than 2 and we can therefore limit
ourselves to small $d$. For the real device experiments in the main text, the
smallest separation $R$ considered corresponds to the largest coupling value
of $\gamma=1.55$.

Figure 7: Ground state energy of one-dimensional QDO pair as a function of
coupling constant $\gamma$ for various values of $d$; the dimensionality of
the truncated Fock space for each QDO. Analytic refers to exact analytic
calculation of the ground state.

##  Appendix C Decomposition and measurement of Hamiltonian

###  C.1 Scaling in the total number of Pauli operators

######  Theorem 1.

For a $d$-level harmonic QDO, there are
$\left(\frac{1}{2}d\log_{2}d\right)^{2}$ non-zero Pauli terms in the binary
encoding of the $\hat{x}\otimes\hat{x}$ interaction.

###### Proof.

The following single qubit operators can be represented in the Pauli basis

$\begin{split}\ket{0}\bra{0}&=\frac{1}{2}(\mathds{1}+Z),\\\
\ket{1}\bra{0}&=\frac{1}{2}(X-iY),\\\
\end{split}\qquad\begin{split}\ket{1}\bra{1}&=\frac{1}{2}(\mathds{1}-Z),\\\
\ket{0}\bra{1}&=\frac{1}{2}(X+iY).\\\ \end{split}$ (25)

We first consider a $\hat{x}=\frac{1}{\sqrt{2}}(\hat{a}^{\dagger}+\hat{a})$
for a single QDO, where the ladder operators $\hat{a}^{\dagger}$ and $\hat{a}$
for the $d$-level Fock space are given by

$\begin{split}\hat{a}^{\dagger}&=\sum_{n=0}^{d-1}\sqrt{n+1}\ket{\underline{n+1}}\bra{\underline{n}},\\\
\hat{a}^{\dagger}&=\sum_{n=0}^{d-1}\sqrt{n+1}\ket{\underline{n}}\bra{\underline{n+1}},\end{split}$
(26)

which follow the usual bosonic commutation relations. We begin by considering
a single $\ket{\underline{n+1}}\bra{\underline{n}}$ which given in terms of
$k$ single qubit operators is

$\begin{split}\ket{\underline{n+1}}\bra{\underline{n}}=&\left(\bigotimes_{i=k+1}^{m-1}\ket{\operatorname*{bin}(n)_{i}}\bra{\operatorname*{bin}(n)_{i}}\right)\\\
&\otimes\ket{1}\bra{0}\otimes\left(\ket{0}\bra{1}\right)^{\otimes
k},\end{split}$ (27)

where $\operatorname*{bin}(n)_{i}\in\set{0,1}$ is the binary digit at position
$i$ of the binary representation of $n$ and
$k(n)=\max\set{i\mid\operatorname*{bin}(n)_{i}=0}$ is the position of the
right-most $0$. Using the the definitions in (25) this is equal to

$\begin{split}\ket{\underline{n+1}}\bra{\underline{n}}=&\left(\bigotimes_{i=k+1}^{m-1}(\mathds{1}+(-1)^{\operatorname*{bin}(n)_{i}}Z)\right)\\\
&\otimes(X+iY)\otimes\left(X-iY\right)^{\otimes k}.\end{split}$ (28)

Expanding (28) will give $2^{m}$ unique Pauli terms of the form
$\hat{P}=c\sigma_{i}^{\otimes m}$ where $\sigma_{i}^{\otimes
m}\in\set{\mathds{1},X,Y,Z}^{m}$ and $c\in\set{\pm 1,\pm i}$. Each $\hat{P}$
will either be hermitian or anti-hermitian, depending on whether $c=\pm 1$ or
$c=\pm i$.

For each term in $\hat{x}$ like
$(\ket{\underline{n+1}}\bra{\underline{n}}+\ket{\underline{n}}\bra{\underline{n+1}})$,
we note that
$\ket{\underline{n+1}}\bra{\underline{n}}=\left(\ket{\underline{n}}\bra{\underline{n+1}}\right)^{\dagger}$
and their sum is Hermitian. Therefore the Pauli terms within this summed term
are simply two times the Hermitian terms in the expansion of (28) (anti-
Hermitian terms cancel). Of the $2^{m}$ different $\hat{P}$ within (28), half
will be Hermitian ($c=\pm 1$) and half will be anti-Hermitian ($c=\pm i$).
Therefore, there are a total of $\frac{1}{2}2^{m}$ unique Pauli operators in
the compact representation of
$(\ket{\underline{n+1}}\bra{\underline{n}}+\ket{\underline{n}}\bra{\underline{n+1}})$.

Finally, in the summation

$\hat{x}=\sum_{n=0}^{d-1}\sqrt{n+1}(\ket{\underline{n+1}}\bra{\underline{n}}+\ket{\underline{n}}\bra{\underline{n+1}}),$
(29)

there are $\log_{2}d=m$ unique values of $k(n)$, each of which will contribute
$\frac{2^{m}}{2}$ unique Pauli operators to the binary encoding of $\hat{x}$.
Thus there are a total of $\left(\frac{1}{2}d\log_{2}d\right)^{2}$ Pauli
operators within $\hat{x}\otimes\hat{x}$, where we have used $m=\log_{2}d$. ∎

###  C.2 Scaling of Number of Circuits Required

######  Lemma 2.

The $\left(\frac{1}{2}d\log_{2}d\right)^{2}$ Pauli operators for the
$\hat{x}\otimes\hat{x}$ interaction in the binary encoding can be measured
using $(d-1)^{2}$ circuits.

###### Proof.

We define the cover of an $M$ qubit Pauli operator
$\hat{P}\in\set{\mathds{1},X,Y,Z}^{M}$ as

$\begin{split}\operatorname*{Cover}(\hat{P})=\big{\\{}\hat{B}&\in\set{\mathds{1},X,Y,Z}^{\otimes
M}\\\ &|\quad\hat{B}^{(i)}=\hat{P}^{(i)}\text{ when
}\hat{P}^{(i)}\neq\mathds{1}\big{\\}},\end{split}$ (30)

where $\hat{P}^{(i)}$ is Pauli operator acting on the $i$th qubit. For a
collection of Pauli operators
$\set{\hat{P}_{1},\hat{P}_{2},\hat{P}_{3},\dots}$ that qubit-wise commute,
there exists at least one basis $\hat{B}$ which covers all of the operators:

$\hat{B}\in\operatorname*{Cover}(\hat{P}_{1})\cup\operatorname*{Cover}(\hat{P}_{2})\cup\operatorname*{Cover}(\hat{P}_{3})\cup\dots.$
(31)

The expectation values of every operator in
$\set{\hat{P}_{1},\hat{P}_{2},\hat{P}_{3},\dots}$ can be measured using a
single circuit in which qubit $i$ is measured in the $\hat{B}^{(i)}$ basis,
where the $\pm 1$ measurement value is used if $\hat{P}^{(i)}\neq\mathds{1}$,
and discarded otherwise. Here we called bases independent if they do not qubit
wise commute.

We show that there are $(d-1)$ independent covers for terms in $\hat{x}$ and
therefore $(d-1)^{2}$ covers needed to measure $\hat{x}\otimes\hat{x}$.

We begin by considering the expansion of
$\ket{\underline{n+1}}\bra{\underline{n}}$ given in (28). Expanding out the
$\bigotimes_{i=k+1}^{m-1}(\mathds{1}+(-1)^{\operatorname*{bin}(n)_{i}}Z)$ term
on the first $(m-k-1)$ qubits will give tensor products of only $\mathds{1}$
and $Z$. The covers for this $2^{(m-k-1)}$ operator subspace will always
contain the basis $\bigotimes_{i=k+1}^{m-1}Z^{(i)}$, which can be used to
measure them all. Then, expanding the $(X+iY)\otimes\left(X-iY\right)^{\otimes
k}$ term on the remaining $(m-k+1)$ qubits, gives $2^{(k+1)}$ terms that do
not qubit-wise commute and so each correspond to a unique basis.

Thus, there are $2^{(k+1)}$ independent bases for the terms in
$\ket{\underline{n+1}}\bra{\underline{n}}$. Using the argument that
$\ket{\underline{n+1}}\bra{\underline{n}}=\ket{\underline{n}}\bra{\underline{n+1}}^{\dagger}$
as in Theorem 4, only the Hermitian terms remain when summing the two. This
halves the number of unique covers and so $2^{k}$ bases are required for
$\ket{\underline{n+1}}\bra{\underline{n}}+\ket{\underline{n}}\bra{\underline{n+1}}$.

Again, for the sum of these terms given in (29), there are $m$ unique values
for $k(n)\in\set{0,1,\dots,m-1}$ when written in the form of (25). For each
unique value of $k(n)$, the value of $\operatorname*{bin}(n)$ does not affect
the coverings required (since the first $(m-k-1)$ qubits are always covered
with $Z^{\otimes(m-k-1)}$), so each unique $k(n)$ contributes a unique set of
$2^{k}$ coverings to the summation. Therefore, the total number of coverings
needed for the Pauli terms in $\hat{x}$ are $\sum_{k=0}^{m-1}2^{k}=2^{m}-1$,
where we have used a standard result for the summation.

Finally, the number of unique coverings for the $\hat{x}\otimes\hat{x}$
interaction is then $\left(2^{m}-1\right)^{2}=(d-1)^{2}$. ∎

######  Definition 1 (Hypergraphs Voloshin (2009)).

A hypergraph is a generalisation of a graph in which edges can connect
arbitrary number (i.e., not just a pair) of vertices. Formally, a hypergraph
$G=(V,E)$ consists of $V$ a finite set of vertices and $E$ a set of non-empty
subsets of $V$ called hyperedges.

######  Definition 2 (Complete Hypergraphs Voloshin (2009)).

For $0\leq r\leq N$, a complete hypergraph $K_{N}^{r}=(V,E)$ is a hypergraph
with vertices $V$ and set of hyperedges $E$ such that $|V|=N$ and $E$ contains
all $r$-subsets of $V$.

######  Lemma 3 (Baranyai’s Theorem Baranyai (1974)).

For a complete hypergraph $K_{N}^{r}$ where $2\leq r<N$ and $r$ divides $N$,
the $N$ vertices can be partitioned into subsets (called 1-factorisation) of
$r$ vertices in ${{N-1}\choose{r-1}}$ different ways, such that each
$r$-element subset (hyper edge) appears in exactly one partition.

###### Proof.

See (Van Lint  _et al._ , 2001, Chapter 38).∎

For example, the complete hypergraph $K_{4}^{2}$ contains four vertices, each
connected by two-element hyperedges (i.e, traditional edges). Baranyai’s
theorem states that the vertices can be partitioned into
${{4-1}\choose{2-1}}=3$ partitions. $K_{4}^{2}$ and its corresponding
partitions are shown in Fig. 8.

Figure 8: Fully connected graph $K^{2}_{4}$ along with the three corresponding
partitions of edges.

######  Theorem 4.

The number of circuits required for measurement of the Hamiltonian of the form

$\hat{H}=\sum_{i=1}^{N}(\hat{x}_{i}^{2}+\hat{p}_{i}^{2})+\sum_{j>i}\gamma_{i,j}\hat{x}_{i}\hat{x}_{j},$
(32)

for $N\geq 2$ $d$-level QDOs in the binary encoding is upper bounded by

$\displaystyle(d-1)^{2}(N-1)+1;\quad$ $\displaystyle N\text{ even}$ (33)
$\displaystyle(d-1)^{2}N+1;$ $\displaystyle N\text{ odd},$

which is $O(d^{2}N)$.

###### Proof.

All
$(\hat{x}^{2}+p^{2})=\sum_{n=0}^{d-1}n\ket{\underline{n}}\bra{\underline{n}}$
terms within the Hamiltonian can be decomposed into Pauli operators given in
(4) of the main text. All of these terms commute within each QDO register as
well between QDO registers and can thus be measured using the single covering
$Z^{\otimes Nm}$.

For the pair-wise coupling terms, from Theorem 2, a single
$\hat{x}\otimes\hat{x}$ term between two QDOs can be measured with $(d-1)^{2}$
circuits. The bases required to measure pairwise coupling terms between
disjoint pairs of QDOs are guaranteed to commute since they act on different
qubit registers. Therefore, by measuring each of the bases concurrently for
collections of mutually disjoint interactions, the number of measurement bases
can be reduced in a predictable way.

We represent interaction terms as a graph in which each vertex corresponds to
one of the $N$ QDOs and edges between vertices correspond represent a non-zero
pairwise coupling. For all-to-all interactions this will be a fully connected
graph. For even $N$ the number of partitions needed to separate interaction
edges into disjoint sets is equal by Baranyai’s theorem (Lemma 3) for the
graph $K_{N}^{2}$. For odd $N$ the edges cannot be partitioned exactly - there
will always be one unpaired vertex when collecting disjoint edges. Considering
a dummy vertex where an edge between that and an QDO vertex gives no
contribution to the Hamiltonian, the number of partitions is given by
Baranyai’s theorem for $K_{N+1}^{2}$.

Using these collections, the number of bases that will need to be measured of
the pairwise interaction terms is

$\displaystyle(d-1)^{2}{N-1\choose 1}$ $\displaystyle=(d-1)^{2}(N-1);\quad$
$\displaystyle N\text{ even},$ (34) $\displaystyle(d-1)^{2}{N\choose 1}$
$\displaystyle=(d-1)^{2}N;\quad$ $\displaystyle N\text{ odd}.$

Combining this with the single basis needed to measure the non-interaction
term gives an upper bound for the number of bases that must be measured and
thus the result. ∎

Figure 9: Number of circuits required to measure Hamiltonian as (a) and (b): a
function of $N$ for two fixed values of $d$ as well as (c) and (d): a function
of $d$ for two fixed values of $N$. Red dashed: upper bound scaling given by
Eqn. (33), black solid: using largest-degree first coloring for grouping
Paulis. Trivially, the number of circuits for $N=1$ is always 1.

This upper bound is plotted as a function of $N$ and $d$ in Fig. 9 as well as
the number of circuits required when using a heuristic sorting algorithm,
namely largest-degree first colouring Verteletskyi  _et al._ (2020). We see
that for $N=1$ both curves agree as there are no interaction terms and that
the heuristic sorted saturates the upper bound for $N=2$ QDOs. For $N>2$ we
see that the sorting significantly reduces the number of circuits below the
upper bound due to additional qubit wise commutation between non-disjoint
interaction terms.

######  Corollary 4.1.

The number of circuits required for measurement of the Hamiltonian for $N\geq
r$ $d$-level QDOs with $r$-order coupling in the binary encoding is
$O(d^{2}N^{r-1})$.

###### Proof.

The proof is the same as that for Theorem 4 except using $r>2$ for the
$\hat{x}^{\otimes r}$ interaction in Lemma 2 and the corresponding hypergraph
$K^{r}_{N}$ in Lemma 3. This gives an upper bound on the total number of
groups of

$(d-1)^{r}{N-1\choose r-1}+1,$ (35)

assuming $r$ divides $N$, which is $O(d^{r}N^{(r-1)})$. If $r$ does not divide
$N$, the $O(d^{r}N^{(r-1)})$ scaling will remain as dummy vertices can be used
to give the number of vertices that is divisible for $r$.

∎

###  C.3 Number of circuits for three-dimensional QDOs

######  Theorem 5.

The number of circuits required for measurement of the Hamiltonian of the form

$\begin{split}H=\sum_{i=1}^{N}&\left(\frac{\bf\hat{p}_{i}^{2}}{2\mu}+2\mu\omega\mathbf{\hat{x}}_{i}^{2}\right)\\\
&+q^{2}\sum_{j>i}\frac{3(\mathbf{\hat{x}}_{i}\cdot\mathbf{R}_{i,j})(\mathbf{\hat{x}}_{j}\cdot\mathbf{R}_{i,j})-\mathbf{\hat{x}}_{i}\cdot\mathbf{\hat{x}}_{j}R_{i,j}^{2}}{R_{i,j}^{5}},\end{split}$
(36)

for $N$ three-dimensional QDOs with $d$ levels per spatial dimension in the
binary encoding is upper bounded by

$\displaystyle(d-1)^{2}(N-1)+1;\quad$ $\displaystyle N\text{ even},$ (37)
$\displaystyle(d-1)^{2}N+1;$ $\displaystyle N\text{ odd}.$

###### Proof.

The uncoupled terms within (36) will correspond to operators
$\sum_{n=0}^{2^{m}-1}n\ket{\underline{n}}\bra{\underline{n}}$ for each of the
three spatial dimensions for each QDO. As before, these require just $Z$
measurements on each of the qubits and can thus be measured using a single
circuit.

Written in terms of operators for each one-dimensional component, the
interaction part of (36) will contain all nine terms of the form
$\set{\hat{x}_{i},\hat{y}_{i},\hat{z}_{i}}\otimes\set{\hat{x}_{j},\hat{y}_{j},\hat{z}_{j}}$,
where across each $3m$ qubit register for a single QDO, the operators
comprising $\mathbf{\hat{x}}_{i}$ are given by
$\hat{x}_{i}=\hat{x}\otimes\mathds{1}^{\otimes m}\otimes\mathds{1}^{\otimes
m}$, $\hat{y}_{i}=\mathds{1}^{\otimes
m}\otimes\hat{x}\otimes\mathds{1}^{\otimes m}$ and
$\hat{z}_{i}=\mathds{1}^{\otimes m}\otimes\mathds{1}^{\otimes
m}\otimes\hat{x}$, following a similar pattern for $\mathbf{\hat{x}}_{j}$,
$\bf\hat{p}_{i}$ and $\bf\hat{p}_{j}$. Since every single term of the form
$\hat{x}\hat{x}$, $\hat{x}\hat{y}$, $\hat{y}\hat{z}$ etc. in the coupling
interactions will consist of the same decomposition of Pauli operators either
acting on the same or different QDO subregisters. These all qubit wise commute
and each term each of the $(\frac{d}{2}\log_{2}d)^{2}$ Pauli terms for all
oscillator-oscillator interaction terms can be measured together. Thus for
each pair of coupled three-dimensional QDOs, only $(\frac{d}{2}\log_{2}d)^{2}$
circuits are needed; the same number required for a pairwise interaction
between one-dimensional QDOs.

Finally, using Lemma 3 and following the same logic as Theorem 4, an upper
bound for the number of circuits needed to measure in (36) is
$(N-1)\left(\frac{d}{2}\log_{2}d\right)^{2}+1$ or
$N\left(\frac{d}{2}\log_{2}d\right)^{2}+1$, for even and odd $N$ respectively.
∎

##  Appendix D Shot requirements

In this section, we derive an upper bound for an estimate of the number of
shots required to achieve a given statistical error on the measurement of
$\langle\hat{H}\rangle$ for a system of 1D coupled QDOs using the grouping of
measurements described in Theorem 4.

We write a Hamiltonian in the form

$\hat{H}=\sum_{g\in\mathcal{G}}\sum_{i\in g}a_{i}\hat{P}_{i},$ (38)

where $\mathcal{G}=\set{g}{[\hat{P}_{i},\hat{P}_{j}]=0,\forall i,j\in g}$ is
the set of all collections for the Pauli operatos $\hat{P}_{i}$ in which all
operators in a collection qubit wise commute. A standard method Wecker  _et
al._ (2015) using Lagrange multipliers to optimally distribute measurements
between circuits gives a total number experimental shots of

$S=\frac{1}{\epsilon^{2}}\left[\sum_{g\in\mathcal{G}}\sqrt{\sum_{i\in
g}|a_{i}|^{2}\operatorname*{Var}{\hat{P}_{i}}}\right]^{2}$ (39)

to achieve an expected measurement error of $\epsilon$, assuming that
$\operatorname*{Cov}[\hat{P}_{i},\hat{P}_{j}]=0$ if $i\neq j$. The variance of
a Pauli operator is given by
$\operatorname*{Var}{\hat{P}_{i}}=1-\langle\hat{P}_{i}\rangle$. In the absence
of knowledge about the state that the expectation value is take with respect
to, we take the expectation of this variance with respect to a spherical
measure (See (Watrous, 2018, Chapter 7) and following a similar method to Ref.
Crawford  _et al._ (2021)):

$\mathbb{E}[\operatorname*{Var}\hat{P}_{i}]=\begin{cases}0;\quad&\hat{P}_{i}=\mathds{1}\\\
1-\frac{1}{2^{n_{i}}+1};\quad&\text{otherwise}\end{cases}$ (40)

where $n_{i}$ is the number of non-identity Pauli operators in $\hat{P}_{i}$.
Additionally, in the uniform spherical measure,
$\mathbb{E}[\operatorname*{Cov}(\hat{P}_{i},\hat{P}_{j})]=0$ if $i\neq j$.
Thus the expected number of shots $S$ to achieve error $\epsilon$ is upper
bounded by

$\mathbb{E}[S]=\frac{1}{\epsilon^{2}}\left[\sum_{g\in\mathcal{G}}\sqrt{\sum_{i\in
g}|a_{i}|^{2}\left(1-\frac{1}{2^{n_{i}}+1}\right)}\right]^{2}.$ (41)

######  Theorem 6 (Expected shots over uniform measure when grouping
measurements).

Given a Hamiltonian

$\hat{H}=\sum_{i=1}^{N}(\hat{x}_{i}^{2}+\hat{p}_{i}^{2})+\sum_{j>i}\gamma_{i,j}\hat{x}_{i}\hat{x}_{j},$
(42)

for $N\geq 2$ $d$-level QDOs in the binary encoding and measuring comuting
operators according to the procedure described in Theorem 4, the expected
number of shots $S$ required to achieve a constant measurement error
$\epsilon$ is upper bounded by

$S\leq\frac{1}{\epsilon^{2}}\left[\frac{\sqrt{2}}{3}\sqrt{N}\sqrt{d^{2}-1}+\left(\sum_{F\in\operatorname*{Part}(G)}\sqrt{\sum_{\langle
i,j\rangle\in
F}\gamma_{i,j}^{2}}\right)\left(2d-\left[4+2\log_{2}d\right]d^{2}+\left[2+2\log_{2}d+\frac{1}{2}(\log_{2}d)^{2}\right]d^{3}\right)\right]^{2},$
(43)

when taken over the spherically symmetric distribution of states, where
$\operatorname*{Part}(G)$ is the set of coverings for a given 1-factorisation
of the interaction graph $G$. Assuming $G$ is fully connected and that values
of $\gamma_{i,j}$ are independent of both $N$ and $d$ and randomly distributed
between different $F\in\operatorname*{Part}(G)$, $S$ is
$O\left(\frac{\gamma^{2}N^{3}(\log_{2}d)^{4}d^{6}}{\epsilon^{2}}\right)$,
where $\gamma$ is a typical scale for $\gamma_{i,j}$.

###### Proof.

We consider the groupings described in Theorem 4. Since interacting and non-
interacting parts of the Hamiltonian have different forms and groupings, we
consider each one’s contribution to the total shots separately:

$\epsilon^{2}S=\left(\epsilon\tilde{S}_{\text{non-
int}}+\epsilon\tilde{S}_{\text{int}}\right)^{2}.$ (44)

Starting with the uncoupled part of the Hamiltonian:

$\sum_{i=1}^{N}(\hat{x}_{i}^{2}+\hat{p}_{i}^{2})=\sum_{i=1}^{N}2^{m}\mathds{1}^{\otimes
mN}-\sum_{j=0}^{m-1}2^{j}Z^{(im+j)},$ (45)

where $Z^{(im+j)}$ is the Z operator acting on the $j$th qubit of QDO $i$, all
terms form a single collection which gives

$\epsilon\tilde{S}_{\text{non-
int}}=\sqrt{N\sum_{j=0}^{m-1}2^{2j}\operatorname*{Var}Z}=\sqrt{\frac{2}{9}N\left(2^{2m}-1\right)},$
(46)

where we have used that $\mathbb{E}[\operatorname*{Var}Z]=\frac{2}{3}$ from
(40).

For the interacting terms, the collection to which a single Pauli operators
given by (28) (i.e., which $g\in\mathcal{G}$ it is sorted into) is determined
by three things:

  1. 1.

Which of the $(N-1)$ (even $N$) or $N$ (odd $N$) 1-factorisations of the
interaction graph the corresponding interaction term belongs to.

  2. 2.

The value of $k(n)$, the position of the rightmost zero in the binary
representation of integer $n$ labelling the Fock state, for _each of the two_
the QDOs involved in an interaction term.

  3. 3.

For a given pair of $k(n)$, which of the possible $2^{k+1}$ of $X$ and $Y$
operators is applied to the right most $(k+1)$ qubits.

As before, the anti-hermitian terms in (28) are not needed in the sum
$\ket{\underline{n+1}}\bra{\underline{n}}+\ket{\underline{n}}\bra{\underline{n+1}}$.

Separating the terms within $\hat{H}$ according to these criteria (their
contributions to the equation labelled with curly braces), the expected number
of shots required is

$\begin{split}\epsilon\tilde{S}_{\text{int}}=\frac{1}{2}&\underbrace{\sum_{k_{1},k_{2}=0}^{m-1}}_{1.}\;\underbrace{\sum_{F\in\operatorname*{Part}(G)}}_{2.}\;\\\
&\underbrace{\sum_{\begin{subarray}{c}\hat{P}_{1}^{(:k_{1}+1)}\in\set{X,Y}^{\otimes(k_{1}+1)}\\\
\hat{P}_{2}^{(:k_{2}+1)}\in\set{X,Y}^{\otimes(k_{2}+1)}\end{subarray}}}_{3.}\Gamma(F,k_{1},k_{2}),\end{split}$
(47)

where $\operatorname*{Part}(G)$ is the set of the interaction graph $G$. From
Lemma 3, if $G$ is fully connected then $|G|=(N-1)$ or $N$ (for $N$ even or
odd). We have defined
$\mathcal{N}_{k_{1}}=\set{n}{k(n)=k_{1},n\in\mathbb{Z},0\leq n<2^{m}}$ as the
set of integers satisfying a given $k(n)$ (similarly for
$\mathcal{N}_{k_{2}}$) and $\hat{P}_{1}^{(:k_{1}+1)}$ is a $(k_{1}+1)$ qubit
Pauli operator (similarly for $\hat{P}_{2}^{(:k_{2}+1)}$). The $1/2$ prefactor
results from the cancellation of anti-hermitian terms which make up half of
the operators. The value of $\Gamma$ is given by

$\displaystyle\Gamma(F,k_{1},k_{2})=$ $\displaystyle\Bigg{[}\sum_{\langle
i,j\rangle\in F}\sum_{\begin{subarray}{c}n_{1}\in\mathcal{N}_{k_{1}}\\\
n_{2}\in\mathcal{N}_{k_{2}}\end{subarray}}\;\sum_{\begin{subarray}{c}\hat{P}^{(k_{1}+1:m)}_{1}\in\set{\mathds{1},Z}^{\otimes(m-k_{1}-1)}\\\
\hat{P}_{2}^{(k_{2}+1:m)}\in\set{\mathds{1},Z}^{\otimes(m-k_{2}-1)}\end{subarray}}$
$\displaystyle\left(2\gamma_{i,j}\sqrt{n_{1}+1}\sqrt{n_{2}+1}\right)^{2}\mathbb{E}[\operatorname*{Var}{\hat{P}_{1}\otimes\hat{P}_{2}}]\Bigg{]}^{\frac{1}{2}},$
(48)

where, $\hat{P}^{(k_{1}+1:m)}_{1}$ is an $(m-(k_{1}+1))$ qubit Pauli operator
and $\hat{P}_{1}=\hat{P}_{1}^{(:k_{1}+1)}\otimes\hat{P}^{(k_{1}+1:m)}_{1}$
(similarly for $\hat{P}^{(k_{2}+1:m)}_{2}$). This is the equivalent of the
square root term in (41); the contribution to the variance from a single
collection of commuting operators. In the above, $\langle i,j\rangle$
correspond to edges within 1-factorisation $F$, $n_{1}$ and $n_{2}$ correspond
to Fock terms satisfying the requirement of $k(n)$; the position of the right-
most X or Y operator. Using the expectation on the spherical measure from (40)
gives

$\mathbb{E}[\operatorname*{Var}{\hat{P}_{1}\otimes\hat{P}_{2}}]=\left(1-\frac{1}{2^{k_{1}+k_{2}+\alpha+\beta}+1}\right)$
(49)

where integers $0<\alpha_{1}<m-k_{1}$ and $0<\alpha_{2}<m-k_{2}$ count the
number of Z operators in $\hat{P}_{1}^{(k_{1}+1:m)}$ and
$\hat{P}_{2}^{(k_{2}+1:m)}$ respectively.

By recognising that when using the spherical measure, many of the expected
variances are equivalent and only depend on the total number of X,Y or Z
operators, the expected number of shots simplifies to

$\epsilon\tilde{S}_{\text{int}}=\frac{1}{2}\sum_{k_{1},k_{2}=0}^{m-1}\;\sum_{F\in\operatorname*{Part}(K_{N}^{2})}\;2^{k_{1}+1}2^{k_{2}+1}\Gamma(F,k_{1},k_{2})$
(50)

with

$\begin{split}\Gamma(F,k_{1},k_{2})=\Bigg{[}&\sum_{\langle i,j\rangle\in
F}\sum_{\begin{subarray}{c}n_{1}\in\mathcal{N}_{k_{1}}\\\
n_{2}\in\mathcal{N}_{k_{2}}\end{subarray}}\;\sum_{\alpha=0}^{m-k_{1}-1}\sum_{\beta=0}^{m-k_{2}-1}\;\\\
&{m-k_{1}-1\choose\alpha}{m-k_{2}-1\choose\beta}\\\
&\times\left(2\gamma_{i,j}\sqrt{n_{1}+1}\sqrt{n_{2}+1}\right)^{2}\\\
&\times\left(1-\frac{1}{2^{k_{1}+k_{2}+\alpha+\beta}+1}\right)\Bigg{]}^{\frac{1}{2}}.\end{split}$
(51)

Now we aim to find a simple upper bound for $S$. We start by grouping terms

$\displaystyle\epsilon\tilde{S}_{\text{int}}=\left(\sum_{F\in\operatorname*{Part}(G)}\sqrt{\sum_{\langle
i,j\rangle\in F}\gamma_{i,j}^{2}}\right)\sum_{k_{1},k_{2}=1}^{m}$
$\displaystyle
2^{k_{1}+1}2^{k_{2}+1}\sqrt{\sum_{\begin{subarray}{c}n_{1}\in\mathcal{N}_{k_{1}}\\\
n_{2}\in\mathcal{N}_{k_{2}}\end{subarray}}(n_{1}+1)(n_{2}+1)}$
$\displaystyle\sqrt{\sum_{\alpha=0}^{m-k_{1}-1}\sum_{\beta=0}^{m-k_{2}-1}\;{m-k_{1}-1\choose\alpha}{m-k_{2}-1\choose\beta}\left(1-\frac{1}{2^{k_{1}+k_{2}+\alpha+\beta}+1}\right)}.$
(52)

We note that

$\begin{split}&\sum_{\begin{subarray}{c}n_{1}\in\mathcal{N}_{k_{1}}\\\
n_{2}\in\mathcal{N}_{k_{2}}\end{subarray}}(n_{1}+1)(n_{2}+1)\\\
&=\left(2^{k_{1}}\sum_{\gamma=0}^{2^{m-k_{1}-1}}(2\gamma+1)\right)\left(2^{k_{2}}\sum_{\delta=0}^{2^{m-k_{2}-1}}(2\delta+1)\right)\\\
&=2^{-k_{1}-2}(2^{k_{1}+1}+2^{m})^{2}\cdot
2^{-k_{2}-2}(2^{k_{2}+1}+2^{m})^{2}\end{split}$ (53)

and

$\begin{split}&\sum_{\alpha=0}^{m-k_{1}-1}\sum_{\beta=0}^{m-k_{2}-1}\;{m-k_{1}-1\choose\alpha}{m-k_{2}-1\choose\beta}\\\
&\hphantom{\sum_{\alpha=0}^{m-k_{1}-1}\sum_{\beta=0}^{m-k_{2}-1}\;}\times\left(1-\frac{1}{2^{k_{1}+k_{2}+\alpha+\beta}+1}\right)\\\
&\leq\sum_{\alpha=0}^{m-k_{1}-1}{m-k_{1}-1\choose\alpha}\sum_{\beta=0}^{m-k_{2}-1}{m-k_{1}-1\choose\beta}\\\
&\hphantom{\sum\leq}=2^{m-k_{1}-1}\cdot 2^{m-k_{2}-1},\end{split}$ (54)

where we have used a standard result for the summation of binomial
coefficients. Substituting (53) and (54) into (D) gives

$\begin{split}\epsilon
S_{\text{int}}\leq&\left(\sum_{F\in\operatorname*{Part}(G)}\sqrt{\sum_{\langle
i,j\rangle\in
F}\gamma_{i,j}^{2}}\right)\sum_{k_{1},k_{2}=0}^{m-1}2^{m-1}\left(4\cdot
2^{k_{1}+k_{2}}+2^{m+1}\cdot 2^{k_{1}}+2^{m+1}\cdot 2^{k_{2}}+2^{2m}\right)\\\
&=\left(\sum_{F\in(G)}\sqrt{\sum_{\langle i,j\rangle\in
F}\gamma_{i,j}^{2}}\right)\left(2d-\left[4+2\log_{2}d\right]d^{2}+\left[2+2\log
d+\frac{1}{2}(\log_{2}d)^{2}\right]d^{3}\right),\end{split}$ (55)

using standard results for the summations and substitution $d=2^{m}$.
Combining this with (46) gives the result. ∎

######  Theorem 7 (Expected shots over uniform measure when not grouping
measurements.).

Given a Hamiltonian

$\hat{H}=\sum_{i=1}^{N}(\hat{x}_{i}^{2}+\hat{p}_{i}^{2})+\sum_{j>i}\gamma_{i,j}\hat{x}_{i}\hat{x}_{j},$
(56)

for $N\geq 2$ $d$-level QDOs in the binary encoding and measuring each
corresponding Pauli operator separately, the expected number of shots required
to achieve a constant measurement error $\epsilon$ is given by the upper bound

$\begin{split}S\leq\frac{1}{\epsilon^{2}}\Bigg{[}&N\frac{\sqrt{6}}{3}(d-1)\\\
&+\left(\sum_{\langle i,j\rangle\in
G}|\gamma_{i,j}|\right)\left(\left[3+2\sqrt{2}\right]d^{2}-\left[2+\sqrt{2}\right]d^{\frac{5}{2}}-\left[\frac{7}{2}+3\sqrt{2}\right]d^{3}+\left[1+\sqrt{2}\right]d^{\frac{7}{2}}+\left[\frac{3}{2}+\sqrt{2}\right]d^{4}\right)\Bigg{]}^{2},\end{split}$
(57)

when taken over the spherically symmetric distribution of states, where $G$ is
the interaction graph. Assuming $G$ is fully connected and, that values of
$\gamma_{i,j}$ are independent of both $N$ and $d$, $S$ is
$O\left(\frac{\gamma^{2}N^{4}d^{8}}{\epsilon^{2}}\right)$ where $\gamma$ is a
typical scale for $\gamma_{i,j}$.

###### Proof.

As in Theorem 6, we consider contributions from non-interaction and
interaction terms separately. Starting with the uncoupled part of the
Hamiltonian

$\epsilon\tilde{S}_{\text{non-
int}}=N\sum_{j=0}^{m-1}\sqrt{2^{2j}\operatorname*{Var}Z}=N\frac{\sqrt{6}}{3}(2^{m}-1).$
(58)

For the interaction part of the Hamiltonian

$\begin{split}\epsilon\tilde{S}_{\text{int}}&=\frac{1}{2}\sum_{\langle
i,j\rangle\in
G}\sum_{k_{1},k_{2}=0}^{m-1}\;\sum_{\begin{subarray}{c}\hat{P}_{1}\in\set{X,Y}^{\otimes
k_{1}+1}\otimes\set{\mathds{1},Z}^{\otimes{m-k_{1}-1}}\\\
\hat{P}_{2}\in\set{X,Y}^{\otimes
k_{2}+1}\otimes\set{\mathds{1},Z}^{\otimes{m-k_{2}-1}}\end{subarray}}\;\left[\sum_{\begin{subarray}{c}n_{1}\in\mathcal{N}_{k_{1}}\\\
n_{2}\in\mathcal{N}_{k_{2}}\end{subarray}}4\gamma_{i,j}^{2}(n_{1}+1)(n_{2}+1)\mathbb{E}(\operatorname*{Var}\hat{P}_{1}\otimes\hat{P}_{2})\right]^{\frac{1}{2}}\\\
&=\sum_{\langle i,j\rangle\in
G}|\gamma_{i,j}|\sum_{k_{1},k_{2}=0}^{m-1}\;2^{k_{1}+1}2^{k_{2}+1}\sqrt{\sum_{\begin{subarray}{c}n_{1}\in\mathcal{N}_{k_{1}}\\\
n_{2}\in\mathcal{N}_{k_{2}}\end{subarray}}(n_{1}+1)(n_{2}+1)}\\\
&\hphantom{=\sum_{\langle i,j\rangle\in
G}|\gamma_{i,j}|\sum_{k_{1},k_{2}=0}^{m-1}\;}\sum_{\alpha=0}^{m-k_{1}-1}\sum_{\beta=0}^{m-k_{2}-1}\;{m-k_{1}-1\choose\alpha}{m-k_{2}-1\choose\beta}\sqrt{\left(1-\frac{1}{2^{k_{1}+k_{2}+\alpha+\beta}+1}\right)}.\end{split}$
(59)

Using (54) and (53) again and writing $d=2^{m}$, we get the inequality

$\begin{split}\epsilon\tilde{S}_{\text{int}}\leq&\sum_{\langle i,j\rangle\in
G}|\gamma_{i,j}|\sum_{k_{1},k_{2}=0}^{m-1}\;\left(2^{k_{1}+1}\cdot
2^{k_{2}+1}\right)\left(2^{m-k_{1}-1}\cdot
2^{m-k_{2}-1}\right)\left(2^{-\frac{k_{1}}{2}-1}(2^{k_{1}+1}+2^{m})\cdot
2^{-\frac{k_{2}}{2}-1}(2^{k_{2}+1}+2^{m})\right)\\\ &=\left(\sum_{\langle
i,j\rangle\in
G}|\gamma_{i,j}|\right)\left(\left[3+2\sqrt{2}\right]d^{2}-\left[2+\sqrt{2}\right]d^{\frac{5}{2}}-\left[\frac{7}{2}+3\sqrt{2}\right]d^{3}+\left[1+\sqrt{2}\right]d^{\frac{7}{2}}+\left[\frac{3}{2}+\sqrt{2}\right]d^{4}\right).\end{split}$
(60)

Combining this with (58) gives the result. ∎

The preceeding results assume no knowledge about the states that are being
measured and use an expectation over all of Hilbert space respect to the
spherically symmetric measure. Here we make a second estimate for the number
of shots required using the pure state $\ket{\psi_{0}}=\ket{0}^{\otimes mN}$.
$\ket{\psi_{0}}$ corresponds to the product state of all QDOs in their lowest
Fock state and is assumed to be reasonably close to the states that we seek
using the VQE procedure.

When not using the spherical measure, contributions to the measurement error
from the covarience of terms within the same collection is non-zero may be
non-zero. In this case, the expected number of shots is given by

$S=\frac{1}{\epsilon^{2}}\left[\sum_{g\in\mathcal{G}}\sqrt{\sum_{i,j\in
g}a_{i}a_{j}\operatorname*{Cov}[\hat{P}_{i},\hat{P}_{j}]}\right]^{2}.$ (61)

where
$\operatorname*{Cov}[\hat{P}_{i},\hat{P}_{j}]=\langle\hat{P}_{i}\hat{P}_{j}\rangle_{\ket{\psi_{0}}}-\langle\hat{P}_{i}\rangle_{\ket{\psi_{0}}}\langle\hat{P}_{j}\rangle_{\ket{\psi_{0}}}$
with subscripts indicating the state that expectation values are taken with
respect to.

######  Theorem 8 (Expected shots for uncoupled state when grouping
measurements.).

Given a Hamiltonian

$\hat{H}=\sum_{i=1}^{N}(\hat{x}_{i}^{2}+\hat{p}_{i}^{2})+\sum_{j>i}\gamma_{i,j}\hat{x}_{i}\hat{x}_{j},$
(62)

for $N\geq 2$ $d$-level QDOs in the binary encoding and measuring commuting
operators according to the procedure described in Theorem 4, the expected
number of shots $S$ required to achieve a constant measurement error
$\epsilon$ for state $\ket{0}^{\otimes mN}$ is upper bounded by

$S=\frac{1}{\epsilon^{2}}\frac{1}{16}\left(\sum_{F\in\operatorname*{Part}(G)}\;\sqrt{\sum_{\langle
i,j\rangle\in F}\gamma_{i,j}^{2}}\right)^{2}d^{4},$ (63)

where $\operatorname*{Part}(G)$ is the set of coverings for a given
1-factorisation of the interaction graph $G$. Assuming $G$ is fully connected
and that values of $\gamma_{i,j}$ are independent of both $N$ and $d$ and
randomly distributed between different $F\in\operatorname*{Part}(G)$, $S$ is
$O\left(\frac{\gamma^{2}N^{3}d^{4}}{\epsilon^{2}}\right)$, where $\gamma$ is a
typical scale for $\gamma_{i,j}$.

###### Proof.

Again we consider non-interacting and interacting contributions separately.
For the non interacting contribution, we note that for
$\ket{\psi}=\ket{0}^{\otimes mN}$, $\langle Z^{(i)}\rangle=\langle
Z^{(j)}\rangle=\langle Z^{(i)}Z^{(j)}\rangle=1\forall i,j$. Therefore,

$\epsilon\tilde{S}_{\text{non-
int}}=\sqrt{N\sum_{i,j=0}^{m-1}2^{i}2^{j}\operatorname*{Cov}[Z^{(i)},Z^{(j)}]}=0.$
(64)

The contribution from the coupling terms is

$\begin{split}\epsilon\tilde{S}_{\text{int}}=&\frac{1}{2}\sum_{k_{1},k_{2}=0}^{m-1}\;\sum_{F\in\operatorname*{Part}(K_{N}^{2})}\\\
&\sum_{\begin{subarray}{c}\hat{P}_{1}^{(:k_{1})}\in\set{X,Y}^{\otimes(k_{1}+1)}\\\
\hat{P}_{2}^{(:k_{2})}\in\set{X,Y}^{\otimes(k_{2}+1)}\end{subarray}}\Gamma(F,k_{1},k_{2}),\end{split}$
(65)

where $\operatorname*{Part}(K_{N}^{2})$ is the set of the $N-1$ or $N$ (for
$N$ even or odd) partitions of the interaction graph $K_{N}^{2}$,
$\mathcal{N}_{k_{1}}=\set{n}{k(n)=k_{1},n\in\mathbb{Z},0\leq n<2^{m}}$ is the
set of integers satisfying a given $k(n)$ (similarly for
$\mathcal{N}_{k_{2}}$) and $\hat{P}_{1}^{(:k_{1}+1)}$ is a $(k_{1}+1)$ qubit
Pauli operator (similarly for $\hat{P}_{2}^{(:k_{2}+1)}$). $\Gamma$ is given
by

$\begin{split}&\Gamma(G,k_{1},k_{2})\\\ &=\Bigg{[}\sum_{\langle i,j\rangle\in
F}\sum_{\langle i^{\prime},j^{\prime}\rangle\in
F}\sum_{\begin{subarray}{c}n_{1}\in\mathcal{N}_{k_{1}}\\\
n_{2}\in\mathcal{N}_{k_{2}}\end{subarray}}\;\sum_{\begin{subarray}{c}n_{1}^{\prime}\in\mathcal{N}_{k_{1}}\\\
n_{2}^{\prime}\in\mathcal{N}_{k_{2}}\end{subarray}}\\\
&\hphantom{=}2^{2}\gamma_{i,j}\gamma_{i^{\prime},j^{\prime}}\sqrt{n_{1}+1}\sqrt{n_{2}+1}\sqrt{n_{1}^{\prime}+1}\sqrt{n_{2}^{\prime}+1}\\\
&\hphantom{=}\times\operatorname*{Cov}\left[\hat{A}_{\langle
i,j\rangle}(k_{1},k_{2},n_{1},n_{2}),\hat{A}_{\langle
i^{\prime},j^{\prime}\rangle}(k_{1},k_{2},n_{1}^{\prime},n_{2}^{\prime})\right]\Bigg{]}^{\frac{1}{2}}\end{split}$
(66)

where

$\begin{split}&\hat{A}_{\langle i,j\rangle}(k_{1},k_{2},n_{1},n_{2})\\\
&=\left(\hat{P}_{1}^{(:k_{1})}\otimes\bigotimes_{\alpha=k_{1}+1}^{m-1}(\mathds{1}+(-1)^{\operatorname*{bin}(n_{1})_{\alpha}}Z)\right)_{\text{QDO
}i}\\\
&\hphantom{=}\otimes\left(\hat{P}_{2}^{(:k_{2})}\otimes\bigotimes_{\beta=k_{2}+1}^{m-1}(\mathds{1}+(-1)^{\operatorname*{bin}(n_{2})_{\beta}}Z)\right)_{\text{QDO
}j}\end{split}$ (67)

acts on the sub registers corresponding to QDOs $i$ and $j$.

We note that $\hat{P}_{1}^{(:k_{1})}\in\set{\pm
1}\times\set{X,Y}^{\otimes(k_{1}+1)}$ (recalling only Hermitian terms remain
in $\hat{H}$) and so
$\langle\hat{P}_{1}^{(:k_{1})}\rangle_{\ket{\psi_{0}}}=0\;\forall k_{1}$ and
similarly for $\hat{P}_{2}^{(:k_{2})}$, from which it follows

$\begin{split}\langle i&,j\rangle\neq\langle
i^{\prime},j^{\prime}\rangle\implies\\\
&\operatorname*{Cov}\left[\hat{A}_{\langle
i,j\rangle}(k_{1},k_{2},n_{1},n_{2}),\hat{A}_{\langle
i^{\prime},j^{\prime}\rangle}(k_{1},k_{2},n_{1}^{\prime},n_{2}^{\prime})\right]=0\end{split}$
(68)

and

$\begin{split}&\operatorname*{Cov}\left[\hat{A}_{\langle
i,j\rangle}(k_{1},k_{2},n_{1},n_{2}),\hat{A}_{\langle
i,j\rangle}(k_{1},k_{2},n_{1}^{\prime},n_{2}^{\prime})\right]\\\
&=\left\langle\hat{A}_{\langle
i,j\rangle}(k_{1},k_{2},n_{1},n_{2})\hat{A}_{\langle
i,j\rangle}(k_{1},k_{2},n_{1}^{\prime},n_{2}^{\prime})\right\rangle_{\ket{\psi_{0}}}\\\
&=\left\langle\hat{P}_{1}^{(:k_{1})}\hat{P}_{1}^{(:k_{1})}\right\rangle_{\ket{\psi_{0}}}\left\langle\hat{P}_{2}^{(:k_{2})}\hat{P}_{2}^{(:k_{2})}\right\rangle_{\ket{\psi_{0}}}\times\\\
&\hphantom{=}\left\langle\bigotimes_{\alpha=k_{1}+1}^{m-1}(\mathds{1}+(-1)^{\operatorname*{bin}(n_{1})_{\alpha}}Z)(\mathds{1}+(-1)^{\operatorname*{bin}(n_{1}^{\prime})_{\alpha}}Z)\right\rangle\\\
&\hphantom{=}\left\langle\bigotimes_{\alpha=k_{1}+1}^{m-1}(\mathds{1}+(-1)^{\operatorname*{bin}(n_{2})_{\alpha}}Z)(\mathds{1}+(-1)^{\operatorname*{bin}(n_{2}^{\prime})_{\alpha}}Z)\right\rangle,\end{split}$
(69)

where all expectation values are with respect to $\ket{\psi_{0}}$. Trivially
$\langle\hat{P}_{1}^{(:k_{1})}\hat{P}_{1}^{(:k_{1})}\rangle_{\ket{\psi_{0}}}=\langle\hat{P}_{2}^{(:k_{2})}\hat{P}_{2}^{(:k_{2})}\rangle_{\ket{\psi_{0}}}=1$.
It can be seen that

$\begin{split}&\left\langle\bigotimes_{\alpha=k_{1}+1}^{m-1}(\mathds{1}+(-1)^{\operatorname*{bin}(n_{1})_{\alpha}}Z)(\mathds{1}+(-1)^{\operatorname*{bin}(n_{1}^{\prime})_{\alpha}}Z)\right\rangle_{\ket{\psi_{0}}}\\\
&\qquad=\;\begin{cases}4^{m-k_{1}-1},&\mbox{if }n_{1}=n_{1}^{\prime}=0,\\\
0,&\mbox{otherwise},\end{cases}\end{split}$ (70)

and similarly for $n_{2}$ and $n_{2}^{\prime}$. Therefore

$\begin{split}&\operatorname*{Cov}\left[\hat{A}_{\langle
i,j\rangle}(k_{1},k_{2},n_{1},n_{2}),\hat{A}_{\langle
i^{\prime},j^{\prime}\rangle}(k_{1},k_{2},n_{1}^{\prime},n_{2}^{\prime})\right]\\\
&=2^{4m-2k_{1}-2k_{2}-4}\delta_{i,i^{\prime}}\delta_{j,j^{\prime}}\delta_{n_{1},0}\delta_{n_{2},0}\delta_{n_{1}^{\prime},0}\delta_{n_{2}^{\prime},0},\end{split}$
(71)

where $\delta$ is the Kronecker delta. We note that
$\delta_{n_{1},0}=\delta_{n_{1},0}\delta_{k_{1},0}\;\forall
k_{1}:n_{1}\in\mathcal{N}_{k_{1}}$ and similarly for $\delta_{n_{2},0}$. Then,
substituting this into (66) gives

$\Gamma(G,k_{1},k_{2})=\sum_{\langle
i,j\rangle}2^{2}\gamma_{i,j}^{2}2^{4m-k_{1}-k_{2}-4}\delta_{k_{1},0}\delta_{k_{2},0},$
(72)

which substituting into (65) gives

$\epsilon\tilde{S}_{\text{int}}=\frac{1}{4}\sum_{F\in\operatorname*{Part}{G}}\sqrt{\sum_{\langle
i,j\rangle\in F}\gamma_{i,j}^{2}d^{4}},$ (73)

from which the result follows. ∎

######  Theorem 9 (Expected shots for uncoupled state when not grouping
measurements.).

Given a Hamiltonian

$\hat{H}=\sum_{i=1}^{N}(\hat{x}_{i}^{2}+\hat{p}_{i}^{2})+\sum_{j>i}\gamma_{i,j}\hat{x}_{i}\hat{x}_{j},$
(74)

for $N\geq 2$ $d$-level QDOs in the binary encoding and measuring each
corresponding Pauli operator separately, the expected number of shots $S$
required to achieve a constant measurement error $\epsilon$ for state
$\ket{0}^{\otimes mN}$ is upper bounded by

$S=\frac{1}{\epsilon^{2}}\frac{1}{16}\left(\sum_{\langle i,j\rangle\in
G}|\gamma_{i,j}|\right)^{2}d^{4},$ (75)

where $G$ is the interaction graph. Assuming $G$ is fully connected and that
values of $\gamma_{i,j}$ are independent of both $N$ and $d$ and randomly
distributed between different $F\in\operatorname*{Part}(G)$, $S$ is
$O\left(\frac{\gamma^{2}N^{4}d^{4}}{\epsilon^{2}}\right)$, where $\gamma$ is a
typical scale for $\gamma_{i,j}$.

###### Proof.

Again we consider non-interacting and interacting contributions separately. As
in theorem 8, the uncoupled contribution is zero.

The contribution from the coupling terms is

$\begin{split}\epsilon\tilde{S}_{\text{int}}=&\frac{1}{2}\sum_{k_{1},k_{2}=0}^{m-1}\sum_{\langle
i,j\rangle\in
G}\sum_{\begin{subarray}{c}\hat{P}_{1}^{(:k_{1})}\in\set{X,Y}^{\otimes(k_{1}+1)}\\\
\hat{P}_{2}^{(:k_{2})}\in\set{X,Y}^{\otimes(k_{2}+1)}\end{subarray}}\\\
&\sum_{\begin{subarray}{c}n_{1}\in\mathcal{N}_{k_{1}}\\\
n_{2}\in\mathcal{N}_{k_{2}}\end{subarray}}\sqrt{2^{2}\gamma_{i,j}^{2}\operatorname*{Var}[\hat{A}_{\langle
i,j\rangle}(k_{1},k_{2},n_{1},n_{2})}]\end{split}$ (76)

where

$\displaystyle\hat{A}_{\langle i,j\rangle}(k_{1},k_{2},n_{1},n_{2})$
$\displaystyle=\left(\hat{P}_{1}^{(:k_{1})}\otimes\bigotimes_{\alpha=k_{1}+1}^{m-1}(\mathds{1}+(-1)^{\operatorname*{bin}(n_{1})_{\alpha}}Z)\right)_{\text{QDO
}i}$
$\displaystyle\hphantom{=}\otimes\left(\hat{P}_{2}^{(:k_{2})}\otimes\bigotimes_{\beta=k_{2}+1}^{m-1}(\mathds{1}+(-1)^{\operatorname*{bin}(n_{2})_{\beta}}Z)\right)_{\text{QDO
}j}$ (77)

acts on the sub registers corresponding to QDOs $i$ and $j$. Using (71), this
simplifies to

$\epsilon\tilde{S}_{\text{int}}=\frac{1}{4}\left(\sum_{\langle i,j\rangle\in
G}|\gamma_{i,j}|\right)d^{2},$ (78)

from which the result follows. ∎

##  Appendix E IBM Quantum Hardware

For practical reasons related to device availability and queue-times we used
three different quantum processors to evaluate our algorithm. All three
processors were provided by IBM via the Qiskit Anis  _et al._ (2021) open
source framework. We used ibmq_montreal (27 qubits, Processor type Falcon r4),
ibm_lagos (7 qubits, Processor type Falcon r5.11H) and ibmq_santiago (5
qubits, Processor type Falcon r4L). The reported quantum volume of the
aforementioned devices is 128, 32 and 32 respectively. Our experiments use 4
qubits on every processor, namely qubits 13,14, 16 and 19 in ibmq_montreal,
qubits 1,3,4 and 5 in ibm_lagos and qubits 1,2,3 and 4 in ibmq_santiago
following the coupling map show in Fig. 10. The experiments were performed in
a timeframe of 2 months and the noise parameters (and qubit properties) were
varied along the period of the experiments.

0 1 4 7 10 12 15 18 21 23 6 17 2 13 24 3 5 8 9 11 14 16 19 22 25 26 20 (a)

0 1 2 3 5 4 6 0 1 2 3 4 (b)(c)

Figure 10: Coupling map of the IBM Quantum devices used to produce the results
of the paper. The devices used are (a) ibmq_montreal, (b) ibm_lagos and (c)
ibmq_santiago with 27, 7 and 5 qubits respectively.

## References

  * Saito (2020) Kazuya Saito,  _Chemical Physics of Molecular Condensed Matter_ (Springer, 2020).
  * Reilly and Tkatchenko (2015) Anthony M Reilly and Alexandre Tkatchenko, “van der Waals dispersion interactions in molecular materials: beyond pairwise additivity,” Chemical Science 6, 3289–3301 (2015).
  * Shtogun and Woods (2010) Yaroslav V Shtogun and Lilia M Woods, “Many-body van der Waals interactions between graphitic nanostructures,” The Journal of Physical Chemistry Letters 1, 1356–1362 (2010).
  * Pastorczak  _et al._ (2017) Ewa Pastorczak, Jun Shen, Michał Hapka, Piotr Piecuch, and Katarzyna Pernal, “Intricacies of van der Waals interactions in systems with elongated bonds revealed by electron-groups embedding and high-level coupled-cluster approaches,” Journal of chemical theory and computation 13, 5404–5419 (2017).
  * Hermann  _et al._ (2017) Jan Hermann, Robert A DiStasio Jr, and Alexandre Tkatchenko, “First-principles models for van der Waals interactions in molecules and materials: Concepts, theory, and applications,” Chemical Reviews 117, 4714–4758 (2017).
  * Hellman (1937) Hans Hellman, “Einführung in die quantenchemie,” Franz Deuticke, Leipzig 285 (1937).
  * Feynman (1939) Richard Phillips Feynman, “Forces in molecules,” Physical review 56, 340 (1939).
  * Born and Oppenheimer (1927) Max Born and Robert Oppenheimer, “Zur quantentheorie der molekeln,” Annalen der physik 389, 457–484 (1927).
  * French  _et al._ (2010) Roger H. French, V. Adrian Parsegian, Rudolf Podgornik, Rick F. Rajter, Anand Jagota, Jian Luo, Dilip Asthagiri, Manoj K. Chaudhury, Yet-ming Chiang, Steve Granick, Sergei Kalinin, Mehran Kardar, Roland Kjellander, David C. Langreth, Jennifer Lewis, Steve Lustig, David Wesolowski, John S. Wettlaufer, Wai-Yim Ching, Mike Finnis, Frank Houlihan, O. Anatole von Lilienfeld, Carel Jan van Oss, and Thomas Zemb, “Long range interactions in nanoscale science,” Rev. Mod. Phys. 82, 1887–1944 (2010).
  * von Lilienfeld  _et al._ (2004) O. Anatole von Lilienfeld, Ivano Tavernelli, Ursula Rothlisberger, and Daniel Sebastiani, “Optimization of effective atom centered potentials for london dispersion forces in density functional theory,” Phys. Rev. Lett. 93, 153004 (2004).
  * Tkatchenko and Scheffler (2009) Alexandre Tkatchenko and Matthias Scheffler, “Accurate molecular van der waals interactions from ground-state electron density and free-atom reference data,” Phys. Rev. Lett. 102, 073005 (2009).
  * Grimme (2011) Stefan Grimme, “Density functional theory with london dispersion corrections,” Wiley Interdisciplinary Reviews: Computational Molecular Science 1, 211–228 (2011).
  * Slater (1929) John C Slater, “The theory of complex spectra,” Physical Review 34, 1293 (1929).
  * Boys (1950) S Francis Boys, “Electronic wave functions-i. a general method of calculation for the stationary states of any molecular system,” Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences 200, 542–554 (1950).
  * Szalay  _et al._ (2012) Peter G Szalay, Thomas Muller, Gergely Gidofalvi, Hans Lischka, and Ron Shepard, “Multiconfiguration self-consistent field and multireference configuration interaction methods and applications,” Chemical reviews 112, 108–181 (2012).
  * Čížek (1966) Jiří Čížek, “On the correlation problem in atomic and molecular systems. calculation of wavefunction components in ursell-type expansion using quantum-field theoretical methods,” The Journal of Chemical Physics 45, 4256–4266 (1966).
  * Bartlett and Musiał (2007) Rodney J Bartlett and Monika Musiał, “Coupled-cluster theory in quantum chemistry,” Reviews of Modern Physics 79, 291 (2007).
  * Pieniazek  _et al._ (2007) Piotr A Pieniazek, Stephen A Arnstein, Stephen E Bradforth, Anna I Krylov, and C David Sherrill, “Benchmark full configuration interaction and equation-of-motion coupled-cluster model with single and double substitutions for ionized systems results for prototypical charge transfer systems: Noncovalent ionized dimers,” The Journal of chemical physics 127, 164110 (2007).
  * Oliveira and Kraka (2017) Vytor Oliveira and Elfi Kraka, “Systematic coupled cluster study of noncovalent interactions involving halogens, chalcogens, and pnicogens,” The Journal of Physical Chemistry A 121, 9544–9556 (2017).
  * Bistoni (2020) Giovanni Bistoni, “Finding chemical concepts in the hilbert space: Coupled cluster analyses of noncovalent interactions,” Wiley Interdisciplinary Reviews: Computational Molecular Science 10, e1442 (2020).
  * Preskill (2018) John Preskill, “Quantum computing in the NISQ era and beyond,” Quantum 2, 79 (2018).
  * Cerezo  _et al._ (2021) Marco Cerezo, Andrew Arrasmith, Ryan Babbush, Simon C Benjamin, Suguru Endo, Keisuke Fujii, Jarrod R McClean, Kosuke Mitarai, Xiao Yuan, Lukasz Cincio, _et al._ , “Variational quantum algorithms,” Nature Reviews Physics , 1–20 (2021).
  * McArdle  _et al._ (2020) Sam McArdle, Suguru Endo, Alán Aspuru-Guzik, Simon C Benjamin, and Xiao Yuan, “Quantum computational chemistry,” Reviews of Modern Physics 92, 015003 (2020).
  * Peruzzo  _et al._ (2014) Alberto Peruzzo, Jarrod McClean, Peter Shadbolt, Man-Hong Yung, Xiao-Qi Zhou, Peter J Love, Alán Aspuru-Guzik, and Jeremy L O’brien, “A variational eigenvalue solver on a photonic quantum processor,” Nature communications 5, 1–7 (2014).
  * Farhi  _et al._ (2014) Edward Farhi, Jeffrey Goldstone, and Sam Gutmann, “A quantum approximate optimization algorithm,” arXiv preprint arXiv:1411.4028 (2014).
  * Farhi and Harrow (2016) Edward Farhi and Aram W Harrow, “Quantum supremacy through the quantum approximate optimization algorithm,” arXiv preprint arXiv:1602.07674 (2016).
  * Bauer  _et al._ (2016) Bela Bauer, Dave Wecker, Andrew J Millis, Matthew B Hastings, and Matthias Troyer, “Hybrid quantum-classical approach to correlated materials,” Physical Review X 6, 031045 (2016).
  * Kreula  _et al._ (2016a) Juha M Kreula, Laura García-Álvarez, Lucas Lamata, Stephen R Clark, Enrique Solano, and Dieter Jaksch, “Few-qubit quantum-classical simulation of strongly correlated lattice fermions,” EPJ Quantum Technology 3, 1–19 (2016a).
  * Kreula  _et al._ (2016b) JM Kreula, Stephen R Clark, and D Jaksch, “Non-linear quantum-classical scheme to simulate non-equilibrium strongly correlated fermionic many-body dynamics,” Scientific reports 6, 1–7 (2016b).
  * Lubasch  _et al._ (2020) Michael Lubasch, Jaewoo Joo, Pierre Moinier, Martin Kiffner, and Dieter Jaksch, “Variational quantum algorithms for nonlinear problems,” Physical Review A 101, 010301 (2020).
  * Lloyd  _et al._ (2020) Seth Lloyd, Giacomo De Palma, Can Gokler, Bobak Kiani, Zi-Wen Liu, Milad Marvian, Felix Tennie, and Tim Palmer, “Quantum algorithm for nonlinear differential equations,” arXiv preprint arXiv:2011.06571 (2020).
  * Ollitrault  _et al._ (2020a) Pauline J Ollitrault, Abhinav Kandala, Chun-Fu Chen, Panagiotis Kl Barkoutsos, Antonio Mezzacapo, Marco Pistoia, Sarah Sheldon, Stefan Woerner, Jay M Gambetta, and Ivano Tavernelli, “Quantum equation of motion for computing molecular excitation energies on a noisy quantum processor,” Physical Review Research 2, 043140 (2020a).
  * Grimsley  _et al._ (2019) Harper R Grimsley, Sophia E Economou, Edwin Barnes, and Nicholas J Mayhall, “An adaptive variational algorithm for exact molecular simulations on a quantum computer,” Nature communications 10, 1–9 (2019).
  * Kandala  _et al._ (2017) Abhinav Kandala, Antonio Mezzacapo, Kristan Temme, Maika Takita, Markus Brink, Jerry M Chow, and Jay M Gambetta, “Hardware-efficient variational quantum eigensolver for small molecules and quantum magnets,” Nature 549, 242–246 (2017).
  * Nam  _et al._ (2020) Yunseong Nam, Jwo-Sy Chen, Neal C Pisenti, Kenneth Wright, Conor Delaney, Dmitri Maslov, Kenneth R Brown, Stewart Allen, Jason M Amini, Joel Apisdorf, _et al._ , “Ground-state energy estimation of the water molecule on a trapped-ion quantum computer,” npj Quantum Information 6, 1–6 (2020).
  * Arute  _et al._ (2020) Frank Arute, Kunal Arya, Ryan Babbush, Dave Bacon, Joseph C Bardin, Rami Barends, Sergio Boixo, Michael Broughton, Bob B Buckley, David A Buell, _et al._ , “Hartree-Fock on a superconducting qubit quantum computer,” Science 369, 1084–1089 (2020).
  * McArdle  _et al._ (2019) Sam McArdle, Alexander Mayorov, Xiao Shan, Simon Benjamin, and Xiao Yuan, “Digital quantum simulation of molecular vibrations,” Chemical science 10, 5725–5735 (2019).
  * Ollitrault  _et al._ (2020b) Pauline J Ollitrault, Alberto Baiardi, Markus Reiher, and Ivano Tavernelli, “Hardware efficient quantum algorithms for vibrational structure calculations,” Chemical science 11, 6842–6855 (2020b).
  * Rungger  _et al._ (2019) Ivan Rungger, Nathan Fitzpatrick, Honxiang Chen, CH Alderete, Harriett Apel, Alexander Cowtan, Andrew Patterson, D Munoz Ramo, Yingyue Zhu, Nhung Hong Nguyen, _et al._ , “Dynamical mean field theory algorithm and experiment on quantum computers,” arXiv preprint arXiv:1910.04735 (2019).
  * Jaderberg  _et al._ (2020) Ben Jaderberg, Abhishek Agarwal, Karsten Leonhardt, Martin Kiffner, and Dieter Jaksch, “Minimum hardware requirements for hybrid quantum–classical dmft,” Quantum Science and Technology 5, 034015 (2020).
  * Helgaker  _et al._ (2014) Trygve Helgaker, Poul Jorgensen, and Jeppe Olsen,  _Molecular electronic-structure theory_ (John Wiley & Sons, 2014).
  * Sherrill (1996) C David Sherrill, “Computational scaling of the configuration interaction method with system size,” Quantum Chemistry Lecture Notes, School of Chem. and Biochem., Georgia Inst. of Technol. (1996).
  * Sherrill (2017) C David Sherrill, “Wavefunction theory approaches to noncovalent interactions,” in  _Non-Covalent Interactions in Quantum Chemistry and Physics_ (Elsevier, 2017) pp. 137–168.
  * Cipcigan  _et al._ (2019) FS Cipcigan, J Crain, VP Sokhan, and GJ Martyna, “Electronic coarse graining: Predictive atomistic modeling of condensed matter,” Reviews of Modern Physics 91, 025003 (2019).
  * Cipcigan  _et al._ (2018) Flaviu Cipcigan, Vlad Sokhan, Glenn Martyna, and Jason Crain, “Structure and hydrogen bonding at the limits of liquid water stability,” Scientific reports 8, 1–8 (2018).
  * Sokhan  _et al._ (2015a) Vlad P Sokhan, Andrew P Jones, Flaviu S Cipcigan, Jason Crain, and Glenn J Martyna, “Signature properties of water: Their molecular electronic origins,” Proceedings of the National Academy of Sciences 112, 6341–6346 (2015a).
  * Sokhan  _et al._ (2015b) VP Sokhan, A Jones, FS Cipcigan, J Crain, and GJ Martyna, “Molecular-scale remnants of the liquid-gas transition in supercritical polar fluids,” Physical review letters 115, 117801 (2015b).
  * Cipcigan  _et al._ (2016) Flaviu S Cipcigan, Vlad P Sokhan, Jason Crain, and Glenn J Martyna, “Electronic coarse graining enhances the predictive power of molecular simulation allowing challenges in water physics to be addressed,” Journal of Computational Physics 326, 222–233 (2016).
  * Gokhale  _et al._ (2020) Pranav Gokhale, Olivia Angiuli, Yongshan Ding, Kaiwen Gui, Teague Tomesh, Martin Suchara, Margaret Martonosi, and Frederic T Chong, “$O(N^{3})$ measurement cost for variational quantum eigensolver on molecular hamiltonians,” IEEE Transactions on Quantum Engineering 1, 1–24 (2020).
  * Westbroek  _et al._ (2018) Marise JE Westbroek, Peter R King, Dimitri D Vvedensky, and Stephan Dürr, “User’s guide to monte carlo methods for evaluating path integrals,” American Journal of Physics 86, 293–304 (2018).
  * Jones  _et al._ (2013a) A Jones, F Cipcigan, VP Sokhan, J Crain, and GJ Martyna, “Electronically coarse-grained model for water,” Physical review letters 110, 227801 (2013a).
  * Jones  _et al._ (2013b) Andrew P Jones, Jason Crain, Vlad P Sokhan, Troy W Whitfield, and Glenn J Martyna, “Quantum Drude oscillator model of atoms and molecules: Many-body polarization and dispersion interactions for atomistic simulation,” Physical Review B 87, 144103 (2013b).
  * Sawaya  _et al._ (2020) Nicolas PD Sawaya, Tim Menke, Thi Ha Kyaw, Sonika Johri, Alán Aspuru-Guzik, and Gian Giacomo Guerreschi, “Resource-efficient digital quantum simulation of d-level systems for photonic, vibrational, and spin-s hamiltonians,” npj Quantum Information 6, 1–13 (2020).
  * Dumitrescu  _et al._ (2018) Eugene F Dumitrescu, Alex J McCaskey, Gaute Hagen, Gustav R Jansen, Titus D Morris, T Papenbrock, Raphael C Pooser, David Jarvis Dean, and Pavel Lougovski, “Cloud quantum computing of an atomic nucleus,” Physical review letters 120, 210501 (2018).
  * Maroulis  _et al._ (1997) George Maroulis, Constantinos Makris, Uwe Hohm, and Dirk Goebel, “Electrooptical properties and molecular polarization of iodine, i2,” The Journal of Physical Chemistry A 101, 953–956 (1997).
  * Yamamoto  _et al._ (2006) Sekika Yamamoto, Yasuhiko Ishibashi, Yasuhiro Inamura, Yoshinori Katayama, Tomobumi Mishina, and Jun’ichiro Nakahara, “Pressure dependence of local structure in liquid carbon disulfide,” The journal of chemical physics 124, 144511 (2006).
  * Schropp and Tavan (2008) Bernhard Schropp and Paul Tavan, “The polarizability of point-polarizable water models: density functional theory/molecular mechanics results,” The Journal of Physical Chemistry B 112, 6233–6240 (2008).
  * Huang  _et al._ (2013) Yongli Huang, Xi Zhang, Zengsheng Ma, Wen Li, Yichun Zhou, Ji Zhou, Weitao Zheng, and Chang Q Sun, “Size, separation, structural order and mass density of molecules packing in water and ice,” Scientific reports 3, 1–5 (2013).
  * Jena  _et al._ (2019) Andrew Jena, Scott Genin, and Michele Mosca, “Pauli partitioning with respect to gate sets,” arXiv preprint arXiv:1907.07859 (2019).
  * Crawford  _et al._ (2021) Ophelia Crawford, Barnaby van Straaten, Daochen Wang, Thomas Parks, Earl Campbell, and Stephen Brierley, “Efficient quantum measurement of Pauli operators in the presence of finite sampling error,” Quantum 5, 385 (2021).
  * Verteletskyi  _et al._ (2020) Vladyslav Verteletskyi, Tzu-Ching Yen, and Artur F Izmaylov, “Measurement optimization in the variational quantum eigensolver using a minimum clique cover,” The Journal of chemical physics 152, 124114 (2020).
  * Izmaylov  _et al._ (2019) Artur F Izmaylov, Tzu-Ching Yen, Robert A Lang, and Vladyslav Verteletskyi, “Unitary partitioning approach to the measurement problem in the variational quantum eigensolver method,” Journal of chemical theory and computation 16, 190–195 (2019).
  * Zhao  _et al._ (2020) Andrew Zhao, Andrew Tranter, William M Kirby, Shu Fay Ung, Akimasa Miyake, and Peter J Love, “Measurement reduction in variational quantum algorithms,” Physical Review A 101, 062322 (2020).
  * Yen  _et al._ (2020) Tzu-Ching Yen, Vladyslav Verteletskyi, and Artur F Izmaylov, “Measuring all compatible operators in one series of single-qubit measurements using unitary transformations,” Journal of chemical theory and computation 16, 2400–2409 (2020).
  * Watrous (2018) John Watrous,  _The theory of quantum information_ (Cambridge university press, 2018).
  * Vatan and Williams (2004) Farrokh Vatan and Colin Williams, “Optimal quantum circuits for general two-qubit gates,” Physical Review A 69, 032315 (2004).
  * Kendrick  _et al._ (2012) John Kendrick, Frank JJ Leusen, and Marcus A Neumann, “Empirical van der Waals corrections to solid-state density functional theory: Iodine and phosphorous containing molecular crystals,” Journal of computational chemistry 33, 1615–1622 (2012).
  * Mitarai  _et al._ (2018) Kosuke Mitarai, Makoto Negoro, Masahiro Kitagawa, and Keisuke Fujii, “Quantum circuit learning,” Physical Review A 98, 032309 (2018).
  * Kingma and Ba (2014) Diederik P Kingma and Jimmy Ba, “Adam: A method for stochastic optimization,” arXiv preprint arXiv:1412.6980 (2014).
  * Jones (2010) Andrew Jones,  _Quantum Drude oscillators for accurate many-body intermolecular forces_ , Ph.D. thesis, University of Edinburgh (2010).
  * Boyd (2020) Robert W Boyd,  _Nonlinear optics_ (Academic press, 2020).
  * Sawaya and Huh (2019) Nicolas PD Sawaya and Joonsuk Huh, “Quantum algorithm for calculating molecular vibronic spectra,” The journal of physical chemistry letters 10, 3586–3591 (2019).
  * Fingerhuth  _et al._ (2018) Mark Fingerhuth, Tomáš Babej, _et al._ , “A quantum alternating operator ansatz with hard and soft constraints for lattice protein folding,” arXiv preprint arXiv:1810.13411 (2018).
  * Robert  _et al._ (2021) Anton Robert, Panagiotis Kl Barkoutsos, Stefan Woerner, and Ivano Tavernelli, “Resource-efficient quantum algorithm for protein folding,” npj Quantum Information 7, 1–5 (2021).
  * (75) https://doi.org/10.5281/zenodo.22558.
  * London (1942) F London, “On centers of van der Waals attraction.” The Journal of Physical Chemistry 46, 305–316 (1942).
  * Hirschfelder  _et al._ (1964) Joseph O Hirschfelder, Charles F Curtiss, Robert Byron Bird, and Maria Goeppert Mayer,  _Molecular theory of gases and liquids_ , Vol. 165 (Wiley New York, 1964).
  * Voloshin (2009) Vitaly I Voloshin,  _Introduction to graph and hypergraph theory_ (Nova Science Publ., 2009).
  * Baranyai (1974) Zsolt Baranyai, “On the factorization of the complete uniform hypergraphs,” Infinite and finite sets (1974).
  * Van Lint  _et al._ (2001) Jacobus Hendricus Van Lint, Richard Michael Wilson, and Richard Michael Wilson,  _A course in combinatorics_ (Cambridge university press, 2001).
  * Wecker  _et al._ (2015) Dave Wecker, Matthew B Hastings, and Matthias Troyer, “Progress towards practical quantum variational algorithms,” Physical Review A 92, 042303 (2015).
  * Anis  _et al._ (2021) Md Sajid Anis _et al._ , “Qiskit: An open-source framework for quantum computing,”  (2021).