# Leveraging Analog Quantum Computing with Neutral Atoms for Solvent
Configuration Prediction in Drug Discovery

Mauro D’Arcangelo  Pasqal, 7 Rue Léonard de Vinci, 91300 Massy, France  these
authors contributed equally to this work  Daniele Loco  Qubit Pharmaceuticals,
Advanced Research Department, 24 rue du Faubourg Saint-Jacques, 75014 Paris,
France  these authors contributed equally to this work  Fresnel team  Pasqal,
7 Rue Léonard de Vinci, 91300 Massy, France  Nicolaï Gouraud  Qubit
Pharmaceuticals, Advanced Research Department, 24 rue du Faubourg Saint-
Jacques, 75014 Paris, France  Sorbonne Université, Laboratoire Jacques-Louis
Lions, UMR 7589 CNRS, 75005, Paris, France  Sorbonne Université, Laboratoire
de Chimie Théorique, UMR 7616 CNRS, 75005, Paris, France  Stanislas Angebault
Qubit Pharmaceuticals, Advanced Research Department, 24 rue du Faubourg Saint-
Jacques, 75014 Paris, France  Jules Sueiro  Qubit Pharmaceuticals, Advanced
Research Department, 24 rue du Faubourg Saint-Jacques, 75014 Paris, France
Pierre Monmarché  Sorbonne Université, Laboratoire Jacques-Louis Lions, UMR
7589 CNRS, 75005, Paris, France  Jérôme Forêt  Qubit Pharmaceuticals, Advanced
Research Department, 24 rue du Faubourg Saint-Jacques, 75014 Paris, France
Louis-Paul Henry  Pasqal, 7 Rue Léonard de Vinci, 91300 Massy, France  Loïc
Henriet  Pasqal, 7 Rue Léonard de Vinci, 91300 Massy, France
<EMAIL_ADDRESS><EMAIL_ADDRESS>Jean-
Philip Piquemal  Qubit Pharmaceuticals, Advanced Research Department, 24 rue
du Faubourg Saint-Jacques, 75014 Paris, France  Sorbonne Université,
Laboratoire de Chimie Théorique, UMR 7616 CNRS, 75005, Paris, France
<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

We introduce quantum algorithms able to sample equilibrium water solvent
molecules configurations within proteins thanks to analog quantum computing.
To do so, we combine a quantum placement strategy to the 3D Reference
Interaction Site Model (3D-RISM), an approach capable of predicting continuous
solvent distributions. The intrinsic quantum nature of such coupling
guarantees molecules not to be placed too close to each other, a constraint
usually imposed by hand in classical approaches. We present first a full
quantum adiabatic evolution model that uses a local Rydberg Hamiltonian to
cast the general problem into an anti-ferromagnetic Ising model. Its solution,
an NP-hard problem in classical computing, is embodied into a Rydberg atom
array Quantum Processing Unit (QPU). Following a classical emulator
implementation, a QPU portage allows to experimentally validate the algorithm
performances on an actual quantum computer. As a perspective of use on next
generation devices, we emulate a second hybrid quantum-classical version of
the algorithm. Such a variational quantum approach (VQA) uses a classical
Bayesian minimization routine to find the optimal laser parameters. Overall,
these Quantum-3D-RISM (Q-3D-RISM) algorithms open a new route towards the
application of analog quantum computing in molecular modelling and drug
design.

###### keywords:

Keyword1, Keyword2, Keyword3

The direct manipulation of quantum systems to perform Quantum Computations has
become an intense field of interdisciplinary research aiming to apply quantum
computing to various fields ranging from theoretical chemistry and many-body
physics to material sciences and drug design. For example, such technologies
are believed to be ultimately able to improve the accuracy of Quantum
Chemistry methods [1, 2, 3, 4] through algorithms exponentially faster than
classical ones. [1, 2, 5]. If no widely accepted proof of exponential
advantage for quantum computing calculations has been yet provided [6], in
practice, a polynomial advantage could still provide more accurate simulations
at hand in near-term quantum devices [1, 6, 5, 7]. More pragmatically, a lower
energy consumption of quantum devices is expected to be more rapidly
envisioned.

Reasons for this present situation are linked to near-term quantum
technologies limitations in term of qubit count and to the absence of fully
error-corrected gate-based machines leading to high levels of noise [2, 7, 8].
In that context, analog quantum computing is an alternative promising strategy
[9, 10].

Interestingly, on the software side, many groups have contributed to develop
algorithms that are expected to scale efficiently with the size of chemical
systems, once that more advanced hardware will become available [1, 2]. Such
techniques include the Variational Quantum Algorithms (VQAs) [5, 7], a family
of algorithms allowing to exploit present hardware through coupling with
classical optimization. In a nutshell, if one is able to formulate a problem
in terms of a mathematical cost function to be minimized/maximized with
respect to the parameters of the quantum machine, the search for the best
parameters can be outsourced to a classical hardware, the quantum computer
being used to compute the cost function, a task usually too expensive to
compute classically. The solution to the problem is then given in terms of the
best parameters for the quantum evolution. Concerning the application of
quantum computing to drug design, a necessary tool for accurate drug-target
affinity predictions in structure-based Drug Discovery (DD) is related to
methods capable of addressing the problem of predicting the hydration sites in
protein structures and in any biomolecular systems to a larger extend.
Predicting the solvent structure inside important protein pockets thanks to
user-friendly and computationally-efficient methods is a determining task to
improve the predicting power of numerical simulations. Indeed, it requires to
deal with the computational effort of sampling a highly dimensional phase
space, generated by the daunting number of degrees of freedom of the combined
protein-environment complex system. There are methods such as 3D Reference
Interaction Site Model (3D-RISM) that can avoid this effort through the use of
continuous distributions [11, 12, 13, 14, 15, 16] based on a statistical
mechanics integral equation formalism. However they lack information about the
exact locations of molecules. In this work, we present the first algorithms,
exploiting the analog quantum computing paradigm, to analyze the hydration
structures within proteins’ binding sites thanks to the 3D-RISM continuous
distribution. Indeed, a poor treatment of the protein-hydration problem is
often a limiting factor preventing the access to accurate protein- ligand
affinity predictions, a crucial step for practitioners in drug discovery. Our
strategy is applied to small molecules of interest in drug design for which
the effect of solvent water molecules is substantial. The main advantage of
this proposed approach is its potential to efficiently sample the distribution
of water molecules inside protein cavities. Therefore, this work presents the
design of novel quantum algorithmic strategies in that direction. To do so, we
implemented the approach on both PASQAL’s Pulser emulator and real neutral
atom quantum computer, which, thanks to its flexible connectivity and ease of
implementation of analog operations, has demonstrated the potential to solve
combinatorial problems. The sampling of water molecules in biological systems
distribution can also be seen as a combinatorial problem, since it corresponds
to finding the most probable disposition of solvent molecules given the
constraints imposed by both the protein structure and the rest of the solvent
atomistic structure.

The importance of water placement. More precisely, the distribution of water
molecules impacts the protein’s structure and determines its overall shape.
Water molecules also mediate the interaction between proteins and other small
molecules, the so-called ligands. The presence of water generally influences
the binding capacity of a ligand to a specific protein site [17, 18], a
problem of particular interest for the discovery of new drugs. In practice,
water molecules can occupy ligand-protein interaction sites, and clusters of
water molecules can deeply affect the thermodynamics of ligand binding modes.
Consequently, accounting for protein solvation effects is crucial in
structure-based drug design, to reliably model the ligand activity through a
rational design of the ligand properties [17, 18]. To a certain extent, the
presence of water molecules can be experimentally established using X-ray
crystallography, or other techniques. However, experiments can present
accuracy limitations due to poor resolution [19] and it is then critical to
refine the structure using accurate molecular modeling tools.

Analog quantum computing with neutral atoms. To address the solvent sampling
problem, we propose an approach based on a neutral atom Quantum Processing
Unit (QPU) coupled to a 3D-RISM computation. The machine employs arrays of
Rubidium atoms arranged in arbitrary 2D configurations defined by a layout of
optical traps of the type shown in Fig. 1. Each atom is described by its
position $\vec{q}_{i}$ and its internal state, which can be in general an
arbitrary superposition of the ground state $\ket{0}$ and a highly excited
Rydberg state $\ket{1}$. The quantum state of the collective $N$-qubit array
at time $t$ will be denoted $\ket{\psi_{t}}$.

A laser system is tuned close to the resonant frequency of the $\ket{0}$ to
$\ket{1}$ transition, so that each atom is effectively a qubit living in a
two-dimensional Hilbert space. The programmable time-dependent control fields
of the driving laser are the Rabi frequency and detuning, denoted $\Omega(t)$
and $\Delta(t)$ respectively. After evolving for a period of time $T$, the
system can be measured in the computational basis by a fluorescence imaging
process that will find certain atoms in $\ket{0}$, while the rest are inferred
to be in $\ket{1}$ [9].

The term ground state will be used in two different contexts. The ground state
of a single isolated atom is denoted $\ket{0}$, and the all-ground state of a
non-interacting $M$-qubit system refers to the state
$\ket{\psi_{0}}:=\ket{0}\otimes\ldots\otimes\ket{0}=\ket{0}^{\otimes M}$ which
is often the initial state of any quantum algorithm. The ground state of the
system coupled with interactions and a driving laser, however, corresponds to
the state of minimal energy of the Hamiltonian, which is not necessarily the
all-ground state. For instance, the Hamiltonian at $\Omega=0$ and $\Delta>0$
can be related to a known graph-theoretical problem. By interpreting each
neutral atom as a vertex in a graph, and connecting with an edge vertices
situated at a distance closer than a certain threshold dictated by the physics
of interactions, the ground state is a solution to the Maximum Independent Set
(MIS) problem [20], i.e. the problem of selecting the largest number of
vertices in a graph that are not connected to each other directly by an edge.

A quantum Ising model formulation of the placement problem. To study the water
placement problem from a quantum perspective, we place ourselves in the
formalim of 3D-RISM. Indeed, this algorithm produces a continuous map of the
oxygen atom density inside the protein cavity. We want to translate this
continuous information into a discrete set of water molecule positions. In
order to do so, we assume that each water molecule randomly oscillates around
a definite stable position in the protein cavity. If one tracks over time the
oscillating motion of a water molecule and assigns to each unit of space the
probability of finding the molecule there, the resulting probability map
should resemble a Gaussian of a certain width centered around the stable
position. In this ansatz, the 3D-RISM density is interpreted as the sum of all
the Gaussians associated to each water molecule in the cavity. Therefore, the
joint density is decomposed into several single-particle densities expressing
the uncertainty in the position of each water molecule, and the problem of
locating the position of water molecules is transposed to the problem of
finding the best Gaussian mixture that approximates the 3D-RISM density.

Formally, a protein cavity can be modeled as a connected subset
$\mathcal{C}\subset\mathbb{R}^{3}$ and a 3D-RISM density is a scalar function
$g:\mathcal{C}\rightarrow\mathbb{R}$ normalized such that
$\int_{\mathcal{C}}g(\mathbf{r})d\mathbf{r}=1$. Denote by
$\mathcal{G}({\bf\it\mu},\sigma^{2})$ a normalized isotropic Gaussian with
mean ${\bf\it\mu}\in\mathbb{R}^{3}$ and variance $\sigma^{2}$. For practical
purposes, we restrict the possible values of the mean of the Gaussians to a
finite discrete lattice
$Q:=\\{\mathbf{q}_{1},\ldots,\mathbf{q}_{M}\\}\subset\mathcal{C}$ rather than
the whole continuous space. Assigning a binary variable $n_{i}\in\\{0,1\\}$ to
each point of $Q$, an arbitrary sum of Gaussians can be written:

$\sum_{i=1}^{M}\mathcal{G}(\mathbf{q}_{i},\sigma^{2})n_{i}.$ (1)

The value of each $n_{i}$ will then act as a switch that indicates whether or
not a Gaussian is placed in position $\mathbf{q}_{i}$. We want to find the
optimal assignment of $\\{n_{i}\\}_{i=1,\ldots,M}$ such that the following
$L^{2}$ norm is minimized:

$I^{2}:=\int_{\mathcal{C}}\left(g(\mathbf{r})-\sum_{i=1}^{M}\mathcal{G}(\mathbf{q}_{i},\sigma^{2})(\mathbf{r})\
n_{i}\right)^{2}d\mathbf{r}.$ (2)

Expanding the square, one sees that $I^{2}$ defines the energy of an Ising
model:

$I^{2}=K-\sum_{i=1}^{M}\Gamma_{i}n_{i}+\sum_{i\neq j=1}^{M}V_{ij}n_{i}n_{j}$
(3)

where $K$ is an unimportant constant, and the coefficients of the linear and
quadratic terms are given by:

$\Gamma_{i}:=2\int_{\mathcal{C}}g(\mathbf{r})\
\mathcal{G}(\mathbf{q}_{i},\sigma^{2})(\mathbf{r})\
d\mathbf{r}-\int_{\mathcal{C}}\Big{(}\mathcal{G}(\mathbf{q}_{i},\sigma^{2})(\mathbf{r})\Big{)}^{2}\
d\mathbf{r}$ (4)

and

$V_{ij}:=\int_{\mathcal{C}}\mathcal{G}(\mathbf{q}_{i},\sigma^{2})(\mathbf{r})\
\mathcal{G}(\mathbf{q}_{j},\sigma^{2})(\mathbf{r})\ d\mathbf{r}.$ (5)

For fixed $\sigma^{2}$, the coefficients can be computed numerically, but it
is helpful to gain further insight into the interaction term $V_{ij}$. When
the Gaussians are concentrated far enough from the boundary of $\mathcal{C}$,
it is possible to approximate the integral by extending the integration to the
whole of $\mathbb{R}^{3}$. By performing the change of variables
$\mathbf{r}\rightarrow\mathbf{r}+\mathbf{q}_{i}$, the term becomes the
convolution of two Gaussians centered in zero:

$V_{ij}\simeq\int_{\mathbb{R}^{3}}\mathcal{G}(\mathbf{0},\sigma^{2})(\mathbf{r})\mathcal{G}(\mathbf{r}_{ij},\sigma^{2})(\mathbf{r})d\mathbf{r}=\left[\mathcal{G}(\mathbf{0},\sigma^{2})*\mathcal{G}(\mathbf{0},\sigma^{2})\right](\mathbf{r}_{ij})$
(6)

where $\mathbf{r}_{ij}=\mathbf{q}_{j}-\mathbf{q}_{i}$. It is a known result
that the convolution of two Gaussians is itself a Gaussian, therefore one has:

$V_{ij}\sim\exp\left(-\alpha|\mathbf{r}_{ij}|^{2}\right).$ (7)

Equation (3) is therefore a classical Ising model with exponentially decaying
interactions. In addition, the water placement problem requires two water
molecules to be placed at a minimal physical distance from each other. This
translates to an extra constraint on the Ising model (3) in the form of
$n_{i}=n_{j}=1\implies|\mathbf{r}_{i}-\mathbf{r}_{j}|>R$ for some $R>0$. If
$\\{n_{1}^{*},\ldots,n_{M}^{*}\\}$ is the ground state of the constrained
Ising model, the solution to the water placement problem is then defined as:

$\displaystyle\mathcal{W}$ $\displaystyle:=\\{\mathbf{q}_{i}\in Q\mid
n_{i}^{*}=1\\}$ (8) $\displaystyle\mathcal{N}$
$\displaystyle:=\sum_{i=1}^{M}n_{i}^{*}$ (9)

with $\mathcal{N}$ the number of placed water molecules and $\mathcal{W}$
their positions. This formulation has far-reaching applications to general
finite Gaussian mixture models. See Supplementary Information for more
details.

The quantum version of an Ising model [21, 22] is obtained by replacing the
binary variable $n_{i}$ with the number operator $\hat{n}_{i}$, whose spectrum
is $\\{0,1\\}$. In terms of Pauli matrices, the number operator is
$\hat{n}_{i}=(\hat{\sigma}_{i}^{z}+\mathds{1})/2$. Performing the substitution
$n_{i}\rightarrow\hat{n}_{i}$, we can establish a direct identification
between the classical Ising model (3) and the (diagonal) Hamiltonian of a
system of interacting spins located in $\mathbf{q}_{1},\ldots,\mathbf{q}_{M}$:

$I^{2}\rightarrow\hat{I}^{2}:=-\sum_{i=1}^{N}\Gamma_{i}\hat{n}_{i}+\sum_{i\neq
j=1}^{N}V_{ij}\hat{n}_{i}\hat{n}_{j}$ (10)

which we will be referring to as the problem Hamiltonian. If $\mathcal{H}$
denotes the Hilbert space of the quantum system and $\mathcal{B}(\mathcal{H})$
denotes its computational basis (the basis in which the operator in eq. (10)
is diagonal), the ground state of (10) corresponds to the bitstring that
minimizes (3), but it is now seen as the computational basis vector
$\ket{e^{*}}\in\mathcal{B}(\mathcal{H})$ that minimizes the expectation value
of $\hat{I}^{2}$:

$\ket{e^{*}}:=\operatorname*{arg\,min}_{\ket{e}\in\mathcal{B}(\mathcal{H})}\bra{e}\hat{I}^{2}\ket{e}.$
(11)

The position and number of the water molecules are then given by:

$\displaystyle\mathcal{W}$ $\displaystyle:=\\{\mathbf{q}_{i}\in Q\mid\
\bra{e^{*}}\hat{n}_{i}\ket{e^{*}}=1\\}$ (12) $\displaystyle\mathcal{N}$
$\displaystyle:=\sum_{i=1}^{M}\bra{e^{*}}\hat{n}_{i}\ket{e^{*}}.$

Solving the quantum Ising problem with a local Rydberg Hamiltonian. A system
of neutral atoms coupled to an optical laser can be crafted in such a way as
to evolve according to the following time-dependent Hamiltonian:

$\hat{H}(t)=\sum_{i=1}^{M}\Omega_{i}(t)\hat{\sigma}^{x}_{i}-\sum_{i=1}^{M}\Delta_{i}(t)\hat{n}_{i}+\sum_{i<j=1}^{M}U_{ij}\hat{n}_{i}\hat{n}_{j}$
(13)

where $\Omega_{i}(t)$ is the Rabi frequency of the driving laser on qubit $i$,
$\Delta_{i}(t)$ is the detuning of the laser, $U_{ij}$ is the interaction
coefficient between Rydberg excitations:

$U_{ij}=\frac{C_{6}}{|\mathbf{r}_{ij}|^{6}}$ (14)

and $C_{6}$ is a physical constant. Before describing each term more in
detail, one can immediately see the similarity between (10) and (13), and
consequently why systems of neutral atoms might be particularly suited for
solving such a problem. The biggest difference lies in the coefficient of the
two-body term: in the Rydberg Hamiltonian $U_{ij}$ decays as a power-law,
while in the problem Hamiltonian $V_{ij}$ decays exponentially. However, one
should keep in mind that (10) does not directly encode the proximity
constraint between neighbouring excitations, which is taken care of precisely
by the $r^{-6}$ part of the Rydberg Hamiltonian.

Having control over the Rydberg Hamiltonian (13), the solution to the water
placement problem can be found by applying the adiabatic theorem [23]. In a
nutshell, the strategy is to identify a path $\hat{H}(s)$ in the space of
admissible Hamiltonians for some normalized time parameter $s\in[0,1]$ such
that the system is initialized in the ground state of $\hat{H}(0)$ (which
therefore has to be known and easy to prepare), while the only known property
of the ground state of $\hat{H}(1)$ is that it encodes the solution to a hard
combinatorial optimization problem. If the evolution is performed slowly
enough, then the system remains at all times in the instantaneous ground state
of $\hat{H}(s)$ for all $s$, and therefore measurements of the system at $s=1$
give the solution to the optimization problem.

We can exploit the quantum adiabatic evolution (QAE) algorithm by choosing
$\hat{H}(0)$ to be (13) with $\Omega_{i}(0)=0$, $\Delta_{i}(0)=-c$ for some
large positive constant $c$, so that the ground state of $\hat{H}(0)$ is, to a
very good approximation, the all-ground state $\ket{0}^{\otimes M}$. The final
Hamiltonian $\hat{H}(1)$ instead is chosen in such a way as to maximize the
overlap between its low-energy spectrum and the one of the problem Hamiltonian
$\hat{I}^{2}$. This can be achieved by mapping the one-body terms $\Gamma_{i}$
in (3) to the final detunings $\Delta_{i}(1)$ in (13). A schematic example of
$\Omega_{i}(s)$ and $\Delta_{i}(s)$ is shown in Fig. 2.

QAE, however, cannot be implemented on the current generation of neutral atom
QPUs because of the lack of local addressing. In other words, instead of the
local terms $\Omega_{i}(t)$ and $\Delta_{i}(t)$, one only has access to a
global control field $\Omega(t)$ and $\Delta(t)$. Nevertheless, the QAE
algorithm, described in Algorithm 1, can be still tested using the Pulser
classical emulator [24].

Algorithm 1 Quantum Adiabatic Evolution (QAE) algorithm that solves the
corresponding Ising problem in eq. (10) using local lasers.
inputs: reference 2D density $g(\vec{r}):\mathcal{C}\rightarrow\mathbb{R}$ and
set of qubits position $Q:=\\{\vec{q}_{i}\\}\subset\mathcal{C}$, with
$\mathcal{C}\subset\mathbb{R}^{2}$. $\Omega$ is set according to the known
adiabatic protocol, fitting the three points $[0,\Omega_{Max},0]$ during the
whole pulse duration $[0,T]$ and $\Delta$ is computed according to eq. (4).
outputs: water molecule positions $\mathcal{W}$ and their number
$\mathcal{N}$.

1:

2:procedure QAEAlgo($g(\vec{r}),\\{\vec{q}_{i}\\}$)

3:

4: $\Delta_{i}(T)\leftarrow$
FinalDetuningPerQubit($g(\vec{r}),\\{\vec{q}_{i}\\}$)

5: $\Omega(t),\Delta_{i}(t)\leftarrow$ BuildAdiabaticPulse($\Delta_{i}(T)$)

6: $\ket{e^{*}}\leftarrow$
QuantumAdiabaticSelection($\\{\vec{q}_{i}\\},\Omega(t),\Delta_{i}(t)$)
$\triangleright$ Quantum evolution according to the built laser pulse and
identification of the best basis state for the given $g(\vec{r})$

7:

8: $(\mathcal{W},\mathcal{N})\leftarrow$ PositionsFromState($\ket{e^{*}}$)
$\triangleright$ $(\mathcal{W},\mathcal{N})$ are defined in eq. (12)

9: return $(\mathcal{W},\mathcal{N})$ $\triangleright$ Best placement of water
molecules in, e.g., target protein cavity

10:

11:end procedure

We first tested this approach using synthetic densities. Results can be found
in SI, Section LABEL:sec:num-sim-supportQAE.

Local algorithm emulation using 3D-RISM densities. To test the algorithm just
presented, we compute the 3D-RISM solvent density within a real protein. We
chose the major urinary protein (MUP-I) pocket, where a small ligand, the
2-sec-butyl-4,5-dihydrothiazole, is binded to the protein (see Fig. 9). This
choice is motivated by the fact that the protein-ligand complex structure has
been co-crystallized, its complete atomic structure being available from the
Protein Data Bank (PDB). Interestingly, the crystal structure presents two
structural (i.e. stable) water molecules present in the vicinity of the
ligand. Such water molecules are clearly visible in the 3D-RISM density which
exhibits high density spots in the same region as the crystal structure (see
Fig. 9). The two water molecules are included in a diffused volume of 3D-RISM
density that is identified by a density isovalue surfaces of medium value (see
Fig. 9 caption) whereas maximal density values are spotted close to the
experimental crystal structures water molecules positions (see orange solid
surface in Fig. 9, b).

Despite the algorithm being well-defined in any dimension, current generation
neutral atom quantum computers can only operate on two-dimensional qubit
registers. For this reason we cut the 3D-RISM density (which, by construction,
is defined in three dimensions) into two-dimensional slices, and we apply the
algorithm separately to each slice. For simplicity, we restrict ourselves to a
small region of the protein around the two crystal water molecules. The plane
corresponding to the first slice is defined by the axis connecting the two
crystallographic water molecules and a third random direction. A series of
slices is then produced from this first one in the direction normal to this
plane, separated by a spacing of 0.5 Å, for a total of 6 different slices. To
smooth the 3D-RISM density in each slice, we apply a Laplacian of Gaussian
(LoG) filter as provided in the SciPy Python library [25], with a $\sigma=8$
(a.u.). The high-density regions of each smoothed-out slice are then covered
with a qubit register. The process is depicted schematically in Fig. 3.
Overall, this test implementation requires 4 to 14 qubits, which remains in
the range of qubit count that can be classically emulated with high precision
by a state vector solver with a realistic noise model. The density slices with
associated qubit registers are shown in Fig. 4. In order to mimic a physical
implementation, we fit the registers in the same trap layout of Fig. 1, as
shown in Fig. 5.

For each register, we design an adiabatic pulse to find the low energy states
of the Rydberg Hamiltonian (13). The detunings $\Delta_{i}$ are computed from
the $\Gamma_{i}$ coefficients of each problem Hamiltonian by using the
following heuristic mapping:

$\tilde{\Delta}_{i}=\Gamma_{i}-\frac{1}{|N(i)|}\sum_{j\in
N(i)}\Gamma_{j},\quad\quad\Delta_{i}=\tilde{\Delta}_{i}\frac{\Delta_{\text{max}}}{\max_{i}|\tilde{\Delta}_{i}|}$
(15)

where $N(i)$ denotes nodes situated within a Rydberg blockade distance from
node $i$, and $\Delta_{\max}$ is a fixed value. The mapping (15) can be
understood as a two-step process. First, each $\Gamma_{i}$ is shifted so that
it is centered around the mean value of neighbouring nodes, obtaining the
intermediate detunings $\tilde{\Delta}_{i}$. Then, the $\tilde{\Delta}_{i}$
are rescaled so that the maximum one corresponds to a certain fixed value
$\Delta_{\max}$, obtaining the final value $\Delta_{i}$.

We simulate the resulting adiabatic dynamics with Pulser, and sample the final
state 1000 times. For each computational basis state sampled, we calculate its
cost according to the problem Hamiltonian, and select as a solution the one
with the minimal cost. The resulting sampling, where the basis states are
represented as bitstrings, is shown in Fig. 6 for all 6 slices. The orange bar
in each plot identifies the bitstring corresponding to the best solution for
that slice. The fact the bitstring with the lowest cost is not always the most
sampled one is a consequence of the differences between the low energy
spectrum of problem and Rydberg Hamiltonians.

Experimental measurement on QPU exploiting the local algorithm results. The
emulations presented in the previous section correspond to the capabilities of
a neutral atom QPU capable of local addressing, which will be available in the
next generation of neutral atom devices. An earlier prototype, with limited
functionality, of such machines was tested here to confirm that solutions to
the water placement problem belong to the space of physically allowed
configurations of the machine.

In this simplified setting we use the same 3D-RISM densities and qubit
registers presented in Fig. 5 and Fig. 5. Each register was associated to the
solution found by simulating the local algorithm on a classical computer, as
reported in Fig. 6. The qubit register is then reproduced in the QPU, and its
quantum state is measured after evolving under a constant global laser pulse.
Tunable parameters of the driving laser are its Rabi frequency $\Omega$,
detuning $\Delta$ and duration. If $b_{i}$ indicates the bitstring encoding
the solution to the water placement problem for the $i$-th register, and
$\ket{b_{i}}$ is the associated computational basis state, we evaluate the
expectation value of its projector $\ket{b_{i}}\bra{b_{i}}$ on the final
quantum state, in a two-dimensional parameter space spanned by detuning and
duration of the driving laser, while the Rabi frequency is fixed to a value
that would enforce a Rydberg blockade between nearest-neighbours.

The experimental measurements of the projector for all four registers are
shown in Fig. 7, together with theoretical expectations given a reasonable
estimate for false positive and false negative error rates on the machine.
Each point in the 2D plane is calculated as an average over roughly 500
measurements. On a qualitative level, the experiments can be seen to be
compatible with theoretical predictions obtained from emulation. These results
indicate that meaningful solutions to the water placement problem belong to
the space of physically allowed configurations of the QPU. To our knowledge,
such use of Rydberg physics had not been explored before in a real quantum
physics experiment. Furthermore, the experiment shows that probability
landscapes in this simplified setting can be accurately resolved by the
machine, even in the absence of error correction or error mitigation
techniques. This is a key result knowing the small magnitude of the involved
probabilities that are capped at only a few percent in certain systems.

In order to better understand and quantitatively validate the experimental
data produced by the set-up, we considered a more sophisticated error model.
Calibration of the control devices can only be achieved with finite precision,
resulting in static uncertainties in global spacing of the atomic array
($\approx 1\%$) or in the spatial homogeneity of $\Omega$ on this array
($\approx 4\%$). In addition, fluctuations of laser intensity induces a shot-
to-shot variation of $\Omega$ ($\approx 5\%$). While the laser frequency can
be set with high precision, variation in $\Omega$ indirectly alters $\Delta$,
resulting in small detuning shifts of the order of $2\pi\times 0.06$ MHz.
Decay processes are also taken into account by solving the Master equation
with an effective decay rate $\Gamma_{eff}/2\pi$. Finally, the measurement
phase is inherently flawed by several physical processes like atomic losses
due to background-gas collisions or Rydberg state finite lifetime, whose
effects can all be encompassed as first approximation into two detection error
terms, $\varepsilon$ and $\varepsilon^{\prime}$. The various values of
experimental and noise parameters are usually fitted by comparing the
expectation value of easy-to-access observables between simulated and
experimentally acquired data. For instance, emulating the dynamics described
by the Hamiltonian (16) in presence or absence of noise processes and
measuring the occupation $\hat{n}_{i}$ at each site enables to estimate for
the magnitude of the various error sources. An example of fitting curves at a
given detuning $\Delta$ for the first register is shown in Fig. 8. While the
errors on calibration parameters are directly measured on the experiment, the
remaining effective noise parameter are found to be $\Gamma_{eff}/2\pi=0.05$
MHz, $\varepsilon=2\%$ and $\varepsilon^{\prime}=18\%$.

Overall, given the experimental errors sources that were discussed, the
experimental data appear fully compatible with the emulated ones (see Fig. 7
and 8) confirming the viability of the proposed approach on an actual cold
atom device. We want to stress here that thanks to the quantum nature of the
algorithm, two water molecules are guaranteed to never be placed too close to
each other, a constraint that is imposed by hand in classical approaches such
as Placevent [26] or GAsol [27] with the risk of incurring in suboptimal local
solutions or non-ergodicity. This constitutes one of the main advantage of the
algorithm.

Beyond a Local Algorithm: Variational Algorithm using a global Rydberg
Hamiltonian. As we discussed, technical limitations of the present hardware
prevent us from using a local Hamiltonian as in equation (13). However,
through emulation, we can still plan for an algorithm that will be able to run
on more short-term devices using a global Rydberg Hamiltonian:

$\hat{H}(t)=\Omega(t)\sum_{i=1}^{M}\hat{\sigma}^{x}_{i}-\Delta(t)\sum_{i=1}^{M}\hat{n}_{i}+\sum_{i<j=1}^{M}\frac{C_{6}}{r_{ij}^{6}}\hat{n}_{i}\hat{n}_{j}.$
(16)

Such formulation allows to further approximate the problem Hamiltonian in eq.
(10), since the linear term in the qubits excitation is no longer local. A
variational procedure is then established using, as cost function to be
minimized the problem Hamiltonian in eq. (10), so to solve

$(\Omega(t),\Delta(t))^{*}=\displaystyle\operatorname*{arg\,min}_{\Omega(t),\Delta(t)}\bra{\Psi^{\Omega,\Delta}}\hat{I}^{2}\ket{\Psi^{\Omega,\Delta}}.$
(17)

This algorithm is therefore part of the Variational Quantum Algorithm (VQA)
family The cost function
$\bra{\Psi^{\Omega,\Delta}}\hat{I}^{2}\ket{\Psi^{\Omega,\Delta}}$ is evaluated
summing up the contributions of each basis state sampled from a trial wave
function, obtained from the laser parameters $(\Omega(t),\Delta(t))$ produced
during each optimization step. The full VQA developed for this scope is
described in Algorithm 2.

Algorithm 2 Hybrid quantum-classical Variational Quantum Algorithm (VQA) using
global lasers, with cost function issued from the Ising model, namely
$\bra{\Psi^{\Omega,\Delta}}\hat{I}^{2}\ket{\Psi^{\Omega,\Delta}}$.
inputs: reference 2D density $g(\vec{r}):\mathcal{C}\rightarrow\mathbb{R}$ and
set of qubits position $Q:=\\{\vec{q}_{i}\\}\subset\mathcal{C}$, with
$\mathcal{C}\subset\mathbb{R}^{2}$. $n_{c}$ and $n_{r}$ are parameters for the
Bayesian optimization procedure, being, respectively, the total number of
cycles performed to generate a single state $\ket{\psi}_{k=1,\dots,n_{c}}$ and
$n_{r}<n_{c}$ is the number of cycles initiated with randomized values of
laser parameters $\Omega$ and $\Delta$.
outputs: water molecule positions $\mathcal{W}$ and their number
$\mathcal{N}$.

1:

2:procedure VQAgo($g(\vec{r}),\\{\vec{q}_{i}\\},n_{r},n_{c}$)

3:

4: $k\leftarrow 0$

5: while $k<n_{c}$ do

6: if $k\leq n_{r}$ then$\triangleright$ The first $n_{r}$ iterations are
uniformly randomized

7: $\Omega_{k}$, $\Delta_{k}\leftarrow$ Random($\mathcal{U}$)

8: else

9: $\Omega_{k},\Delta_{k}\leftarrow$
BayesianMinimization($\\{\Omega_{l},\Delta_{l}\\}_{0\leq l<k}$,
$\bra{\Psi^{\Omega,\Delta}}\hat{I}^{2}\ket{\Psi^{\Omega,\Delta}}$)

10: end if

11: $\ket{\Psi_{k}}\leftarrow$
QuantumEvolution($\\{\vec{q}_{i}\\},\Omega_{k},\Delta_{k}$)

12: $J_{k}\leftarrow\bra{\Psi_{k}}\hat{I}^{2}\ket{\Psi_{k}}$ $\triangleright$
Collection of $J_{k}\forall k$ to evaluate later best laser parameters

13: end while

14: $\Omega^{*},\Delta^{*}\leftarrow$ ArgMin($\\{J_{k}\\}$) $\triangleright$
Select the $\Omega(t)$ and $\Delta(t)$ that minimize the cost function
according to eq. (17)

15: $\ket{\Psi^{*}}\leftarrow$
QuantumEvolution($\\{\vec{q}_{i}\\},\Omega^{*},\Delta^{*}$)

16: $\ket{e^{*}}\leftarrow$ BestBasisState($\ket{\Psi^{*}}$) $\triangleright$
Basis state whose GMA best represents $g(\vec{r})$ sampled from the optimized
$N$-qubits wavefunction

17:

18: $(\mathcal{W},\mathcal{N})\leftarrow$ PositionsFromState($\ket{e^{*}}$)
$\triangleright$ $(\mathcal{W},\mathcal{N})$ are defined in eq. (12)

19: return $(\mathcal{W},\mathcal{N})$ $\triangleright$ Best placement of
water molecules in, e.g., target protein cavity

20:

21:end procedure

$\ket{\Psi^{\Omega,\Delta}}$ in Eq. (17) is a short notation for a quantum
state obtained from the evolution dictated by the Hamiltonian in eq. (16) for
a specific set of laser parameter $(\Omega(t),\Delta(t))$.

The objective of the numerical procedure is to maximize the probability of
sampling the basis state (11), from the optimized $N$-qubits wavefunction.
This corresponds to finding one or more configurations of excited qubits best
representing the 3D-RISM density distribution as a sum of gaussian
distributions. The final output is still described by the quantities in Eq.
(12).

To test the performances of the algorithm we use again simple synthetic
densities to limit the number of qubits to employ, so that a cycle of
optimization, emulating multiple times the quantum evolution of the system on
the CPU, is performed for each test case. The results are reported in SI (see
section LABEL:num-sim-support). They show that the algorithm is able to give
the correct positions of the Gaussian distributions centroids. Due to the high
computational cost involved in the emulation of a numerical solution of the
time-dependent Shrödinger equation, we could not extend our tests towards
larger systems.

Furthermore, We used as an additional test case one of the slices obtained
from the 3D-RISM density, and already presented for the QAE approach. The
sliced density, together with the qubit register, are shown in Fig. 9 (panel
c). We perform 50 cycles of Bayesian optimization, using 200 samples to
represent each wavefunction produced by the corresponding set of laser
parameters. The gaussians used to dress the excited qubits, so to assign
scores to the corresponding bitstring, have amplitudes proportional to the
local value of the 3D-RISM density and an uniform variance $\sigma^{2}=5$
(a.u.). From the resulting final wavefunction, we obtain the 3D coordinates of
the oxygen atoms of the placed water molecules. In Fig. 10 we report the best
water molecules configuration found by the algorithm (cyan atoms),
corresponding to the best bitstring $\ket{101001}$ found in the optimization.
Also two excited qubits states, as $\ket{010100}$ and $\ket{010001}$ show good
sampling probability (see corresponding histogram).

Additional emulated experiments are performed on all the other slices and the
results are reported in the SI. The algorithm exhibits equivalent performances
in all these cases. Again, we stress that this variational quantum algorithm,
like the version using the local lasers, prevents two water molecules from
being placed too close to each other, a constraint that is imposed by hand in
classical approaches.

Conclusion and perspectives. In this work we presented two new quantum
algorithms, able to sample equilibrium solvent configurations within proteins.
We proposed a first fully local quantum adiabatic evolution version to be used
on next generation devices whereas, a second version, belonging to the
variational quantum algorithm family of algorithms has also been introduced as
a viable short-term alternative. This new class of algorithms correspond to
quantum versions of the 3D Reference Interaction Site Model (Quantum-3D-RISM
or Q-3D-RISM) since they use continuous solvent distributions and are able to
efficiently locate density distributions’ maxima using a neutral atom quantum
computer. As a proof of concept, both algorithms have been shown to
successfully be able to locate density maxima in non-trivial densities. In the
case of the second VQA algorithm, a classical optimization is performed to
find the best set of laser pulses, capable of producing the expected
distribution of maxima. A more simplified version was implemented on a real
neutral atom QPU, which represents to the best of our knowledge the first
application of an analog quantum algorithm to the protein solvation problem.
Presently, we limited ourselves to a qubit count of 14. This number was
constrained by the present machine layout (see Fig. 1 and Fig. 5 ), by the
state vector emulation capabilities (roughly 16 qubits can be emulated easily
with Pulser) and the need of coupling Pulser and the actual QPU for the
quantum adiabatic evolution model implementation. In that connection,
concerning the local algorithm, future availability of time-dependent pulses
and local addressing will totally remove the need of the Pulser emulation to
prepare QPU data. In the same line, only time-dependent pulses are required
for the VQA version of Quantum-3D-RISM: it should enable us to perform larger
simulations at short-term since the upgrade of the PASQAL’s machine is
ongoing. Moreover, since the quantum versions of 3D-RISM comes with the native
advantage of preventing water molecules from being placed on top of each
other, it will be interesting to compare the performances of the classical and
quantum versions of 3D-RISM since the next QPU implementation will provide us
with the possibility to test at large scale the accuracy of such techniques.
To do so, a careful study of the machine noise will be necessary but one key
advantage of analog computing is to exhibit relatively constant noise levels
with increasing system size making us optimistic about the prospect of QPU
simulations encompassing a high number of qubits. At that stage, it will be
then possible to prepare any molecular system with such algorithms in order to
couple them to state-of-the-art molecular dynamics engines[28, 29] for further
properties evaluations. Overall, this Quantum-3D-RISM (Q-3D-RISM) family of
algorithms demonstrates promises in predicting the solvation structure within
biomolecular systems of interest for drug discovery applications, providing
concrete use cases for the application of analog quantum computing in life
sciences.

## Appendix A: Quantum computational resources

For the quantum computing sequences, we make use of Fresnel, an industrial
neutral atom QPU made of single 87Rb atoms trapped in arrays of optical
tweezers, conceived and manufactured by PASQAL. We operate the QPU in the
Ground-Rydberg qubit basis with global analog control [9]. The qubits are
encoded into the ground state $\ket{0}=\ket{5S_{1/2},F=2,m_{F}=2}$ and a
Rydberg state $\ket{1}=\ket{60S_{1/2},m_{J}=1/2}$ . This effective two-level
system is addressed with a two-photon laser excitation through an intermediate
state $6P_{3/2}$. The first (respectively second) photon excitation is
generated by a 420-nm (1013-nm) laser beam. Details about the Pulser control
software used to program the experiment can be found in [24].

## Acknowledgements

This work was made possible thanks to the Pack Quantique grant from région Ile
de France and GENCI, project ACQMED (convention N°20012758). Funding from the
European Research Council (ERC) under the European Union’s Horizon 2020
research and innovation program, project EMC2 (grant agreement N°810367), is
also acknowledged (J.-P. P) as funding from PEPR Epiq (ANR-22-PETQ-0007) and
HQI programs.

### Authors Contribution Statement

M.D, D. L performed simulations and contributed new code;
M.D, D.L, N. G, P. M, S. A, J. S, J. F, L.-P. H, J.-P. P contributed new
methodology (theory);
The Fresnel team, L. H designed and performed the experiments;
M.D, D.L, the Fresnel team, N. G., L.-P.H, L. H, J.-P. P analyzed data;
M.D, D.L, part of the Fresnel team, L.-P. H, L. H., J.-P. P wrote the paper
with the input of all authors;
L. H and J–P. P designed the research.

The Fresnel team (Pasqal): J. ARMOUGOM, D. BENVENUTTI, L. BEGUIN, L.
BOURACHOT, J. BRIAND, C. BRIOSNE FREJAVILLE, N. CARREZ, T. CARTRY, A.
CHARPENTIER, D. CLAVEAU, L. COLIN, G. COURNEZ, L. COUTURIER, J. DE HOND, S.
DESIRE, A. DUMAS, S. DUTARTRE, P. FAVIER, G. FIRENZE, D. KACZOR , C. HAMOT, G.
HERCE, J. HEURTEBIZE, V. HULLY, B. LABARRE, L. LASSABLIERE, H. LE BARS , L.
LECLERC, A. LINDBERG, G. MERIAUX, F. NAMBI, T. PANSIOT, G. PARIENTE, J.
PELLEGRINO, L. PONSOT, S. ROCHE, H. SILVERIO, G. VILLARET, J.-M. WIPFF.

## References

  * [1] Cao, Y. _et al._ Quantum chemistry in the age of quantum computing.  _Chemical Reviews_ 119, 10856–10915, DOI: 10.1021/acs.chemrev.8b00803 (2019).  PMID: 31469277, https://doi.org/10.1021/acs.chemrev.8b00803.
  * [2] McArdle, S., Endo, S., Aspuru-Guzik, A., Benjamin, S. C. & Yuan, X.  Quantum computational chemistry.  _Rev. Mod. Phys._ 92, 015003, DOI: 10.1103/RevModPhys.92.015003 (2020).
  * [3] Feniou, C. _et al._ Overlap-ADAPT-VQE: Practical quantum chemistry on quantum computers via overlap-guided compact ansätze.  _Communications Physics_ 6, 192 (2023).
  * [4] Feniou, C. _et al._ Greedy gradient-free adaptive variational quantum algorithms on a noisy intermediate scale quantum computer.  _arXiv preprint arXiv:2306.17159_ (2023).
  * [5] Bharti, K. _et al._ Noisy intermediate-scale quantum algorithms.  _Rev. Mod. Phys._ 94, 015004, DOI: 10.1103/RevModPhys.94.015004 (2022).
  * [6] Lee, S. _et al._ Is there evidence for exponential quantum advantage in quantum chemistry?, DOI: 10.48550/ARXIV.2208.02199 (2022).
  * [7] Cerezo, M. _et al._ Variational quantum algorithms.  _Nature Reviews Physics_ 3, 625–644, DOI: 10.1038/s42254-021-00348-9 (2021).
  * [8] Sennane, W., Piquemal, J.-P. & Rančić, M. J.  Calculating the ground-state energy of benzene under spatial deformations with noisy quantum computing.  _Phys. Rev. A_ 107, 012416, DOI: 10.1103/PhysRevA.107.012416 (2023).
  * [9] Henriet, L. _et al._ Quantum computing with neutral atoms.  _Quantum_ 4, 327 (2020).
  * [10] Scholl, P. _et al._ Quantum simulation of 2d antiferromagnets with hundreds of rydberg atoms.  _Nature_ 595, 233–238 (2021).
  * [11] Chandler, D. & Andersen, H. C.  Optimized Cluster Expansions for Classical Fluids. II. Theory of Molecular Liquids.  _The Journal of Chemical Physics_ 57, 1930–1937, DOI: 10.1063/1.1678513 (1972).
  * [12] Chandler, D.  Derivation of an integral equation for pair correlation functions in molecular fluids.  _The Journal of Chemical Physics_ 59, 2742–2746, DOI: 10.1063/1.1680393 (1973).
  * [13] Chandler, D.  Equilibrium structure and molecular motion in liquids.  _Accounts of Chemical Research_ 7, 246–251, DOI: 10.1021/ar50080a002 (1974).
  * [14] Beglov, D. & Roux, B.  An Integral Equation To Describe the Solvation of Polar Molecules in Liquid Water.  _The Journal of Physical Chemistry B_ 101, 7821–7826, DOI: 10.1021/jp971083h (1997).
  * [15] Kovalenko, A.  Three-dimensional Rism Theory for Molecular Liquids and Solid-Liquid Interfaces.  In Hirata, F. (ed.) _Molecular Theory of Solvation_ , Understanding Chemical Reactivity, 169–275, DOI: 10.1007/1-4020-2590-4_4 (Springer Netherlands, Dordrecht, 2003).
  * [16] Roy, D. & Kovalenko, A.  Biomolecular Simulations with the Three-Dimensional Reference Interaction Site Model with the Kovalenko-Hirata Closure Molecular Solvation Theory.  _International Journal of Molecular Sciences_ 22, 5061, DOI: 10.3390/ijms22105061 (2021).
  * [17] Samways, M. L., Taylor, R. D., Bruce Macdonald, H. E. & Essex, J. W.  Water molecules at protein–drug interfaces: computational prediction and analysis methods.  _Chem. Soc. Rev._ 50, 9104–9120, DOI: 10.1039/D0CS00151A (2021).
  * [18] Bucher, D., Stouten, P. & Triballeau, N.  Shedding light on important waters for drug design: Simulations versus grid-based methods.  _Journal of Chemical Information and Modeling_ 58, 692–699, DOI: 10.1021/acs.jcim.7b00642 (2018).  PMID: 29489352, https://doi.org/10.1021/acs.jcim.7b00642.
  * [19] Wlodawer, A., Minor, W., Dauter, Z. & Jaskolski, M.  Protein crystallography for non-crystallographers, or how to get the best (but not more) from published macromolecular structures.  _The FEBS Journal_ 275, 1–21, DOI: https://doi.org/10.1111/j.1742-4658.2007.06178.x (2008).  https://febs.onlinelibrary.wiley.com/doi/pdf/10.1111/j.1742-4658.2007.06178.x.
  * [20] Bourgeois, N., Escoffier, B., Paschos, V. & van Rooij, J. M. M.  Fast algorithms for max independent set.  _HAL open science_ (2013).
  * [21] Lucas, A.  Ising formulations of many np problems.  _Frontiers in Physics_ 2, DOI: 10.3389/fphy.2014.00005 (2014).
  * [22] Das, A. & Chakrabarti, B. K.  Colloquium: Quantum annealing and analog quantum computation.  _Rev. Mod. Phys._ 80, 1061–1081, DOI: 10.1103/RevModPhys.80.1061 (2008).
  * [23] Farhi, E. _et al._ A quantum adiabatic evolution algorithm applied to random instances of an np-complete problem.  _Science_ 292, 472–475, DOI: 10.1126/science.1057726 (2001).  https://www.science.org/doi/pdf/10.1126/science.1057726.
  * [24] Silvé rio, H. _et al._ Pulser: An open-source package for the design of pulse sequences in programmable neutral-atom arrays.  _Quantum_ 6, 629, DOI: 10.22331/q-2022-01-24-629 (2022).
  * [25] Virtanen, P. _et al._ SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python.  _Nature Methods_ 17, 261–272, DOI: 10.1038/s41592-019-0686-2 (2020).
  * [26] Sindhikara, D. J., Yoshida, N. & Hirata, F.  Placevent: an algorithm for prediction of explicit solvent atom distribution-application to HIV-1 protease and F-ATP synthase.  _J Comput Chem_ 33, 1536–1543 (2012).
  * [27] Fusani, L., Wall, I., Palmer, D. & Cortes, A.  Optimal water networks in protein cavities with GAsol and 3D-RISM.  _Bioinformatics_ 34, 1947–1948, DOI: 10.1093/bioinformatics/bty024 (2018).
  * [28] Lagardère, L. _et al._ Tinker-hp: a massively parallel molecular dynamics package for multiscale simulations of large complex systems with advanced point dipole polarizable force fields.  _Chem. Sci._ 9, 956–972, DOI: 10.1039/C7SC04531J (2018).
  * [29] Adjoua, O. _et al._ Tinker-hp: Accelerating molecular dynamics simulations of large complex systems with advanced point dipole polarizable force fields using gpus and multi-gpu systems.  _Journal of Chemical Theory and Computation_ 17, 2034–2053, DOI: 10.1021/acs.jctc.0c01164 (2021).  PMID: 33755446, https://doi.org/10.1021/acs.jctc.0c01164.

Figure 1: An example of layout of optical traps in a neutral atom QPU. Each
circle represents a site were a Rubidium atom can be trapped and used for
quantum computations.

---

---

Figure 2: Schematic example of an adiabatic path parametrized by $\Omega_{i}(s)$ and $\Delta_{i}(s)$ in the Rydberg Hamiltonian (13) for three qubits. The Rabi frequency is kept the same for all qubits, so $\Omega_{i}(s)=\Omega(s)$ for all $i$, and it vanishes at the extremal points of the path. The Detunings $\Delta_{i}(s)$ are ramped up from a negative value to some final value related to the one-body terms of the problem Hamiltonian (3), and therefore they will be different for each qubit. Figure 3: a) Example representation of a 3D-RISM density as a closed isosurface, with a define, constant, density value, as obtained starting from the AmberTools output, discretized on a given chosen grid of point. b) projection of the 3D-RISM density on 2D planes, namely slices, passing through the cross-section of the density; each slice is a 2D density map, where the 3D-RISM density value discretized on the grid is projected on the plane according to a proximity criterion. c) Process of selection of the qubits array on each 2D density slice: starting from the slice, a regular mesh of traps is disposed uniformly on the density, to then select a limited number of locations where to place qubits according to the local density value; a threshold on the density is set, so that only traps close to a defined density or higher will be occupied. |
---|---
|
|
Figure 4: Two-dimensional 3D-RISM density slices and qubit registers used in the local algorithm emulation. The registers vary in size from 4 to 14 qubits. The units displayed in the $x$ and $y$ directions correspond to the 3D-RISM discretization grid, and they are not representative of neither the size of the protein cavity in Å, nor the size of the qubit registers in $\mu m$. To get an idea of the physical size of the registers, consider the spacing between neighbouring qubits to be fixed at 5 $\mu m$.  |
---|---
|
|
Figure 5: Registers of Fig. 4, reshaped to fit in the trap layout available on the QPU (see Fig. 1). The units correspond to the physical dimension of the registers. |
---|---
|
|
Figure 6: Histogram of 1000 samples of the quantum state obtained with the
local algorithm based on adiabatic evolution. Each plot corresponds to one of
the registers shown in Fig. 4. The orange bar corresponds to the bitstring
with the minimal cost among the ones that were sampled, and it is taken as the
solution of the water placement problem for that register.

---

---

---

---

Figure 7: Expectation value of the projector operator $\ket{b_{i}}\bra{b_{i}}$
associated to the solutions $b_{i}$, $i=1,2,3,4$ of the four of the 3D-RSM 2D
slices . The area scanned is a two-dimensional region parametrized by pulse
duration (horizontal axis) and detuning (vertical axis) of the driving laser.
The plot on the left represents the experimental values, while the plot on the
right represents the values obtained from a classical simulation of the
Rydberg Hamiltonian with measurement errors $\epsilon=0.02$ and
$\epsilon^{\prime}=0.15$ related to the probability of false positive and
false negative detection. Figure 8: (a) A laser pulse of constant $\Omega$ and
$\Delta$ is applied for various times on the first atomic register comprising
of 4 qubits (colored). (b) For a fixed detuning $\Delta/2\pi=-0.87$ MHz, the
expectation value of the occupation of each qubit site, $\langle
n_{i}\rangle$, is plotted for noiseless simulation (dashed), noisy simulation
(line) and experimental data (dot). The error bars and uncertainty regions are
the standard deviation computed over $N_{shots}=500$ samples,
$\sigma_{n_{i}}=\sqrt{\langle n_{i}\rangle(1-\langle
n_{i}\rangle)/N_{shots}}$. Figure 9: a) View from inside the MUP-I protein
pocket, in a complex with the 2-sec-butyl-4,5-dihydrothiazole, in licorice
style (PDB entry: 1i06.pdb). The 3D-RISM isodensity with isovalue 5.5 is
represented in wireframe style, in cyan; the oxygen of the crystal water
molecules placed inside the protein cavity are represented in red. The VMD
software is used for the visualization. b) Zoom on the crystal water molecules
position: the cyan wireframe 3D-RISM isosurface (isovalue $=5$) is compared
with the orange solid one (isovalue $=8.5$). c) Representation of the smoothed
3D-RISM density slice, with 6 qubits. d) Histograms reporting the 6-qubit
quantum state composition, in terms of basis stats, as obtained from the VQA
emulated through Pulser on CPU, without noise, with most sampled bitstring
highlighted, corresponding to basis states 101001 (first most sampled) and
010100 (second most sampled). Figure 10: MUP-I protein complexed with the
2-sec-butyl-4,5-dihydrothiazole, inside view of the protein binding pocket.
Red atoms are oxygen atoms of the crystal water molecules, and the cyan one
are the three water molecules placed by the HQC algorithm . The 3D-RISM
density isosurface (isovalue 5.5) is represented in wireframe style, in cyan.