# The variational quantum eigensolver self-consistent field method within a
polarizable embedded framework

Erik Rosendahl Kjellgren<EMAIL_ADDRESS>Department of Physics, Chemistry
and Pharmacy, University of Southern Denmark, DK-5230 Odense M, Denmark  Peter
Reinholdt  Department of Physics, Chemistry and Pharmacy, University of
Southern Denmark, DK-5230 Odense M, Denmark  Aaron Fitzpatrick  Algorithmiq
Ltd, Kanavakatu 3C, FI-00160 Helsinki, Finland  Walter N. Talarico
Algorithmiq Ltd, Kanavakatu 3C, FI-00160 Helsinki, Finland  Department of
Applied Physics, QTF Centre of Excellence, Center for Quantum Engineering,
Aalto University School of Science, FIN-00076 Aalto, Finland  Phillip W. K.
Jensen  Department of Chemistry, University of Copenhagen, DK-2100 Copenhagen
Ø, Denmark  Stephan P. A. Sauer  Department of Chemistry, University of
Copenhagen, DK-2100 Copenhagen Ø, Denmark  Sonia Coriani  DTU Chemistry -
Department of Chemistry, Technical University of Denmark, DK-2800 Kongens
Lyngby, Denmark  Stefan Knecht  Algorithmiq Ltd, Kanavakatu 3C, FI-00160
Helsinki, Finland  ETH Zürich, Department of Chemistry and Applied Life
Sciences, Vladimir-Prelog-Weg 1-5/10, CH-8093 Zürich, Switzerland  Jacob
Kongsted  Department of Physics, Chemistry and Pharmacy, University of
Southern Denmark, DK-5230 Odense M, Denmark

###### Abstract

We formulate and implement the Variational Quantum Eigensolver Self Consistent
Field (VQE-SCF) algorithm in combination with polarizable embedding (PE),
thereby extending PE to the regime of quantum computing. We test the resulting
algorithm, PE-VQE-SCF, on quantum simulators and demonstrate that the
computational stress on the quantum device is only slightly increased in terms
of gate counts compared to regular VQE-SCF. On the other hand, no increase in
shot noise was observed. We illustrate how PE-VQE-SCF may lead to the modeling
of real chemical systems using a simulation of the reaction barrier of the
Diels-Alder reaction between furan and ethene as an example.

##  I Introduction

Quantum computers are expected to enable the solution of many computational
problems that are intractable on classical high-performance computers. The
probably most often mentioned example is quantum chemistry, i.e., the solution
of the electronic Schrödinger equation. Hoefler  _et al._ (2023) Quantum
chemists have always been very receptive to advances in computer technology to
improve the performance of quantum chemistry programs, devising theories and
algorithms that could take advantage of the new hardware. Gordon  _et al._
(2020) Not surprisingly, since the first (proof-of-principle) application of a
quantum computing algorithm for computational chemistry problems in 2005,
Aspuru-Guzik  _et al._ (2005) the interest of the quantum chemistry community
in the development of quantum chemistry methods on quantum processing units
has been steadily growing, and almost exploded during the last five years in
concomitance with the advances on the hardware site. Setting out from
different angles and perspectives, a number of recent comprehensive reviews
discuss the current state of affairs concerning the development of theories,
methods, and algorithms that aim at performing quantum chemical calculations
on both near-term noisy intermediate-scale quantum (NISQ) devices and
early/far-future, fault-tolerant ones. Cao  _et al._ (2019); Bauer  _et al._
(2020); Elfving  _et al._ (2020); McArdle  _et al._ (2020); Motta and Rice
(2022); Aulicino  _et al._ (2022); Bharti  _et al._ (2022); Liu  _et al._
(2022) As of today, potential quantum advantage has only been demonstrated for
model systems, though. Liu  _et al._ (2022)

The current paradigm for calculations on NISQ devices is to employ a hybrid
quantum-classical approach, McClean  _et al._ (2016) where only parts of the
quantum chemical calculation are actually carried out on the quantum computer,
while the remaining parts are run on classical CPU- or GPU-based computers.
The calculations carried out on the QPUs mostly involve evaluating expectation
values over a Hamiltonian or other quantum mechanical operators. Often a
Variational Quantum Eigensolver (VQE) algorithm Yung  _et al._ (2014); Cao
_et al._ (2019); Romero  _et al._ (2019); Wang  _et al._ (2019); Fedorov  _et
al._ (2022) is employed in this step.

Within the realm of quantum chemical methods for NISQ devices, the majority of
the published work concerning wavefunction ansätze focuses on some variants of
unitary coupled cluster (UCC) theory. McClean  _et al._ (2016); Romero  _et
al._ (2019); Fedorov  _et al._ (2022) It is believed Romero  _et al._ (2019)
that UCC methods intrinsically are, or can be made, more easily multi-
configurational than the traditional coupled cluster ones. Helgaker  _et al._
(2013); Christiansen (2006); Bartlett (2012) This is a very important aspect
as several interesting chemical problems, especially within life sciences,
pertain to large and complicated multiconfigurational systems, the properties
and chemistry of which stand little chance of being accurately described by
classical computers. Elfving  _et al._ (2020) Approaches like complete-active
space self-consistent field (CASSCF), Siegbahn  _et al._ (1980); Roos  _et
al._ (1980); Siegbahn  _et al._ (1981) complete-active space configuration
interaction (CASCI) Olsen  _et al._ (1988) or density matrix renormalization
group (DMRG), Baiardi and Reiher (2020) which all involve a full configuration
interaction type wavefunction within the space of the chemical relevant
orbitals, are typically used to treat molecular systems with
multiconfigurational character on classical computers, but quickly hit the
ceiling of what can be handled in terms of size of the relevant active space.
Therefore, exploring the possibility of performing CASSCF/CI calculations
using quantum processors is natural. Parrish  _et al._ (2019); Bauman  _et
al._ (2021); Fitzpatrick  _et al._ (2022)

However, most of chemistry and all of biochemistry happen in solution.
Furthermore, the current most accurate approach to model the reactions
catalyzed by enzymes is to treat the active site with some quantum chemical
method and describe the protein environment using some molecular mechanics
model, i.e., the hybrid QM/MM approach. For quantum chemical simulations on
classical computers, many methods for treating solvent or environment effects
have been developed and implemented, ranging from simple continuum solvent
models like the polarizable continuum model (PCM) Mennucci (2012) to an
explicit treatment of solvent molecules or the fragments of molecules in the
environment via the Polarizable Embedding (PE) model. Olsen  _et al._ (2010);
Olsen and Kongsted (2011) For calculations on quantum computers, on the other
hand, only simple approaches such as the PCM method Castaldo  _et al._ (2022)
and point-charge embeddingsHohenstein  _et al._ (2023) have so far been
implemented in combination with VQE calculations of ground state energies.
More sophisticated quantum embedding schemesRossmannek  _et al._ (2023) could
potentially be applied to include solvation effects, although a practical
demonstration has yet to be reported. In this work, we present, therefore, an
implementation of the PE model for quantum computers in combination with the
Adaptive Derivative-Assembled Problem-Tailored Ansatz Variational Quantum
Eigensolver self-consistent field approach (ADAPT-VQE-SCF), Fitzpatrick  _et
al._ (2022) and illustrate its performance for ammonia in water and for the
Diels-Alder reaction of furan with ethene in water.

##  II Theoretical summary

In the following, we outline the theoretical background for our novel
developments. This theoretical summary is split into three parts: first, we
introduce the polarizable embedding model; second, the working equations from
the variational quantum eigensolver self-consistent field method are briefly
discussed; and third, we show how to integrate the PE model with the ADAPT-
VQE-SCF scheme.

###  II.1 Polarizable Embedding

The Polarizable EmbeddingOlsen  _et al._ (2010); Olsen and Kongsted (2011)
(PE) model divides a total molecular system into a smaller quantum region and
an environment. The quantum region is described using an accurate quantum
chemistry model, whereas the environment is represented semi-classically
through multipoles and polarizabilities. The parameters describing the
embedding potential are generated by dividing the environment into small
fragments and carrying out individual quantum-chemical calculations for each
of them to obtain atomic-based multipoles and polarizabilities using a multi-
center multipole expansion. The multipoles and polarizabilities represent the
permanent and induced charge density of the fragments in the environment,
respectively. The partitioning of the system in QM and polarizable MM regions
is graphically illustrated in Figure 1.

Figure 1: The QM/MM partitioning in the PE model. A complete system is divided
into QM and MM subsystems, where the QM part is described through an accurate
quantum chemical method. The environment is divided into smaller fragments,
described by an atomic multi-center multipole expansion and atomic
polarizabilities.

The effects of the environment on the quantum region are included by
constructing an effective Hamiltonian

$\hat{H}^{\mathrm{eff}}=\hat{H}^{\mathrm{vac}}+\hat{v}^{\mathrm{PE}}$ (1)

where $\hat{H}^{\mathrm{vac}}$ represents the regular (isolated) Hamiltonian
for the quantum region without the presence of the environment. The operator
$\hat{v}^{\mathrm{PE}}$ is the PE operator, which describes the potential from
the environment. The PE operator is further divided into an electrostatic and
an induction (polarization) term

$\hat{v}^{\mathrm{PE}}=\hat{v}^{\mathrm{es}}+\hat{v}^{\mathrm{ind}}$ (2)

where the electrostatic operator, $\hat{v}^{\mathrm{es}}$, describes the
potential from the permanent charge distributions of the $N$ environment sites

$\hat{v}^{\mathrm{es}}=\sum_{s=1}^{N}\sum_{\left|k\right|=0}^{K}\frac{\left(-1\right)^{\left|k\right|}}{k!}M_{s}^{\left(k\right)}\hat{V}_{s,\mathrm{el}}^{\left(k\right)}$
(3)

written using a multi-index notation,Olsen (2012) i.e., $k$ runs over the
Cartesian components of the multipole moment tensors. Here,
$M_{s}^{\left(k\right)}$ are the multipoles on expansion site $s$, i.e., it is
the $k$’th components of the atomic charge, dipole and quadrupole moment of
site $s$. The operator $\hat{V}_{s,\mathrm{el}}^{\left(k\right)}$ gives a
derivative of the electric potential at site $s$, and is defined as

$\hat{V}_{s,\mathrm{el}}^{(k)}=\sum_{pq}t^{(k)}_{pq}(\mathbf{R}_{s})\hat{E}_{pq}$
(4)

where
$\hat{E}_{pq}=\hat{a}^{\dagger}_{p\alpha}\hat{a}_{q\alpha}+\hat{a}^{\dagger}_{p\beta}\hat{a}_{q\beta}$
is the singlet one-electron excitation operator,Helgaker  _et al._ (2013) and
the potential derivative integrals are

$t^{(k)}_{pq}(\mathbf{R}_{s})=\int\phi^{*}_{p}(\mathbf{r})\nabla^{k}\left(\frac{1}{\left|\mathbf{r}-\mathbf{R}_{s}\right|}\right)\phi_{q}(\mathbf{r})\mathrm{d}\mathbf{r}~{}.$
(5)

The induction operator, $\hat{v}^{\mathrm{ind}}$, originates from the
polarized charge distribution of the environment, described by induced
dipoles, and it is defined as

$\hat{v}^{\mathrm{ind}}=-\sum_{s=1}^{N}\boldsymbol{\mu}_{s}^{\mathrm{ind}}\left({\mathbf{F}}_{\mathrm{tot}}\right)\hat{V}^{(1)}_{s,\mathrm{el}}\
.$ (6)

The induced dipoles $\boldsymbol{\mu}_{s}^{\mathrm{ind}}$ are generated in
response to the total electric field ($\mathbf{F}_{\mathrm{tot}}$) at the
polarizable site $s$, which is the sum of fields from the electrons and nuclei
in the quantum region, the permanent multipole moments in the environment, as
well as other induced dipoles in the environment. The operator
$\hat{{V}}^{(1)}_{s,\mathrm{el}}$ gives the electric field from the electrons
at a polarizable site $s$. The induced dipoles entering the expression for the
induction operator are derived from

$\boldsymbol{\mathbf{\mu}}_{s}^{\mathrm{ind}}\left(\mathbf{F}_{\mathrm{tot}}\right)=\boldsymbol{\alpha}_{s}\mathbf{F}_{\mathrm{tot}}(\mathbf{R}_{s})=\boldsymbol{\alpha}_{s}\left(\mathbf{F}(\mathbf{R}_{s})+\sum_{s^{\prime}\neq
s}\mathbf{T}_{ss^{\prime}}^{(2)}\boldsymbol{\mathbf{\mu}}_{s^{\prime}}^{\mathrm{ind}}\right)\
,$ (7)

where $\mathbf{F}(\mathbf{R}_{s})$ is the electric field at site $s$ from the
nuclei, electrons, and permanent multipole moments, and
$\mathbf{T}_{ss^{\prime}}^{(2)}$ is a dipole–dipole interaction tensor.Stone
(2013) Equation (7) leads to a set of coupled equations which can be
formulated as a linear system of equations by introducing a column vector
containing the induced dipoles
$\boldsymbol{\mu}^{\mathrm{ind}}=\left(\boldsymbol{\mu}_{1}^{\mathrm{ind}},\boldsymbol{\mu}_{2}^{\mathrm{ind}},\ldots,\boldsymbol{\mu}_{N}^{\mathrm{ind}}\right)^{T}$,
and one containing the electric fields
$\mathbf{F}=\left(\mathbf{F}(\mathbf{R}_{1}),\mathbf{F}(\mathbf{R}_{2}),\ldots,\mathbf{F}(\mathbf{R}_{N})\right)^{T}$.
The induced dipole moments may formally be obtained as the solution of

$\boldsymbol{\mu}^{\mathrm{ind}}=\mathbf{BF}\ ,$ (8)

where $\mathbf{B}$ is the ($3N\times 3N$) classical response matrix

$\mathbf{B}=\begin{pmatrix}\boldsymbol{\alpha}_{1}^{-1}&-\mathbf{T}_{12}^{(2)}&\ldots&-\mathbf{T}_{1N}^{(2)}\\\
-\mathbf{T}_{21}^{(2)}&\boldsymbol{\alpha}_{2}^{-1}&\ddots&\vdots\\\
\vdots&\ddots&\ddots&-\mathbf{T}_{(N-1)N}^{(2)}\\\
-\mathbf{T}_{N1}^{(2)}&\dots&-\mathbf{T}_{N(N-1)}^{(2)}&\boldsymbol{\alpha}_{N}^{-1}\end{pmatrix}^{-1}\
.$ (9)

In practice, the induced dipoles can be obtained by either explicitly
constructing the response matrix and directly solving Equation (8) by
conventional means or by iterative schemes when an explicit construction is
not computationally feasible due to the potentially large dimensions of the
environment. The induced dipoles depend on the fields exerted by the electron
density of the quantum region, which in turn depends on the induced dipoles
through the induction operator, and it is thus required to use a self-
consistent scheme in which the induced dipoles are updated in every quantum
mechanically self-consistent-field cycle. Such a scheme finally leads to a
mutual relaxation of the quantum wavefunction/density and the induced dipoles
of the environment.

###  II.2 Variational Quantum Eigensolver Self-Consistent Field

The following account of the variational quantum eigensolver self-consistent
field (VQE-SCF) presents the most important expressions from Fitzpatrick et
al.Fitzpatrick  _et al._ (2022) The molecular spin-free full (non-
relativistic) electronic Hamiltonian, $\hat{H}$, is the starting point for the
VQE-SCF approach

$\hat{H}=\sum_{pq}h_{pq}\hat{E}_{pq}+\frac{1}{2}\sum_{pqrs}g_{pqrs}\left(\hat{E}_{pq}\hat{E}_{rs}-\delta_{qr}\hat{E}_{ps}\right)$
(10)

with the one- and two-electron integrals, respectively,

$h_{pq}=\int\phi^{*}_{p}\left(\mathbf{r}\right)\hat{h}\
\phi_{q}\left(\mathbf{r}\right)\mathrm{d}\mathbf{r}$ (11)
$g_{pqrs}=\int\phi^{*}_{p}\left(\mathbf{r}_{1}\right)\phi^{*}_{r}\left(\mathbf{r}_{2}\right)r^{-1}_{12}\phi_{q}\left(\mathbf{r}_{1}\right)\phi_{s}\left(\mathbf{r}_{2}\right)\mathrm{d}\mathbf{r}_{1}\mathrm{d}\mathbf{r}_{2}~{}.$
(12)

The set of functions $\phi_{p}$ represents an orthonormal molecular orbital
(MO) basis, i.e., $\left<\phi_{p}\left|\phi_{q}\right.\right>=\delta_{pq}$. In
the framework of VQE-SCF, the wavefunction is parameterized by orbital
rotation coefficients, $\kappa_{pq}$, and unitary qubit rotation angles,
$\theta_{i}$. The parameterized wavefunction takes the form,

$\left|\Psi\left(\boldsymbol{\kappa},\boldsymbol{\theta}\right)\right>=U\left(\boldsymbol{\kappa}\right)U\left(\boldsymbol{\theta}\right)\left|0\right>$
(13)

with $\left|0\right>$ being the reference state. The orbital rotation
parameterization takes the same form as in the conventional multi-
configurational self-consistent field (MCSCF) methods, Helgaker  _et al._
(2013)

$U\left(\boldsymbol{\kappa}\right)=\exp\left({\hat{\kappa}}\right)$ (14)

with the anti-hermitian orbital rotation operator being

${\hat{\kappa}}=\sum_{p>q}\kappa_{pq}\left(\hat{E}_{pq}-\hat{E}_{qp}\right)~{}.$
(15)

The parameterization of the unitary qubit rotation angles is where the VQE-SCF
differs from the conventional MCSCF. In conventional MCSCF, the wavefunction
is usually parameterized with configuration interaction (CI) linear-
combination coefficients or, in unitary form, by the state-transfer operator.
Instead, to target a quantum computer, the parameterization uses unitary qubit
rotation angles. This parameterization takes the form

$U\left(\boldsymbol{\theta}\right)=\prod_{i}\exp\left(i\theta_{i}\hat{P}_{i}\right)$
(16)

with $\hat{P}_{i}$ being strings of Pauli operators. The ground state energy
is determined by minimization with respect to both the orbital rotation
coefficients and the unitary qubit rotation angles,

$E=\min_{\boldsymbol{\theta},\boldsymbol{\kappa}}\left<0\left|U^{\dagger}\left(\boldsymbol{\theta}\right)U^{\dagger}\left(\boldsymbol{\kappa}\right)\hat{H}U\left(\boldsymbol{\kappa}\right)U\left(\boldsymbol{\theta}\right)\right|0\right>~{}.$
(17)

The choice of specific unitary qubit rotation parameters in the ansatz follows
the ADAPT approach and can be either ADAPT-VQE Grimsley  _et al._ (2019) or
qubit-ADAPT-VQE.Tang  _et al._ (2021)

###  II.3 PE-VQE-SCF

The following briefly outlines how to integrate the VQE approach with PE. In
the PE-VQE-SCF model, instead of minimizing the energy with respect to the
molecular Hamiltonian, as in Eq. (10), we minimize it with respect to the
free-energy in solution

$E_{\text{PE}}=\min_{\boldsymbol{\theta},\boldsymbol{\kappa}}\left<0\left|U^{\dagger}\left(\boldsymbol{\theta}\right)U^{\dagger}\left(\boldsymbol{\kappa}\right)\hat{H}^{\text{PE}}U\left(\boldsymbol{\kappa}\right)U\left(\boldsymbol{\theta}\right)\right|0\right>$
(18)

where the Hamiltonian is defined as

$\hat{H}^{\text{PE}}=\hat{H}^{\mathrm{eff}}-\frac{1}{2}\hat{v}^{\mathrm{ind}}.$
(19)

For the VQE-SCF algorithm, the effective Hamiltonian needs to be decomposed
into Pauli strings. The underlying fermionic operators in the operators
$\hat{V}_{s,\mathrm{el}}^{\left(k\right)}$ and
$\hat{\mathbf{F}}_{s,\mathrm{el}}$, Eq. (3) and Eq. (6), respectively, are the
singlet single excitation operators $\hat{E}_{pq}$. In the gas-phase, this
corresponds to the following decomposition of the Hamiltonian

$\hat{H}^{\text{vac}}(\boldsymbol{\kappa})=\sum_{i}c_{i}(\boldsymbol{\kappa})\hat{P}_{i}~{},$
(20)

where $c_{i}(\kappa)$ is the weight of the Pauli strings from the
decomposition and $\hat{P}_{i}$ is a unique Pauli string. After introducing
the PE environment, the decomposition of the Hamiltonian becomes

$\hat{H}^{\text{PE}}\left(\boldsymbol{\kappa},\mathbf{F}_{\text{tot}}\right)=\sum_{i}\tilde{c}_{i}(\boldsymbol{\kappa},\mathbf{F}_{\text{tot}})\hat{P}_{i},$
(21)

which is similar to eq. (20), but with modified expansion coefficients
$\tilde{c}_{i}$. Crucially, these weights depend on electric fields from the
polarizable environment instead of just the orbital rotation parameters.

Due to the induced dipoles in the environment, the PE Hamiltonian is non-
linear in the wavefunction parametrization. This leads to a double-SCF
algorithm where the wavefunction parameters and the induced dipoles are
iteratively minimized, thereby leading to a computational overhead within PE
compared to calculating the corresponding molecule in a vacuum.

Figure 2: Flow chart of the PE-VQE-SCF algorithm. In this flow chart $D$ is
the one-electron reduced density matrix.

The algorithmic approach underlying the PE-VQE-SCF methods is illustrated in
Figure 2. In this figure, the overall structure of the PE-VQ-SCF algorithm can
be seen. The blue box shows the SCF loop that holds together the PE model and
the VQE-SCF procedure. The yellow box shows what the PE library handles, with
input and output. Finally, the green box illustrates the VQE-SCF algorithm. It
can be seen that the PE library provides the PE contribution to the
Hamiltonian, which is then used to solve the VQE-SCF problem. The VQE-SCF
procedure provides a density matrix to the PE library that is needed to solve
for the PE contribution (the induced dipoles). This cycle is continued until
convergence is achieved.

##  III Computational Details

The VQE-SCF calculations were run using AURORA,Algorithmiq Ltd. (2023) with
PySCFSun  _et al._ (2017, 2020) as the backend for the orbital optimization
and CPPEScheurer  _et al._ (2019) as the induced dipole moment solver for the
PE model. For the ammonia systems (as described below), ADAPT-VQE was used as
the ansatz for the VQE-SCF calculations. The VQE-SCF calculations were
optimized according to the algorithm described in Fitzpatrick et
al.Fitzpatrick  _et al._ (2022) Furthermore, to simulate the probabilistic
nature of the measurement outcome of a quantum computer, we used Adaptive
Informationally complete generalized Measurements (AIM) implemented as
positive operator valued measurements (IC-POVMs)García-Pérez  _et al._ (2021);
Glos  _et al._ (2022) from the AURORA package. As we use informationally
complete POVM for the estimation of physical observables $\bar{O}$ with a
finite set of measurement outcomes (very much like Monte Carlo sampling), the
estimation of these observables comes with a statistical error originating
from the finite number of samples $S$ considered and has the expected scaling
$\sigma_{\bar{O}}\sim\frac{1}{\sqrt{S}}$. This scheme not only allows us to
simulate the shot noise of the measurements but also to convert the AIM
outcomes into estimates of different observablesFitzpatrick  _et al._ (2022);
Nykänen  _et al._ (2022); Filippov  _et al._ (2023) needed for the proposed
PE-VQE-SCF method.

All calculations performed on the ammonia systems utilized the aug-cc-
pVTZDunning (1989); Kendall  _et al._ (1992) basis set for the quantum region
(ammonia), and the loprop-aug-cc-pVTZGagliardi  _et al._ (2004) basis set for
the calculation of the atomic multipole moments and atomic polarizabilities
within the PE region. The active space comprised 6 electrons in 6 spatial
orbitals, and was picked on the basis of MP2 natural orbital occupation
numbers.

To illustrate the PE-VQE-SCF approach further, the reaction profile for the
Diels-Alder reaction between furan and ethene (reactant) to form
(1R,4S)-7-oxabicyclo[2.2.1]hept-2-ene (product) was simulated using the nudged
elastic band (NEB) method.Henkelman  _et al._ (2000); Goumans  _et al._ (2009)
The reactant and product structures were placed in a sphere of 40 water
molecules using the PackmolMartínez  _et al._ (2009) program. To ensure that
the two structures were as similar as possible for the following NEB
calculation, the geometry of the two structures was minimized using
CAM-B3LYPYanai  _et al._ (2004)-D3BJGrimme  _et al._ (2010,
2011)/6-31G*Ditchfield  _et al._ (1971); Hehre  _et al._ (1972); Hariharan and
Pople (1973) using TerachemUfimtsev and Martinez (2009); Titov  _et al._
(2013); Seritan  _et al._ (2021) (version 1.95A) for energies and gradients
with shared solvent coordinates, i.e., we minimize the sum of the energies of
the solvated reactant and product structures with respect to the coordinates
of the reactant, the coordinates of the product, and the coordinates of the
water (shared between the solvated reactant/product structures). The
minimization was done in cartesian coordinates with a custom Python script
using SLSQPKraft (1988) from SciPy.Virtanen  _et al._ (2020) The NEB path was
generated from these starting structures with Terachem with
CAM-B3LYP-D3BJ/6-31G*. We used a ”free-end” NEB with 25 NEB images and an NEB
spring constant of 0.01 a.u.. The (PE-)VQE-SCF calculations of the reaction
profile were performed using the 6-31G* basis set and an active space
consisting of 6 electrons in 6 spatial orbitals, and was picked on the basis
of MP2 natural orbital occupation numbers. The ADAPT approach was used in all
molecular examples to find a suitable parametrized quantum circuit describing
the multiconfigurational wavefunction within the considered active spaces.

##  IV Results

The result section is structured into two parts. In part A, we aim to test the
implemented polarizable embedding strategy in relation to regular VQE-SCF. Our
focus in this respect will be on important parameters regarding the efficiency
of quantum devices, e.g., gate count and simulated shot noise. For this
investigation, we will use a chemical model system consisting of ammonia
solvated in water, and we will monitor the above-mentioned parameters as a
function of the number of water molecules included in the calculation. All
water molecules will be described using the polarizable embedding potential,
whereas ammonia will be considered at the QM level of theory. In part B, we
illustrate the PE-VQE-SCF approach for the simulation of a reaction coordinate
profile of a Diels-Alder reaction, taking directly into account the effects of
solvation. We include all the solvating water molecules in the polarizable
embedding description. We note that we use a limited number of snapshots in
the PE calculations in these examples. In general, statistical sampling of the
solvent configurations, for example, sampled from a molecular dynamics
trajectory, is needed for realistic computational applications.

###  IV.1 Ammonia

In this section, we will consider ammonia using PE-VQE-SCF as a test system to
access key performance parameters related to the VQE-SCF algorithm.

Figure 3: Relative CNOT count with respect to ammonia in isolation (no
polarizable embedding) as a function of the number of water molecules in the
polarizable embedding potential. The different colors represent different
snapshots of water configurations.

In Figure 3, the number of required CNOT gates as a function of the size of
the PE region (increasing number of water molecules) relative to the
corresponding energy calculation of ammonia in isolation is shown. Five
different snapshots have been considered in order to avoid any bias towards
just a single solvent configuration. From Figure 3, we observe that adding the
PE model to the system, on average, will slightly increase the gate count of
the calculation compared to the corresponding calculation of ammonia in
isolation. On average, the gate count is increased by 16% when using the PE
model for this system. Using ADAPT-VQE, it is to be expected that a larger
gate count will appear when including the environment in the system since the
wave function will be more complicated (e.g., lack of point-group symmetry)
than the wave function describing the system in isolation. However, as further
observed from Figure 3, we do not observe a systematic increase in the CNOT
count when increasing the size of the embedding part of the system, and the
increase in the number of water molecules does not lead, on this basis, to a
systematic upward trend regarding the complexity of the underlying circuit.

Figure 4: Mean energy and shot noise indicated as plus/minus the standard
error. Dashed curves represent the energy, while the full lines represent the
error interval. Orange color refers to VQE-SCF while blue color represents PE-
VQE-SCF. In total 46 water molecules were included in the PE region.

In Figure 4, we inspect the mean energy and the simulated shot noise as a
function of the number of sampling shots for either VQE-SCF or PE-VQE-SCF. It
is immediately clear that the energy of the PE-VQE-SCF calculations is lower
than the corresponding energy of the VQE-SCF calculation. This is expected
since the PE model will add terms to the Hamiltonian that will decrease the
total energy of the system, i.e., the solvating water will stabilize the
ammonia molecule. As expected, the sampling error reduces with an increasing
number of shots for either VQE-SCF or PE-VQE-SCF, and the solutions converge
towards the classically computed CASSCF and PE-CASSCF results. It is further
observed that the noise decreases roughly at the same rate for VQE-SCF and PE-
VQE-SCF upon increasing the number of shots. On this basis, it is concluded
that adding the embedding term to the molecular Hamiltonian does not
negatively affect how many shots are needed to obtain reliable energies within
the PE-VQE-SCF framework.

###  IV.2 Reaction curve

Next, we consider in this section the ability of the PE-VQE-SCF approach to
model real applications showcased based on the Diels-Alder reaction between
furan and ethene.

Figure 5: Calculated reaction profile of the Diels-Alder reaction between
furan and ethene based on either VQE-SCF (blue line) or PE-VQE-SCF (orange
line).

In Figure 5, the reaction profile of the furan-ethene Diels-Alder reaction is
calculated based on either VQE-SCF or PE-VQE-SCF and is shown with respect to
the reactant energies obtained in either VQE-SCF or PE-VQE-SCF, i.e.,
$E-E_{0}$. In these calculations, furan and ethene are included in the QM part
of the system, while all water molecules are described using the embedding
potential. From Figure 5, we observe a decrease in the calculated barrier upon
inclusion of the water solvent, and further, we observe a stabilization of the
products. Both these observations, i.e., reduction of the reaction barrier and
stabilization of the products, are caused mainly by dipolar stabilization of
the chemical reactive part of the system due to the presence of the polar
(water) environment. The specific reduction of the reaction barrier is around
0.015 au, corresponding to 9.4 kcal/mol, and the presence of the water
environment thereby represents an important part of the chemical system.

##  V Conclusion

In this work, we have formulated and implemented the VQE-SCF algorithm in
combination with polarizable embedding, thereby enabling hybrid quantum-
classical simulations on quantum devices. The resulting algorithm, PE-VQE-SCF,
has been tested on simulators, where we show that the computational stress on
the quantum device is not significantly increased, neither in terms of gate
count nor in terms of additional shot noise, upon extending VQE-SCF to PE-VQE-
SCF. We have further showcased our PE-VQE-SCF implementation by calculating
the Diels-Alder reaction barrier between furan and ethene. We conclude that
PE-VQE-SCF represents a robust computational quantum computing strategy for
modeling larger and extended chemical systems. Further work will address the
ability to use PE-VQE-SCF to tackle biological challenges and extend the
computational framework to excited states and molecular property calculations.

###### Acknowledgements.

We acknowledge the support of the Novo Nordisk Foundation (NNF) for the
focused research project “Hybrid Quantum Chemistry on Hybrid Quantum
Computers” (grant number NNFSA220080996).

## DATA AVAILABILITY

The data that support the findings of this study are openly available in
Zenodo,at https://zenodo.org/records/10279684, reference number
10.5281/zenodo.10279684.

## References

  * Hoefler  _et al._ (2023) T. Hoefler, T. Häner, and M. Troyer, Communication of the ACM 66, 83 (2023).
  * Gordon  _et al._ (2020) M. S. Gordon, G. Barca, S. S. Leang, D. Poole, A. P. Rendell, J. L. Galvez Vallejo, and B. Westheimer, J. Phys. Chem. A 124, 4557 (2020).
  * Aspuru-Guzik  _et al._ (2005) A. Aspuru-Guzik, A. D. Dutoi, P. J. Love, and M. Head-Gordon, Science 309, 1704 (2005).
  * Cao  _et al._ (2019) Y. Cao, J. Romero, J. P. Olson, M. Degroote, P. D. Johnson, M. Kieferová, I. D. Kivlichan, T. Menke, B. Peropadre, N. P. D. Sawaya, S. Sim, L. Veis, and A. Aspuru-Guzik, Chem. Rev. 119, 10856 (2019).
  * Bauer  _et al._ (2020) B. Bauer, S. Bravyi, M. Motta, and G. Kin-Lic Chan, Chem. Rev. 120, 12685 (2020).
  * Elfving  _et al._ (2020) V. E. Elfving, B. W. Broer, M. Webber, J. Gavartin, M. D. Halls, K. P. Lorton, and A. Bochevarov, arxiv , 1 (2020), arXiv:2009.12472 .
  * McArdle  _et al._ (2020) S. McArdle, S. Endo, A. Aspuru-Guzik, S. C. Benjamin, and X. Yuan, Rev. Mod. Phys. 92, 015003 (2020).
  * Motta and Rice (2022) M. Motta and J. E. Rice, WIREs Comput. Mol. Sci. 12, e1580 (2022).
  * Aulicino  _et al._ (2022) J. C. Aulicino, T. Keen, and B. Peng, Int. J. Quant. Chem. 122 (2022), 10.1002/qua.26853.
  * Bharti  _et al._ (2022) K. Bharti, A. Cervera-Lierta, T. H. Kyaw, T. Haug, S. Alperin-Lea, A. Anand, M. Degroote, H. Heimonen, J. S. Kottmann, T. Menke, W.-K. Mok, S. Sim, L.-C. Kwek, and A. Aspuru-Guzik, Rev. Mod. Phys. 94, 015004 (2022).
  * Liu  _et al._ (2022) H. Liu, G. H. Low, D. S. Steiger, T. Häner, M. Reiher, and M. Troyer, Mater. theory 6 (2022), 10.1186/s41313-021-00039-z, arXiv:2102.10081 .
  * McClean  _et al._ (2016) J. R. McClean, J. Romero, R. Babbush, and A. Aspuru-Guzik, New J. Phys. 18, 023023 (2016).
  * Yung  _et al._ (2014) M.-H. Yung, J. Casanova, A. Mezzacapo, J. McClean, L. Lamata, A. Aspuru-Guzik, and E. Solano, Sci. Rep. 4, 3589 (2014).
  * Romero  _et al._ (2019) J. Romero, R. Babbush, J. R. McClean, C. Hempel, P. J. Love, and A. Aspuru-Guzik, Quantum Sci. Technol. 4, 014008 (2019).
  * Wang  _et al._ (2019) D. Wang, O. Higgott, and S. Brierley, Phys. Rev. Lett. 122, 140504 (2019).
  * Fedorov  _et al._ (2022) D. A. Fedorov, B. Peng, N. Govind, and Y. Alexeev, Mater. theory 6 (2022), 10.1186/s41313-021-00032-6, arXiv:2103.08505 .
  * Helgaker  _et al._ (2013) T. Helgaker, P. Jørgensen, and J. Olsen,  _Molecular electronic-structure theory_ (John Wiley & Sons, Nashville, TN, 2013).
  * Christiansen (2006) O. Christiansen, Theor. Chem. Acc. 116, 106 (2006).
  * Bartlett (2012) R. J. Bartlett, WIREs Comput. Mol. Sci. 2, 126 (2012).
  * Siegbahn  _et al._ (1980) P. Siegbahn, A. Heiberg, B. Roos, and B. Levy, Phys. Scr. 21, 323 (1980).
  * Roos  _et al._ (1980) B. O. Roos, P. R. Taylor, and P. E. Sigbahn, Chem. Phys. 48, 157 (1980).
  * Siegbahn  _et al._ (1981) P. E. M. Siegbahn, J. Almlöf, A. Heiberg, and B. O. Roos, J. Chem. Phys. 74, 2384 (1981).
  * Olsen  _et al._ (1988) J. Olsen, B. O. Roos, P. Jo/rgensen, and H. J. A. Jensen, J. Chem. Phys. 89, 2185–2192 (1988).
  * Baiardi and Reiher (2020) A. Baiardi and M. Reiher, J. Chem. Phys. 152, 040903 (2020).
  * Parrish  _et al._ (2019) R. M. Parrish, E. G. Hohenstein, P. L. McMahon, and T. J. Martínez, Phys. Rev. Lett. 122, 230401 (2019).
  * Bauman  _et al._ (2021) N. P. Bauman, H. Liu, E. J. Bylaska, S. Krishnamoorthy, G. H. Low, C. E. Granade, N. Wiebe, N. A. Baker, B. Peng, M. Roetteler, M. Troyer, and K. Kowalski, J. Chem. Theory Comput. 17, 201 (2021).
  * Fitzpatrick  _et al._ (2022) A. Fitzpatrick, A. Nykänen, N. W. Talarico, A. Lunghi, S. Maniscalco, G. García-Pérez, and S. Knecht, arxiv (2022), arXiv:2212.11405 [quant-ph] .
  * Mennucci (2012) B. Mennucci, WIREs Comput. Mol. Sci. 2, 386 (2012).
  * Olsen  _et al._ (2010) J. M. H. Olsen, K. Aidas, and J. Kongsted, J. Chem. Theory Comput. 6, 3721 (2010).
  * Olsen and Kongsted (2011) J. M. H. Olsen and J. Kongsted, Adv. Quantum Chem. 61, 107 (2011).
  * Castaldo  _et al._ (2022) D. Castaldo, S. Jahangiri, A. Delgado, and S. Corni, J. Chem. Theory Comput. 18, 7457 (2022).
  * Hohenstein  _et al._ (2023) E. G. Hohenstein, O. Oumarou, R. Al-Saadon, G.-L. R. Anselmetti, M. Scheurer, C. Gogolin, and R. M. Parrish, J. Chem. Phys. 158 (2023).
  * Rossmannek  _et al._ (2023) M. Rossmannek, F. Pavosevic, A. Rubio, and I. Tavernelli, J. Phys. Chem. Lett. 14, 3491 (2023).
  * Olsen (2012) J. M. H. Olsen, “Development of Quantum Chemical Methods Towards Rationalization and Optimal Design of Optically Active Proteins. Ph.D. thesis,” (2012), University of Southern Denmark, Faculty of Science, https://doi.org/10.6084/m9.figshare.156852.
  * Stone (2013) A. Stone,  _The Theory of Intermolecular Forces_ , 2nd ed. (Oxford University Press, 2013) iSBN: 0199672393.
  * Grimsley  _et al._ (2019) H. R. Grimsley, S. E. Economou, E. Barnes, and N. J. Mayhall, Nat. Commun. 10, 3007 (2019).
  * Tang  _et al._ (2021) H. L. Tang, V. Shkolnikov, G. S. Barron, H. R. Grimsley, N. J. Mayhall, E. Barnes, and S. E. Economou, PRX Quantum 2, 020310 (2021).
  * Algorithmiq Ltd. (2023) Algorithmiq Ltd., “Aurora,” (2023).
  * Sun  _et al._ (2017) Q. Sun, T. C. Berkelbach, N. S. Blunt, G. H. Booth, S. Guo, Z. Li, J. Liu, J. D. McClain, E. R. Sayfutyarova, S. Sharma, S. Wouters, and G. K.-L. Chan, WIREs Comput. Mol. Sci. 8 (2017), 10.1002/wcms.1340.
  * Sun  _et al._ (2020) Q. Sun, X. Zhang, S. Banerjee, P. Bao, M. Barbry, N. S. Blunt, N. A. Bogdanov, G. H. Booth, J. Chen, Z.-H. Cui, J. J. Eriksen, Y. Gao, S. Guo, J. Hermann, M. R. Hermes, K. Koh, P. Koval, S. Lehtola, Z. Li, J. Liu, N. Mardirossian, J. D. McClain, M. Motta, B. Mussard, H. Q. Pham, A. Pulkin, W. Purwanto, P. J. Robinson, E. Ronca, E. R. Sayfutyarova, M. Scheurer, H. F. Schurkus, J. E. T. Smith, C. Sun, S.-N. Sun, S. Upadhyay, L. K. Wagner, X. Wang, A. White, J. D. Whitfield, M. J. Williamson, S. Wouters, J. Yang, J. M. Yu, T. Zhu, T. C. Berkelbach, S. Sharma, A. Y. Sokolov, and G. K.-L. Chan, J. Chem. Phys. 153, 024109 (2020).
  * Scheurer  _et al._ (2019) M. Scheurer, P. Reinholdt, E. R. Kjellgren, J. M. Haugaard Olsen, A. Dreuw, and J. Kongsted, J. Chem. Theory Comput. 15, 6154 (2019).
  * García-Pérez  _et al._ (2021) G. García-Pérez, M. A. Rossi, B. Sokolov, F. Tacchino, P. K. Barkoutsos, G. Mazzola, I. Tavernelli, and S. Maniscalco, PRX Quantum 2, 040342 (2021).
  * Glos  _et al._ (2022) A. Glos, A. Nykänen, E.-M. Borrelli, S. Maniscalco, M. A. Rossi, Z. Zimborás, and G. García-Pérez, “Adaptive povm implementations and measurement error mitigation strategies for near-term quantum devices,”  (2022).
  * Nykänen  _et al._ (2022) A. Nykänen, M. A. C. Rossi, E.-M. Borrelli, S. Maniscalco, and G. García-Pérez, “Mitigating the measurement overhead of ADAPT-VQE with optimised informationally complete generalised measurements,”  (2022), arXiv: 2212.09719.
  * Filippov  _et al._ (2023) S. Filippov, M. Leahy, M. A. C. Rossi, and G. García-Pérez, “Scalable tensor-network error mitigation for near-term quantum computing,” (2023), arXiv:2307.11740 [quant-ph] .
  * Dunning (1989) T. H. Dunning, J. Chem. Phys. 90, 1007 (1989).
  * Kendall  _et al._ (1992) R. A. Kendall, T. H. Dunning, and R. J. Harrison, J. Chem. Phys. 96, 6796 (1992).
  * Gagliardi  _et al._ (2004) L. Gagliardi, R. Lindh, and G. Karlström, J. Chem. Phys. 121, 4494 (2004).
  * Henkelman  _et al._ (2000) G. Henkelman, B. P. Uberuaga, and H. Jónsson, J. Chem. Phys. 113, 9901 (2000).
  * Goumans  _et al._ (2009) T. Goumans, C. R. A. Catlow, W. A. Brown, J. Kästner, and P. Sherwood, Phys. Chem. Chem. Phys. 11, 5431 (2009).
  * Martínez  _et al._ (2009) L. Martínez, R. Andrade, E. G. Birgin, and J. M. Martínez, J. Comput. Chem. 30, 2157 (2009).
  * Yanai  _et al._ (2004) T. Yanai, D. P. Tew, and N. C. Handy, Chem. Phys. Lett. 393, 51 (2004).
  * Grimme  _et al._ (2010) S. Grimme, J. Antony, S. Ehrlich, and H. Krieg, J. Chem. Phys. 132, 154104 (2010).
  * Grimme  _et al._ (2011) S. Grimme, S. Ehrlich, and L. Goerigk, J. Comput. Chem. 32, 1456 (2011).
  * Ditchfield  _et al._ (1971) R. Ditchfield, W. J. Hehre, and J. A. Pople, J. Chem. Phys. 54, 724 (1971).
  * Hehre  _et al._ (1972) W. J. Hehre, R. Ditchfield, and J. A. Pople, J. Chem. Phys. 56, 2257 (1972).
  * Hariharan and Pople (1973) P. C. Hariharan and J. A. Pople, Theor. Chim. Acta 28, 213 (1973).
  * Ufimtsev and Martinez (2009) I. S. Ufimtsev and T. J. Martinez, J. Chem. Theory Comput. 5, 2619 (2009).
  * Titov  _et al._ (2013) A. V. Titov, I. S. Ufimtsev, N. Luehr, and T. J. Martinez, J. Chem. Theory Comput. 9, 213 (2013).
  * Seritan  _et al._ (2021) S. Seritan, C. Bannwarth, B. S. Fales, E. G. Hohenstein, C. M. Isborn, S. I. Kokkila-Schumacher, X. Li, F. Liu, N. Luehr, J. W. Snyder Jr, _et al._ , WIREs Comput. Mol. Sci. 11, e1494 (2021).
  * Kraft (1988) D. Kraft, Forschungsbericht- Deutsche Forschungs- und Versuchsanstalt fur Luft- und Raumfahrt (1988).
  * Virtanen  _et al._ (2020) P. Virtanen, R. Gommers, T. E. Oliphant, M. Haberland, T. Reddy, D. Cournapeau, E. Burovski, P. Peterson, W. Weckesser, J. Bright, _et al._ , Nature methods 17, 261 (2020).