# Adaptive Belief Discretization for POMDP Planning

Divya Grover Department of Computer Science and Engineering
Chalmers University, Sweden  Christos Dimitrakakis  Department of Computer
Science and Engineering
Chalmers University, Sweden  Department of Informatics
University of Oslo, Norway

###### Abstract

Partially Observable Markov Decision Processes (POMDP) is a widely used model
to represent the interaction of an environment and an agent, under state
uncertainty. Since the agent does not observe the environment state, its
uncertainty is typically represented through a probabilistic belief. While the
set of possible beliefs is infinite, making exact planning intractable, the
belief space’s complexity (and hence planning complexity) is characterised by
its covering number. Many POMDP solvers uniformly discretize the belief space
and give the planning error in terms of the (typically unknown) covering
number. We instead propose an adaptive belief discretization scheme, and give
its associated planning error. We furthermore characterise the covering number
with respect to the POMDP parameters. This allows us to specify the exact
memory requirements on the planner, needed to bound the value function error.
We then propose a novel, computationally efficient solver using this scheme.
We demonstrate that our algorithm is highly competitive with the state of the
art in a variety of scenarios.

##  1 Introduction

We are interested in sequential decision making problems, where an agents is
interacting with a partially observable environment, meaning that it is
unaware of the actual state of the environment. This interaction can be
formalised through a Partially Observable Markov Decision Process (POMDP)
[Smallwood and Sondik, 1973]. Since the agent does not know the true state of
the environment, it must instead maintain a subjective belief representing its
uncertainty about the current environment state. This is typically expressed
as a probability distribution over the states.

Online POMDP algorithms plan ahead by starting from the current belief and
building a tree-like structure in the space of possible future beliefs, by
enumerating different environment responses to the agent’s actions. However,
exact solutions to some lookahead horizon require an exponential size tree
with respect to that horizon. Typical approximations include action selection
heuristics that only consider promising actions and leaf node approximations,
that truncate the tree and replace the value of the leaf nodes with
heuristics, or upper or lower bounds. While both these techniques are
effective when building shallow trees, we instead propose an online POMDP
algorithm that builds deep trees. In particular, we propose an adaptive
discretization scheme that generates a cover with a small size, thus reducing
memory and computational complexity.

###  1.1 Related work

Current state-of-the-art online POMDP solvers include POMCP [Silver and
Veness, 2010] and DESPOT Ye et al. [2017]. Like many older [Ross et al., 2007]
and newer [Kurniawati and Yadav, 2016] works they are based on inference using
particle filtering. In building the lookahead tree, they keep its size in
check by only considering ‘promising’ sequence of actions, by maintaining a
statistical estimate of the action values (known as Upper Confidence Bound).
DESPOT in addition, uses importance sampling to not miss important beliefs.
But neither algorithm gives a structured way to check the approximate
inference’s effect on the planning error. While the belief discretization we
propose results in a form of approximate inference as well, the focus here is
on how the approximations made in the inference and in the planning process
_interact_.

A different method for approximate inference in POMDPs includes belief
representation methods like [Roy et al., 2005]. These implicitly define a
kernel function over the belief space, but give no bound on the resultant
solver’s planning error.

One concept we use in our theoretical analysis is the covering number of the
belief space. It is a measure of the density of a set of points in any
continuous space. Only recently has it been used as a complexity measure for
the POMDP. In particular, Zhang et al. [2012] performs an experimental
estimation of the covering number for various environments, by using breadth
first search (BFS). However, unlike this paper, the authors do not provide a
POMDP solver or derive any dependence between the underlying POMDP parameters
and the covering number of the belief space.

Covering number was also used in the analysis of SARSOP [Kurniawati et al.,
2008] by Lee et al. [2007]. In particular, it shows how the planner’s error
depends on the (unknown) covering number. Such solvers still depends
exponentially on the lookahead, as noted by [Zhang et al., 2012] (refer Lemma
1). Finally, Zhang et al. [2014] also propose PGVI, which like SARSOP, uses
depth first search (DFS) to estimate the optimal action. Analytically, they
are the same, since PGVI’s complexity is identical to that of SARSOP.111Verify
that the tuned parameters are set to identical values in both works.

#### Contributions:

  1. 1.

We propose explicit lookahead dependent covers for the belief space. This
significantly reduces the memory requirement on the planner.

  2. 2.

We give the exact dependence of the planning error on this adaptive
discretization scheme.

  3. 3.

We give an upper bound on the memory requirement of the planner by giving
bounds on the covering number of the belief space.

  4. 4.

We propose a novel, computationally efficient planner with direct control on
the planning error with respect to its tuning parameters.

#### Organization.

In Section 2, we define the POMDP model, state the used equations and
introduce the concept of cover sets. In Section 3, we follow up with analysis,
giving a proof of sketch first, then stating our main result, and finally the
proofs. In Section 4, we describe our proposed algorithm and Section 5 is
dedicated to experiments.

##  2 Setting

We divide the section into three parts: First giving the definitions of POMDP
model, second giving the equations used to compute various quantities of
interest in a POMDP and third part introducing the concept of
$\epsilon$-cover.

###  2.1 Model Definition

######  Definition 1 (MDP).

A Markov Decision Process (MDP) is a discrete-time stochastic process that
provides a formal framework for decision making agents. An MDP
$\mu\mathrel{\triangleq}(S,A,T,R,\gamma)$ is composed of a state space $S$, an
action space $A$, a reward distribution $R$ and a transition function $T$. The
transition function
$T\mathrel{\triangleq}\mathop{\mathbb{P}}\nolimits(s_{t+1}|s_{t},a_{t})$
dictates the distribution over next states $s_{t+1}$ given the present state-
action pair $(s_{t},a_{t})$. The reward distribution
$R\mathrel{\triangleq}\mathop{\mathbb{P}}\nolimits(r_{t+1}|s_{t},a_{t})$
dictates the obtained reward with $r\in[0,\mathcal{R}_{\max}]$.

The discount $\gamma\in[0,1)$ signifies the effective horizon
$H\mathrel{\triangleq}\frac{1}{1-\gamma}$ of the agent’s interaction with its
environment. _Henceforth, we consider only finite state, action and
observation spaces_.

In many interesting settings, the agent can only indirectly estimate the state
$s_{t}$ via a set of observations $o_{t}\in O$ and their (state-dependent)
emission probability. Formally,

######  Definition 2 (POMDP).

A Partially Observable MDP, or POMDP, denoted by $(S,A,T,R,,\gamma,O,Z)$ tuple
is an MDP model augmented with an observation $o_{t}\in O$ set and their
corresponding emission probabilities
$Z\mathrel{\triangleq}\mathop{\mathbb{P}}\nolimits(o_{t}|s_{t})$.

In such cases, the best an agent could hope to achieve is to act optimally in
an expected sense (averaging over its estimation of the current state).
Formally, the objective is to find action maximizing

$V^{*}(b_{0})\mathrel{\triangleq}\lim_{H\rightarrow\infty}V_{H}(b_{0})$ (1)

where $V_{H}(b_{0})$, called the belief value function, is defined as

$\displaystyle
V_{H}(b_{0})=\max_{a}\mathbb{E}(\sum_{k=0}^{H-1}\gamma^{k}r_{k+1}|b_{0},a)$
(2) $\displaystyle b_{0}\in\Delta^{|S|}$
$\displaystyle\Delta^{|S|}\mathrel{\triangleq}\\{b\in{[0,1]}^{|S|}\;|\;\sum_{s}b(s)=1\\}$
(3)

Where $b_{0}$ denotes its current state estimation, known as belief. Formally,
belief is a probability distribution over the state space whose element $b(s)$
gives the probability that the system’s true state is $s$. It is a compact
representation of the agent’s complete history of interaction with its
environment, i.e, the action-observation sequence it followed.

Figure 1 shows the interaction of an agent with its environment. At each time
step, it evaluates the best possible action by anticipating the effect of a
sequence of actions on future belief and reward. This process of anticipation
is called planning and the corresponding tree-like structure is called the
belief tree. Although in Figure 1, we only calculate effect of a fixed action,
you can still note the exponential growth of the number of nodes with depth.

###  2.2 Equations

Here we describe the standard computational equations of the POMDP model. To
compute a new belief $b_{k+1}$, given the current belief $b_{k}$, action $a$
and next observation $o_{k+1}$, one uses the following set of equations:

$\displaystyle\mathop{\mathbb{P}}\nolimits(s_{k+1})=\sum_{s^{\prime}}T(s_{k+1}|s^{\prime},a)b_{k}(s^{\prime})$
(4) $\displaystyle
b_{k+1}(s_{k+1})=\frac{1}{\eta}\;Z(o_{k+1}|s_{k+1})\mathop{\mathbb{P}}\nolimits(s_{k+1})$
(5)
$\displaystyle\eta\mathrel{\triangleq}\sum_{s_{k+1}}Z(o_{k+1}|s_{k+1})\mathop{\mathbb{P}}\nolimits(s_{k+1})$

collectively referred to as the _Bayes filter_. Equation 4 is state marginal
distribution after transition has occurred under the given action, while
equation 5 gives the Bayesian posterior (called the next belief) of the
states.
Solving a POMDP involves first generating the belief tree, then doing backward
induction [Bellman, 1952] on it to calculate the optimal action at the root
node. The corresponding Bellman equation for POMDP is:

$V(b_{k})=\max_{a}\mathbb{E}_{\mathop{\mathbb{P}}\nolimits(r_{k+1},o_{k+1}|a,b_{k})}(\;r_{k+1}+V(b_{k+1})\;)$
(6)

Where

$\displaystyle\mathop{\mathbb{P}}\nolimits(r_{k+1},o_{k+1}|a,b_{k})$
$\displaystyle=\mathop{\mathbb{P}}\nolimits(o_{k+1}|a,b_{k})\mathop{\mathbb{P}}\nolimits(r_{k+1}|a,b_{k})$
$\displaystyle\mathop{\mathbb{P}}\nolimits(r_{k+1}|a,b_{k})$
$\displaystyle=\sum_{s^{\prime}}R(r_{k+1}|s^{\prime},a)b_{k}(s^{\prime})$
$\displaystyle\mathop{\mathbb{P}}\nolimits(o_{k+1}|a,b_{k})$
$\displaystyle=\eta$

Note that $b_{k+1}$ is uniquely determined by $o_{k+1}$ using Bayes filter.

$b_{t}$$b_{t+1}$$a$$a$$b_{o_{2}}$$b_{o_{1}}$$b_{o_{2}}$$b_{o_{1}}$$a$$a$$b_{o_{4}}$$b_{o_{3}}$$a$$a$$b_{o_{4}}$$b_{o_{3}}$$t$$d=1$$d=2$
Figure 1: Online POMDP planning

###  2.3 Cover sets and covering number

For a set $\mathcal{X}\subset\mathbb{R}^{n}$, define metric $d_{p}:x\times
y\rightarrow\mathbb{R}\;|\;x,y\in\mathcal{X}$:

$d_{p}(x,y)\mathrel{\triangleq}{\left(\sum^{n}_{i=1}|x_{i}-y_{i}|^{p}\right)}^{\frac{1}{p}}$

######  Definition 3.

Given a metric space $(\mathcal{X},d_{p})$. Define $\epsilon$-ball of a point
$x_{0}\in\mathcal{X}$ as

$c(x_{0},\epsilon)\mathrel{\triangleq}\\{x\in\mathcal{X}\;|\;d_{p}(x_{0},x)<\epsilon\\}$

Then for a given set $B\subset\mathcal{X}$, define an $\epsilon$-cover of $B$
as

$C_{B}(\epsilon)\mathrel{\triangleq}\\{x\in
B\;|\;B\subset\cup\;c(x,\epsilon)\\}$

Similarly an improper cover of $B$ is defined as

$\displaystyle\hat{C}_{B}(\epsilon)\mathrel{\triangleq}\\{x\in\mathcal{X}\;|\;B\subset\cup\;c(x,\epsilon)\\}$

Therefore the elements of an improper cover do not necessarily have to come
from set $B$ itself.

The covering number is defined as

$\mathcal{N}(\epsilon,B,d_{p})\mathrel{\triangleq}\min|C_{B}(\epsilon)|$

i.e, the minimum cardinality sets amongst such sets.

An improper covering number is similarly defined

$\hat{\mathcal{N}}(\epsilon,B,d_{p})\mathrel{\triangleq}\min|\hat{C}_{B}(\epsilon)|$

A standard result used here is:

$\mathcal{N}(2\epsilon,B,d_{p})\leq\hat{\mathcal{N}}(\epsilon,B,d_{p})\leq\mathcal{N}(\epsilon,B,d_{p})$
(7)

Note that we refer to an $\epsilon$-cover as a cover set when the value
$\epsilon$ is either implied or not important in the context.

##  3 Analysis

Our analysis relies on creating a sequence of covers for the belief, so that
covers further away in the future are sparser. We show that we can uniformly
bound the value function error by an appropriate increase in the sparsity that
depends on the depth and discount factor.

Our main theorem, the _planning error theorem_ , gives a bound on the number
of representative beliefs we need to keep track of, to bound the value
function error at the root. This is accompanied by a lemma bounding the memory
required by our discretization scheme and algorithm.

In particular, for exact $H$-step lookahead planning, the corresponding belief
tree would have order $O({|A|}^{H}{|O|}^{H})$ nodes. This is prohibitive to
store in memory even for a moderate horizon. A remedy to this problem comes
from the fact that belief values functions are Lipschitz continuous. While
previous work [c.f. Lee et al., 2007] used this property to obtain a uniform
cover, in our work we make the cover sparser as we go deeper in the tree.

The Lipschitz property allows us to control tree size by expanding only a
certain representative node’s subtree and substituting its value for all
beliefs represented by that particular belief node. We can then relate the
approximation errors at different depths through discounting [c.f. Kearns et
al., 2002]. Combined together these ideas gives the required planning error
theorem.

###  3.1 Main Result

######  Definition 4.

A $\epsilon$-discretization of a belief space $\Delta^{|S|}$ is
defined222Ignoring the normalization constant. as:

$A(\epsilon)=\\{\sum_{i=1}^{S}\alpha_{i}{\mathbb{1}}_{i}|\;\alpha_{i}\in\\{0,\epsilon,2\epsilon,\ldots,1-\epsilon,1\\}\\}$

i.e, a vector space with ${\mathbb{1}}_{i}$ as basis and coefficients
$\alpha_{i}$. ${\mathbb{1}}_{i}$ denotes a unit vector with
${\mathbb{1}}_{i}(j)=0\;|j\neq i$.

The error between root belief’s value function to its true value function is
given by the planning error theorem:

######  Theorem 1.

For any belief $b_{0}$, and target error $\epsilon_{v}$, there exists a belief
tree $\mathcal{T}$ of height $h_{0}$ rooted at $b_{0}$, with separate
$\epsilon_{d}$-discretization of the belief spaces at each depth $d$, given by
the sequence
$\epsilon_{d}\mathrel{\triangleq}\frac{\epsilon_{b_{0}}}{\gamma^{d}}$, such
that $b_{0}$’s value $\hat{V}(b_{0})$ at the root satisfies

$|\hat{V}(b_{0})-V^{*}(b_{o})|\leq\epsilon_{v}$

for parameter values:

$\displaystyle
h_{0}=\log_{\gamma}\frac{{(1-\gamma)}\epsilon_{v}}{2\mathcal{R}_{\max}}$
$\displaystyle\epsilon_{b_{0}}=\frac{{(1-\gamma)}\epsilon_{v}}{2\gamma\mathcal{R}_{\max}h_{0}}$

Contrast $\epsilon_{b_{0}}$ with the value of parameter
$\delta(=\frac{{{(1-\gamma)}^{2}}\epsilon_{v}}{2\gamma\mathcal{R}_{\max}})$ in
[Lee et al., 2007]. Both denote the level of discretization required. While
our scheme demands a finer initial discretization ($\epsilon_{b_{0}}$) at the
root (by a factor of $\frac{1}{h_{0}(1-\gamma)}$), it quickly (because
$\epsilon_{d}\mathrel{\triangleq}\frac{\epsilon_{b_{0}}}{\gamma^{h}}$) becomes
coarser (better) as we expand the tree deeper.
Now we state our Lemma giving total memory required to build the discretized
belief tree $\mathcal{T}$.

######  Lemma 1.

For a tree $\mathcal{T}$ of total height $h_{0}$ with separate
$\epsilon_{d}$\- discretized belief spaces at each level $d$, as defined in
Theorem 1. The total memory required $M_{h_{0}}$ is bounded by:

$M_{h_{0}}\leq{\lceil\frac{2}{\epsilon_{b_{0}}(1-\gamma)}\rceil}^{S}$

Note that the memory requirement of SARSOP and PGVI is the covering number of
the minimal $\delta$-cover for their specified tolerance value $\delta$ [Lee
et al., 2007]. But they don’t give an upper bound on its value, because the
corresponding cover sets are only implicit.

###  3.2 Proofs

Since most proofs depend on backward induction equation 6, we use the
following indexing convention: For a belief tree of total height $h_{0}$, we
define its height by index $i$, hence $i=0$ denotes the leaf nodes level,
while $i=h_{0}$ denotes the root level. Also for depth $d$, the relation
$d=h_{0}-i$ holds. This convention is also followed in the previous
literature.
We first give a Lemma on Lipschitz continuity of value function at different
heights:

######  Lemma 2 (Lipschitz continuity).

For any two belief nodes $b$ and $b^{\prime}$ at height $i$, we have

${||V_{i}(b)-V_{i}(b^{\prime})||}_{1}\leq L_{i}\;\delta_{i}$

where

$\displaystyle
L_{i}\leq\mathcal{R}_{\max}\frac{1-\gamma^{i}}{1-\gamma}+\gamma^{i}L_{0}$
$\displaystyle\delta_{i}={||b-b^{\prime}||}_{1}$

Proof is a simple consequence of H$\ddot{o}$lder’s inequality.

$\displaystyle{||V_{i}(b)-V_{i}(b^{\prime})||}_{1}$
$\displaystyle\mathrel{\triangleq}$
$\displaystyle=\sum_{s}V_{i}(s)(b(s)-b^{\prime}(s))$
$\displaystyle\leq{||V_{i}(s)||}_{\infty}\times{||b(s)-b^{\prime}(s)||}_{1}$
$\displaystyle=\max_{s}V_{i}(s)\times\delta_{i}$
$\displaystyle=L_{i}\delta_{i}$

where second line follows from H$\ddot{o}$lder’s inequality.
Hence by using equation 6 we get

$\displaystyle L_{i}$ $\displaystyle\mathrel{\triangleq}\max_{s}V_{i}$
$\displaystyle\max_{s}V_{i}$ $\displaystyle=\max_{s}\;(r+\gamma V_{i-1})$
$\displaystyle\max_{s}V_{i}$
$\displaystyle=\mathcal{R}_{\max}+\gamma\max_{s}V_{i-1}$

which applied recursively gives the result $\square$

We next give a Lemma relating value function error across different heights:

######  Lemma 3.

Maximum error for any node at height $h$, denoted by $\epsilon_{h}$ is given
by:

$\epsilon_{h}\leq\sum^{h-1}_{i=0}\gamma^{h-i}L_{i}\delta_{i}+\gamma^{h}\epsilon_{0}$

where $\epsilon_{0}$ denotes the error at leaf nodes.

Proof. For a belief $b$ at level $i=0$, consider its error in value function
$\epsilon(b)$ as:

$\displaystyle\epsilon(b)$ $\displaystyle=|\hat{V}(b)-V^{*}(b)|$
$\displaystyle\leq|\hat{V}(b)-V(b^{\prime})|+|V(b^{\prime})-V^{*}(b)|$
$\displaystyle\leq L_{0}\delta_{0}+\epsilon_{0}$

Where line 2 comes from the fact that it maybe the case that value of $b$ is
substituted by value of another belief $b^{\prime}$, constituting the
$\delta_{0}$-cover at level $i=0$. Moreover, for any parent of $b$ (say
$b_{i}$), at level $i=1$, its error $\epsilon_{1}\leq\gamma\epsilon(b)$
[Kearns et al., 2002]. Hence by recursion, we get:

$\displaystyle\epsilon_{1}$
$\displaystyle\leq\gamma(L_{0}\delta_{0}+\epsilon_{0})$
$\displaystyle\epsilon_{2}$ $\displaystyle\leq\gamma(L_{1}\delta_{1}+\gamma
L_{0}\delta_{0}+\gamma\epsilon_{0})$ $\displaystyle\ldots$
$\displaystyle\epsilon_{h}$
$\displaystyle\leq\sum^{h-1}_{i=0}\gamma^{h-i}L_{i}\delta_{i}+\gamma^{h}\epsilon_{0}$

Where $\epsilon_{i}$ denotes the error in value function for any belief at
level $i$.

#### Proof of Theorem 1:

Equating root error in Lemma 3 to target error $\epsilon_{v}$ gives:

$\sum^{h-1}_{i=0}\gamma^{h-i}L_{i}\delta_{i}+\gamma^{h}\epsilon_{0}=\epsilon_{v}$

Combining this with Lemma 2 gives:

$\sum^{h-1}_{i=0}\gamma^{h-i}\left(\mathcal{R}_{\max}\frac{1-\gamma^{i}}{1-\gamma}+\gamma^{i}L_{0}\right)\delta_{i}+\gamma^{h}\epsilon_{0}=\epsilon_{v}$

Assuming $L_{0}=\frac{\mathcal{R}_{\max}}{1-\gamma}$, and leaf error
$\epsilon_{0}=\frac{\mathcal{R}_{\max}}{1-\gamma}$, we get

$\displaystyle\sum^{h-1}_{i=0}\gamma^{h-i}\epsilon_{0}(1-\gamma^{i}+\gamma^{i})\delta_{i}+\gamma^{h}\epsilon_{0}$
$\displaystyle=\epsilon_{v}$
$\displaystyle\gamma^{h}\epsilon_{0}\left(\sum^{h-1}_{i=0}\frac{\delta_{i}}{\gamma^{i}}+1\right)$
$\displaystyle=\epsilon_{v}$
$\displaystyle\left(\sum^{h-1}_{i=0}\frac{\epsilon_{b_{0}}}{\gamma^{h-i-1}\gamma^{i}}+1\right)$
$\displaystyle=\frac{\epsilon_{v}}{\gamma^{h}\epsilon_{0}}$
$\displaystyle\frac{h\epsilon_{b_{0}}}{\gamma^{h-1}}+1=\frac{\epsilon_{v}}{\gamma^{h}\epsilon_{0}}$

Where we use the fact that

$\displaystyle\delta_{i}$
$\displaystyle\mathrel{\triangleq}\epsilon_{d}=\frac{\epsilon_{b_{0}}}{\gamma^{d}}=\frac{\epsilon_{b_{0}}}{\gamma^{h-i-1}}$

Finally, substituting

$h=h_{0}=\log_{\gamma}\frac{\epsilon_{v}}{2\epsilon_{0}}=\log_{\gamma}\frac{\epsilon_{v}(1-\gamma)}{2\mathcal{R}_{\max}}$

gives $\epsilon_{b_{0}}=\frac{\gamma^{h_{0}-1}}{h_{0}}\;\square$

We now calculate the total memory required to store a discretized tree
$\mathcal{T}$ by proving Lemma 1. Note that $A(\epsilon)$ acts as
$\epsilon$-cover for the set ${[0,1]}^{|S|}$. Since belief space
$\Delta^{|S|}\subset{[0,1]}^{|S|}$, $A(\epsilon)$ also forms an improper cover
for it. We bound the covering number for each depth $d$, by bounding the
cardinality of the corresponding sets $A(\epsilon_{d})$. We consider the
$\epsilon_{d}\mathrel{\triangleq}\frac{\epsilon_{b}}{\gamma^{d}}$-discretization
of any arbitrary sequence of ${[0,1]}^{|S|}$ spaces.

#### Proof of Lemma 1:

For any $h_{0}$ length sequence of discretized spaces, total memory required:

$\displaystyle\hat{M}_{h_{0}}$
$\displaystyle\mathrel{\triangleq}\sum_{d}\mathcal{N}(\epsilon_{d},A(\epsilon_{d}),d_{1})$
$\displaystyle=\sum_{d}|A(\frac{\epsilon_{b}}{\gamma^{d}})|$ $\displaystyle
log\sum_{d}|A(\frac{\epsilon_{b}}{\gamma^{d}})|$
$\displaystyle\leq|S|\sum_{d}\log(1+\lfloor\frac{\gamma^{d}}{\epsilon_{b}}\rfloor)$
$\displaystyle\leq|S|\sum_{d=0}^{h_{0}}\lfloor\frac{\gamma^{d}}{\epsilon_{b}}\rfloor\>\>|\;\gamma^{h_{0}}\geq\epsilon_{b}$
$\displaystyle\leq|S|\sum_{d=0}^{h_{0}}\frac{\gamma^{d}}{\epsilon_{b}}$
$\displaystyle\hat{M}_{h_{0}}$
$\displaystyle\leq{\lceil\frac{1}{\epsilon_{b}(1-\gamma)}\rceil}^{|S|}$

From equation 7

$\mathcal{N}(\epsilon,\Delta^{|S|},d_{1})\leq\mathcal{N}(\epsilon/2,{[0,1]}^{|S|},d_{1})$

Combing the two results proves the Lemma $\square$

##  4 Finite Memory Planner

We now introduce our main algorithm, Finite Memory Planner (FMP). This is an
approximate online POMDP algorithm, with the following two main steps:

  1. 1.

Building an approximation $\mathcal{T}$ to the full belief tree, achieved here
by recursively calling Algorithm 1.

  2. 2.

Calculating the optimal action at the root and its value function using
equation 8 on $\mathcal{T}$.

The lookahead algorithm for building the belief tree $\mathcal{T}$ is
described in Algorithm 1.

Algorithm 1 FMP (Finite Memory Planner)
Parameters: Horizon $h_{0}$, discretization levels $\epsilon_{d}$

1: Input: Nodes with possible beliefs $b_{d}\in B_{d}$ at depth $d$,tolerances
$\epsilon_{d}$, current height $h$

2: if $h=h_{0}$ then

3: return []

4: end if

5: children = []

6: for all $b_{d}\in B_{d}$ do

7: for all actions a do

8: for all next observations $o_{d+1}$ do

9: $b_{d+1},o_{prob},r=$ BayesFilter$(b_{d},o_{d+1},a)$

10: $\hat{b}_{d+1}=NN(b_{d+1},A(\epsilon_{d+1}))$

11: children.append$([\hat{b}_{d+1},o_{prob},a,r])$

12: end for

13: end for

14: end for

15: FMP(children,$\epsilon_{d+1}$,$h+1$)

#### Building the tree.

The first stage is building a tree to a fixed depth, as follows.

  1. 1.

For each belief $b_{d}$ at depth $d$, action $a$ and observation $o$, we call
‘BayesFilter’ function at line 9, which uses equations 5 and 4 to compute its
true posterior and other attributes. They are later needed to compute the root
value and optimal action using backward induction (equation 8).

  2. 2.

At line 10, we replace the true posterior $b_{d+1}$ by its nearest neighbour
in the $A(\epsilon_{d})$ set which serves as its $\epsilon_{d}$-cover. This is
computationally cheap because of uniform discretization (refer Definition 4).

We remark that the computational cost of Bayesian inference (line 9) is a
shared burden for all POMDP solvers, and is not any more expensive for FMP,
than for other solvers. One possible advantage of FMP is that inference can be
quite flexible, e.g., one may use exact, variational inference or sampling
based inference at line 9.333Contrast this to POMCP, which is solely based on
Monte Carlo (particle filter) inference. Finally, we make a note that even
$\epsilon$-cover based algorithms SARSOP and PGVI are significantly different.
In contrast to FMP, both perform DFS to build the tree, thereby forgoing any
chance of building a minimal $\epsilon$-cover at each depth. Computationally,
FMP is cheaper as it does not need to cluster points at any level, while both
SARSOP and PGVI do.

#### Calculating the value.

In our scheme, the subsequent beliefs $b_{k+1}$ of any belief $b_{k}$ are
replaced by some approximations $\hat{b}_{k+1}$ in the $\epsilon_{k+1}$-cover.
Hence, the backwards induction equation is now

$V(b_{k})=\max_{a}\mathbb{E}_{\mathop{\mathbb{P}}\nolimits(r_{k+1},o_{k+1}|a,b_{k})}(\;r_{k+1}+V(\hat{b}_{k+1})\;)$
(8)

When the observation space is large we approximate the above by sampling
$o_{k+1}\sim\mathop{\mathbb{P}}\nolimits(o_{k+1}|a,b_{k})$ and
$r_{k+1}\sim\mathop{\mathbb{P}}\nolimits(r_{k+1}|a,b_{k})$.444In particular,
we sample $(o_{k+1},r_{k+1})$ in _RandomPOMDP_ environment but only
$(r_{k+1})$ in _RockSample_ environment since it has binary observation space.
Note that the algorithm returns the complete approximate value function for
every depth $k+1$ to the previous depth $k$.

##  5 Experiment

We compare with State-of-the-art planner POMCP [Silver and Veness, 2010]
mainly, on RandomPOMDP and RockSample environments. We also compare with PBVI
[Pineau et al., 2003] on RandomPOMDP and DESPOT [Ye et al., 2017] on
RockSample environment. The environment description is as follows:

  1. 1.

RandomPOMDP ($RP[s,sp]$): A random POMDP of given state, action and
observation space. The reward, transition and observation distributions is
generated using a sparsity parameter $sp\in\\{0.3,0.6,0.9\\}$ with
$\mathcal{R}_{\max}=1$. We experiment with different state space sizes $|S|$
and sparsity parameters. We fix $|Z|=|S|$ and $|A|=4$. Denote a particular
configuration by $RP[s,sp]$.

  2. 2.

RockSample ($RS[n,k]$): Environment proposed in [Smith and Simmons, 2012a],
its state space size ranges from 12 thousand to 7 million. Although the belief
space is very sparse, consisting of only 8 to 15 partially observed
quantities. RS[n,k] denotes a map of size $n\times n$, searched by a agent,
looking to collect good rocks between k possible total rocks. Its has a binary
observation space, with distance dependent observation likelihood, giving more
accurate observations as the agent moves closer to the rocks. Its goal is to
collect the good rocks and exit the map. The agent knows its location and that
of rocks, but is unaware of the state (good/bad) of each rock.

###  5.1 setup

We take discount $\gamma=0.95$ (standard) and run each experiment for $90$
steps555Inline with POMCP and DESPOT paper for RockSample, and sufficient for
RandomPOMDP (since $0.95^{90}\approx 0.01$). We perform 50 experiments for
each (Planner,POMDP) combination. We use POMDPy library [Emami et al., 2015]
for RockSample simulator and PyPOMDP library [Lu, 2018] to solve our generated
RandomPOMDP environment with POMCP and PBVI. While for the RockSample
environment, we simply report the published (and tuned) results for DESPOT and
POMCP (Table.1 in [Ye et al., 2017]). Like [Smith and Simmons, 2012b], we mask
the known states in RockSample environment. Finally, we also use GNU parallel
[Tange, 2018] to run parallel experiments.

#### Parameter tuning:

It is worth noting that there are minimal tuning parameters to begin with,
just $h_{0}$ and the sequences $\\{\epsilon_{d}\\}$. The solution quality is
fairly linear with fineness of discretization, and inversely proportional to
the error tolerance levels. In the experiments, the values are selected such
that computation remains within satisfactory online planning limit (5 second
for all environments up to RS[11,11], while 30 seconds for RS[15,15]). In
particular, the following values were selected:
$\epsilon_{d}=\\{\epsilon_{0}=\frac{1}{5},\ldots,\epsilon_{H}=\frac{1}{5}\\}$
or
$\\{\epsilon_{0}=\frac{1}{5},\epsilon_{1}=\frac{1}{4},\ldots,\epsilon_{H}=\frac{1}{2}\\}$
(depending on planning limit) and tree depth $H=5$ or $H=4$ (depending on
planning limit). Note that the $\epsilon_{d}$-sequence and $H$ chosen are
arbitrary and gave sufficiently good results without finetuning. For
RandomPOMDP, we the allocate maximum possible budget (simulation time) to
POMCP666We tuned for exploration constant $c\in\\{3,10,100\\}$. and PBVI.

The source code is available online and with supplementary material. The
general FMP planner is made available in a pseudo-code style implementation in
Python language.

###  5.2 Results and discussion

RS[n,k] | RS[7,8] | RS[11,11] | RS[15,15]
---|---|---|---
|S| | 12,544 | 247,808 | 7,372,800
|A| | 13 | 16 | 20
FMP | 20.24$\pm$0.85 | 21.49$\pm$0.63 | 15.67$\pm$0.37
POMCP∗ | 20.71$\pm$0.21 | 20.01$\pm$0.23 | 15.32$\pm$0.28
DESPOT∗ | 20.93$\pm$0.30 | 21.74$\pm$0.30 | 18.64$\pm$0.28
Table 1: Total Discounted Reward. ∗Reported values. RP[s,sp] | FMP | POMCP | PBVI
---|---|---|---
RP[100,0.9] | 6.02 | 5.22 | 5.82
RP[100,0.6] | 4.38 | 3.99 | 4.65
RP[100,0.3] | 2.42 | 1.57 | 2.40
RP[60,0.9] | 6.04 | 5.90 | 5.01
RP[60,0.6] | 4.65 | 4.94 | 4.23
RP[60,0.3] | 2.42 | 2.67 | 2.27
RP[30,0.9] | 6.36 | 6.37 | 6.26
RP[30,0.6] | 4.12 | 5.33 | 4.26
RP[30,0.3] | 3.40 | 3.53 | 3.74
Table 2: Total Discounted Reward.

The results are given in Table 1 and 2 for RockSample and RandomPOMDP777For
RandomPOMDP, the variance was low and we do not report it simply for clarity.
respectively. It is clear that the performance of FMP is competitive in all
settings, even with relatively lax parameters. For RockSample, FMP’s python
implementation is significantly handicapped compared to reported values of
POMCP and DESPOT, obtained by parameter tuning and with much faster C++
implementation. We believe FMP further would benefit from such speed
optimization. The RandomPOMDP performance gives some hints to FMP behaviour.
It is evident that FMP performs best when there are many beliefs to begin
with, i.e, it performs best in low sparsity setting. For high sparsity, it
performed competitively to the rest (for larger state sizes), and we guess it
misses out in other cases simple due to less lookahead (only $H=5$, again due
to planning limit and Python language implementation).

Finally, we contrast FMP with other planners by noting their tuning and
implementation complexity. For example, POMCP requires depth of tree, number
of particles used for inference, simulation budget, UCB constant, as
parameters. DESPOT in addition to all that, requires maintaining upper/lower
bounds on belief nodes. Most other solvers apply such heuristics as well.

##  6 Conclusion

We have a presented a simple algorithm for online planning in POMDPs that is
highly competitive with the state-of-the-art despite its apparent simplicity.
The main idea behind FMP was to design an algorithm where the inference and
planning are tightly coupled to minimise the estimation error. In particular,
the algorithm performs approximate inference with adaptive belief
discretization so as to upper bound the planning error, with the coarseness of
the approximation carefully tuned so as not to overly increase the size of the
planning tree. This makes the algorithm both theoretically attractive and
practical.

In future work, there is possibility to extend in multiple directions. we
would like to extend our method to continuous spaces, which can be done as the
proof strategy only needs kernel function over beliefs. For example, kernels
for the sufficient statistics of the belief over continuous state space could
be explored. The analysis could also be extended to get regret bounds in
Reinforcement Learning setting. Both directions may lead to potential solvers
that can handle a wide variety of problems.

## References

  * Bellman [1952] Richard Bellman.  On the theory of dynamic programming.  _Proceedings of the National Academy of Sciences of the United States of America_ , 38(8):716, 1952.
  * Emami et al. [2015] Patrick Emami, Alan J Hamlet, and Carl D Crane.  Pomdpy: An extensible framework for implementing partially-observable markov decision processes in python.  2015\.
  * Kearns et al. [2002] Michael Kearns, Yishay Mansour, and Andrew Y Ng.  A sparse sampling algorithm for near-optimal planning in large markov decision processes.  _Machine learning_ , 49(2-3):193–208, 2002.
  * Kurniawati and Yadav [2016] Hanna Kurniawati and Vinay Yadav.  An online pomdp solver for uncertainty planning in dynamic environment.  In _Robotics Research_ , pages 611–629. Springer, 2016.
  * Kurniawati et al. [2008] Hanna Kurniawati, David Hsu, and Wee Sun Lee.  Sarsop: Efficient point-based pomdp planning by approximating optimally reachable belief spaces.  In _Robotics: Science and systems_ , volume 2008. Zurich, Switzerland., 2008.
  * Lee et al. [2007] Wee Lee, Nan Rong, and David Hsu.  What makes some pomdp problems easy to approximate?  _Advances in neural information processing systems_ , 20:689–696, 2007.
  * Lu [2018] Di Lu.  Pypomdp.  https://github.com/namoshizun/PyPOMDP, 2018.
  * Pineau et al. [2003] Joelle Pineau, Geoff Gordon, Sebastian Thrun, et al.  Point-based value iteration: An anytime algorithm for pomdps.  In _IJCAI_ , volume 3, pages 1025–1032. Citeseer, 2003.
  * Ross et al. [2007] Stephane Ross, Brahim Chaib-draa, and Joelle Pineau.  Bayes-adaptive pomdps.  In _NIPS_ , pages 1225–1232, 2007.
  * Roy et al. [2005] Nicholas Roy, Geoffrey Gordon, and Sebastian Thrun.  Finding approximate pomdp solutions through belief compression.  _Journal of artificial intelligence research_ , 23:1–40, 2005.
  * Silver and Veness [2010] David Silver and Joel Veness.  Monte-carlo planning in large pomdps.  In _Advances in neural information processing systems_ , pages 2164–2172, 2010.
  * Smallwood and Sondik [1973] Richard D Smallwood and Edward J Sondik.  The optimal control of partially observable markov processes over a finite horizon.  _Operations research_ , 21(5):1071–1088, 1973\.
  * Smith and Simmons [2012a] Trey Smith and Reid Simmons.  Heuristic search value iteration for pomdps.  _arXiv preprint arXiv:1207.4166_ , 2012a.
  * Smith and Simmons [2012b] Trey Smith and Reid Simmons.  Point-based pomdp algorithms: Improved analysis and implementation.  _arXiv preprint arXiv:1207.1412_ , 2012b.
  * Tange [2018] Ole Tange.  Gnu parallel.  _DOI: https://doi. org/10.5281/zenodo_ , 1146014, 2018.
  * Ye et al. [2017] Nan Ye, Adhiraj Somani, David Hsu, and Wee Sun Lee.  Despot: Online pomdp planning with regularization.  _Journal of Artificial Intelligence Research_ , 58:231–266, 2017.
  * Zhang et al. [2012] Zongzhang Zhang, Michael Littman, and Xiaoping Chen.  Covering number as a complexity measure for pomdp planning and learning.  In _Twenty-Sixth AAAI Conference on Artificial Intelligence_ , 2012\.
  * Zhang et al. [2014] Zongzhang Zhang, David Hsu, and Wee Sun Lee.  Covering number for efficient heuristic-based pomdp planning.  In _International conference on machine learning_ , pages 28–36, 2014\.


# Modeling Nuclear Quantum Effects on Long Range Electrostatics in Nonuniform
Fluids

Richard C. Remsing  Department of Chemistry and Chemical Biology, Rutgers
University, Piscataway, NJ 08854

###### Abstract

Nuclear quantum effects play critical roles in a variety of molecular
processes, especially in systems that contain hydrogen and other light nuclei,
such as water. For water at ambient conditions, nuclear quantum effects are
often interpreted as local effects resulting from a smearing of the hydrogen
atom distribution. However, the orientational structure of water at interfaces
determines long range effects like electrostatics through the O-H bond
ordering that is impacted by nuclear quantum effects. In this work, I examine
nuclear quantum effects on long range electrostatics of water confined between
hydrophobic walls using path integral simulations. To do so, I combine
concepts from local molecular field (LMF) theory with path integral methods at
varying levels of approximation to develop an efficient and physically
intuitive approaches for describing long range electrostatics in nonuniform
quantum systems. Using these approaches, I show that quantum water requires
larger electrostatic forces to achieve the same level of interfacial screening
as the corresponding classical system. This work highlights subtleties of
electrostatics in nonuniform classical and quantum molecular systems, and the
methods presented here are expected to be of use to efficiently model nuclear
quantum effects in large systems.

##  I Introduction

The quantum mechanical behavior of atomic nuclei can have important effects on
molecular processes. Chemical kinetics, isotope effects, and the structure of
molecular systems is influenced by nuclear quantum effects (NQEs). Therefore,
NQEs must be adequately described in theoretical or computational models of
molecular systems, especially when light nuclei or low temperatures are
involved, because NQEs are prevalent at these conditions.

In water at ambient conditions, NQEs primarily impact hydrogen atoms and
affect the structure, dynamics, and chemical reactivity of liquid water
Markland and Ceriotti (2018); Ceriotti  _et al._ (2016); Marsalek and Markland
(2017); Wang, Ceriotti, and Markland (2014); Marx, Tuckerman, and Parrinello
(2000); Tuckerman, Marx, and Parrinello (2002); Xu  _et al._ (2021); Raugei
and Klein (2003); Tang  _et al._ (2021); Morrone and Car (2008); Ceriotti  _et
al._ (2013); Marx, Tuckerman, and Parrinello (2000); Zhuang  _et al._ (2019);
Paesani  _et al._ (2010); Spura  _et al._ (2015); Ojha, Henao, and Kühne
(2018); Litman  _et al._ (2018). Quantum mechanical treatment of water is
often considered to smear out the positions of the hydrogens, although non-
trivial electronic quantum effects can couple to the nuclear fluctuations in
ab initio treatments. In all of these situations, NQEs are considered mainly
local effects, impacting the hydrogen bond structure and dynamics of water.
However, any changes in the hydrogen bond structure can also impact non-local
properties of water, especially in nonuniform systems where changes in liquid
structure are magnified by interfaces. In particular, water preferentially
adopts specific orientations at interfaces Lee, McCammon, and Rossky (1984);
Lee and Rossky (1994); Lum, Chandler, and Weeks (1999); Chandler (2005);
Remsing and Weeks (2013); Jamadagni, Godawat, and Garde (2011);
Giovambattista, Rossky, and Debenedetti (2012); Rego and Patel (2022). These
orientational preferences control interfacial electrostatics and are dictated
by the water hydrogen bond network. This hydrogen bond network, in turn, is
sensitive to NQEs. Ultimately, this suggests that NQEs can indirectly impact
orientational properties of water in nonuniform systems and the consequent
electrostatic properties of the system. Here, I investigate NQEs on the long
range electrostatic properties of water in nonuniform systems using molecular
simulations and quantum statistical mechanics combined with liquid state
theory.

A standard approach to modeling NQEs is path integral molecular dynamics
(PIMD) simulations Parrinello and Rahman (1984); Tuckerman  _et al._ (1996);
Ceperley (1995), as well as the closely related ring polymer molecular
dynamics (RPMD) simulation approach Habershon  _et al._ (2013). In these
approaches, each quantum particle is replaced by an isomorphic classical _ring
polymer_ composed of $P$ beads (or monomers). For large enough $P$, one can
determine exact statistical averages of static quantities. Using RPMD, one can
model approximate quantum dynamics by evolving the system in time using
classical dynamics, but in the extended phase space of the ring polymer system
Habershon  _et al._ (2013). Modeling the ring polymer system amounts to
simulating $P$ coupled copies or replicas of the classical system, enabling
the use and extension of standard computational algorithms.

Replicating the system $P$ times leads to significant computational cost,
especially when $P$ is large. This is especially true for long range
interactions, like electrostatics, which are already expensive in the
classical system ($P=1$). However, by judiciously separating electrostatic
interactions into short and long range components, the number of beads needed
to evaluate the long range part can be significantly reduced, even down to the
$P=1$ limit Markland and Manolopoulos (2008a, b). This approach — ring polymer
contraction (RPC) — evaluates electrostatic interactions on a subset of beads
and distributes the result over all the beads of the polymer, reducing the
number of evaluations of long range interactions from $P$ to $P^{\prime}<P$.
Here, I focus on the limit where long range interactions are evaluated only at
the centroid of each ring polymer $P^{\prime}=1$.

While RPC can significantly reduce the computational cost of PIMD and RPMD
simulations, this approach is still plagued with the usual problems associated
with evaluating long range electrostatics — Ewald summations are expensive and
conceptually difficult. The conceptual issues are particularly problematic,
because Ewald sums and other lattice summation techniques can often lead to
significant geometry-dependent finite size effects, as well as other artifacts
associated with periodic boundary conditions (PBCs) Figueirido, Del Buono, and
Levy (1995); Hub  _et al._ (2014); Pan, Yi, and Hu (2017). An appealing
alternative to lattice sums is local molecular field (LMF) theory and related
developments Rodgers and Weeks (2008a); Remsing, Liu, and Weeks (2016); Hu
(2014); Gao, Remsing, and Weeks (2020, 2023). LMF theory relies on a
physically intuitive separation of electrostatic interactions into short and
long range components, much like RPC, and replaces all electrostatic
interactions by the short range part. The averaged effects of the long range
interactions are then modeled through an effective external field that is
chosen to produce the same structure as the full system with long range
interactions. Accurate structure, thermodynamics, and even dynamics can be
efficiently and accurately predicted by LMF simulations, making it a useful
alternative to lattice summations.

In this work, I combine LMF theory and the path integral isomorphism to
develop RPC-based approximations and model nuclear quantum effects in water
confined between model nonpolar walls. After a brief review of RPC and LMF
theory, I discuss strategies for combining the two to obtain approaches for
the efficient and accurate predictions of water structure between model
hydrophobic walls. These approaches can aid in reducing the cost of PIMD and
RPMD simulations while also helping to build physical intuition regarding the
effects of long range interactions in heterogeneous quantum systems.

##  II Theory

###  II.1 Ring Polymer Contraction

In standard PIMD and RPMD simulation approaches, Feynman’s path integral
formulation of quantum mechanics is used to model a system of $N$
distinguishable particles with Hamiltonian
$\mathcal{H}({\textbf{P}},{\textbf{R}})$, where P and R represent the momenta
and positions of all $N$ particles at a single point in phase space Feynman
and Hibbs (2005); Chandler and Wolynes (1981); Habershon  _et al._ (2013). The
partition function, $\mathcal{Z}$, of this system can be approximated by that
of an isomorphic system composed of classical ring polymers, each composed of
$P$ beads,

$\displaystyle\mathcal{Z}$
$\displaystyle\approx\mathcal{Z}_{P}=\frac{1}{(2\pi\hbar)^{3NP}}\int
d{\textbf{P}}\int
d{\textbf{R}}e^{-\beta\mathcal{H}_{P}({\textbf{P}},{\textbf{R}})/P},$ (1)

which becomes an equality in the $P\rightarrow\infty$ limit. Here,
$\beta=(k_{\rm B}T)^{-1}$ and the ring polymer Hamiltonian is

$\displaystyle\mathcal{H}_{P}({\textbf{P}},{\textbf{R}})$
$\displaystyle=\sum_{i=1}^{N}\sum_{\alpha=1}^{P}\left(\frac{\left|\mathbf{p}_{i}^{(\alpha)}\right|^{2}}{2m_{i}}+\frac{1}{2}m_{i}\omega_{P}^{2}\left|\mathbf{r}_{i}^{(\alpha)}-\mathbf{r}_{i}^{(\alpha-1)}\right|^{2}\right)$
$\displaystyle+\sum_{\alpha=1}^{P}V({\textbf{R}}^{(\alpha)}),$ (2)

where the last term is the sum of the potential energy over all beads,
$\mathbf{r}_{i}^{(\alpha)}$ and $\mathbf{p}_{i}^{(\alpha)}$ are the position
and momentum of bead $\alpha$ of site $i$, $\mathbf{R}^{(\alpha)}$ represents
the position vector for bead $\alpha$ of all $N$ particles in a single
configuration, $m_{i}$ is the mass of site $i$, and $\omega_{P}=P/\beta\hbar$
is the spring constant. In general, this Hamiltonian can involve many-body
interactions, but here I focus on one- and two-body interactions,

$V({\textbf{R}}^{(\alpha)})=\sum_{i=1}^{N}\phi(\mathbf{r}_{i}^{(\alpha)})+\sum_{i=1}^{N-1}\sum_{j=i+1}^{N}w\left(\left|\mathbf{r}_{i}^{(\alpha)}-\mathbf{r}_{j}^{(\alpha)}\right|\right),$
(3)

where the one-body potential $\phi(\mathbf{r})$ arises from an external field.
The two-body potential $w(r)$ involves Lennard-Jones (LJ) and electrostatic
interactions,

$w(r_{ij})=u_{\rm LJ}(r_{ij})+q_{i}q_{j}v(r_{ij}),$ (4)

where $r_{ij}=\left|\mathbf{r}_{i}-\mathbf{r}_{j}\right|$, $u_{\rm LJ}(r)$ is
the LJ potential, and $v(r)=1/r$.

Lennard-Jones interactions are typically truncated at some distance and their
effects beyond the cutoff accounted for with a correction. In contrast,
electrostatic interactions are long-ranged and are typically evaluated using
lattice summation techniques like Ewald summation. Lattice summations are
generally expensive and significantly increase the cost of PI simulations when
the system is replicated $P$ times to construct the ring polymers. Ring
polymer contraction (RPC) can be used to reduce this cost by lowering the
number of beads on which long range interactions need to be evaluated Markland
and Manolopoulos (2008a, b).

RPC splits the interparticle interactions into two components, $V(r)=V_{\rm
S}(r)+V_{\rm L}(r)$, where $V_{\rm S}(r)$ is the short range part of the
potential and $V_{\rm L}(r)$ is the long range part of the potential. The
splitting is chosen such that $V_{\rm L}(r)$ varies slowly on the length scale
of the largest ring polymer, estimated by the ensemble averaged radius of
gyration. For water at ambient conditions, this length scale is close to the
free particle limit, $\lambda_{\rm free}=\sqrt{\beta\hbar^{2}/4m}$. When the
potential varies slowly on the scale of $\lambda_{\rm free}$, the total
interaction between two ring polymers can be approximated by $P$ times the
interaction between their centroids,

$\sum_{\alpha=1}^{P}V_{\rm
L}\left(\left|\mathbf{r}_{i}^{(\alpha)}-\mathbf{r}_{j}^{(\alpha)}\right|\right)\approx
PV_{\rm L}(\left|\bar{\mathbf{r}}_{i}-\bar{\mathbf{r}}_{j}\right|)$ (5)

where $\bar{\mathbf{r}}_{i}$ is the centroid of ring polymer $i$. Forces can
be readily evaluated following previous work Markland and Manolopoulos (2008a,
b). The centroid RPC approximation, Eq. 5, significantly reduces the cost of
evaluating long range interactions without sacrificing accuracy Markland and
Manolopoulos (2008a, b). Using local molecular field theory, summarized in the
next section, we can further reduce the cost of evaluating long range
interactions.

###  II.2 Local Molecular Field Theory

LMF theory accounts for the averaged effects of long range electrostatics with
a renormalized or effective electrostatic potential Rodgers and Weeks (2008a);
Remsing, Liu, and Weeks (2016). The first step in determining this potential
is to separate intermolecular Coulomb interactions into short and long range
components. For LMF theory to be valid, long range interactions must vary
slowly over typical nearest neighbor distances. As such, LMF theory separates
the $1/r$ portion of the Coulomb potential according to

$\displaystyle v(r)$ $\displaystyle=\frac{1}{r}$ (6)
$\displaystyle=\frac{\mathop{\rm erfc}(r/\sigma)}{r}+\frac{\mathop{\rm
erf}(r/\sigma)}{r}$ (7) $\displaystyle\equiv v_{0}(r)+v_{1}(r),$ (8)

where $\sigma$ is the LMF smoothing length that is on the order of
intermolecular correlations, $v_{0}(r)$ is the short-range component of the
electrostatic interactions, and $v_{1}(r)$ is the long-range component. For
water at ambient conditions, previous work has shown that $\sigma\geq 3$ Å
Rodgers, Hu, and Weeks (2011), and here I use a conservative value of
$\sigma=4.5$ Å.

In LMF theory, the full model is replaced by its Gaussian-truncated (GT)
counterpart, in which $v(r)$ is replaced by $v_{0}(r)$ for all sites in the
system. The averaged effects of fluid-fluid long range electrostatics are
taken into account through the renormalized electrostatic potential

$\displaystyle\mathcal{V}_{\rm R}(\mathbf{r})$
$\displaystyle=\mathcal{V}(\mathbf{r})+\int
d\mathbf{r}^{\prime}\rho^{q}(\mathbf{r}^{\prime})v_{1}(\left|\mathbf{r}-\mathbf{r}^{\prime}\right|)$
(9) $\displaystyle\equiv\mathcal{V}(\mathbf{r})+\mathcal{V}_{\rm
S}(\mathbf{r}),$ (10)

where $\mathcal{V}(\mathbf{r})$ is the external electrostatic potential
present in the full system, equal to zero for the systems studied here,
$\rho^{q}(\mathbf{r})$ is the ensemble averaged singlet charge density of the
system, given by

$\rho^{q}(\mathbf{r})=\left<\rho^{q}(\mathbf{r};{\textbf{R}})\right>,$ (11)

and $\rho^{q}(\mathbf{r};{\textbf{R}})$ is the charge density operator
evaluated in configuration R. The charge density operator is defined in the
next section. In general, the external field contains can be split into short
and long range parts,
$\mathcal{V}(\mathbf{r})=\mathcal{V}_{0}(\mathbf{r})+\mathcal{V}_{1}(\mathbf{r})$,
such that the LMF can be rewritten in a way that isolates the long range
interactions as $\mathcal{V}_{\rm
R}(\mathbf{r})=\mathcal{V}_{0}(\mathbf{r})+\mathcal{V}_{\rm R1}(\mathbf{r})$.

Equation 10 is self-consistent and can be solved through brute-force
simulations or using methods based in linear response theory Hu and Weeks
(2010). Simulating the GT model in the presence of $\mathcal{V}_{\rm
R}(\mathbf{r})$ yields structures in agreement with the full system Rodgers
and Weeks (2008b); Baker III, Rodgers, and Weeks (2020); Cox (2020), and
thermodynamics can be obtained with previously-derived corrections Rodgers and
Weeks (2009); Remsing, Rodgers, and Weeks (2011); Remsing, Liu, and Weeks
(2016); Remsing and Weeks (2018); Gao, Remsing, and Weeks (2020).

###  II.3 Solving the Local Molecular Field Equation for Quantum Systems

The LMF potential can be obtained by writing the ring polymer expression for
the charge density operator of the system in configuration R as Craig and
Manolopoulos (2004); Habershon  _et al._ (2013)

$\displaystyle\rho^{q}(\mathbf{r};{\textbf{R}})$
$\displaystyle=\frac{1}{P}\sum_{\alpha=1}^{P}\sum_{i=1}^{N}q_{i}^{(\alpha)}\delta\left(\mathbf{r}-\mathbf{r}_{i}^{(\alpha)}({\textbf{R}})\right)$
(12)
$\displaystyle=\frac{1}{P}\sum_{\alpha=1}^{P}\rho^{q\alpha}(\mathbf{r};{\textbf{R}}).$
(13)

Using this expression for the charge density, the LMF potential is given by

$\displaystyle\mathcal{V}_{\rm R}(\mathbf{r})$
$\displaystyle=\mathcal{V}(\mathbf{r})+\frac{1}{P}\sum_{\alpha=1}^{P}\int
d\mathbf{r}^{\prime}\rho^{q\alpha}(\mathbf{r}^{\prime})v_{1}(\left|\mathbf{r}-\mathbf{r}^{\prime}\right|)$
(14)
$\displaystyle\equiv\mathcal{V}(\mathbf{r})+\frac{1}{P}\sum_{\alpha=1}^{P}\mathcal{V}_{\rm
S}^{(\alpha)}(\mathbf{r})$ (15)

I will refer to this as the path integral local molecular field (PI-LMF).

The PI-LMF Equation 15 must be solved self-consistently. A self-consistent
solution for $\mathcal{V}_{\rm R}(\mathbf{r})$ can be found by iterating with
molecular simulations, but this is expensive, especially for quantum systems.
This can be circumvented by iterating to self-consistency using linear
response theory (LRT) instead of simulations to predict the density induced by
each field, as described by Hu and Weeks for classical fluids Hu and Weeks
(2010). I now describe how to extend this LRT approach to solving the LMF
equation to path integral models.

In a system of quantum particles described within the path integral formalism,
the Hamiltonian is replaced by the corresponding approximation involving ring
polymers composed of $P$ beads each. Ignoring the momenta — we are only
concerned with configurational averages here — the path integral Hamiltonian
can be written as

$\displaystyle\mathcal{H}_{P}({\textbf{R}})$
$\displaystyle=\frac{1}{P}\sum_{\alpha=1}^{P}\left[U_{0}^{(\alpha)}({\textbf{R}})+\Phi_{0}^{(\alpha)}({\textbf{R}})+\Phi_{\rm
R1}^{(\alpha)}({\textbf{R}})\right]$ (16)
$\displaystyle=\mathcal{H}_{P,0}({\textbf{R}})+\frac{1}{P}\sum_{\alpha=1}^{P}\Phi_{\rm
R1}^{(\alpha)}({\textbf{R}}),$ (17)

where the bond potentials between neighboring beads are included in
$U_{0}^{(\alpha)}({\textbf{R}})$. The Hamiltonian
$\mathcal{H}_{P,0}({\textbf{R}})=\frac{1}{P}\sum_{\alpha=1}^{P}\left[U_{0}^{(\alpha)}({\textbf{R}})+\Phi_{0}^{(\alpha)}({\textbf{R}})\right]$
represents the purely short ranged (reference) system,
$U_{0}^{(\alpha)}({\textbf{R}})$ is the total potential energy of the short
range pair interactions for bead $\alpha$, and
$\Phi_{0}^{(\alpha)}({\textbf{R}})$ is the corresponding total short range
one-body potential energy. The total potential energy of the long range
interactions for each bead is contained within $\Phi_{\rm
R1}^{(\alpha)}({\textbf{R}})=\int
d\mathbf{r}\rho^{q\alpha}(\mathbf{r};{\textbf{R}})\mathcal{V}_{\rm
R1}(\mathbf{r})$.

Using this separation of the Hamiltonian into short and long range components,
the average charge density can be written as an ensemble average in the short
range system according to

$\displaystyle\rho^{q}(\mathbf{r})$
$\displaystyle=\left<\frac{1}{P}\sum_{\alpha=1}^{P}\rho^{q(\alpha)}(\mathbf{r};{\textbf{R}})\right>$
$\displaystyle=\frac{\left<\frac{1}{P}\sum_{\alpha=1}^{P}\rho^{q(\alpha)}(\mathbf{r};{\textbf{R}})e^{-\frac{\beta}{P}\sum_{\gamma=1}^{P}\Phi_{\rm
R1}^{(\gamma)}({\textbf{R}})}\right>_{0}}{\left<e^{-\frac{\beta}{P}\sum_{\alpha=1}^{P}\Phi_{\rm
R1}^{(\alpha)}({\textbf{R}})}\right>_{0}}.$ (18)

Now, noting that the instantaneous bead-averaged field energy is

${\Phi}_{\rm R1}({\textbf{R}})=\frac{1}{P}\sum_{\alpha=1}^{P}\Phi_{\rm
R1}^{(\alpha)}({\textbf{R}}),$ (19)

we can rewrite the charge density as an average over configurations in the
short range system,

$\rho^{q}(\mathbf{r})=\frac{\left<\rho^{q}(\mathbf{r};{\textbf{R}})e^{-\beta{\Phi}_{\rm
R1}({\textbf{R}})}\right>_{0}}{\left<e^{-\beta{\Phi}_{\rm
R1}({\textbf{R}})}\right>_{0}}.$ (20)

We can then linearize this expression for the charge density of quantum
particles to obtain the linear response approximation

$\rho^{q}(\mathbf{r})\approx\left<\rho^{q}(\mathbf{r};{\textbf{R}})\right>_{0}-\beta\left<\delta\rho^{q}(\mathbf{r};{\textbf{R}})\delta{\Phi}_{\rm
R1}({\textbf{R}})\right>_{0},$ (21)

where
$\delta\rho^{q}(\mathbf{r};{\textbf{R}})=\rho^{q}(\mathbf{r};{\textbf{R}})-\left<\rho^{q}(\mathbf{r};{\textbf{R}})\right>_{0}$
and $\delta{\Phi}_{\rm R1}({\textbf{R}})=\Phi_{\rm
R1}({\textbf{R}})-\left<\Phi_{\rm R1}({\textbf{R}})\right>_{0}$. Equation 21
is analogous to the classical result of Hu and Weeks Hu and Weeks (2010),
except the classical operators are replaced by their bead-averaged
counterparts Craig and Manolopoulos (2004); Habershon  _et al._ (2013).
Equations 15 and 21 are the main results of this section and are used to
obtain a self-consistent solution to the LMF equation through iteration.

The benefits of using Eq. 21 instead of the more traditional form of the
linear response approximation are magnified in path integral treatments of
quantum systems. In this case, the more traditional expression

$\rho^{q}(\mathbf{r})\approx\left<\rho^{q}(\mathbf{r};{\textbf{R}})\right>_{0}-\frac{\beta}{P^{2}}\sum_{\alpha=1}^{P}\sum_{\gamma=1}^{P}\int
d\mathbf{r}^{\prime}\phi_{\rm
R1}(\mathbf{r}^{\prime})\chi^{qq}_{\alpha\gamma}(\mathbf{r},\mathbf{r}^{\prime})$
(22)

involves pair correlations between all beads in the system, including those on
different slices of imaginary time ($\alpha$ and $\gamma$), through the
quantum charge-charge linear response function Chandler and Wolynes (1981);
Shinoda, Miura, and Okazaki (2001); Sesé (2002)

$\chi^{qq}_{\alpha\gamma}(\mathbf{r},\mathbf{r}^{\prime})=\left<\delta\rho^{q(\alpha)}(\mathbf{r};{\textbf{R}})\delta\rho^{q(\gamma)}(\mathbf{r}^{\prime};{\textbf{R}})\right>.$
(23)

In addition to the difficulties of evaluating a six-dimensional correlation
function in a nonuniform system, the need to evaluate correlations between
different points in imaginary time further increases the expense of using Eq.
22. Because of these difficulties, the much more efficient Eq. 21 is preferred
to solve the LMF equation for path integral models.

###  II.4 Combining Local Molecular Field Theory and Ring Polymer Contraction

While the solution to the LMF equation in the previous section can be obtained
from simulation results and linear response theory, the slowly-varying nature
of the long range potentials suggests that simpler approximations can be
exploited to more efficiently solve the LMF equation. One approach is to
combine LMF theory with RPC.

A fundamental concept in both LMF theory and RPC is the separation of
interaction potentials into short and long range components based on physical
principles. In RPC, electrostatic interactions are separated so that the long
range component is slowly-varying over the size of the ring polymer Markland
and Manolopoulos (2008a, b). In LMF theory, electrostatic interactions are
separated so that the long range component is uniformly slowly-varying over
typical nearest-neighbor distances (or a correlation length) in the liquid
Rodgers and Weeks (2008a); Remsing, Liu, and Weeks (2016). In liquids like
water at ambient conditions, the typical separation length scales are similar,
and I will indicate this by the LMF smoothing length $\sigma$.

I follow these principles and use the typical LMF smoothing length of
$\sigma=4.5$ Å to separate the potential into short and long range components,
$v_{0}(r)$ and $v_{1}(r)$, respectively, as described in the previous section.
In RPC, the electrostatic (pair) potential $v_{1}(r)$ is evaluated between
centroid positions. By combining RPC and LMF, long range pair interactions are
completely removed, $V_{\rm L}=0$. The averaged effects of long range
interactions are instead accounted for via the effective one-body
electrostatic (LMF) potential $\mathcal{V}_{\rm R}(\mathbf{r})$, and RPC can
be used to evaluate the long range part of the LMF potential at centroid
positions only;

$\sum_{i=1}^{N-1}\sum_{j>i}^{N}PV_{\rm
L}(\left|\bar{\mathbf{r}}_{i}-\bar{\mathbf{r}}_{j}\right|)\rightarrow\sum_{i=1}^{N}P\mathcal{V}_{\rm
R1}(\bar{\mathbf{r}}_{i}).$ (24)

This strategy results in the LMF-RPC scheme for evaluating long range
interactions.

To determine the effective field $\mathcal{V}_{\rm R}(\mathbf{r})$ within the
LMF-RPC scheme using Eq. 21, the long range potential energy is evaluated only
at the location of the centroid. Defining the centroid charge density
operator,

$\displaystyle\bar{\rho}^{q}(\mathbf{r};{\textbf{R}})=\sum_{i=1}^{N}q_{i}\delta\left(\mathbf{r}-\bar{\mathbf{r}}_{i}({\textbf{R}})\right),$
(25)

where $\bar{\mathbf{r}}_{i}$ is the position of the centroid of particle $i$,
the total long range potential energy within the LMF-RPC approximation is

$\Phi_{\rm R1}^{\rm RPC}({\textbf{R}})=\int
d\mathbf{r}\bar{\rho}^{q}(\mathbf{r};{\textbf{R}})\mathcal{V}_{\rm
R1}(\mathbf{r}).$ (26)

This corresponds to evaluating the field $\mathcal{V}_{\rm R1}(\mathbf{r})$,
determined using all beads, at the location of the centroids only. As a
result, the linear response approximation for the LMF-RPC charge density is

$\rho^{q}(\mathbf{r})\approx\left<\rho^{q}(\mathbf{r};{\textbf{R}})\right>_{0}-\beta\left<\delta\rho^{q}(\mathbf{r};{\textbf{R}})\delta{\Phi}^{\rm
RPC}_{\rm R1}({\textbf{R}})\right>_{0},$ (27)

with an analogous expression for the centroid charge density. The converged
charge density and $\mathcal{V}_{\rm R}(\mathbf{r})$ are obtained within the
LMF-RPC approximation by iterating Eqs. 15 and 27 to self consistency.

###  II.5 Centroid Approximation

The LMF-RPC approach involves evaluating the charge density and centroid
charge density, as well as the full $\mathcal{V}_{\rm R}(\mathbf{r})$. Such a
complicated approach might not be necessary. Because $v_{1}(r)$ essentially
smears the charge distribution over the length scale $\sigma$, one might
anticipate that a completely centroid approximation to the long range
interactions should be reasonable when $\sigma$ is larger than the typical
size of the ring polymers in the system. To determine the effective field via
Eq. 15 within the centroid approximation, the charge density operator
appearing in Eq. 21 is approximated everywhere by the centroid charge density
operator,
$\rho^{q}(\mathbf{r};{\textbf{R}})\approx\bar{\rho}^{q}(\mathbf{r};{\textbf{R}})$.
In this case, the centroid approximation to the LMF potential is
$\mathcal{V}_{\rm R}(\mathbf{r})\approx\bar{\mathcal{V}}_{\rm R}(\mathbf{r})$,
where

$\displaystyle\bar{\mathcal{V}}_{\rm R}(\mathbf{r})$
$\displaystyle=\mathcal{V}(\mathbf{r})+\int
d\mathbf{r}^{\prime}\bar{\rho}^{q}(\mathbf{r};{\textbf{R}})v_{1}(\left|\mathbf{r}-\mathbf{r}^{\prime}\right|)$
(28) $\displaystyle\equiv\mathcal{V}(\mathbf{r})+\bar{\mathcal{V}}_{\rm
S}(\mathbf{r}).$ (29)

The centroid approximation is then inserted into Eq. 21 to iterate the LMF
equation to self-consistency in conjunction with the linear response
approximation for the centroid charge density

$\bar{\rho}^{q}(\mathbf{r})=\left<\bar{\rho}^{q}(\mathbf{r};{\textbf{R}})\right>_{0}-\beta\left<\delta\bar{\rho}^{q}(\mathbf{r};{\textbf{R}})\delta\bar{\Phi}_{\rm
R1}({\textbf{R}})\right>_{0},$ (30)

where $\bar{\Phi}_{\rm R1}({\textbf{R}})=\int
d\mathbf{r}\bar{\rho}^{q}(\mathbf{r};{\textbf{R}})\bar{\mathcal{V}}_{\rm
R1}(\mathbf{r})$ is the instantaneous energy from the centroid approximation
to the long range field evaluated at the centroid positions.

###  II.6 Feynman-Kleinert Approximation

The LMF-RPC and centroid approaches reduce the number of sites needed to
evaluate $\mathcal{V}_{\rm R}(\mathbf{r})$, but both still require a path
integral simulation in the purely short range GT system to evaluate ensemble
averages. Instead, we could first model a classical ($P=1$) GT system and use
the Feynman-Kleinert (FK) procedure to estimate the quantum LMF from its
classical counterpart Feynman and Kleinert (1986). This second approximation,
here called the FK approximation, in essence, corresponds to approximating a
quantum observable by a Gaussian smoothing of its classical counterpart
Feynman and Kleinert (1986). This can be used to determine the LMF for the
quantum system.

First, I determine the LMF potential for the classical system,
$\mathcal{V}_{\rm R}^{\rm cl}(\mathbf{r})$, by self-consistent iteration using
the classical linear response approximation Hu and Weeks (2010). The linear
response approximation is then used to predict the oxygen and hydrogen site
densities in the LMF system, $\rho_{\rm O}(\mathbf{r})$ and $\rho_{\rm
H}(\mathbf{r})$, respectively. I then smooth these densities over the lengths
$l_{\rm O}$ and $l_{\rm H}$ to convert the classical charge density into an
approximation of the quantum charge density,

$\displaystyle\rho^{ql}(\mathbf{r})$ $\displaystyle=\int
d\mathbf{r}^{\prime}\big{[}q_{\rm O}\rho_{\rm
O}(\mathbf{r}^{\prime})\rho_{G}(\left|\mathbf{r}-\mathbf{r}^{\prime}\right|;l_{\rm
O})$ $\displaystyle+q_{\rm H}\rho_{\rm
H}(\mathbf{r}^{\prime})\rho_{G}(\left|\mathbf{r}-\mathbf{r}^{\prime}\right|;l_{\rm
H})\big{]},$ (31)

where

$\rho_{G}(r;l)=\frac{1}{l^{3}\pi^{3/2}}e^{-r^{2}/l^{2}}$ (32)

is a spherical Gaussian of width $l$. Physically, $l$ corresponds to the
average size of a ring polymer, quantified by its radius of gyration, for
example. Because of the different masses of oxygen and hydrogen, and
consequently different spreads, we need different smoothing lengths for each,
$l_{\rm O}$ and $l_{\rm H}$, respectively, and we need to separate the charge
density into its components from oxygen and hydrogen sites. Here, I
approximate the size of the ring polymers by their free particle values,
$l_{\rm O}\approx\lambda_{\rm free,O}\approx 0.05$ Å and $l_{\rm
H}\approx\lambda_{\rm free,H}\approx 0.2$ Å. This crude approximation is
reasonable for water at ambient conditions, where the average radius of
gyration of the ring polymers is close to their free particle values Markland
and Manolopoulos (2008b).

After the densities are smoothed to account for nuclear quantum effects within
the FK approximation, I then perform a second smoothing over the length scale
$\sigma$ by convoluting the quantum charge density with $v_{1}$. The resulting
FK approximation to the LMF potential is $\mathcal{V}_{\rm
R}(\mathbf{r})\approx\mathcal{V}_{\rm R}^{\rm FK}(\mathbf{r})$, where

$\displaystyle\mathcal{V}_{\rm R}^{\rm FK}(\mathbf{r})$
$\displaystyle=\mathcal{V}(\mathbf{r})+\int
d\mathbf{r}^{\prime}\rho^{ql}(\mathbf{r}^{\prime})v_{1}(\left|\mathbf{r}-\mathbf{r}^{\prime}\right|)$
(33) $\displaystyle\equiv\mathcal{V}(\mathbf{r})+\mathcal{V}_{\rm S}^{\rm
FK}(\mathbf{r}).$ (34)

In water at room temperature, the long range electrostatic interactions vary
much more slowly than the ring polymers for the nuclei, $\sigma\gg\l_{\rm
H}>l_{\rm O}$. Therefore, one may anticipate that any deficiencies in the FK
approximation — deviations of the ring polymer from a spherical Gaussian —
will be washed out by smoothing over $\sigma$, and the FK approximation to the
LMF potential will be reasonably accurate at these conditions.

Figure 1: (a) Charge density, (b) Electrostatic potential, $\Phi(z)$, and (c)
Gaussian smoothed charge density for the full, Gaussian-truncated (GT), and
PI-LMF mimic systems for $P=32$ beads. Results are shown for the left wall
only.

##  III Results and Discussion

I demonstrate the utility of the LMF-RPC scheme using what has become the
canonical system for examining LMF theory-based methods — water confined
between idealized smooth hydrophobic walls. Near the walls, the effects of
long range interactions do not cancel, as they do in bulk, but instead create
forces near the water-wall interface. The averaged effects of long range
electrostatics provides a torque on interfacial water molecules, resulting in
dipolar screening and the proper orientational ordering of water molecules
near the wall Rodgers and Weeks (2008b); Baker III, Rodgers, and Weeks (2020);
Cox (2020). Neglecting long range electrostatics in this system results in
over-orientation of interfacial water molecules and a non-zero electric field
in the bulk due to the absence of dielectric screening in the purely short
range system. Here, I show that the same general physics arises in path
integral representations of water confined between hydrophobic walls and that
LMF theory adapted to path integral methods can account for the averaged
effects of long range electrostatics.

###  III.1 Ring Polymer Local Molecular Field Theory can Account for Long
Range Electrostatics

The role of long range electrostatics in determining the structure of confined
water can be observed by comparing simulation results obtained with the Full
and truncated models, Fig. 1. The charge density of water is similar in the
Full, GT, and PI-LMF systems. However, previous work has shown that
differences relevant to long range electrostatics are often hidden under
large, atomic-scale fluctuations in the charge density. Rodgers and Weeks
(2008b, a); Remsing and Weeks (2015a) These differences influence the
orientational preferences of interfacial water. In the GT systems, O-H bonds
point toward the wall more than in the full systems, consistent with
expectations from classical simulations. The orientation of interfacial water
molecules alters the electrostatic potential, given by

$\Phi(z)=-\int_{-\infty}^{z}dz^{\prime}\int_{-\infty}^{z^{\prime}}dz^{\prime\prime}\rho^{q}(z^{\prime\prime}).$
(35)

The resulting electrostatic potential determined from GT configurations does
not plateau in the bulk region. Including long range electrostatics, through
Ewald summation (Full) or through the PI-LMF potential, corrects this
behavior, resulting in less orientational ordering of water at the interface
and the expected plateau of the electrostatic potential in the bulk region.

Figure 2: (a) Charge density, (b) Electrostatic potential, $\Phi(z)$, (c)
Gaussian smoothed charge density, and (d) gradient of the LMF potential
(negative of force) for the various methods of solving the LMF equation: the
“exact” PI-LMF, LMF-RPC, the centroid approximation, and the Feynman-Kleinert
approximation (FK-LMF), all for $P=32$ beads. Results are shown for the left
wall only. Vertical dashed line shows the location of the wall, $z=0$.

The long range part of the LMF potential satisfies a Poisson equation,

$\nabla^{2}\mathcal{V}_{\rm
R1}(\mathbf{r})=-\frac{4\pi}{\varepsilon}\rho^{q\sigma}(\mathbf{r}),$ (36)

involving the Gaussian-smoothed charge density Rodgers and Weeks (2008a)

$\rho^{q\sigma}(\mathbf{r})=\int
d\mathbf{r}^{\prime}\rho^{q}(\mathbf{r}^{\prime})\rho_{G}(\left|\mathbf{r}-\mathbf{r}^{\prime}\right|;\sigma),$
(37)

which is shown in Fig. 1c for the Full, GT, and LMF systems. The Gaussian-
smoothed charge density therefore represents the portion of the charge density
that is relevant for the long range response of water, and one can think of
$\rho^{q\sigma}(\mathbf{r})$ as a macroscopic charge density Robinson (1973);
Zangwill (2013); Remsing and Weeks (2015b). The Full $\rho^{q\sigma}(z)$
displays a (macroscopic) dipole layer at the interface. In contrast, GT water
overpolarizes and produces a large positive peak in $\rho^{q\sigma}(z)$ at the
interface. Moreover, $\rho^{q\sigma}(z)$ does not go to zero in the bulk,
again due to overpolarization in the GT system. The PI-LMF potential corrects
this overpolarization and results in a $\rho^{q\sigma}(z)$ consistent with the
dipole layer produced in the full system. This indicates that the PI-LMF
system reproduces the long range behavior of confined quantum water.

In addition to the PI-LMF solution to the LMF equation, I described three
approximate solutions to the LMF equation: the LMF-RPC approximation, the
centroid approximation, and the Feynman-Kleinert approximation (FK-LMF). The
results obtained using these approximations are compared to the PI-LMF results
in Fig. 2. The charge densities agree for all systems but the FK
approximation, which produces slightly smaller first peaks. The discrepancy in
the FK results stems from the inability of the linear response approximation
to sufficiently shift the first peaks in the atomic densities; direct
simulation of GT water in the presence of $\mathcal{V}_{\rm R}^{\rm FK}(z)$
produces a charge density in agreement with the others. Despite this
difference in the results of the FK approximation, the electrostatic
potentials and smoothed charge densities all agree, indicating that these
approximations for the long range interactions are reasonable for water at
ambient conditions. This good agreement among the various approximations is a
consequence of washing out molecular-scale details in the coarse-graining
inherent within the LMF approach. However, one might expect that these
approximations could break down when the thermal radius of the quantum
particles is comparable to or greater than the smoothing length, $\sigma$, as
is the case for light nuclei at low temperatures Ceperley (1995); Myung,
Hirshberg, and Parrinello (2022), or for light particles like electron or hole
quasiparticles Remsing and Bates (2020); Park, Obliger, and Limmer (2022). In
cases like this, approximating a quantum particle by its centroid could prove
inaccurate, especially if the ring polymers are aspherical.

###  III.2 Quantum Water Works Harder for Dipolar Screening

The impact of NQEs on long range electrostatics can be assessed by examining
the extent to which the quantum and classical GT systems deviate from our
expectations for the full system. The charge densities for the quantum and
classical GT systems, shown in Fig. 3a, are qualitatively similar but display
small differences near the wall. In particular, the magnitudes of the first
two peaks are slightly larger for the quantum GT model. These small
differences between the charge densities lead to large differences when
integrated to obtain the polarization potential, as shown in Fig. 3b. The
potential in the bulk region of the quantum system is significantly larger
than that of the classical one. However, the differences in the electric
fields (derivative of the potential) are localized near the interface, and the
fields are similar in the bulk for both quantum and classical systems. The
overpolarization of interfacial water results from the lack of long range
forces that reorient molecules near the wall. The LMF potential provides these
long range forces and corrects interfacial water structure.

Figure 3: (a) The charge density and (b) polarization potential for the
classical ($P=1$) and quantum ($P=32$) GT systems. The corresponding (c) LMF
potential and (d) its gradient (negative of force) for the same two systems.
Results are shown for the left wall only.

Because quantum GT water overpolarizes more than classical GT water, the
forces needed to reorient quantum GT water should be larger. Indeed, the
corresponding LMF potentials and forces are larger in quantum GT water due to
NQEs, Fig. 3c,d. The LMF, $\mathcal{V}_{\rm R}(z)$, exhibits a larger change
across the interface for the quantum systems at all levels of approximation.
Moreover, the LMF force, $-\nabla\mathcal{V}_{\rm R}(\mathbf{r})$, is larger
in magnitude for the quantum systems. This suggests that the forces required
to achieve proper dielectric screening are larger in quantum systems than
their classical counterparts, by several $k_{\rm B}T$ for water at room
temperature, which may be anticipated from the larger zero point energy of the
quantized nuclei.

###  III.3 Computational Efficiency

Most of the results above were obtained using linear response theory and
simulations of a short range system. However, simulations can be performed
with the presence of the LMF field to also obtain accurate predictions of the
structure of nonuniform fluids. Compared to typical lattice summation
techniques for evaluating long range electrostatics, the LMF-RPC scheme
reduces the cost significantly in two ways. First, RPC reduces the number of
sites at which long range interactions need to be evaluated — here, from 32 to
1. Second, LMF theory replaces $N^{2}$ two-body interactions with $N$ one-body
interactions. This reduction in scaling is beneficial when simulating large
numbers of molecules characteristic of biological and materials systems.

To illustrate the increased efficiency of the LMF-RPC approach in comparison
to Ewald summation-based approaches for evaluating long range electrostatic
interactions in path integral models, I evaluated the time required to perform
a PIMD time step as a function of the number of water molecules in the system,
shown in Fig. 4. The number of water molecules was varied by replicating the
simulation box in the lateral directions ($x$ and $y$). The computational time
was evaluated for the LMF-RPC, particle-particle-particle mesh (PPPM) Ewald,
and GT systems with $P=32$ using a single 2.3 GHz Intel Xeon W core of a 2017
iMac Pro. For the small system size used in the majority of the text (1024
molecules), the simulation time is similar for all approaches. However, as the
size of the system grows, the increased efficiency of the LMF-based approaches
becomes apparent. Moreover, the timings for the GT system are nearly identical
to the LMF systems, indicating that the evaluation of the (one-body) LMF
potential requires minimal overhead and the calculation is dominated by the
evaluation of short range pairwise interactions. Of course, this means that
there is negligible speedup gained by using RPC with LMF theory for the system
sizes studied here, but differences may appear for systems with large numbers
of beads. In contrast to the LMF results, the PPPM timings are slowed by the
evaluation of the long range electrostatic interactions. This suggests that
LMF-based approaches can drastically reduce the computational cost of PIMD
calculations in large-scale molecular systems.

Figure 4: Computation time per MD simulation time step as a function of the
number of water molecules in the system. Lines are linear fits to the data.
Error bars are smaller than the symbol size.

##  IV Conclusions

In this work, I have examined nuclear quantum effects on long range
electrostatic properties of confined water. To do so, I demonstrated that LMF
theory can be used to efficiently and accurately account for long range
electrostatics in path integral simulations of nonuniform liquids. Moreover, a
RPC-based approximations were introduced that leverage the complementary ideas
that underlie the separation of short and long range interactions in both LMF
theory and RPC — long range interactions are slowly varying over relevant
molecular length scales. I expect that the LMF-RPC scheme will be useful for
modeling NQEs in large systems with many light particles at low temperatures
(many beads).

The LMF-RPC scheme can be readily combined with developments in LMF theory to
evaluate NQEs on free energy differences Remsing, Liu, and Weeks (2016);
Remsing and Weeks (2018). The general ideas presented here may also be
valuable for modeling NQEs with the symmetry preserving mean field (SPMF)
theory, which replaces pairwise long range electrostatics with a symmetry-
dependent effective field in each configuration Hu (2014); Yi  _et al._
(2017). The LMF-RPC approach may be particularly powerful when combined with
the short solvent model of Gao  _et al._ Gao, Remsing, and Weeks (2020, 2023)
for molecular assembly. In the short solvent model, long range solvent-solvent
and solute-solvent interactions are truncated everywhere, and the averaged
effects of all long range interactions are accounted for with effective
solute-solute interactions. Therefore, the only possible long range
interactions are between solutes, greatly reducing the number of charged sites
and the associated cost of evaluating long range interactions. Combining the
short solvent model with RPC will result in a model where long range
interactions only need to be evaluated for the solute centroids. All other
interactions are short range. The resulting SSM-RPC scheme could be of great
use for modeling NQEs in self-assembly processes.

The results reported here use empirical force fields to represent intra- and
intermolecular interactions, which neglect coupling of electronic polarization
to nuclear quantum fluctuations. These effects can be taken into account with
ab initio simulations. Previous work has shown that RPC can significantly
speed up ab initio simulations by using a cheap approach to evaluate
interactions on all beads, e.g. density functional tight binding (DFTB)
methods, while higher level ab initio methods are evaluated only on the
centroid and the resulting bead-bead interactions are obtained from the
differences between these two levels of theory Marsalek and Markland (2016,
2017). Such an approach does not readily lend itself to treatment with LMF
theory. However, if LMF theory or similar approaches can be extended to ab
initio models, this would facilitate the use of the LMF-RPC scheme in ab
initio path integral simulations.

An alternative to costly ab initio simulations is to use machine learning
approaches to develop neural network potentials that can produce ab initio
accuracy with much smaller cost Zhang  _et al._ (2018); Behler (2021); Batzner
_et al._ (2022). Traditional neural network potentials lack a good description
of long range interactions Yue  _et al._ (2021), but recent developments
include some description of long electrostatics Ko  _et al._ (2021); Grisafi
and Ceriotti (2019); Yue  _et al._ (2021); Niblett, Galib, and Limmer (2021);
Gao and Remsing (2022); Dhattarwal, Gao, and Remsing (2023). Of particular
interest are neural network potentials that are informed LMF ideas Niblett,
Galib, and Limmer (2021); Gao and Remsing (2022); Dhattarwal, Gao, and Remsing
(2023), like the self-consistent field neural network (SCFNN) Gao and Remsing
(2022); Dhattarwal, Gao, and Remsing (2023). These networks focus on training
short range GT interactions separately from long range interactions. I
anticipate that many of these neural network potentials could be combined with
the LMF-RPC scheme to treat the averaged effects of long range interactions in
path integral and ring polymer MD simulations with ab initio accuracy.

##  V Simulation Details

Path integral molecular dynamics (PIMD) simulations were performed using the
i-PI package Kapil  _et al._ (2019) interfaced with the LAMMPS software
package Plimpton (1995), modified to include the truncated potential
$v_{0}(r)$ and the LMF potential $\mathcal{V}_{\rm R}(z)$. Equations of motion
were integrated in the normal mode representation using the Cayley integration
scheme with a timestep of 0.5 fs Korol, Bou-Rabee, and Miller III (2019). All
simulations were performed in the canonical (NVT) ensemble at a constant
temperature of 298 K maintained using a stochastic velocity rescaling
thermostat Bussi, Donadio, and Parrinello (2007), with 1024 water molecules in
a $27.72\times 27.72\times 150.00$ Å3 simulation cell. All simulations
employed the flexible quantum SPC water model of Voth and coworkers Paesani
_et al._ (2006). Lennard-Jones and short range Coulomb interactions we
truncated at a distance of 9.8 Å. Idealized hydrophobic walls were each
represented with a 9-3 Lennard-Jones potential,

$U_{\rm w}(z)=\varepsilon_{\rm w}\left[\frac{2}{15}\left(\frac{\sigma_{\rm
w}}{z}\right)^{9}-\left(\frac{\sigma_{\rm w}}{z}\right)^{3}\right],$ (38)

with $\sigma_{\rm w}=3.461$ Å and $\varepsilon_{\rm w}=0.43875491$ kcal/mol.
Wall-water interactions are cut off at a distance of 30 Å. Walls are
positioned at $z=0$ and $z=43.06$ Å, in accord with previous work Rodgers and
Weeks (2008b). All ring polymers used $P=32$ beads for each particle in the
system, which has been shown to be sufficiently converged Markland and
Manolopoulos (2008b). The modified LAMMPS source code is available at
github.com/remsing-group/lmf-rpc/.

###### Acknowledgements.

This work is supported by the National Aeronautics and Space Administration
under grant number 80NSSC20K0609, issued through the NASA Exobiology Program.
I acknowledge the Office of Advanced Research Computing (OARC) at Rutgers, The
State University of New Jersey for providing access to the Amarel cluster and
associated research computing resources that have contributed to some of the
results reported here. I thank Atul Thakur for helpful comments on the
manuscript and D. Rodman for inspiration on the figures.

## References

  * Markland and Ceriotti (2018) T. E. Markland and M. Ceriotti, “Nuclear quantum effects enter the mainstream,” Nature Reviews Chemistry 2, 1–14 (2018).
  * Ceriotti  _et al._ (2016) M. Ceriotti, W. Fang, P. G. Kusalik, R. H. McKenzie, A. Michaelides, M. A. Morales, and T. E. Markland, “Nuclear quantum effects in water and aqueous systems: Experiment, theory, and current challenges,” Chemical reviews 116, 7529–7550 (2016).
  * Marsalek and Markland (2017) O. Marsalek and T. E. Markland, “Quantum dynamics and spectroscopy of ab initio liquid water: The interplay of nuclear and electronic quantum effects,” The journal of physical chemistry letters 8, 1545–1551 (2017).
  * Wang, Ceriotti, and Markland (2014) L. Wang, M. Ceriotti, and T. E. Markland, “Quantum fluctuations and isotope effects in ab initio descriptions of water,” The Journal of chemical physics 141, 104502 (2014).
  * Marx, Tuckerman, and Parrinello (2000) D. Marx, M. E. Tuckerman, and M. Parrinello, “Solvated excess protons in water: quantum effects on the hydration structure,” Journal of Physics: Condensed Matter 12, A153 (2000).
  * Tuckerman, Marx, and Parrinello (2002) M. E. Tuckerman, D. Marx, and M. Parrinello, “The nature and transport mechanism of hydrated hydroxide ions in aqueous solution,” Nature 417, 925–929 (2002).
  * Xu  _et al._ (2021) J. Xu, Z. Sun, C. Zhang, M. DelloStritto, D. Lu, M. L. Klein, and X. Wu, “Importance of nuclear quantum effects on the hydration of chloride ion,” Physical Review Materials 5, L012801 (2021).
  * Raugei and Klein (2003) S. Raugei and M. L. Klein, “Nuclear quantum effects and hydrogen bonding in liquids,” J. Am. Chem. Soc. 125, 8992–8993 (2003).
  * Tang  _et al._ (2021) F. Tang, J. Xu, D. Y. Qiu, and X. Wu, “Nuclear quantum effects on the quasiparticle properties of the chloride anion aqueous solution within the gw approximation,” Physical Review B 104, 035117 (2021).
  * Morrone and Car (2008) J. A. Morrone and R. Car, “Nuclear quantum effects in water,” Physical review letters 101, 017801 (2008).
  * Ceriotti  _et al._ (2013) M. Ceriotti, J. Cuny, M. Parrinello, and D. E. Manolopoulos, “Nuclear quantum effects and hydrogen bond fluctuations in water,” Proceedings of the National Academy of Sciences 110, 15591–15596 (2013).
  * Zhuang  _et al._ (2019) D. Zhuang, M. Riera, G. K. Schenter, J. L. Fulton, and F. Paesani, “Many-body effects determine the local hydration structure of cs+ in solution,” The journal of physical chemistry letters 10, 406–412 (2019).
  * Paesani  _et al._ (2010) F. Paesani, S. Yoo, H. J. Bakker, and S. S. Xantheas, “Nuclear quantum effects in the reorientation of water,” The Journal of Physical Chemistry Letters 1, 2316–2321 (2010).
  * Spura  _et al._ (2015) T. Spura, C. John, S. Habershon, and T. D. Kühne, “Nuclear quantum effects in liquid water from path-integral simulations using an ab initio force-matching approach,” Molecular Physics 113, 808–822 (2015).
  * Ojha, Henao, and Kühne (2018) D. Ojha, A. Henao, and T. D. Kühne, “Nuclear quantum effects on the vibrational dynamics of liquid water,” The Journal of chemical physics 148, 102328 (2018).
  * Litman  _et al._ (2018) Y. Litman, D. Donadio, M. Ceriotti, and M. Rossi, “Decisive role of nuclear quantum effects on surface mediated water dissociation at finite temperature,” The Journal of chemical physics 148, 102320 (2018).
  * Lee, McCammon, and Rossky (1984) C.-Y. Lee, J. A. McCammon, and P. Rossky, “The structure of liquid water at an extended hydrophobic surface,” The Journal of chemical physics 80, 4448–4455 (1984).
  * Lee and Rossky (1994) S. H. Lee and P. J. Rossky, “A comparison of the structure and dynamics of liquid water at hydrophobic and hydrophilic surfaces—a molecular dynamics simulation study,” The Journal of chemical physics 100, 3334–3345 (1994).
  * Lum, Chandler, and Weeks (1999) K. Lum, D. Chandler, and J. D. Weeks, “Hydrophobicity at small and large length scales,” J. Phys. Chem. B 103, 4570–4577 (1999).
  * Chandler (2005) D. Chandler, “Interfaces and the driving force of hydrophobic assembly,” Nature 437, 640–7 (2005).
  * Remsing and Weeks (2013) R. C. Remsing and J. D. Weeks, “Dissecting hydrophobic hydration and association,” J. Phys. Chem. B 117, 15479–15491 (2013).
  * Jamadagni, Godawat, and Garde (2011) S. N. Jamadagni, R. Godawat, and S. Garde, “Hydrophobicity of proteins and interfaces: Insights from density fluctuations,” Annual review of chemical and biomolecular engineering 2, 147–171 (2011).
  * Giovambattista, Rossky, and Debenedetti (2012) N. Giovambattista, P. Rossky, and P. Debenedetti, “Computational studies of pressure, temperature, and surface effects on the structure and thermodynamics of confined water,” Annual Review of Physical Chemistry 63, 179–200 (2012).
  * Rego and Patel (2022) N. B. Rego and A. J. Patel, “Understanding hydrophobic effects: Insights from water density fluctuations,” Annual Review of Condensed Matter Physics 13, 303–324 (2022).
  * Parrinello and Rahman (1984) M. Parrinello and A. Rahman, “Study of an f center in molten kcl,” J. Chem. Phys. 80, 860–867 (1984).
  * Tuckerman  _et al._ (1996) M. E. Tuckerman, D. Marx, M. L. Klein, and M. Parrinello, “Efficient and general algorithms for path integral car–parrinello molecular dynamics,” The Journal of Chemical Physics 104, 5579–5588 (1996).
  * Ceperley (1995) D. M. Ceperley, “Path integrals in the theory of condensed helium,” Rev. Mod. Phys. 67, 279 (1995).
  * Habershon  _et al._ (2013) S. Habershon, D. E. Manolopoulos, T. E. Markland, and T. F. Miller, 3rd, “Ring-polymer molecular dynamics: quantum effects in chemical dynamics from classical trajectories in an extended phase space,” Annu. Rev. Phys. Chem. 64, 387–413 (2013).
  * Markland and Manolopoulos (2008a) T. E. Markland and D. E. Manolopoulos, “A refined ring polymer contraction scheme for systems with electrostatic interactions,” Chem. Phys. Lett. 464, 256 – 261 (2008a).
  * Markland and Manolopoulos (2008b) T. E. Markland and D. E. Manolopoulos, “An efficient ring polymer contraction scheme for imaginary time path integral simulations,” J. Chem. Phys. 129, 024105 (2008b).
  * Figueirido, Del Buono, and Levy (1995) F. Figueirido, G. S. Del Buono, and R. M. Levy, “On finite-size effects in computer simulations using the ewald potential,” The Journal of chemical physics 103, 6133–6142 (1995).
  * Hub  _et al._ (2014) J. S. Hub, B. L. De Groot, H. Grubmüller, and G. Groenhof, “Quantifying artifacts in ewald simulations of inhomogeneous systems with a net charge,” Journal of chemical theory and computation 10, 381–390 (2014).
  * Pan, Yi, and Hu (2017) C. Pan, S. Yi, and Z. Hu, “The effect of electrostatic boundaries in molecular simulations: symmetry matters,” Phys. Chem. Chem. Phys. 19, 4861–4876 (2017).
  * Rodgers and Weeks (2008a) J. M. Rodgers and J. D. Weeks, “Local molecular field theory for the treatment of electrostatics,” Journal of Physics: Condensed Matter 20, 494206 (2008a).
  * Remsing, Liu, and Weeks (2016) R. C. Remsing, S. Liu, and J. D. Weeks, “Long-ranged contributions to solvation free energies from theory and short-ranged models,” Proc. Natl. Acad. Sci. USA 113, 2819–2826 (2016).
  * Hu (2014) Z. Hu, “Symmetry-preserving mean field theory for electrostatics at interfaces,” Chemical Communications 50, 14397–14400 (2014).
  * Gao, Remsing, and Weeks (2020) A. Gao, R. C. Remsing, and J. D. Weeks, “Short solvent model for ion correlations and hydrophobic association,” Proc. Natl. Acad. Sci. U.S.A. 117, 1293–1302 (2020).
  * Gao, Remsing, and Weeks (2023) A. Gao, R. C. Remsing, and J. D. Weeks, “Local molecular field theory for coulomb interactions in aqueous solutions,” The Journal of Physical Chemistry B 127, 809–821 (2023).
  * Feynman and Hibbs (2005) R. P. Feynman and A. R. Hibbs,  _Quantum Mechanics and Path Integrals_ (Dover Publications, 2005).
  * Chandler and Wolynes (1981) D. Chandler and P. G. Wolynes, “Exploiting the isomorphism between quantum theory and classical statistical mechanics of polyatomic fluids,” J. Chem. Phys. 74, 4078–4095 (1981).
  * Rodgers, Hu, and Weeks (2011) J. M. Rodgers, Z. Hu, and J. D. Weeks, “On the efficient and accurate short-ranged simulations of uniform polar molecular liquids,” Molecular Physics 109, 1195–1211 (2011).
  * Hu and Weeks (2010) Z. Hu and J. D. Weeks, “Efficient solutions of self-consistent mean field equations for dewetting and electrostatics in nonuniform liquids,” Physical review letters 105, 140602 (2010).
  * Rodgers and Weeks (2008b) J. M. Rodgers and J. D. Weeks, “Interplay of local hydrogen-bonding and long-ranged dipolar forces in simulations of confined water,” Proceedings of the National Academy of Sciences 105, 19136–19141 (2008b), https://www.pnas.org/content/105/49/19136.full.pdf .
  * Baker III, Rodgers, and Weeks (2020) E. B. Baker III, J. M. Rodgers, and J. D. Weeks, “Local molecular field theory for nonequilibrium systems,” The Journal of Physical Chemistry B 124, 5676–5684 (2020).
  * Cox (2020) S. J. Cox, “Dielectric response with short-ranged electrostatics,” Proceedings of the National Academy of Sciences 117, 19746–19752 (2020).
  * Rodgers and Weeks (2009) J. M. Rodgers and J. D. Weeks, “Accurate thermodynamics for short-ranged truncations of coulomb interactions in site-site molecular models,” The Journal of chemical physics 131, 244108 (2009).
  * Remsing, Rodgers, and Weeks (2011) R. C. Remsing, J. M. Rodgers, and J. D. Weeks, “Deconstructing classical water models at interfaces and in bulk,” J. Stat. Phys. 145, 313–334 (2011).
  * Remsing and Weeks (2018) R. C. Remsing and J. D. Weeks, “Alchemical free energy calculations and umbrella sampling with local molecular field theory,” J. Theor. Comput. Chem. 17, 1840003 (2018).
  * Craig and Manolopoulos (2004) I. R. Craig and D. E. Manolopoulos, “Quantum statistics and classical mechanics: Real time correlation functions from ring polymer molecular dynamics,” J. Chem. Phys. 121, 3368–3373 (2004).
  * Shinoda, Miura, and Okazaki (2001) K. Shinoda, S. Miura, and S. Okazaki, “A molecular approach to quantum fluids based on a generalized ornstein–zernike integral equation,” The Journal of Chemical Physics 114, 7497–7505 (2001).
  * Sesé (2002) L. M. Sesé, “Properties of the path-integral quantum hard-sphere fluid in k space,” The Journal of Chemical Physics 116, 8492 (2002).
  * Feynman and Kleinert (1986) R. Feynman and H. Kleinert, “Effective classical partition functions,” Physical Review A 34, 5080 (1986).
  * Remsing and Weeks (2015a) R. C. Remsing and J. D. Weeks, “Hydrophobicity scaling of aqueous interfaces by an electrostatic mapping,” The Journal of Physical Chemistry B 119, 9268–9277 (2015a).
  * Robinson (1973) F. Robinson,  _Macroscopic Electromagnetism_, International series of monographs in natural philosophy (Pergamon Press, 1973).
  * Zangwill (2013) A. Zangwill,  _Modern Electrodynamics_ (Cambridge University Press, 2013).
  * Remsing and Weeks (2015b) R. C. Remsing and J. D. Weeks, “Hydrophobicity scaling of aqueous interfaces by an electrostatic mapping,” J. Phys. Chem. B 119, 9268–9277 (2015b).
  * Myung, Hirshberg, and Parrinello (2022) C. W. Myung, B. Hirshberg, and M. Parrinello, “Prediction of a supersolid phase in high-pressure deuterium,” Physical Review Letters 128, 045301 (2022).
  * Remsing and Bates (2020) R. C. Remsing and J. E. Bates, “Effective mass path integral simulations of quasiparticles in condensed phases,” J. Chem. Phys. 153, 121104 (2020).
  * Park, Obliger, and Limmer (2022) Y. Park, A. Obliger, and D. T. Limmer, “Nonlocal screening dictates the radiative lifetimes of excitations in lead halide perovskites,” Nano Letters 22, 2398–2404 (2022).
  * Yi  _et al._ (2017) S. Yi, C. Pan, L. Hu, and Z. Hu, “On the connections and differences among three mean-field approximations: A stringent test,” Phys. Chem. Chem. Phys. 19, 18514–18518 (2017).
  * Marsalek and Markland (2016) O. Marsalek and T. E. Markland, “Ab initio molecular dynamics with nuclear quantum effects at classical cost: Ring polymer contraction for density functional theory,” The Journal of chemical physics 144, 054112 (2016).
  * Zhang  _et al._ (2018) L. Zhang, J. Han, H. Wang, R. Car, and W. E, “Deep potential molecular dynamics: A scalable model with the accuracy of quantum mechanics,” Phys. Rev. Lett. 120, 143001 (2018).
  * Behler (2021) J. Behler, “Four generations of high-dimensional neural network potentials,” Chemical Reviews 121, 10037–10072 (2021).
  * Batzner  _et al._ (2022) S. Batzner, A. Musaelian, L. Sun, M. Geiger, J. P. Mailoa, M. Kornbluth, N. Molinari, T. E. Smidt, and B. Kozinsky, “E (3)-equivariant graph neural networks for data-efficient and accurate interatomic potentials,” Nature communications 13, 2453 (2022).
  * Yue  _et al._ (2021) S. Yue, M. C. Muniz, M. F. Calegari Andrade, L. Zhang, R. Car, and A. Z. Panagiotopoulos, “When do short-range atomistic machine-learning models fall short?” The Journal of Chemical Physics 154, 034111 (2021).
  * Ko  _et al._ (2021) T. W. Ko, J. A. Finkler, S. Goedecker, and J. Behler, “A fourth-generation high-dimensional neural network potential with accurate electrostatics including non-local charge transfer,” Nature communications 12, 1–11 (2021).
  * Grisafi and Ceriotti (2019) A. Grisafi and M. Ceriotti, “Incorporating long-range physics in atomic-scale machine learning,” The Journal of chemical physics 151, 204105 (2019).
  * Niblett, Galib, and Limmer (2021) S. P. Niblett, M. Galib, and D. T. Limmer, “Learning intermolecular forces at liquid-vapor interfaces,” arXiv preprint arXiv:2107.06208 (2021).
  * Gao and Remsing (2022) A. Gao and R. C. Remsing, “Self-consistent determination of long-range electrostatics in neural network potentials,” Nature communications 13, 1572 (2022).
  * Dhattarwal, Gao, and Remsing (2023) H. S. Dhattarwal, A. Gao, and R. C. Remsing, “Dielectric saturation in water from a long-range machine learning model,” The Journal of Physical Chemistry B 127, 3663–3671 (2023).
  * Kapil  _et al._ (2019) V. Kapil, M. Rossi, O. Marsalek, R. Petraglia, Y. Litman, T. Spura, B. Cheng, A. Cuzzocrea, R. H. Meißner, D. M. Wilkins, _et al._ , “i-pi 2.0: A universal force engine for advanced molecular simulations,” Computer Physics Communications 236, 214–223 (2019).
  * Plimpton (1995) S. Plimpton, “Fast parallel algorithms for short-range molecular dynamics,” J. Comp. Phys. 117, 1–19 (1995).
  * Korol, Bou-Rabee, and Miller III (2019) R. Korol, N. Bou-Rabee, and T. F. Miller III, “Cayley modification for strongly stable path-integral and ring-polymer molecular dynamics,” The Journal of chemical physics 151, 124103 (2019).
  * Bussi, Donadio, and Parrinello (2007) G. Bussi, D. Donadio, and M. Parrinello, “Canonical sampling through velocity rescaling,” The Journal of chemical physics 126, 014101 (2007).
  * Paesani  _et al._ (2006) F. Paesani, W. Zhang, D. A. Case, T. E. Cheatham III, and G. A. Voth, “An accurate and simple quantum model for liquid water,” The Journal of chemical physics 125, 184507 (2006).