# Coded Relaying in LoRa Sensor Networks

Siddhartha S. Borkotoky, Pavan Datta Abbineni, Vatsalya Chaubey, and Sonu
Rathi  This work was supported by the Science and Engineering Research Board
(SERB), Government of India, through its Start-up Research Grant (SRG) under
SRG/2020/001491. The authors are with the Indian Institute of Technology
Bhubaneswar, India. (E-mail<EMAIL_ADDRESS>
###### Abstract

To enhance the reliability of LoRa sensor networks, we propose two
opportunistic forwarding protocols in which relays overhear sensor-to-gateway
messages and periodically forward them to a gateway. To reduce the number of
transmissions, the relays forward XOR sums of the overheard messages. One
protocol employs a single relay node, while two relays cooperate in the other.
With a single relay, our proposed protocol requires fewer relay transmissions
to achieve similar loss rates as conventional uncoded forwarding in which
every received message is immediately forwarded. We demonstrate up to 55%
reduction in relay transmission time relative to uncoded forwarding. The
cooperative protocol provides better loss and delay performance than the
single-relay protocol, and we demonstrate up to 16% reduction in the message-
loss rate along with 33% reduction in relay transmission times compared to
conventional forwarding.

##  I Introduction

The LoRa technology has established itself as a major enabler of IoT
connectivity due to its low-power operations, robustness against propagation
loss, and simplicity of implementation. However, interference-induced frame
losses are a major concern in LoRa [1]. Use of acknowledgements (ACKs) and
retransmissions to combat frame loss are often infeasible since scanning for
ACKs after every transmission drains sensor battery and reduces network
lifetime; duty-cycle regulations prevent the gateway from sending frequent
ACKs; and ACKs and retransmissions may themselves cause frame loss [2].
Alternative reliability mechanisms are therefore of interest. In this paper,
we propose and evaluate one such mechanism in which relay nodes are used to
provide additional routes between the sensors and gateway. Relays
opportunistically overhear the sensor messages and forward them to the
gateway, so that a message that did not reach the gateway directly may reach
via a relay. One of our primary motives is to examine the efficacy of coded
forwarding, in which a relay forwards bitwise XORs of sensor messages, as
opposed to uncoded forwarding, in which each message is forwarded as it is.

Relaying in LoRa networks and its potential benefits are reported in [3]–[8].
Uncoded forwarding is employed in [3]–[7]. In the coded relaying scheme of
[8], multiple relays employ random linear network coding to combine received
sensor messages before sending. Relays use feedback from the gateway to decide
when to code and transmit, and the coding is performed in Galois Field of
order 16 (GF(16)).

By contrast, we assume unsophisticated, low-cost relays that receive
periodically (as in [3]) and then transmit XOR sums of received messages. Such
GF(2) operations lead to low computational requirements. We assume that the
gateway does not send any feedback either to the sensors or to the relays. The
sensors transmit each message only once and do not listen to downlink frames
(except infrequent control messages carrying configuration data) so as to
prevent battery drain. No routing protocol is required; the sensors transmit
to the gateway without knowledge of the relay’s presence.

Through simulations and analysis, we demonstrate that the duration for which a
relay stays in the receive mode is a crucial design parameter. With its
appropriate selection, our proposed methods outperform uncoded forwarding in
terms of both loss rate and the average number of transmissions by the relay.

Figure 1: A LoRa sensor network with relays.

##  II Relaying Protocols

We consider up to two relays placed in between a cluster of sensors and the
gateway (see Fig. 1). We assume slotted ALOHA transmissions by the sensors,
whereas the relays transmit over predefined time slots. Time slotting is also
assumed in [4]–[7]. No frame-level synchronization is required at the sensors,
which only need to be aware of the slot boundaries. Regardless, our protocols
can be easily adapted to pure ALOHA transmissions from sensors. We refer to
the two protocols as the single-relay protocol and the cooperative protocol,
depending on whether they use one relay or two relays. In each, the sensors
transmit messages to the gateway, and the relays try to overhear the sensors’
frames.

###  II-A The Single-Relay Protocol

In the single-relay protocol, the relay alternates between two states – a
receive window (RW) having $n_{r}$ slots and a transmit window (TW) having a
single slot. This is illustrated in Fig. 2. In the receive window, the relay
tries to overhear the sensors’ transmissions and stores the received messages.
In the transmit window, the relay transmits a coded frame containing the sum
(bitewise XOR) of the stored messages. The gateway can recover one message
from a coded frame if all other XORed messages have already been received.
Coded frames that fail to deliver a message are discarded.

Figure 2: Relay timing diagram (shown for $n_{r}\\!=\\!3$ and
$n_{s}\\!=\\!2$).

###  II-B The Cooperative Protocol

The cooperative protocol employs two relays, each having an $n_{r}$-slot
receive window and a single-slot transmit window. In addition, the relays also
have a sleep window of $n_{s}$ slots. As in the single-relay protocol, each
relay overhears sensor frames in the receive window, and transmits the sum of
the overheard messages in the following transmit window. In the sleep window,
a relay stays idle to conserve energy. The relay windows, shown in Fig. 2, are
such that $n_{r}\\!=\\!n_{s}\\!+\\!1$, and one relay’s receive window
coincides with the other’s transmit and sleep windows. This ensures that at
any given instant, there is one relay available to receive.

##  III System Model

We denote the number of sensors by $n$ and the slot length by $l_{s}$. A
sensor frame contains $b_{\mathrm{pl}}$ bytes of measurement data and a header
with a $b_{\mathrm{id}}$-byte sensor ID and a $b_{\mathrm{seq}}$-byte sequence
number. Since a coded message is the sum of multiple sensor messages, each
$b_{\mathrm{pl}}$ bytes long, the coded messages are also of length
$b_{\mathrm{pl}}$ bytes. In addition to the sum, the coded frame also contains
the identifiers (originating node IDs and sequence numbers) of the summed
messages as its payload.

We assume the interval between two consecutive transmissions from a sensor to
be exponentially distributed with parameter $\lambda$. Thus $\lambda^{-1}$ is
the mean interval between two consecutive frames, which is assumed to be much
larger than the slot length $l_{s}$. Each transmission begins at the slot
boundary, and a slot is long enough to accommodate one full frame. Such
assumptions are consistent with typical LoRa applications in which EDs send
small amounts of data infrequently.

The distances from the sensors to the gateway are independent and identically
distributed (i.i.d.) with probability density function (pdf) $f_{D}(x)$.
Similarly, the distances from the EDs to the relays are i.i.d. with pdf
$f_{D^{\prime}}(x)$. The power-fading coefficient $A$ on a link remains
constant over the reception of a frame but varies from frame to frame with pdf
$f_{A}(x)$. The cumulative distribution functions (cdfs) of $D$, $D^{\prime}$,
and $A$ are denoted by $F_{D}(x)$, $F_{D^{\prime}}(x)$, and $F_{A}(x)$.

All frames are transmitted with the same power. The received power at the
gateway from an ED at distance $D$ is $R\\!=\\!\gamma AD^{-\alpha}$, where the
constant $\gamma$ depends on antenna characteristics and transmit power. The
cdf of $R$ is

$\displaystyle F_{R}(x)$ $\displaystyle=P(\gamma
AD^{-\alpha}<x)=\int_{u}F_{A}(\gamma^{-1}u^{\alpha}x)f_{D}(u)\>\mathrm{d}{u}.$
(1)

Similarly, the received power $R^{\prime}$ at the relay has the cdf

$\displaystyle F_{R^{\prime}}(x)$
$\displaystyle=\int_{u}F_{A}(\gamma^{-1}u^{\alpha}x)f_{D^{\prime}}(u)\>\mathrm{d}{u}.$
(2)

For simplicity, the distance between the gateway and the relay(s) is assumed
to be fixed at $d_{r}$. Consequently, the received power at the gateway from
the relays has the cdf

$\displaystyle F_{R^{\prime\prime}}(x)$
$\displaystyle=F_{A}(\gamma^{-1}d_{r}^{\alpha}x).$ (3)

We denote the receiver sensitivity by $\zeta$. An incoming frame goes
undetected if its power at the receiver is below $\zeta$.

Frame duration and receiver sensitivity both depend on the spreading factor
(SF), which is an integer between 7 and 12. Larger SFs result in longer frames
and lower sensitivity. The SFs are quasi-orthogonal, therefore the
interference experienced by a frame from signals of different SFs is rather
small. We denote by $L(s,b)$ the duration of a frame with $b$ bytes payload
and SF $s$. The frame length computation is described in [9]. We assume
sensors and relays to use different SFs and neglect interference between relay
and sensor transmissions.

The LoRa modulation exhibits the capture effect. In order for a frame to be
correctly demodulated in the presence of other same-SF frames, its power must
exceed the strongest interferer’s power by at least an amount equal to the
capture threshold, which we denote by $\xi$ (in dB) [1].

##  IV Performance Analysis

Our performance metrics are the message loss rate (MLR) and the relay duty
cycle (RDC). The MLR is the probability that a sensor message fails to reach
the gateway, that is, neither the sensor-to-gateway direct link nor a relay is
able to deliver the message. The RDC is the fraction of time the relays spend
transmitting frames. With all relay frames having the same transmit power, the
RDC metric captures the energy spent by the relays in forwarding sensor
messages. Since LoRa networks operate in unlicensed frequency bands, the RDC
may also be interpreted as a measure of the interference contributed by the
relays to other devices occutpying the same band.

###  IV-A MLR for the Single-Relay Protocol

The MLR for the single-relay protocol is

$\mathrm{MLR}_{\text{single-relay}}=1-S_{\mathrm{dir}}-S_{\mathrm{rel}},$ (4)

where $S_{\mathrm{dir}}$ is probability that the sensor delivers the message
directly to the gateway (the direct delivery probability), and
$S_{\mathrm{rel}}$ is the probability that the direct transmission fails but
the relay is able to deliver the message to the gateway (the relay delivery
probability). We can express $S_{\mathrm{rel}}$ as

$S_{\mathrm{rel}}=S_{\text{s-}\mathrm{r}}S_{\mathrm{r}\text{-g}}S_{\mathrm{c}}$
(5)

where $S_{\text{s-}\mathrm{r}}$ is the probability that the sensor’s
transmission fails at the gateway but is received by the relay,
$S_{\mathrm{r}\text{-g}}$ is the probability that the relay’s coded frame is
delivered to the gateway, and $S_{\mathrm{c}}$ is the probability that of all
XORed messages in the frame, only one message is missing at the gateway.

Consider a message from an arbitrary sensor located at distance $D$ from the
gateway. We call this message the desired message. We will derive
$S_{\mathrm{dir}}$ and $S_{\mathrm{rel}}$ for the desired message in the
following subsections.

####  IV-A1 Direct Delivery Probability

Suppose that $K$ other nodes also transmitted their messages in the same slot
as the desired message. Denote the received power at the gateway in the frame
carrying the desired message by $R_{0}\\!=\\!\gamma AD^{-\alpha}$. Conditioned
on $D\\!=\\!d_{0}$, $A\\!=\\!a_{0}$, and $K\\!=\\!k$, we write the conditional
probability that this frame is received over the direct link as
$S_{\mathrm{dir}}(d_{0},a_{0},k)$. It is the probability that the frame
survives both fading and interference, and can be written as

$S_{\mathrm{dir}}(d_{0},a_{0},k)=\theta(d_{0},a_{0})S_{I}(a_{0},d_{0},k),$ (6)

where $\theta$ specifies whether the desired frame survives fading. It is
defined as $\theta(d_{0},a_{0})\\!=\\!1$ if $\gamma
a_{0}d_{0}^{-\alpha}\\!>\\!\zeta$, and $\theta(d_{0},a_{0})\\!=\\!0$
otherwise, where $\zeta$ is the receiver sensitivity. The term
$S_{I}(a_{0},d_{0},k)$ is the conditional probability that the frame survives
interference, and is equal to the probability that the received power
$r_{0}\\!=\\!\gamma a_{0}d_{0}^{-\alpha}$ exceeds the maximum of the powers in
the $k$ interferers by at least $\xi$ dB. Let $\hat{R}$ be the strongest
interferer power. Then

$\displaystyle S_{I}(a_{0},d_{0},k)$
$\displaystyle=P(r_{0}/\hat{R}>\tilde{\xi})=\left(F_{{R}}(\tilde{\xi}^{-1}\gamma
a_{0}d_{0}^{-\alpha})\right)^{k},$ (7)

where $\tilde{\xi}\\!=\\!10^{0.1\xi}$. To find the probability that there are
$k$ interferers, note that the probability that a sensor makes a transmission
in any given slot is approximately

$p_{\mathrm{tx}}\approx 1-\mathrm{e}^{-\lambda l_{s}}.$ (8)

To see this, note that the exponential($\lambda$) distribution of the inter-
arrival times result in a Poisson arrival process with rate $\lambda$ for each
sensor. The probability that a sensor produces no measurement within a slot is
$\mathrm{e}^{-\lambda l_{s}}$. Since the probability that a sensor generates
multiple measurements within a slot is negligible,
$1\\!-\\!\mathrm{e}^{-\lambda l_{s}}$ is the approximate probability that it
transmits in a given slot. Thus the number of interferers $K$ is a binomial
$(n\\!-\\!1,p_{\mathrm{tx}})$ random variable, to which we employ the Poisson
approximation

$P_{K}(k)\approx\frac{\nu^{k}\mathrm{e}^{-\nu}}{k!},\quad 0\leq k\leq
n\\!-\\!1,$ (9)

where $\nu\\!=\\!(n\\!-\\!1)p_{\mathrm{tx}}$. Now using (7) in (6) and
deconditioning on $K$, we obtain

$\displaystyle S_{\mathrm{dir}}(a_{0},d_{0})$
$\displaystyle=\sum_{k=0}^{n-1}\frac{\nu^{k}\mathrm{e}^{-\nu}}{k!}\theta(a_{0},d_{0})\left(F_{{R}}(\tilde{\xi}^{-1}\gamma
a_{0}d_{0}^{-\alpha})\right)^{k}$
$\displaystyle=\theta(a_{0},d_{0})\mathrm{e}^{-\nu}\sum_{k=0}^{n-1}\frac{\left(\nu\lambda_{R}(a_{0},d_{0})\right)^{k}}{k!},$
(10)

where

$\displaystyle\lambda_{R}(a_{0},d_{0})=F_{{R}}(\tilde{\xi}^{-1}\gamma
a_{0}d_{0}^{-\alpha}).$ (11)

Approximating the summation by the exponential function,

$\displaystyle S_{\mathrm{dir}}(a_{0},d_{0})$
$\displaystyle\approx\theta(a_{0},d_{0})\mathrm{e}^{-\nu}\mathrm{e}^{\nu\lambda_{R}(a_{0},d_{0})}$
$\displaystyle=\theta(a_{0},d_{0})\mathrm{e}^{-\nu(1-\lambda_{R}(a_{0},d_{0}))}.$
(12)

The final expression for $S_{\mathrm{dir}}$ is given by

$\displaystyle
S_{\mathrm{dir}}=\mathop{\mathbb{E}}_{(A,D)}\left[\theta(A,D)\mathrm{e}^{-\nu(1-\lambda_{R}(A,D))}\right],$
(13)

where

$\displaystyle\mathop{\mathbb{E}}_{(A,D)}[(g(A,D))]=\int_{u}\int_{a}g(a,u)f_{A}(a)f_{D}(u)\>\mathrm{d}{a}\>\mathrm{d}{u},$
(14)

####  IV-A2 Probability of Relay Delivery

A necessary – but not sufficient – condition for the delivery of a frame to
the relay is that the frame must be transmitted within the relay’s receive
window. This occurs with probability

$p_{\mathrm{rw}}=\frac{n_{r}}{n_{r}+1}.$ (15)

Now consider a frame transmitted during the relay’s receive window. Let
$S_{\mathrm{rx}}$ be the probability that it is lost at the gateway but is
received by the relay. To find $S_{\mathrm{rx}}$, we condition on $K\\!=\\!k$
interferers in the slot, $D\\!=\\!d_{0}$ and $D^{\prime}\\!=\\!d_{1}$ being
the sensor-to-gateway and sensor-to-relay separations, respectively, and
$A\\!=\\!a_{0}$ and $A^{\prime}\\!=\\!a_{1}$ being the fading coefficients on
the sensor-to-gateway and sensor-to-relay links, respectively. We then have
the following conditional reception probability

$\displaystyle S_{\mathrm{rx}}^{\mathrm{cond}}$
$\displaystyle=(1-S_{\mathrm{dir}}(a_{0},d_{0},k))S_{r}(a_{1},d_{1},k),$ (16)

where $S_{r}$ is the probability of correct reception over the sensor-to-relay
link during a slot in which the relay is receiving. Its derivation is
identical to that of $S_{\mathrm{dir}}$, except that the distance pdf is
$f_{D^{\prime}}(\cdot)$ instead of $f_{D}(\cdot)$. Therefore,

$\displaystyle S_{\mathrm{rx}}^{\mathrm{cond}}$
$\displaystyle=\left[1-\theta(a_{0},d_{0})(\lambda_{R}(a_{0},d_{0}))^{k}\right]\theta(a_{1},d_{1})(\lambda_{R^{\prime}}(a_{1},d_{1}))^{k}.$
(17)

Removing the conditioning, we obtain

$\displaystyle S_{\text{s-}\mathrm{r}}$
$\displaystyle=p_{\mathrm{rw}}\sum_{k=0}^{n-1}\frac{\nu^{k}\mathrm{e}^{-\nu}}{k!}\mathop{\mathbb{E}}_{(A,A^{\prime},D,D^{\prime})}[S_{\mathrm{rx}}^{\mathrm{cond}}],$
(18)

where,

$\displaystyle\mathop{\mathbb{E}}_{(A,A^{\prime},D,D^{\prime})}$
$\displaystyle[g(A,A^{\prime},D,D^{\prime}]$
$\displaystyle=\int_{a}\int_{a^{\prime}}\int_{u}\int_{v}g(a,a^{\prime},u,v)$
$\displaystyle
f_{A}(a)f_{A}(a^{\prime})f_{D}(u)f_{D^{\prime}}(v).\>\mathrm{d}{a}\>\mathrm{d}{a^{\prime}}\>\mathrm{d}{u}\>\mathrm{d}{v}.$

The expression above simplifies greatly if the sensor locations are such that
the loss due to fading is small, and the sensors are clustered such that each
of them has approximately the same distance $d_{0}$ from the gateway and
$d_{1}$ from the relay. In such cases, we may substitute
$\theta(a_{0},d_{0})\approx\theta(a_{1},d_{1})\approx 1$ in (17) and remove
the conditioning to obtain

$\displaystyle S_{\text{s-}\mathrm{r}}$
$\displaystyle=p_{\mathrm{rw}}\cdot\Bigg{(}\mathop{\mathbb{E}}_{A_{1}}\left[\mathrm{e}^{-\nu(1-\lambda_{R}(A_{1},d_{1}))}\right]$
$\displaystyle+\mathop{\mathbb{E}}_{(A_{0},A_{1})}\left[\mathrm{e}^{-\nu(1-\lambda_{R}(A_{0},d_{0})\lambda_{R}(A_{1},d_{1}))}\right]\Bigg{)},$
(19)

where $A_{0}$ and $A_{1}$ both have the pdf $f_{A}(x)$.

To find $S_{\mathrm{r}\text{-g}}$, the success probability on the relay-to-
gateway link, only frame loss due to fading is of concern since the relay uses
a different spreading factor than the sensors. Hence,

$S_{\mathrm{r}\text{-g}}=1-F_{A}(\gamma^{-1}\zeta\tilde{d}_{r}^{\alpha}),$
(20)

where $\tilde{d}$ is the relay-to-gateway separation.

The only parameter left to find is $S_{c}$, the probability that the gateway
is able to recover the desired message from a coded message. Since recovery is
possible if and only if all XORed messages other than the desired message have
already been delivered to the gateway, we wish to derive the probability that
in the preceding receive window, the relay received $m$ messages (other than
the desired message), each of which was also received by the gateway. Using
the usual conditioning and deconditioning, the probability that a frame
transmitted within the relay’s receive window is received by both the relay
and the gateway is

$\displaystyle
p_{\mathrm{gr}}\\!=\\!\sum_{k=0}^{n-1}\frac{\nu^{k}\mathrm{e}^{-\nu}}{k!}\mathop{\mathbb{E}}_{(A_{0},A_{1},D,D^{\prime})}[S_{\mathrm{dir}}(A_{0},D,k)S_{r}(A_{1},D^{\prime},k)].$
(21)

We can now determine $S_{c}$ as follows

$\displaystyle S_{c}=\sum_{m=0}^{n_{r}-1}{n_{r}-1\choose
m}p_{\mathrm{gr}}^{m}f^{n_{r}-m-1},$ (22)

where $f$ is the probability that no message is successfully received at the
relay over a given slot. It is given by

$\displaystyle
f=1-np_{\mathrm{tx}}\mathop{\mathbb{E}}_{(A,D^{\prime})}\left[\theta(A,D^{\prime})\mathrm{e}^{-\nu(1-\lambda_{R}(A,D^{\prime}))}\right].$
(23)

The product of $p_{\mathrm{tx}}$ and the expectation in (23) is the
probability that an arbitrary sensor transmits a message that is correctly
received by the relay. The multiplier $n$ accounts for the fact that there are
$n$ such sensors.

###  IV-B RDC of the Single-Relay Protocol

Suppose that the relay received $m\\!\geq\\!1$ messages in the receive window
preceding a certain transmit window. Then it must send a coded frame of length
$b_{\mathrm{pl}}$ bytes along with a header containing sender ids and sequence
numbers of the messages that are XORed. Assuming all these to be part of the
relay payload, the payload length in bytes is

$b_{c}(m)=b_{\mathrm{pl}}+m(b_{\mathrm{id}}+b_{\mathrm{seq}}).$ (24)

The average relay frame length is therefore

$l_{\mathrm{av}}=\sum_{m=1}^{n_{r}}{n_{r}\choose
m}(1-f)^{m}f^{n_{r}-m}L(\mathrm{sf}_{r},b_{c}(m)),$ (25)

where $\mathrm{sf}_{r}$ is the SF used by the relay. The average transmission
duration per transmit window is hence $l_{\mathrm{av}}$, and the fraction of
time the relay spends transmitting is

$\displaystyle\mathrm{RDC}_{\text{single-
relay}}=\frac{l_{\mathrm{av}}}{(n_{r}+1)l_{s}}.$ (26)

###  IV-C Analysis of the Cooperative Protocol

The MLR analysis for the cooperative protocol is the same as that of single-
relay case if we assume that the distance from the gateway to both relays are
approximately the same, and so are the pdfs of the distances from the sensors
to the two relays. We may then visualize the two relays as a single full-
duplex relay that is always receiving. The MLR for can then be analyzed using
an approach identical to the one for the single-relay protocol, except that we
must use ${p_{\mathrm{rw}}}\\!=\\!1$.

For the RDC too, the derivation is identical except that now there are two
relays and each sleeps for $n_{s}$ slots in addition to transmitting and
receiving. The resulting expression is

$\displaystyle\mathrm{RDC}_{\mathrm{cooperative}}=\frac{2l_{\mathrm{av}}}{(n_{r}+n_{s}+1)l_{s}}.$
(27)

Since $n_{r}\\!=\\!n_{s}\\!+\\!1$, the above simplifies to

$\displaystyle\mathrm{RDC}_{\mathrm{cooperative}}=\frac{l_{\mathrm{av}}}{(n_{r}+0.5)l_{s}}.$
(28)

##  V Performance Results

We first consider a network of 20 sensors, each having a 1-byte ID, and
transmitting 10-byte messages with a mean inter-frame arrival time of 17.5 s.
Each frame has a 1-byte sequence number. The sensors and relays use SF 8 and
7, respectively. The links have block Rayleigh fading. Our performance graphs
show the analytical results, which we found to be in very close agreement with
simulations.

Fig. 3 compares the MLR for the single-relay protocol with three baselines:
(1) The no-relay protocol does not employ relays. (2) The immediate-forwarding
protocol uses a single relay that forwards every received message immediately
in the following slot. This is essentially the protocol of [6]. (3) The
uncoded-forwarding protocol is same as our proposed single-relay protocol
except that the relay forwards the individual messages received during the
previous receive window instead of transmitting their sum. If all messages do
not fit within the transmit window, some messages are discarded at random.
This approach is identical to that of [3] with one relay. To distinguish our
proposed single-relay protocol from the benchmarks that also employ one relay,
we occasionally refer to it as the single-relay sum-and-forward protocol.

Figure 3: MLR comparison of various protocols.

It is clear from Fig. 3 that any form of relaying provides better MLR than the
system without a relay, even though each relay-based protocol compared here
has only a single relay against 20 sensors. The no-relay and immediate-
forwarding protocols do not have the concept of $n_{r}$, so we plot their MLRs
as straight lines. The uncoded-forwarding protocol – although better than no-
relay – performs poorly relative to the other relay protocols, especially for
larger receive window sizes. A long receive window results in more received
messages than the relay can send in the transmit window slot, resulting in
dropped massages. By contrast, our proposed single-relay protocol performs
significantly better, despite having the same transmit window duration. This
is the benefit of the sum-and-forward approach, which only needs a transmit
window long enough to accommodate one XORed message (whose length is the same
as any of the individual information messages) and the necessary header
information. Like uncoded forwarding, the MLR of the proposed protocol is
sensitive to the receive window size. Having too small a receive window often
leads to the relay not being able to receive any frames during a receive
window, and consequently having to sit idle in the next transmit window.
Having too large a receive window means that the relay’s transmitted frames
are XORs of a large number of messages, and the probability that the gateway
is missing only one of them is small. The optimal receive window for this
scenario is found to be of 11 slots.

Figure 4: RDC comparison of the proposed single-relay protocol and immediate
forwarding.

One might conclude from Fig. 3 that the proposed sum-and-forward approach
hardly provides any benefit over immediate forwarding. But Fig. 4 demonstrates
the advantage of the sum-and-forward approach by comparing the respective
RDCs, which are significantly higher for immediate forwarding. For example,
with $n_{r}\\!=\\!11$, immediate forwarding requires 42% more relay
transmission time than our proposed technique. Thus, single-relay sum-and-
forward would ensure significantly longer operation of battery-powered relays
and also contribute less interference to the network.

(a) MLR.

(b) RDC.

Figure 5: Comparison of the proposed protocols.

Fig. 5(a) compares the MLR for the two proposed protocols. The cooperative
protocol is able to significantly reduce the MLR relative to the single-relay
approach, especially for short receive windows. For $n_{r}\\!=\\!1$, the
cooperative protocol provides the lowest MLR, which is 27% lower than than
that of the single relay protocol. Since $n_{r}\\!=\\!1$ means that a relay
has a maximum of one message to send per transmit window, there is no scope
for coding in this case. By contrast, if only one relay is available, it is
beneficial to forward sums of messages with an appropriately chosen receive
window size.

In addition to the MLR, the receive window size also determines the delay
characteristic of the protocol, since it specifies how long the gateway has to
wait before possibly recovering a lost message via a relay. With an optimal
receive window comprising a single slot, it is clear the cooperative protocol
provides better delay performance than single-relay sum-and-forward at the
cost of an extra relay.

Fig. 5(b) compares the RDCs of the proposed protocols. The RDC for the
cooperative protocol exceeds that of the single-relay protocol for all receive
window sizes. However, for the cooperative protocol, the RDC includes
transmissions by both relays; the RDC per relay in this case is much smaller
than the RDC for the single-relay protocol.

Fig. 6(a) plots the MLR against the number of sensors. For the proposed
protocols, the minimum achievable MLR (with respect to $n_{r}$) is plotted.
For the no-relay and immediate-forwarding protocols, no such minimization is
needed since the MLR is fixed for a given sensor count. For each protocol, the
MLR increases with the number of sensors due to increased interference. The
cooperative protocol provides the lowest MLR, while the immediate-forwarding
and single-relay sum-and-forward protocols provide comparable results.
However, as seen from Fig. 6(b), immediate forwarding requires more
transmission time at the relay (i.e., higher RDC) compared to single-relay
sum-and-forward. For 40 sensors, single-relay sum-and-forward requires about
55% less relay transmission time compared to immediate forwarding. The RDC for
the cooperative protocol is also significantly lower (33% with 40 sensors)
than immediate forwarding, despite using two relays.

The $n_{r}$ value that provides the lowest MLR for the proposed protocols are
shown in Fig. 7. We observe that $n_{r}\\!=\\!1$ is the optimal choice for the
cooperative protocol irrespective of the number of sensors. By contrast, the
optimal receive window size for the former decreases with increasing sensor
count. With large number of nodes, loss rates are high on the links. For
coding to be beneficial, gateway must have exactly one missing message out of
all XORed messages in a coded frame. With greater frame loss rates, coded
frames containing fewer messages are more likely to be successful. A shorter
receive window leads to more such coded frames.

(a) MLR.

(b) RDC.

Figure 6: Protocol performance vs. number of sensors. Figure 7: Optimal
receive window size vs. number of sensors.

##  VI Conclusion

We proposed and evaluated two low-complexity relaying protocols for LoRa
sensor networks. Using either one or two relays, the protocols improve the
message delivery performance of the network. The proposed transmission schemes
achieve similar (with a single relay) or better (with two relays) loss
performance compared to conventional uncoded forwarding, while requiring fewer
transmissions by the relay(s).

##  VII Future Work

In this paper, we restricted attention to a maximum of two relays and 40
sensor nodes. As the sensor count increases, it may be necessary to employ
more relays. By partitioning the sensors to multiple clusters, and assigning
one or two relays to each cluster, our proposed protocols can be scaled to
large networks. Devising such approaches is left as future work.

Adaptation of the protocols to pure ALOHA channel access by sensors is also of
practical importance. With pure ALOHA, the receive window of a relay may end
while the reception of a sensor frame is still in progress. Analysis of
relaying with pure ALOHA transmissions by the sensors is provided in [3] under
the assumption of periodic sensor transmissions and uncoded relaying. Analysis
of coded relaying with possibly aperiodic sensor frames is left as future
work.

## References

  * [1] O. Georgiou and U. Raza, “Low Power Wide Area Network Analysis: Can LoRa Scale?,” in IEEE Wireless Communications Letters, vol. 6, no. 2, pp. 162–165, April 2017.
  * [2] A. Pop, U. Raza, P. Kulkarni, and M. Sooriyabandara, “Does bidirectional traffic do more harm than good in LoRaWAN based LPWA networks?," in Proc. IEEE GLOBECOM, December 2017, pp. 1–6.
  * [3] S. S. Borkotoky, U. Schilcher, and C. Bettstetter, “Cooperative relaying in LoRa sensor networks,” in Proc. IEEE Global Communications Conference, Waikoloa, HI, USA, December 2019.
  * [4] C.-H. Liao, G. Zhu, D. Kuwabara, M. Suzuki, and H. Morikawa, “Multi-hop LoRa networks enabled by concurrent transmission,” IEEE Access, vol. 5, pp. 21430–21446, 2017.
  * [5] M. Bor, J. E. Vidler, and U. Roedig, “LoRa for the Internet of Things,” in Proc. Intern. Conf. on Embedded Wireless Systems and Networks (EWSN), Graz, Austria, Feb. 2016, pp. 361–366.
  * [6] P. Branch and T. Cricenti, “A LoRa based wireless relay network for actuator data,” Proc. ICOIN, 2020, pp. 190–195.
  * [7] H. P. Tran, W. Jung, T. Yoon, D. Yoo, and H. Oh, “A Two-Hop Real-Time LoRa Protocol for Industrial Monitoring and Control Systems,” in IEEE Access, vol. 8, pp. 126239–126252, 2020.
  * [8] D. Tanjung, S. Byeon, D. H. Kim, and J. Deok Kim, “OODC: An Opportunistic and On-Demand Forwarding Mechanism for LPWA Networks,” Proc. ICOIN 2020\.
  * [9] SX1272/3/6/7/8: LoRa Modem Designer’s Guide AN1200.13, Semtech Corporation, July 2013.
  * [10] T. Bouguera, J. F. Diouris, J. J. Chaillout, R. Jaouadi, and G. Andrieux, “Energy Consumption Model for Sensor Nodes Based on LoRa and LoRaWAN,” Sensors, vol. 18, 2018.


# Molecules in Environments: Towards Systematic Quantum Embedding of Electrons
and Drude Oscillators

Matej Ditte<EMAIL_ADDRESS>Department of Physics and Materials Science,
University of Luxembourg, L-1511 Luxembourg City, Luxembourg  Matteo Barborini
Department of Physics and Materials Science, University of Luxembourg, L-1511
Luxembourg City, Luxembourg  Leonardo Medrano Sandonas  Department of Physics
and Materials Science, University of Luxembourg, L-1511 Luxembourg City,
Luxembourg  Alexandre Tkatchenko<EMAIL_ADDRESS>Department of
Physics and Materials Science, University of Luxembourg, L-1511 Luxembourg
City, Luxembourg

(March 28, 2023)

###### Abstract

We develop a quantum embedding method that enables accurate and efficient
treatment of interactions between molecules and an environment, while
explicitly including many-body correlations. The molecule is composed of
classical nuclei and quantum electrons, whereas the environment is modeled via
charged quantum harmonic oscillators. We construct a general Hamiltonian and
introduce a variational ansatz for the correlated ground state of the fully
interacting molecule/environment system. This wavefunction is optimized via
variational Monte Carlo and the ground state energy is subsequently estimated
through diffusion Monte Carlo. The proposed scheme allows an explicit many-
body treatment of electrostatic, polarization, and dispersion interactions
between the molecule and the environment. We study solvation energies and
excitation energies of benzene derivatives, obtaining excellent agreement with
explicit ab initio calculations and experiment.

The challenges of computational electronic-structure methods have evolved from
the precise description of short-range interactions involving few electrons to
the accurate modeling of longer-range many-electron effects for molecules and
materials Sun and Chan (2016); Jones  _et al._ (2020); Hermann  _et al._
(2017). This shift has been enabled by the enormous progress in hardware
capabilities combined with the development of new effective physical models
and numerical techniques.

In particular, embedding approaches, such as the hybrid QM/MM Fersht (2013);
Khare  _et al._ (2007); Torras  _et al._ (2009); Zimmerman  _et al._ (2011);
Nasluzov  _et al._ (2001, 2003); Brändle and Sauer (1998); Guareschi  _et al._
(2016); Senn and Thiel (2009); Acevedo and Jorgensen (2010); Huang  _et al._
(2017); Alves  _et al._ (2015); Silva  _et al._ (2015); Bao  _et al._ (2019);
Marín  _et al._ (2019); Robertazzi and Platts (2006); Gkionis and Platts
(2009); Kupfer  _et al._ (2014); Osoegawa  _et al._ (2019); Loco and Piquemal
(2022); Cho and Lester (2010) or QM/QM strategies Kotliar  _et al._ (2006);
Lan  _et al._ (2015); Knizia and Chan (2012); Wesolowski  _et al._ (2015);
Libisch  _et al._ (2014); Manby  _et al._ (2012); Katin  _et al._ (2017);
Friedrich  _et al._ (2007); Seth  _et al._ (2002); Vogiatzis  _et al._ (2015);
Huo  _et al._ (2016); Lee  _et al._ (2019); Coughtrie  _et al._ (2018);
Neugebauer  _et al._ (2005); Jacob and Visscher (2006); Exner  _et al._
(2012); Exner and Mezey (2004); Bockstedte  _et al._ (2018); Lau  _et al._
(2021); Schäfer  _et al._ (2021); Ma  _et al._ (2021); Petocchi  _et al._
(2020); Yeh  _et al._ (2021); Nowadnick  _et al._ (2015); Chen  _et al._
(2022); Nusspickel and Booth (2022), based on the partitioning of the system
of interest in different subspaces, each treated with a different level of
theory and thus of accuracy Sun and Chan (2016); Jones  _et al._ (2020), have
been devised to study excitation energies, solvation energies or charge
transfer between fragments in a complex environment Bovi  _et al._ (2013);
Torras  _et al._ (2009); Guareschi  _et al._ (2016); Marín  _et al._ (2019);
Wesolowski  _et al._ (2015); Neugebauer  _et al._ (2005); Schäfer  _et al._
(2021); Cho and Lester (2010).

Even though embedding approaches can significantly improve the computational
efficiency, this comes at the cost of approximating collective interactions
within and between the subsystems, thus excluding important many-body effects
in the overall description Sun and Chan (2016); Jones  _et al._ (2020).
Furthermore, systems dominated by non-covalent interactions Benali  _et al._
(2014); Al-Hamdani  _et al._ (2021) or by static correlations Reta Mañeru  _et
al._ (2014); Boyn and Mazziotti (2021); Bao  _et al._ (2016); Varsano  _et
al._ (2017); Barborini  _et al._ (2022); Pozun  _et al._ (2013); Barborini and
Coccia (2015); Abe (2013) remain challenging, since they require high-level
_ab initio_ methods, such as coupled cluster (CC) Bartlett and Musiał (2007)
or quantum Monte Carlo (QMC) Foulkes  _et al._ (2001); Kalos and Whitlock
(2008); Becca and Sorella (2017), for which systematic convergence can be
guaranteed only for relatively small systems Benali  _et al._ (2014); Al-
Hamdani  _et al._ (2021).

In this Letter, we propose a novel approach, based on the embedding of
electronic systems (atoms/molecules) within an environment of charged quantum
harmonic (Drude) oscillators, which can reproduce the quantum fluctuations
responsible for the non-covalent interactions and for the long-range response
of realistic molecules and materials (see Fig. 1) Jones  _et al._ (2009);
Tkatchenko  _et al._ (2012); Jones  _et al._ (2013a, b, c); Odbadrakh  _et
al._ (2015); Odbadrakh and Jordan (2016); Cipcigan  _et al._ (2019);
Tkatchenko and Scheffler (2009); Vaccarelli  _et al._ (2021); Poier  _et al._
(2022); Karimpour  _et al._ (2022a, b). Hence, our approach explicitly
describes many-body correlations between the subsystem and the environment.

Although harmonic oscillators have been used in classical polarizable force
fields Lemkul  _et al._ (2016); Bui  _et al._ (2023); Sommerfeld and Jordan
(2005), the novelty of our method lies in the quantum description of the full
system of atoms and quantum Drude oscillators (QDOs) through one unified
Hamiltonian for which we construct a ground state depending on all degrees of
freedom.

Figure 1: Schematic representation of the embedding model for a single atom
and two Drude oscillators. Strings indicate the harmonic interactions between
the drudon and their corresponding centers, the dashed lines represent the
Coulomb interactions between all other pairs of charged particles (for
simplicity, only selected interactions are explicitly shown).

Furthermore, the proposed approach can be used with arbitrary ab initio
methods, including CC Bartlett and Musiał (2007) and QMC Foulkes  _et al._
(2001); Kalos and Whitlock (2008); Becca and Sorella (2017).

Within the Born-Oppenheimer approximation, each QDO consists of a pair of
charged particles corresponding respectively to the center of the oscillator,
with charge $q$ and position $\textbf{R}^{O}$ and to a distinguishable quantum
particle called drudon of charge -q, mass $\mu$ and position $\textbf{r}^{d}$.
The drudon interacts with its center via a quadratic potential
$V_{\text{QDO}}\left(r\right)=\frac{1}{2}\mu\omega^{2}r^{2}$ that depends on
the mass $\mu$ and on the frequency $\omega$ that determine the slope of the
quadratic well, and on their relative Euclidean distance
$r=|\textbf{r}^{d}-\textbf{R}^{O}|$. For systems containing more QDOs, all
pairs of charged particles (except the drudon and its own center) interact via
the Coulomb potential. In such a model the non-covalent interactions, such as
dispersion, arise from the quantum fluctuations of the charge densities that
are obtained through the coupling of the quantum states of each QDO.
Naturally, each QDO depends on a set of well defined parameters
$\\{\textbf{R}_{i}^{O},q_{i},\mu_{i},\omega_{i}\\}_{i=1}^{N_{d}}$ that are
obtained by reproducing the leading order polarizabilities and dispersion
coefficients of real matter Jones  _et al._ (2013c) from accurate ab initio or
experimental data Margoliash and Meath (1978); Tkatchenko and Scheffler
(2009); Tkatchenko  _et al._ (2012); Vaccarelli  _et al._ (2021); Góger  _et
al._ (2023).

Once the parameters have been chosen, the analytical expression of the
Hamiltonian of a system of $N_{d}$ interacting QDOs, the drudonic Hamiltonian,
is written in atomic units as

$\hat{H}^{d}=\sum_{i=1}^{N_{d}}\hat{h}_{i}^{d}\left(\textbf{r}_{i}^{d}\right)+\sum_{i=1}^{N_{d}}\sum_{j>i}^{N_{d}}\frac{q_{i}q_{j}}{\left|\textbf{R}_{i}^{O}-\textbf{R}_{j}^{O}\right|}+\\\
+\sum_{i=1}^{N_{d}}\sum_{j>i}^{N_{d}}\frac{q_{i}q_{j}}{\left|\textbf{r}_{i}^{d}-\textbf{r}_{j}^{d}\right|},$
(1)

where
$\bar{\textbf{R}}^{O}=\left(\textbf{R}_{1}^{O}\ldots\textbf{R}_{N_{d}}^{O}\right)$
and
$\bar{\textbf{r}}^{d}=\left(\textbf{r}_{1}^{d}\ldots\textbf{r}_{N_{d}}^{d}\right)$
are the vectors of the $3N_{d}$ coordinates respectively of the oscillators’
centers and of the drudons, and
$\hat{h}_{i}^{d}\left(\textbf{r}_{i}^{d}\right)$ are the single-body operators

$\hat{h}_{i}^{d}\left(\textbf{r}_{i}^{d}\right)=-\frac{1}{2\mu_{i}}\nabla^{2}_{\textbf{r}_{i}^{d}}+\frac{1}{2}\mu_{i}\omega_{i}^{2}\left|\textbf{r}_{i}^{d}-\textbf{R}_{i}^{O}\right|^{2}-\\\
-\sum_{j\neq
i}^{N_{d}}\frac{q_{i}q_{j}}{\left|\textbf{r}_{i}^{d}-\textbf{R}_{j}^{O}\right|},$
(2)

that include the kinetic energy of the $i$th drudon and its interaction
energies with the oscillators’ centers.

This drudonic Hamiltonian can be used to model the embedding potential that
acts on a molecular system containing $N_{n}$ atomic nuclei and $N_{e}$
electrons, described respectively by the coordinate vectors
$\bar{\textbf{R}}^{n}$ and $\bar{\textbf{r}}^{e}$. The total Hamiltonian of
electrons and QDOs (El-QDO) takes the form

$\hat{H}^{tot}=\hat{H}^{e}+\hat{H}^{d}+\hat{V}_{int}^{d-e},$ (3)

where $\hat{H}^{e}$ is the standard electronic Hamiltonian describing the
interaction between the electrons and the atomic nuclei (that are parameters
of the system) and $\hat{V}_{int}^{d-e}$ is the interaction between the QDOs
and the atoms

$\hat{V}_{int}^{d-e}=\sum_{i=1}^{N_{e}}\sum_{j=1}^{N_{d}}\left(\frac{q_{j}}{\left|\textbf{r}_{i}^{e}-\textbf{r}_{j}^{d}\right|}-\frac{q_{j}}{\left|\textbf{r}_{i}^{e}-\textbf{R}_{j}^{O}\right|}\right)+\\\
+\sum_{i=1}^{N_{n}}\sum_{j=1}^{N_{d}}\left(\frac{Z_{i}q_{j}}{\left|\textbf{R}_{i}^{n}-\textbf{R}_{j}^{O}\right|}-\frac{Z_{i}q_{j}}{\left|\textbf{R}_{i}^{n}-\textbf{r}_{j}^{d}\right|}\right).$
(4)

In order to describe molecular systems with static multipole moments of the
charge distribution (such as the water molecule) it is necessary to introduce
additional point charges in the drudonic Hamiltonian that modify the external
potential felt by the drudons and by the electrons. For example in the QDO
model of water introduced by Jones et al. Cipcigan  _et al._ (2019) and used
in this work, the authors include in the model three additional point charges
(centered on the two hydrogen atoms and near the oxygen atom), see
Supplemental Material (SM) for details not . These point charges do not
interact with the QDOs of the corresponding molecule but only with the QDOs of
the other water molecules, and in our case also with the electronic subsystem.
A schematic representation of the interactions that are present in the total
Hamiltonian is shown in Fig. 1. As explained below, the QDO model will be
additionally corrected by an exchange-repulsion energy term that dominates the
short-range interaction region Sadhukhan and Manby (2016); Cipcigan  _et al._
(2019); Vaccarelli  _et al._ (2021).

The ground state of the many-body Hamiltonian defined in Eq. (3) is
constructed through an appropriate variational ansatz including the explicit
many-body correlations between all quantum particles. In order to integrate
and optimize such an ansatz over the set of electronic and drudonic
coordinates we choose to work in the framework of QMC methods Foulkes  _et
al._ (2001); Austin  _et al._ (2012); Becca and Sorella (2017), as implemented
in the QMeCha package Barborini , although the approach can be generalized to
other quantum-chemistry methods. QMC methods are stochastic integration
techniques that are able to compute and minimize the energy functional over a
chosen trial wave function (see SM not ). In particular we apply the
variational Monte Carlo (VMC) method Foulkes  _et al._ (2001); Austin  _et
al._ (2012); Becca and Sorella (2017), optimizing the wave function with the
stochastic reconfiguration algorithm Sorella (2001, 2005). Subsequently, to
improve our energy estimations we use the diffusion Monte Carlo (DMC)
projection method Foulkes  _et al._ (2001); Austin  _et al._ (2012); Becca and
Sorella (2017) that, within the fixed-node approximation used to solve the
electronic sign problem, selects the ground state component of the approximate
nodal-surface given by the optimized wave function (see SM not ). Clearly, in
the case of pure QDO systems in which the wave function is always positive and
there is no nodal surface, the DMC algorithm converges to the exact solution.

We proceed to build the trial wave function for the system’s ground state by
expressing it as the product of three contributions

$\Psi_{tot}=\Psi_{e}(\bar{\textbf{r}}^{e})\Psi_{d}(\bar{\textbf{r}}^{d})\mathcal{J}_{e-d}(\bar{\textbf{r}}^{e},\bar{\textbf{r}}^{d}),$
(5)

which are respectively the electronic ansatz $\Psi_{e}(\bar{\textbf{r}}^{e})$,
the drudonic one $\Psi_{d}(\bar{\textbf{r}}_{d})$ and a positive-definite
function $\mathcal{J}_{d-e}(\bar{\textbf{r}}^{e},\bar{\textbf{r}}^{d})$, which
contains the correlation effects between the electronic and drudonic systems..
In this work, the pure electronic part $\Psi_{e}(\bar{\textbf{r}}^{e})$ is
constructed through a Slater determinant of molecular orbitals formed as
linear combinations of Gaussian basis sets, multiplied by an electronic
Jastrow factor used to explicitly describe the many-body correlation effects
between electrons and nuclei (see SM not ).

The approximate form of the drudonic wave function
$\Psi_{d}(\bar{\textbf{r}}_{d})$ is chosen to match the exact solution of a
system of QDOs interacting through the dipole potential, which can be used to
approximate the Coulomb potential at large distances as in the many-body
dispersion MBD method Tkatchenko  _et al._ (2012); Ambrosetti  _et al._
(2014). For this reason, it is written as the exponential function of a
vector-matrix-vector product

$\Psi_{d}(\bar{\textbf{r}}_{d})=\exp\left[-\bar{\textbf{r}}_{dO}^{\top}\textbf{A}\bar{\textbf{r}}_{dO}\right],$
(6)

where $\bar{\textbf{r}}_{dO}=\bar{\textbf{r}}^{d}-\bar{\textbf{R}}^{O}$ is the
$3N_{d}$ dimensional vector containing the distances between each drudon and
its corresponding center, and A is a (symmetric) coupling matrix containing
$3N_{d}\left(3N_{d}+1\right)/2$ independent parameters. The coupling matrix
explicitly correlates the drudons subject to the field of their charged
centers.

Finally, the electron-drudon factor
$\mathcal{J}_{e-d}(\bar{\textbf{r}}^{e},\bar{\textbf{r}}^{d})$ represents the
coupling between the two types of quantum systems. Again we hypothesize that
the interactions between the QDO environment and molecular one are dominated
by dipole interactions, thus the ansatz will resemble the form in Eq. (6)

$\mathcal{J}_{e-d}(\bar{\textbf{r}}^{e},\bar{\textbf{r}}^{d})=\exp\left[\textbf{d}^{\top}\textbf{B}\bar{\textbf{r}}_{dO}\right],$
(7)

coupling the distance vector $\bar{\textbf{r}}_{dO}$ of the QDO’s with the
dipole moment of the molecular part

$\textbf{d}=\sum_{i=1}^{N_{n}}Z_{i}\textbf{R}_{i}^{n}-\sum_{i=1}^{N_{e}}\textbf{r}_{i}^{e},$
(8)

through the rectangular matrix B containing a set of $3\times 3N_{d}$ free
parameters. We remark that the restriction of the QDO wave function to the
dipolar approximation is not a significant limitation, since the variational
optimization of the wave function parameters is followed by a stochastic DMC
calculation of the many-body energy for both electrons and drudons.

Figure 2: Binding energy curve for the Ar dimer obtained using the QDO-QDO
model and electrons-QDO (El-QDO) embedding approach (schematically shown in
the figure) at the DMC level of theory, with the exponential fit of the short-
range repulsion not . The results are compared to the CCSD(T) Patkowski  _et
al._ (2005) reference curve.

Figure 3: Binding energy curve for the water dimer as a function of the
oxygen-oxygen distance in the two energetically lowest geometric
configurations (C1 and C2) Metz  _et al._ (2019) obtained using the QDO-QDO
and El-QDO (with QDO on the right, approximating the acceptor) models at the
DMC level of theory not . The results are compared to the CCSD(T) Metz  _et
al._ (2019) reference curve.

First, we apply our approach to the Ar Patkowski  _et al._ (2005) and the
water dimers Metz  _et al._ (2019). While these systems are small, when taken
together they are challenging benchmark examples, requiring an accurate
description of exchange, dispersion, polarization, and hydrogen bonding. For
both Ar and water, a single QDO represents the response of all 8 valence
electrons. Accordingly, the QDO model allows a significant reduction in the
electronic degrees of freedom. We remark that a bosonic QDO model can be
extended to describe exchange Fedorov  _et al._ (2018); Vaccarelli  _et al._
(2021), but this would require a wavefunction with mixed quantum statistics
for electrons and QDOs. Instead, in our current embedding scheme we describe
exchange via parametrized short-range exponentials (see SM not ).

In Ar2 Patkowski  _et al._ (2005), the interaction energy comes from an
interplay between correlation in the long range and exchange in the short
range. The attractive dispersion energy is correctly described by the QDO-QDO
model and essentially coincides with its counterpart computed from electronic
symmetry-adapted perturbation theory (SAPT) Shirkov and Sladek (2017) (see SM
not ). Noticeably, the electrons-QDO (El-QDO) embedding model reproduces the
full interaction curve upon subtracting the exchange contributions from SAPT,
highlighting the ability of our method to capture electrostatic and
polarization terms. Together with the corresponding repulsive part, both
models match the reference CCSD(T) binding curve Patkowski  _et al._ (2005),
as shown in Fig. 3.

For the water dimer, where intermolecular interactions are strongly
directional, we apply the El-QDO embedding model to study the dissociation
energies of its two lowest-energy geometric configurations, C1 and C2, along
the O-O axis (Fig. 3). Here we use the extended QDO model of Martyna et al.
Cipcigan  _et al._ (2019) that contains three point charges and a QDO located
at the center-of-mass, substituting the acceptor water molecule (always on the
right in Fig. 3). From Fig. 3 one observes that both C1 and C2 binding curves
obtained with the QDO-QDO and El-QDO models are in excellent agreement with
the CCSD(T) references Metz  _et al._ (2019). The models are thus able to
describe the different interactions when varying the molecular orientation.
Moreover, the comparison of the bare El-QDO interaction curve of the C1
conformer with the SAPT decomposition from Ref. 95 (see SM not ), shows that
our model is again close to the full binding energy curve, upon subtracting
the repulsive exchange contributions, thus correctly capturing both dispersion
and polarization effects. From the atomic and molecular examples of Ar and
water dimers, we conclude that our embedding approach is able to accurately
describe the manifold of interactions between an electronic subsystem and an
embedding environment, in which the latter is described through a combination
of QDOs and point charges.

Among the most intriguing applications of embedding approaches are the effects
of solvent on the binding energies between molecules and on their electronic
excitations DiStasio  _et al._ (2007); Carter-Fenk and Herbert (2020). Here,
we study the binding in T-shaped benzene dimer, prototype for C-H$-\pi$
interactions DiStasio  _et al._ (2007); Carter-Fenk and Herbert (2020), and
the Triplet-Singlet excitation energy (T-S) of ortho-benzyne, which is the
non-diradical singlet conformer of benzyne whose excitation energy requires a
proper description of dynamical correlation as well as solvent polarization
Shee  _et al._ (2019). These two systems are studied in cages of water from 4
to 50 molecules (see Fig. 4) with the aim of probing collective polarization
and dispersion effects. In order to avoid empirical treatment of exchange in
these many-body systems, we impose a minimal distance of $\sim$ 3 Å between
the solute and water.

First we compute the solvation energies,
$E_{\text{Solv}}\left(\text{X}\right)=E\left(\text{X}_{\text{Cage}}\right)-E\left(\text{X}_{\text{Vacuum}}\right)-E_{\text{Cage}}$,
where X is the T-shaped benzene dimer or its monomers (Table 3), in a cage
composed of 50 water molecules (structures c) and d) in Fig. 4), or the ortho-
benzyne in its singlet (S) and triplet (T) spin states (Table 3) in a cage of
4 and 30 water molecules (structures a) and b) in Fig. 4). For all systems the
El-QDO DMC results are compared with state-of-the-art DFT calculations using
the PBE0 functional Adamo and Barone (1999) augmented with the Tkatchenko-
Scheffler pairwise dispersion energy (PBE0+TS) Tkatchenko and Scheffler (2009)
or the many-body dispersion method (PBE0+MBD) Tkatchenko  _et al._ (2012);
Ambrosetti  _et al._ (2014); not .

Figure 4: Relative runtimes (see SM not ) of water-embedded El-QDO DMC
calculations with respect to the molecular energies computed with DMC in
vacuum. Runtime for electronic DMC is also shown for ortho-benzyne in 4W cage.
Inset figures a) to d) correspond to the molecular structures Sirianni  _et
al._ (2022); Shee  _et al._ (2019). Table 1: Solvation energies (in kcal/mol)
of the T-shaped benzene dimer and monomers in the 50W water cage Sirianni  _et
al._ (2022) (Molecular structures c) and d) in Fig. 4).

| PBE0 | PBE0+TS | PBE0+MBD | El-QDO DMC
---|---|---|---|---
Monomer1 | -1.5 | -4.1 | -3.8 | -4.0(1)
Monomer2 | -3.4 | -6.6 | -6.1 | -6.4(1)
Dimer | -4.9 | -10.9 | -10.0 | -10.3(2)

Table 2: Solvation energies (in kcal/mol) of ortho-benzyne Shee  _et al._
(2019) in its singlet (S) and triplet (T) spin states in cages of four (4W)
and thirty (30W) water molecules (Molecular structures a) and b) in Fig. 4).

| PBE0 | PBE0+TS | PBE0+MBD | DMC | El-QDO DMC
---|---|---|---|---|---
S (4W) | -1.01 | -1.50 | -1.43 | -1.5(1) | -1.39(5)
T (4W) | -0.77 | -1.24 | -1.18 | -1.0(1) | -1.12(5)
S (30W) | -2.82 | -6.78 | -6.11 | - | -6.20(7)
T (30W) | -2.37 | -6.41 | -5.71 | - | -5.73(7)

Table 3: Singlet-Triplet adiabatic excitation energies (S-T) of the ortho-
benzyne Shee  _et al._ (2019) in vacuum (V), 4-water cage (4W) and 30-water
cage (30W) (in kcal/mol), corresponding to the molecular structures a) and b)
in Fig. 4. Experimental value from Ultraviolet Photoelectron Spectroscopy is
estimated to be 37.5(3) Wenthold  _et al._ (1998).222DMC and El-QDO DMC values
are corrected with respect to the pseudopotential error shown in the SM not

| PBE0 | PBE0+TS | PBE0+MBD | DMC | El-QDO DMC
---|---|---|---|---|---
S-T (V) | 28.23 | 28.23 | 28.26 | 37.23(5) | -
S-T (4W) | 28.46 | 28.48 | 28.51 | 37.7(1) | 37.50(5)
S-T (30W) | 28.68 | 28.61 | 28.66 | - | 37.70(8)

From the results reported in Table 3 and Table 3 we conclude that El-QDO DMC
solvation energies are in excellent agreement with fully electronic DMC
calculations (for the 4W cage for which they are feasible) and with PBE0+MBD
results. It is known that PBE0+MBD binding energies are often quantitatively
accurate for molecular systems, especially when compared to electronic DMC
calculations Al-Hamdani  _et al._ (2021). In fact, we note that many-body
effects are particularly evident for ortho-benzyne, where the PBE0+TS method
overestimates the energies when compared to the El-QDO DMC and PBE0+MBD
methods, which both include, in a different manner, the many-body correlation
effects of the environment on the solute.

The key role of proper quantum embedding achieved with the El-QDO DMC approach
can be observed in the singlet-triplet (S-T) adiabatic excitation energy of
ortho-benzyne (Table 2). Due to delocalization error, DFT severely
underestimates the S-T gap Shee  _et al._ (2019), while electronic DMC yields
accurate results Zhou and Wang (2019) in vacuum in excellent agreement with
the experimental value of 37.5(3) kcal/mol Wenthold  _et al._ (1998). For the
smallest cage of four water molecules, the excitation energies computed with
electronic DMC and the El-QDO DMC approach are in excellent agreement between
each other. Yet, the error bar of El-QDO calculations is much lower compared
to electronic DMC computations for an equivalent sampling trajectory.
Moreover, as Table 2 shows, when increasing the size of the water cage, the
excitation energy changes by 0.2(1) kcal/mol. This effect, captured by the El-
QDO DMC embedding method arises from the higher polarization of the T state
with respect to the S state of benzyne. We expect that future extensions of
our method will allow more realistic embedding geometries where water
molecules can come closer to the solute, yielding more substantial solvation
effects. Remarkably, the computational cost of the El-QDO method is
essentially independent on the number of QDOs in the embedding region (see
Fig. 4) with respect to the cost of the molecule in vacuum. For the 4-water
cage, the cost increases by only $\sim$1%, and by $\sim$16% for the 30-water
cage. This stems from the lack of nodal structure in the QDO wavefunctions and
low noise in their DMC sampling with respect to the electronic wavefunctions.

Our study illustrates the need for fully quantum embedding methods when
treating binding and excitations energies of solvated systems, and introduces
a computationally efficient approach. The El-QDO embedding method is a firm
first step in this direction as we demonstrated on the example of prototypical
dispersion and hydrogen-bonded dimers and for the solvation of benzene
derivates. A future challenge includes the incorporation of a more general
representation of the QDO environment and exchange effects, which could be
achieved through a synergy with machine-learning methods Unke  _et al._
(2021a); Stöhr  _et al._ (2020); Unke  _et al._ (2021b); Hermann  _et al._
(2020).

###### Acknowledgements.

MD and MB thank Jorge Alfonso Charry Martinez for participating in the
development of the QMeCha code. MB acknowledges financial support from the
Luxembourg National Research Fund (INTER/DFG/18/12944860). The DFT
calculations presented in this paper were carried out using the HPC facilities
of the University of Luxembourg Varrette  _et al._ (2014) (see hpc.uni.lu). An
award of computer time was provided by the Innovative and Novel Computational
Impact on Theory and Experiment (INCITE) program. This research used resources
of the Argonne Leadership Computing Facility, which is a DOE Office of Science
User Facility supported under Contract DE-AC02-06CH11357.

## References

  * Sun and Chan (2016) Q. Sun and G. K.-L. Chan, Acc. Chem. Res. 49, 2705 (2016).
  * Jones  _et al._ (2020) L. O. Jones, M. A. Mosquera, G. C. Schatz, and M. A. Ratner, J. Am. Chem. Soc. 142, 3281 (2020).
  * Hermann  _et al._ (2017) J. Hermann, R. A. DiStasio, and A. Tkatchenko, Chem. Rev. 117, 4714 (2017).
  * Fersht (2013) A. R. Fersht, Proc. Nat. Acad. Sci. 110, 19656 (2013).
  * Khare  _et al._ (2007) R. Khare, S. L. Mielke, J. T. Paci, S. Zhang, R. Ballarini, G. C. Schatz, and T. Belytschko, Phys. Rev. B 75, 075412 (2007).
  * Torras  _et al._ (2009) J. Torras, F. Rodríguez-Ropero, O. Bertran, and C. Alemán, J. Phys. Chem. C 113, 3574 (2009).
  * Zimmerman  _et al._ (2011) P. M. Zimmerman, M. Head-Gordon, and A. T. Bell, J. Chem. Theory Comput. 7, 1695 (2011).
  * Nasluzov  _et al._ (2001) V. A. Nasluzov, V. V. Rivanenkov, A. B. Gordienko, K. M. Neyman, U. Birkenheuer, and N. Rösch, J. Chem. Phys. 115, 8157 (2001).
  * Nasluzov  _et al._ (2003) V. A. Nasluzov, E. A. Ivanova, A. M. Shor, G. N. Vayssilov, U. Birkenheuer, and N. Rösch, J. Phys. Chem. B 107, 2228 (2003).
  * Brändle and Sauer (1998) M. Brändle and J. Sauer, J Am. Chem. Soc. 120, 1556 (1998).
  * Guareschi  _et al._ (2016) R. Guareschi, H. Zulfikri, C. Daday, F. M. Floris, C. Amovilli, B. Mennucci, and C. Filippi, J. Chem. Theory Comput. 12, 1674 (2016).
  * Senn and Thiel (2009) H. M. Senn and W. Thiel, Angew. Chem., Int. Ed. Engl. 48, 1198 (2009).
  * Acevedo and Jorgensen (2010) O. Acevedo and W. L. Jorgensen, Acc. Chem. Res. 43, 142 (2010).
  * Huang  _et al._ (2017) M. Huang, T. Dissanayake, E. Kuechler, B. K. Radak, T.-S. Lee, T. J. Giese, and D. M. York, J. Chem. Theory Comput. 13, 3975 (2017).
  * Alves  _et al._ (2015) C. N. Alves, J. R. A. Silva, and A. E. Roitberg, Phys. Chem. Chem. Phys. 17, 17790 (2015).
  * Silva  _et al._ (2015) J. R. A. Silva, A. E. Roitberg, and C. N. Alves, J. Phys. Chem. B 119, 1468 (2015).
  * Bao  _et al._ (2019) L. Bao, W. Liu, Y. Li, X. Wang, F. Xu, Z. Yang, Y. Yue, C. Zuo, Q. Zhang, and W. Wang, Chem. Res. Toxicol. 32, 603 (2019).
  * Marín  _et al._ (2019) M. d. C. Marín, L. De Vico, S. S. Dong, L. Gagliardi, D. G. Truhlar, and M. Olivucci, J. Chem. Theory Comput. 15, 1915 (2019).
  * Robertazzi and Platts (2006) A. Robertazzi and J. A. Platts, J. Phys. Chem. A 110, 3992 (2006).
  * Gkionis and Platts (2009) K. Gkionis and J. Platts, J. Biol. Inorg. Chem. 14, 1165 (2009).
  * Kupfer  _et al._ (2014) S. Kupfer, L. Zedler, J. Guthmuller, S. Bode, M. D. Hager, U. S. Schubert, J. Popp, S. Gräfe, and B. Dietzek, Phys. Chem. Chem. Phys. 16, 12422 (2014).
  * Osoegawa  _et al._ (2019) S. Osoegawa, R. Miyoshi, K. Watanabe, Y. Hirose, T. Fujisawa, M. Ikeuchi, and M. Unno, J. Phys. Chem. B 123, 3242 (2019).
  * Loco and Piquemal (2022) D. Loco and J.-P. Piquemal, in  _Photochemistry: Volume 50_ (The Royal Society of Chemistry, 2022) pp. 386–402.
  * Cho and Lester (2010) H. M. Cho and W. A. Lester, Jr., J. Phys. Chem. Lett. 1, 3376 (2010).
  * Kotliar  _et al._ (2006) G. Kotliar, S. Y. Savrasov, K. Haule, V. S. Oudovenko, O. Parcollet, and C. A. Marianetti, Rev. Mod. Phys. 78, 865 (2006).
  * Lan  _et al._ (2015) T. N. Lan, A. A. Kananenka, and D. Zgid, J. Chem. Phys. 143, 241102 (2015).
  * Knizia and Chan (2012) G. Knizia and G. K.-L. Chan, Phys. Rev. Lett. 109, 186404 (2012).
  * Wesolowski  _et al._ (2015) T. A. Wesolowski, S. Shedge, and X. Zhou, Chem. Rev. 115, 5891 (2015).
  * Libisch  _et al._ (2014) F. Libisch, C. Huang, and E. A. Carter, Acc. Chem. Res. 47, 2768 (2014).
  * Manby  _et al._ (2012) F. R. Manby, M. Stella, J. D. Goodpaster, and T. F. Miller, III, J. Chem. Theory Comput. 8, 2564 (2012).
  * Katin  _et al._ (2017) K. P. Katin, V. S. Prudkovskiy, and M. M. Maslov, Phys. Lett. A 381, 2686 (2017).
  * Friedrich  _et al._ (2007) J. Friedrich, M. Hanrath, and M. Dolg, J. Chem. Phys. 126, 154110 (2007).
  * Seth  _et al._ (2002) M. Seth, P. M. Margl, and T. Ziegler, Macromolecules 35, 7815 (2002).
  * Vogiatzis  _et al._ (2015) K. D. Vogiatzis, W. Klopper, and J. Friedrich, J. Chem. Theory Comput. 11, 1574 (2015).
  * Huo  _et al._ (2016) P. Huo, C. Uyeda, J. D. Goodpaster, J. C. Peters, and T. F. Miller, III, ACS Catal. 6, 6114 (2016).
  * Lee  _et al._ (2019) S. J. R. Lee, M. Welborn, F. R. Manby, and T. F. Miller, III, Acc. Chem. Res. 52, 1359 (2019).
  * Coughtrie  _et al._ (2018) D. J. Coughtrie, R. Giereth, D. Kats, H.-J. Werner, and A. Köhn, J. Chem. Theory Comput. 14, 693 (2018).
  * Neugebauer  _et al._ (2005) J. Neugebauer, M. J. Louwerse, E. J. Baerends, and T. A. Wesolowski, J. Chem. Phys. 122, 094115 (2005).
  * Jacob and Visscher (2006) C. R. Jacob and L. Visscher, J. Chem. Phys. 125, 194104 (2006).
  * Exner  _et al._ (2012) T. E. Exner, A. Frank, I. Onila, and H. M. Möller, J. Chem. Theory Comput. 8, 4818 (2012).
  * Exner and Mezey (2004) T. E. Exner and P. G. Mezey, J. Phys. Chem. A 108, 4301 (2004).
  * Bockstedte  _et al._ (2018) M. Bockstedte, F. Schütz, T. Garratt, V. Ivády, and A. Gali, npj Quantum Mater. 3, 31 (2018).
  * Lau  _et al._ (2021) B. T. G. Lau, G. Knizia, and T. C. Berkelbach, J. Phys. Chem. Lett. 12, 1104 (2021).
  * Schäfer  _et al._ (2021) T. Schäfer, A. Gallo, A. Irmler, F. Hummel, and A. Grüneis, J. Chem. Phys. 155, 244103 (2021).
  * Ma  _et al._ (2021) H. Ma, N. Sheng, M. Govoni, and G. Galli, J. Chem. Theory Comput. 17, 2116 (2021).
  * Petocchi  _et al._ (2020) F. Petocchi, F. Nilsson, F. Aryasetiawan, and P. Werner, Phys. Rev. Res, 2, 013191 (2020).
  * Yeh  _et al._ (2021) C.-N. Yeh, S. Iskakov, D. Zgid, and E. Gull, Phys. Rev. B 103, 195149 (2021).
  * Nowadnick  _et al._ (2015) E. A. Nowadnick, J. P. Ruf, H. Park, P. D. C. King, D. G. Schlom, K. M. Shen, and A. J. Millis, Phys. Rev. B 92, 245109 (2015).
  * Chen  _et al._ (2022) H. Chen, A. Hampel, J. Karp, F. Lechermann, and A. J. Millis, Front. Phys. 10, 835942 (2022).
  * Nusspickel and Booth (2022) M. Nusspickel and G. H. Booth, Phys. Rev. X 12, 011046 (2022).
  * Bovi  _et al._ (2013) D. Bovi, D. Narzi, and L. Guidoni, Angew. Chem., Int. Ed. Engl. 52, 11744 (2013).
  * Benali  _et al._ (2014) A. Benali, L. Shulenburger, N. A. Romero, J. Kim, and O. A. von Lilienfeld, J. Chem. Theory Comput. 10, 3417 (2014).
  * Al-Hamdani  _et al._ (2021) Y. S. Al-Hamdani, P. R. Nagy, A. Zen, D. Barton, M. Kállay, J. G. Brandenburg, and A. Tkatchenko, Nat. Commun. 12, 3927 (2021).
  * Reta Mañeru  _et al._ (2014) D. Reta Mañeru, A. K. Pal, I. d. P. R. Moreira, S. N. Datta, and F. Illas, J. Chem. Theory Comput. 10, 335 (2014).
  * Boyn and Mazziotti (2021) J.-N. Boyn and D. A. Mazziotti, J. Chem. Phys. 154, 134103 (2021).
  * Bao  _et al._ (2016) J. L. Bao, A. Sand, L. Gagliardi, and D. G. Truhlar, J. Chem. Theory Comput. 12, 4274 (2016).
  * Varsano  _et al._ (2017) D. Varsano, S. Sorella, D. Sangalli, M. Barborini, S. Corni, E. Molinari, and M. Rontani, Nat. Commun. 8, 1461 (2017).
  * Barborini  _et al._ (2022) M. Barborini, M. Calandra, F. Mauri, L. Wirtz, and P. Cudazzo, Phys. Rev. B 105, 075122 (2022).
  * Pozun  _et al._ (2013) Z. D. Pozun, X. Su, and K. D. Jordan, J. Am. Chem. Soc. 135, 13862 (2013).
  * Barborini and Coccia (2015) M. Barborini and E. Coccia, J. Chem. Theory Comput. 11, 5696 (2015).
  * Abe (2013) M. Abe, Chem. Rev. 113, 7011 (2013).
  * Bartlett and Musiał (2007) R. J. Bartlett and M. Musiał, Rev. Mod. Phys. 79, 291 (2007).
  * Foulkes  _et al._ (2001) W. M. C. Foulkes, L. Mitas, R. J. Needs, and G. Rajagopal, Rev. Mod. Phys. 73, 33 (2001).
  * Kalos and Whitlock (2008) M. H. Kalos and P. A. Whitlock, “Quantum monte carlo,” in  _Monte Carlo Methods_ (John Wiley & Sons, Ltd, 2008) Chap. 8, pp. 159–178.
  * Becca and Sorella (2017) F. Becca and S. Sorella,  _Quantum Monte Carlo Approaches for Correlated Systems_ (Cambridge University Press, 2017).
  * Jones  _et al._ (2009) A. Jones, A. Thompson, J. Crain, M. H. Müser, and G. J. Martyna, Phys. Rev. B 79, 144119 (2009).
  * Tkatchenko  _et al._ (2012) A. Tkatchenko, R. A. DiStasio, R. Car, and M. Scheffler, Phys. Rev. Lett. 108, 236402 (2012).
  * Jones  _et al._ (2013a) A. Jones, F. Cipcigan, V. P. Sokhan, J. Crain, and G. J. Martyna, Phys. Rev. Lett. 110, 227801 (2013a).
  * Jones  _et al._ (2013b) A. Jones, J. Crain, F. Cipcigan, V. Sokhan, M. Modani, and G. Martyna, Mol. Phys. 111, 3465 (2013b).
  * Jones  _et al._ (2013c) A. P. Jones, J. Crain, V. P. Sokhan, T. W. Whitfield, and G. J. Martyna, Phys. Rev. B 87, 144103 (2013c).
  * Odbadrakh  _et al._ (2015) T. Odbadrakh, V. Voora, and K. Jordan, Chem. Phys. Lett. 630, 76 (2015).
  * Odbadrakh and Jordan (2016) T. T. Odbadrakh and K. D. Jordan, J. Chem. Phys. 144, 034111 (2016).
  * Cipcigan  _et al._ (2019) F. S. Cipcigan, J. Crain, V. P. Sokhan, and G. J. Martyna, Rev. Mod. Phys. 91, 025003 (2019).
  * Tkatchenko and Scheffler (2009) A. Tkatchenko and M. Scheffler, Phys. Rev. Lett. 102, 073005 (2009).
  * Vaccarelli  _et al._ (2021) O. Vaccarelli, D. V. Fedorov, M. Stöhr, and A. Tkatchenko, Phys. Rev. Res 3, 033181 (2021).
  * Poier  _et al._ (2022) P. P. Poier, L. Lagardère, and J.-P. Piquemal, J. Chem. Theory Comput. 18, 1633 (2022).
  * Karimpour  _et al._ (2022a) M. R. Karimpour, D. V. Fedorov, and A. Tkatchenko, Phys. Rev. Res. 4, 013011 (2022a).
  * Karimpour  _et al._ (2022b) M. R. Karimpour, D. V. Fedorov, and A. Tkatchenko, J. Phys. Chem. Lett. 13, 2197 (2022b).
  * Lemkul  _et al._ (2016) J. A. Lemkul, J. Huang, B. Roux, and A. D. MacKerell, Jr., Chem. Rev. 116, 4983 (2016).
  * Bui  _et al._ (2023) A. T. Bui, F. L. Thiemann, A. Michaelides, and S. J. Cox, Nano Lett. 23, 580 (2023).
  * Sommerfeld and Jordan (2005) T. Sommerfeld and K. D. Jordan, J. Phys. Chem. A 109, 11531 (2005).
  * Margoliash and Meath (1978) D. J. Margoliash and W. J. Meath, J. Chem. Phys. 68, 1426 (1978).
  * Góger  _et al._ (2023) S. Góger, A. Khabibrakhmanov, O. Vaccarelli, D. V. Fedorov, and A. Tkatchenko, arXiv , 2212.07305v1 (2023).
  * (84) See supporting information for further details on the drudon-electron model, the quantum Monte Carlo methods and the computational details of the calculations, and the construction and geometries of the water cages.
  * Sadhukhan and Manby (2016) M. Sadhukhan and F. R. Manby, Phys. Rev. B 94, 115106 (2016).
  * Austin  _et al._ (2012) B. M. Austin, D. Y. Zubarev, and W. A. Lester, Jr., Chem. Rev. 112, 263 (2012).
  * (87) M. Barborini, “Quantum Mecha (QMeCha) package (private repository March 2023),” .
  * Sorella (2001) S. Sorella, Phys. Rev. B 64, 024512 (2001).
  * Sorella (2005) S. Sorella, Phys. Rev. B 71, 241103 (2005).
  * Ambrosetti  _et al._ (2014) A. Ambrosetti, A. M. Reilly, R. A. DiStasio, and A. Tkatchenko, J. Chem. Phys. 140, 18A508 (2014).
  * Patkowski  _et al._ (2005) K. Patkowski, G. Murdachaew, C.-M. Fou, and K. Szalewicz, Mol. Phys. 103, 2031 (2005).
  * Metz  _et al._ (2019) M. P. Metz, K. Szalewicz, J. Sarka, R. Tóbiás, A. G. Császár, and E. Mátyus, Phys. Chem. Chem. Phys. 21, 13504 (2019).
  * Fedorov  _et al._ (2018) D. V. Fedorov, M. Sadhukhan, M. Stöhr, and A. Tkatchenko, Phys. Rev. Lett. 121, 183401 (2018).
  * Shirkov and Sladek (2017) L. Shirkov and V. Sladek, J. Chem. Phys. 147, 174103 (2017).
  * Altun  _et al._ (2018) A. Altun, F. Neese, and G. Bistoni, Beilstein J. Org. Chem. 14, 919 (2018).
  * DiStasio  _et al._ (2007) R. A. DiStasio, G. von Helden, R. P. Steele, and M. Head-Gordon, Chem. Phys. Lett. 437, 277 (2007).
  * Carter-Fenk and Herbert (2020) K. Carter-Fenk and J. M. Herbert, Chem. Sci. 11, 6758 (2020).
  * Shee  _et al._ (2019) J. Shee, E. J. Arthur, S. Zhang, D. R. Reichman, and R. A. Friesner, J. Chem. Theory Comput. 15, 4924 (2019).
  * Adamo and Barone (1999) C. Adamo and V. Barone, J. Chem. Phys. 110, 6158 (1999).
  * Sirianni  _et al._ (2022) D. A. Sirianni, X. Zhu, D. F. Sitkoff, D. L. Cheney, and C. D. Sherrill, J. Chem. Phys. 156, 194306 (2022).
  * Wenthold  _et al._ (1998) P. G. Wenthold, R. R. Squires, and W. C. Lineberger, J. Am. Chem. Soc. 120, 5279 (1998).
  * Zhou and Wang (2019) X. Zhou and F. Wang, Phys. Chem. Chem. Phys. 21, 20422 (2019).
  * Unke  _et al._ (2021a) O. T. Unke, S. Chmiela, H. E. Sauceda, M. Gastegger, I. Poltavsky, K. T. Schütt, A. Tkatchenko, and K.-R. Müller, Chem. Rev. 121, 10142 (2021a).
  * Stöhr  _et al._ (2020) M. Stöhr, L. Medrano Sandonas, and A. Tkatchenko, J. Phys. Chem. Lett. 11, 6835 (2020).
  * Unke  _et al._ (2021b) O. Unke, S. Chmiela, M. Gastegger, K. Schütt, H. E. Sauceda, and K.-R. Müller, Nat. Commun. 12, 7273 (2021b).
  * Hermann  _et al._ (2020) J. Hermann, Z. Schätzle, and F. Noé, Nat. Chem. 12, 891 (2020).
  * Varrette  _et al._ (2014) S. Varrette, P. Bouvry, H. Cartiaux, and F. Georgatos, in  _Proc. of the 2014 Intl. Conf. on High Performance Computing & Simulation (HPCS 2014)_ (IEEE, Bologna, Italy, 2014) pp. 959–967.