#  Donnan equilibrium in charged slit-pores from a hybrid nonequilibrium
Molecular Dynamics / Monte Carlo method with ions and solvent exchange

Jeongmin Kim  Department of Energy Engineering, Korea Institute of Energy
Technology (KENTECH), Naju 58330, Republic of Korea  Benjamin Rotenberg
<EMAIL_ADDRESS>Sorbonne Université, CNRS, Physico-
chimie des Électrolytes et Nanosystemès Interfaciaux, PHENIX, F-75005 Paris,
France  Réseau sur le Stockage Electrochimique de l’Energie (RS2E), FR CNRS
3459, 80039 Amiens Cedex, France

###### Abstract

Ion partitioning between different compartments (_e.g._ a porous material and
a bulk solution reservoir), known as Donnan equilibrium, plays a fundamental
role in various contexts such as energy, environment, or water treatment. The
linearized Poisson-Boltzmann (PB) equation, capturing the thermal motion of
the ions with mean-field electrostatic interactions, is practically useful to
understand and predict ion partitioning, despite its limited applicability to
conditions of low salt concentrations and surface charge densities. Here, we
investigate the Donnan equilibrium of coarse-grained dilute electrolytes
confined in charged slit-pores in equilibrium with a reservoir of ions and
solvent. We introduce and use an extension to confined systems of a recently
developed hybrid nonequilibrium molecular dynamics / grand canonical Monte
Carlo simulation method ("H4D"), which enhances the efficiency of solvent and
ion-pair exchange via a fourth spatial dimension. We show that the validity
range of linearized PB theory to predict the Donnan equilibrium of dilute
electrolytes can be extended to highly charged pores, by simply considering
renormalized surface charge densities. We compare with simulations of implicit
solvent models of electrolytes and show that in the low salt concentrations
and thin electric double layer limit considered here, an explicit solvent has
a limited effect on the Donnan equilibrium and that the main limitations of
the analytical predictions are not due to the breakdown of the mean-field
description, but rather to the charge renormalization approximation, because
it only focuses on the behavior far from the surfaces.

††preprint:

##  I Introduction

The composition of an electrolyte near a charged surface or in a charged
porous material differs from that of a bulk solution reservoir with which it
is in equilibrium. This ion partitioning between different compartments, known
as Donnan equilibrium [1], plays an important role in biology [2] or in the
environment [3, 4, 5, 6]; understanding and predicting it is necessary for
various applications _e.g._ water treatment [7], membrane technology [8, 9,
10, 11], energy [12], or electrocatalysis [13, 14]. The classical approach to
this Donnan equilibrium is to consider the equality of the electrochemical
potentials of the various species in the reservoir and in the "phase" or
region of interest. This results in the celebrated quadratic relation between
the salt concentrations in the reservoir and the confined medium, and a so-
called Donnan potential difference between the two regions [1]. This Donnan
model can be refined to include various interactions between the ions or with
the medium [15, 16, 17, 11, 10].

This simple approach neglects the inhomogeneities in the distribution of ions
and the corresponding electrostatic potential inside the pores. In fact, the
balance between electrostatic interactions (energy) and the thermal motion of
ions (entropy) results near charged walls in electric double-layers (EDL),
with an enrichment of counterions and a depletion of co-ions with respect to
the bulk [18, 19, 20, 21]. Poisson-Boltzmann (PB) theory can capture the main
features of such EDL by considering point charges in an implicit solvent, and
treating electrostatic interactions at the mean-field level. For sufficiently
small salt concentrations and surface charge densities, PB theory can further
be linearized, resulting in Debye-Hückel (DH) theory, which predicts an
exponential decay of density and potential profiles as a function of the
distance from the wall, over the so-called Debye length $\lambda_{D}$. The
classical Donnan approach typically applies in the limit where this length is
large compared to the pore size, so that spatial variations can be neglected.
PB or DH theory can be used in simple geometries to predict ion partioning as
a function of pore size, surface charge density and salt concentration in the
reservoir. However its applicability is usually limited to conditions of low
salt concentrations and surface charge densities. For larger surface charge
densities, the large concentration of counterions (overestimated by PB theory,
which neglects in particular the excluded volume interactions due to the
finite size of the ions) is usually lumped into a (Stern) layer of condensed
ions, resulting in the so-called Gouy-Chapman-Stern theory of the interface
[18, 19].

Despite its great merits, this simplified picture of the EDL fails to account
for microscopic features which may play an important role at the interface.
For example, the molecular nature of the solvent results in packing of
hydration effects known to influence the distribution of ions at the surface
[22, 23, 24, 25, 26, 27, 28, 29, 10], or the force between atomically flat
mica surfaces in Surface Force Balance experiments (SFB) [30, 31, 32]. Efforts
have been made to modify PB theory to account for effects such as the finite
size of ions [33, 34, 35, 36], ion-surface interactions [37, 38, 39, 11],
electrostatic correlations [40], solvent polarization [41, 42], or surface
charge regulation [43, 44, 45, 46, 47, 48, 49, 50].

In order to go beyond these theories, Grand-canonical (GC) Monte Carlo (MC)
simulations have been used to investigate the Donnan equilibrium of confined
electrolytes in equilibrium with a reservoir [51, 52, 53, 54, 55]. However,
while GC simulations of electrolytes described as ions in an implicit solvent
can be performed with sufficiently large acceptance rates of trial
insertion/deletion MC moves to allow for a proper sampling of configurations,
such simulations with an explicit solvent are much more challenging, not only
due to the larger computational cost of interactions, but also to the very low
acceptance rate of such moves. Several advanced MC methods to enhance exchange
efficiency have been developed, including a cavity-bias method [56], a
configuration-biased method [57, 58, 52], a Boltzmann-bias method [59], a
continuous fraction component method [60, 61], identity exchange methods [62,
63, 64], the adaptive resolution method [65, 66], and hybrid methods [67, 68,
69, 53, 70, 71, 72, 73, 74, 75, 76]. Recently, Belloni introduced a hybrid
MC/nonequilibrium molecular dynamics (MD) method [75], called "H4D", with
trial moves during which the system adjusts to the gradual addition/deletion
of particles via a non-physical additional dimension. This results in a
significant increase in the acceptance rate of such moves, as demonstrated
i.e. for both ion-pairs and solvent exchange for bulk aqueous NaCl
electrolytes [75, 76].

Here, we extend the H4D method to the case of confined systems and use it to
investigate the Donnan equilibrium of coarse-grained dilute electrolytes in
charged slit-pores in equilibrium with a reservoir of ions and solvent. We
show that the validity range of linearized PB theory to predict the Donnan
equilibrium of dilute electrolytes can be extended to highly charged pores, by
considering renormalized surface charge densities [43, 44, 45, 46, 49, 50],
which can be computed analytically by comparing the full and linearized PB
equations, instead of the bare ones. By comparing with simulations of implicit
solvent models of electrolytes, we find that for the small salt concentration
considered here, an explicit solvent introduces oscillations in the ionic
density profiles, but has a limited effect on the excess salt concentration
inside the pore. In the low concentration and thin electric double-layer limit
considered here, the main limitations of the analytical predictions are not
due to the breakdown of the mean-field description, but rather to the charge
renormalization approximation, which is not sufficient to capture the overall
excess ion density because it focuses on the behavior far from the walls.

The rest of the manuscript is organized as follows. Section II discusses the
classical description of the Donnan equilibrium as well as its prediction from
Poisson-Boltzmann theory and the surface charge renormalization approach.
Section III summarizes the basics of the H4D method for grand-canonical
molecular dynamics (GCMD) simulations introduced in Refs. [75, 76] and extends
it to confined electrolytes. Section III also presents the model system for
which we studied the Donnan equilibrium: Lennard-Jones model electrolytes
confined in a charged slit-like pore. Section IV then discusses the results of
GCMD simulations and their comparison with the linearized PB theory, including
the surface charge renormalization approach. The main conclusions are
summarized in Section V

##  II Donnan equilibrium of electrolytes between a pore and a reservoir

As mentioned above, the composition of an electrolyte near a charged surface
or in a charged porous material differs from that of a bulk solution reservoir
with which it is in equilibrium. Here, we first recall the classical Donnan
description of this equilibrium (Section II.1) before turning to the
predictions obtained for a charged slit pore using the Poisson-Boltzmann
equation and its linearized version (Section II.2). We finally introduce an
effective surface charge density of the pore walls by comparing the full and
linearized Poisson-Boltzmann equations, according to the charge
renormalization approach: a fraction of the counterions is included as part of
the surface charge, thereby reducing the latter sufficiently for the
linearization to be a good approximation far from the surface [77, 44]
(Section II.3).

###  II.1 Classical Donnan approach

At equilibrium, the electrochemical potential of the various species is the
same inside the pore of interest and in the reservoir with which the solution
exchanges particles. For ionic species $i\in\\{+,-\\}$ with charges
$q_{\pm}=\pm e$ with $e$ the elementary charge, the equality of
electrochemical potential between the bulk reservoir $B$, with concentrations
$\rho_{i}^{B}$ and potential $\Psi^{B}$, and the confined region $C$, with
position-dependent concentrations and potential $\rho^{C}_{i}({\bf r})$ and
$\Psi^{C}({\bf r})$, reads:

$\begin{split}\ln(\gamma^{B}_{i}\rho^{B}_{i}/\rho_{\text{ref}})+\beta
q_{i}\Psi^{B}=\\\ \ln[\gamma^{C}_{i}({\bf r})\rho^{C}_{i}({\bf
r})/\rho_{\text{ref}}]+\beta q_{i}\Psi^{C}({\bf r}),\end{split}$ (1)

where $\gamma_{i}$ is the activity coefficient of species $i$,
$\rho_{\text{ref}}$ is a reference concentration and $\beta=1/k_{B}T$ with
$k_{B}$ the Boltzmann constant and $T$ the temperature. Then, the ion
concentrations in the confined region is:

$\begin{split}\rho^{C}_{i}({\bf
r})=\rho^{B}_{i}\frac{\gamma^{B}_{i}}{\gamma^{C}_{i}({\bf r})}\exp[-\beta
q_{i}(\Psi^{C}({\bf r})-\Psi^{B})]\;.\end{split}$ (2)

This can be further simplified by assuming an ideal behavior of the
electrolyte (activities $\gamma_{i}=1$) and no spatial dependence of the
potential and of the ion concentrations inside the pore. Assuming without loss
of generality that the surface is negatively charged, with a corresponding
excess of cations inside the pore with concentration $\rho^{C}_{+,ex}$
counterbalancing the surface charge, one obtains the well-known quadratic
equation for the salt concentration $\rho^{C}_{\text{salt}}$ in the pore:

$\begin{split}(\rho_{\text{salt}})^{2}=\rho^{C}_{+}\rho^{C}_{-}=(\rho^{C}_{\text{salt}}+\rho^{C}_{+,\text{ex}})\rho^{C}_{\text{salt}}\;,\end{split}$
(3)

where $\rho_{\text{salt}}=\rho^{B}_{\pm}=\rho^{B}_{\text{salt}}$, the salt
concentration in the electrically neutral reservoir. Inside a charged pore,
the co-ion (resp. counterion) density is expected to decrease (resp. increase)
with respect to the bulk salt density, as a result of the balance between
electrostatic and osmotic pressures [78]. As in Equation 2, the electric
potential determines the ion distributions in a confined region.

###  II.2 Linearized Poisson-Boltzmann (Debye-Hückel) theory

For a 1:1 electrolyte in a symmetric charged slit-pore, with two walls with
surface charge density $\Sigma_{s}$ located at $z=\pm H/2$, one obtains the
one-dimensional (1D) electrostatic potential profile $\Psi(z)$ by solving the
following 1D Poisson-Boltzmann equation (taking the electrostatic potential in
the reservoir as reference $\Psi_{B}=0$) [78]:

$\frac{d^{2}}{dz^{2}}\Psi(z)=-\frac{\rho_{q}(z)}{\varepsilon_{0}\varepsilon_{s}}=\frac{\kappa^{2}_{D}}{\beta
e}\sinh[\beta e\Psi(z)]$ (4)

where $\rho_{q}(z)=e\left[\rho^{C}_{+}(z)-\rho^{C}_{-}(z)\right]$ is the local
charge density, $\varepsilon_{s}$ is the solvent dielectric constant, and
$\kappa_{D}=\lambda^{-1}_{D}$, with

$\lambda_{D}=\sqrt{\frac{\varepsilon_{0}\varepsilon_{s}}{2\beta\rho_{\text{salt}}e^{2}}}$
(5)

the Debye screening length, with $\rho_{\text{salt}}$ the salt concentration
in the electroneutral reservoir which also sets the potential reference. For
sufficiently low surface charge densities, potential variations across the
slit-pores are small ($|e\beta\Psi(z)|\ll 1$) and Eq. 4 can be linearized,
resulting in:

$\Psi(z)=\Psi_{0}\frac{\cosh(\kappa_{D}z)}{\sinh(\kappa_{D}H/2)},$ (6)

where $\Psi_{0}=\Sigma_{s}/\varepsilon_{s}\kappa_{D}$ is the potential at the
walls. Accordingly, the mean ion concentrations ($\bar{\rho}_{+}$ and
$\bar{\rho}_{-}$) inside the negatively charged pore can be computed from the
linearized Boltzmann distributions
$\rho^{C}_{\pm}(z)=\rho_{\text{salt}}\left[1\mp\beta e\Psi(z)\right]$:

$\bar{\rho}_{+}(H)=\frac{1}{H}\int_{-H/2}^{{H}/2}dz~{}\rho^{C}_{+}(z)=\rho_{\text{salt}}+\frac{|\Sigma_{s}|}{eH},$
(7)

and

$\bar{\rho}_{-}(H)=\frac{1}{H}\int_{-H/2}^{{H}/2}dz~{}\rho^{C}_{-}(z)=\rho_{\text{salt}}-\frac{|\Sigma_{s}|}{eH}.$
(8)

The mean ion densities deviate from their bulk value, depending on the pore
size $H$, with an enrichment in counterions and depletion of co-ions, as
expected. Furthermore, these deviations are equal and opposite, so that the
total ion density ($2\rho_{\text{salt}}$) is independent of the pore size and
the excess ion density vanishes
($\Delta\rho_{\text{ex}}(H)=\bar{\rho}_{+}(H)+\bar{\rho}_{-}(H)-2\rho_{\text{salt}}=0$).

###  II.3 Effective surface charge density

When the distance between the confining walls is much larger than the
thickness of the EDL ($\lambda_{D}\ll H$), the so-called thin-EDL or single-
wall approximation works well, resulting in the following analytical solution
of the non-linear Eq. 4 near each wall:

$\Psi(z_{d})=\frac{4}{\beta e}\tanh^{-1}[\gamma\exp(-\kappa_{D}z_{d})]$ (9)

with $z_{d}$ the distance from the wall and the dimensionless parameter
$\gamma\in\left.\right[0,1\left]\right.$, determined by the boundary condition
(fixed surface charge density $\Sigma_{s}$), is

$\gamma=-l_{GC}\kappa_{D}+\sqrt{(l_{GC}\kappa_{D})^{2}+1}$ (10)

with the Gouy-Chapman length

$l_{GC}=\frac{2\varepsilon_{0}\varepsilon_{s}}{\beta e|\Sigma_{s}|}.$ (11)

At large distances from the charged surface, the electric potential in Eq. 9
decays as $\frac{4\gamma}{\beta e}\exp(-\kappa_{D}z_{d})$. Comparison with the
solution of the linearized solution Eq. 6 then allows to define an effective
surface charge density such that the two solutions coincide at large distances
[77, 79, 80], as $\Sigma_{s,\text{eff}}=a_{\text{eff}}\Sigma_{s}$, with:

$a_{\text{eff}}=2\gamma l_{GC}\kappa_{D}.$ (12)

This effective surface charge accounts for the screening of the electric field
by the ions close to the surface, up to the point where it has sufficiently
decayed so that the linearization becomes accurate. This is closely related to
the concept of charge renormalization [43, 44, 45, 46, 49, 50], widely adopted
in the colloid community. As shown in Fig. 1, $a_{\text{eff}}$ increases from
0 for highly charged surfaces ($l_{GC}\lesssim\lambda_{D}$), reflecting the
fact that the strong attraction of counterions towards the wall screens most
of the electric field at the surface, to 1 for lightly charged surfaces
($l_{GC}\gtrsim\lambda_{D}$). We note that the expression of $a_{\text{eff}}$
in Eq. 12 is still at the mean field level; the contribution of ionic
correlations can further decrease $a_{\text{eff}}$, particularly at high
$\Sigma_{s}$. [79, 80]

Figure 1: Charge renormalization parameter $a_{\text{eff}}$ (see Eq. 12)
characterizing the effective surface charge density, as a function of the
ratio between Gouy-Chapman length $l_{GC}$ and Debye screening length
$\lambda_{D}$. Four red markers indicates the conditions at which simulations
in this study were conducted, with a fixed $\lambda_{D}$ corresponding to
$\rho_{\text{salt}}=0.0144$ in reduced LJ units (see Section III), while
$l_{GC}$ varies according to the surface charge density.

Replacing the non-linear potential profile by its linearized counterpart for
an effective surface charge density $a_{\text{eff}}\Sigma_{s}$ leads to the
following estimates of the average concentrations in the slit-pore:

$\bar{\rho}_{+}(H)=\frac{1}{H}\int_{-H/2}^{H/2}dz~{}\rho^{C}_{+}(z)=\rho_{\text{salt}}+(2-a_{\text{eff}})\frac{|\Sigma_{S}|}{eH},$
(13)

and

$\bar{\rho}_{-}(H)=\frac{1}{H}\int_{-H/2}^{H/2}dz~{}\rho^{C}_{-}(z)=\rho_{\text{salt}}-a_{\text{eff}}\frac{|\Sigma_{S}|}{eH},$
(14)

which reduce to Eqs. 7\- 8 in the limit $l_{GC}/\lambda_{D}\to\infty$, where
$a_{\text{eff}}=1$. While the system remains electroneutral (as it should),
with an average charge concentration in the liquid
$e(\bar{\rho}_{+}-\bar{\rho}_{-})=2|\Sigma_{s}|/H$ balancing the bare charge
of both surfaces, the salt concentration in the pore
($\rho^{C}_{\text{salt}}=\rho^{C}_{-}$) differs from that in the reservoir
($\rho_{\text{salt}}$). This Donnan equilibrium is quantified by the excess
average ion density, which at this level of description reads

$\Delta\rho_{\text{ex}}=\bar{\rho}_{+}(H)+\bar{\rho}_{-}(H)-2\rho_{\text{salt}}=2(1-a_{\text{eff}})\frac{|\Sigma_{s}|}{eH}.$
(15)

We recall that the above results involving the renormalized surface charge
density assume sufficiently large pores, i.e. $H\gg l_{GC},\lambda_{D}$.

##  III Model systems and simulation methods

Here, we introduce the Lennard-Jones (LJ) electrolytes confined in a charged
slit-pore studied in this work (Section III.1), and present the hybrid
nonequilibrium MD / Monte Carlo algorithm with particle exchange via a fourth
spatial dimension, which we extend to the case of confined electrolytes
(Section III.2).

###  III.1 Model Lennard-Jones electrolytes in bulk or confined in a slit-like
pore

LJ electrolytes. The model electrolytes consist of neutral and charged LJ
particles, representing solvent molecules and ions, respectively, all of which
of the same size and same mass $m$ [39]. The electrostatic interactions
between ions are screened by the solvent relative permittivity
$\varepsilon_{s}$. All the LJ interactions $U_{LJ}$ were truncated and shifted
at a cut-off distance $r_{c}$.

$\begin{split}U_{LJ}(r)=&4\epsilon\bigg{[}\bigg{(}\frac{\sigma}{r}\bigg{)}^{12}-\bigg{(}\frac{\sigma}{r}\bigg{)}^{6}-\bigg{(}\frac{\sigma}{r_{c}}\bigg{)}^{12}+\bigg{(}\frac{\sigma}{r_{c}}\bigg{)}^{6}\bigg{]},\end{split}$
(16)

assuming the same LJ energy $\epsilon$, and diameter $\sigma$ for interactions
between all types of particles. The cut-off distance for LJ interactions
between ions and solvent particles and between solvent particles is
$r^{*}_{c}=r_{c}/\sigma=2.5$. Here and in the following, the asterisk
represents a quantity in reduced LJ unit. In order to explore the effect of
packing due to the presence of explicit solvent particles, we also consider
systems where the solvent is only represented by its permittivity, so that
only ions are explicitly simulated. In that case, we either use the same LJ
interactions between ions with $r^{*}_{c}=2.5$ (implicit LJ model) or purely
repulsive Weeks-Chandler-Andersen interactions with $r^{*}_{c}=2^{1/6}$
(implicit WCA model).

The Coulomb interaction ($U_{C}$) between LJ (or WCA) ions is:

$U_{C}(r)=\frac{1}{4\pi\varepsilon_{0}\varepsilon_{s}}\frac{q_{i}q_{j}}{r}=\frac{\epsilon}{4\pi\varepsilon_{0}\varepsilon_{s}}\frac{q^{*}_{i}q^{*}_{j}}{r^{*}},$
(17)

where $\varepsilon_{0}$ is the vacuum permittivity and $\varepsilon_{s}$ the
dielectric constant of the solvent, fixed to unity. LJ ions carry either
$q^{*}_{i}=q_{i}/\sqrt{\sigma\epsilon}=+1$ or -1, while solvent particles
carry no charge. Coulomb interactions are calculated using the particle-
particle and particle-mesh (PPPM) method, with a real-space cut-off at
$r^{*}=3.5$.

Slit-pore. The LJ electrolytes are confined in a slit-pore between two
atomically flat walls modeled as in Ref. 45. Each wall consists of 5 layers of
a cubic lattice with lattice spacing $\sigma$, and the surface is a two-
dimensional square lattice with $n_{x}\times n_{y}=20\times 20$ unit cells, so
that the dimensions of the box in the directions parallel to the walls are
$L^{*}_{x}\times L^{*}_{y}=20\times 20$. The distance $H^{*}=H/\sigma$ between
two inner-most layers (closest to the electrolyte) varies between 12 and 36,
ensuring $H\gg\lambda_{D}$ in all cases studied in this work. There are
$n_{x}\times n_{y}=N_{w,\text{in}}$ atoms in each layer (hence
$10N_{w,\text{in}}$ atoms for the two walls), whose positions are fixed in
space. The total system size, including the walls on both sides, in the $z$
direction is thus $L^{*}_{z}=22,~{}28,~{}34$, and 46 for
$H^{*}=12,~{}18,~{}24$, and 36, respectively. Periodic boundary conditions are
used in the $x$ and $y$ directions only.

Wall atoms interact with the electrolyte via WCA and Coulomb interactions. No
LJ interactions between wall atoms were included as they are fixed in space.
Only the atoms in the inner-most layer of each wall carry partial charges
$q^{*}_{w}$ that determine the surface charge density
$\Sigma^{*}_{s}=N_{w,\text{in}}q^{*}_{w}/L^{*}_{x}L^{*}_{y}$, while all other
wall atoms are neutral. We consider $\Sigma^{*}_{s}$ varying from -0.05 to
-0.5, and accordingly the Gouy-Chapman length $l^{*}_{GC}$
($=\frac{\varepsilon_{s}T^{*}}{2\pi|\Sigma^{*}_{s}|}$) varies from 3.2 to
0.32. With non-zero $q^{*}_{w}$, the wall atoms only in the inner-most layers
interact with the ions via Coulomb interactions. The Coulomb interactions,
including the charged wall atoms and the ionic species, are calculated using
the PPPM method, with the same real-space cut-off as for bulk systems,
$r^{*}=3.5$.

###  III.2 Non-equilibrium MD/MC with H4D

This section discusses a brief overview of the H4D method, and introduces its
extension to confined systems with a bias to alleviate the steric overlaps
with confining walls. Monte Carlo simulations in the grand-canonical ensemble
are challenging for dense systems, in particular when both steric and
electrostatic interactions are present, due to the low acceptance probability
of insertion/deletion moves. Hybrid GCMD simulations, whereby during an
equilibrium MD simulation Monte Carlo moves are proposed using nonequilibrium
MD (NEMD) instead of a "brutal" insertion/deletion, allow to improve the
acceptance probability, at the price of increasing the computational cost of
each trial move. An elegant type of such trial moves, which proved efficient
even for aqueous electrolytes with explicit solvent and ion exchanges, was
proposed by Luc Belloni in Ref. 75 and refined and implemented in the LAMMPS
simulation package in Ref. 76.

Trial moves The essence of the algorithm, H4D, is to insert/delete particles
utilizing a fourth dimension, called the “vertical” axis, orthogonal to all
other physical dimensions, in order to enhance the efficiency of particle
exchange by allowing the system to relax as the considered particles are
introduced/removed. All the interactions, including LJ and Coulomb potentials,
are treated in (3+1)D, instead of 3D. In practice, one introduces a time-
dependent altitude, $w(t)$, of a so-called “flying” particle to exchange,
which is changed with a constant-velocity $v_{f}$ in this work:

$w(t)=v_{f}\cdot(t-t_{i})+w_{max},$ (18)

where $t\in[t_{i},t_{f}]$ in steps of $\delta t$. The sign of $v_{f}$
determines whether a trial move is for addition or removal with proper
boundary conditions of $w(t)$: For a trial insertion, $w(t_{i})=w_{max}$ and
$w(t_{f})=0$, while for a trial deletion, $w(t_{i})=0$ and $w(t_{f})=w_{max}$.
The key parameters of H4D include the maximum altitude $w_{max}$ and vertical
velocity $v_{f}(t)$, both of which are crucial to determine the efficiency of
particle exchange. The altitude schedule (Eq. 18) satisfies the detailed
balance condition, being symmetric and pre-determined for the trial insertion
and deletion moves. In the case of ion-pair exchange, we choose the same
boundary conditions of $w(t)$ and the same $v_{f}$ for both flying ions. In
practice, $v_{f}$ is determined by the number of NEMD steps $N_{f}$:
$v_{f}=w_{max}\cdot(N_{f}\delta t)^{-1}$. For an instantaneous (infinitely
fast) exchange as in a conventional GC MC, $N_{f}=0$.

The trial insertion/deletion moves are obtained by propagating the trajectory
in the NVE ensemble, using a deterministic, time-reversible, and volume-
preserving integrator [68]. In this case, a proposed move is solely determined
by the preparation of all momenta; the momenta of flying particles are
prepared randomly according to the Maxwell-Boltzmann distribution, and the
momenta of non-flying particles are taken from the equilibrium MD. No momentum
in the "vertical" direction needs to be assigned, along which an external
force determines the altitude of flying particles according to their altitude
schedule, and leaves non-flying particles in 3D at zero altitude. Furthermore,
we use the symmetric two-end momentum reversal scheme [70], which satisfies
the detailed balance, since we evolve the systems in the equilibrium phase
using MD as well. In this scheme, for every trial addition or deletion, the
momenta of all particles are reversed with a probability of one-half both at
the beginning and at the end of NEMD.

Generation and acceptance probabilities A proposed trial move via the H4D
method is accepted according to a Metropolis criterion, with a probability
ensuring detailed balance [68, 81, 75]. Here, we consider electrolytes
confined in a slit-like a pore, with impenetrable walls. As commonly done
under such circumstances, we introduce a bias in addition to the interactions
between the liquid and the wall to favor trial moves consistent with the
exclusion of the former from the latter. With this extension of previous work
for bulk electrolytes [75, 76], the Metropolis acceptance probability
$f_{ins}$ for trial ion-pair insertions ($N_{salt}\rightarrow N_{salt}+1$) is:

$\begin{split}f_{ins}(\vec{r}_{a},\vec{r}_{c})&=\min\bigg{[}1,\exp(-\beta\Delta
H_{N_{salt}\rightarrow
N_{salt}+1})\exp(\beta\mu_{salt})\bigg{(}\frac{L_{x}L_{y}H}{\Lambda_{s}^{3}}\frac{1}{N_{salt}+1}\bigg{)}^{2}\frac{B^{del}(\vec{r}_{a}|\vec{r}_{c})}{B^{ins}(\vec{r}_{a}|\vec{r}_{c})}\frac{C^{del}(\vec{r}_{c})}{C^{ins}(\vec{r}_{c})}\bigg{]},\\\
&\equiv\min\bigg{[}1,\exp(-\beta\Delta M+\beta\mu_{salt})\bigg{]},\end{split}$
(19)

where $\Delta M$ is the nonequilibrium work, $\Delta H=\Delta U+\Delta K$ the
total mechanical energy difference (with $U$ and $K$ the potential and kinetic
energies, respectively), $\mu_{salt}$ the chemical potential of the ion pair,
and $\Lambda_{s}(=\sqrt{\Lambda_{+}\Lambda_{-}})$ the geometrical mean of
thermal de Broglie wavelengths of an ion pair. Here, $f_{ins}$ depends
explicitly on the 3D positions of flying ions ($\vec{r}_{a}$ for the anion and
$\vec{r}_{c}$ for the cation) via two biases $B(\vec{r}_{a}|\vec{r}_{c})$ and
$C(\vec{r}_{c})$ to select the cation and anion: the former favors shorter
distances between flying ions (as in the bulk [75, 76]), while the latter
introduces a bias in the position of the flying cation due to the imposed
confinement.

In our previous work [76], $B(\vec{r}_{a}|\vec{r}_{c})$ was found to be
crucial to enhance the efficiency of H4D for bulk electrolytes, in addition to
modulating the interactions via the fourth dimension. On the one hand, for a
trial insertion, the probability to insert the anion at
$\vec{r}_{a}=(x_{a},y_{a},z_{a})$ for a given cation position
$\vec{r}_{c}=(x_{c},y_{c},z_{c})$ is:

$B^{ins}(\vec{r}_{a},\vec{r}_{c})=b^{ins}(x_{a}|x_{c})b^{ins}(y_{a}|y_{c})b^{ins}(z_{a}|z_{c})\cdot
V,$ (20)

where $V=L_{x}L_{y}H$ is the volume of the confined system (with $H$ is the
distance between two innermost layers of the confining walls, while $L_{x}$
and $L_{y}$ are the simulation box dimensions), and $b^{ins}$ applies
separately to the three directions of space. Without such a bias,
$b^{ins}=1/L_{i}$ for $i\in\\{x,y\\}$ and $b^{ins}=1/H$ for the $z$ direction
so that $B^{ins}=1$ (see Section III.3 for the functional form and the
corresponding parameters used in this work). On the other hand, for a trial
deletion, one of the $N_{salt}+1$ anions, knowing the position $\vec{r}_{c}$
of the flying cation, is chosen with a probability:

$B^{del}(\vec{r}_{a},\vec{r}_{c})=\frac{B^{ins}(\vec{r}_{a},\vec{r}_{c})}{\sum_{n=1}^{N_{salt}+1}B^{ins}(\vec{r}_{a,n},\vec{r}_{c})}\cdot(N_{salt}+1).$
(21)

In the present work, we use the same functional form for the bias in both
insertion and deletion trial moves. In the absence of such a bias,
$B^{del}=1$.

For confined electrolytes, we introduce another bias, not considered in bulk
electrolytes, to in favor of 3D positions of flying ions away from the
surfaces. Indeed, a large overlap between a confining wall and a flying ion
with its randomly assigned momenta could lead to numerical instabilities or
require a quite small $\delta t$ for NEMD, significantly lowering the
efficiency of H4D. Even though in principle one could bias the positions of
both ions with respect to the walls separately, here we introduce the bias in
a sequential manner: first the cation is introduced with the bias due to the
confinement on $\vec{r}_{c}$ with $C(\vec{r}_{c})$, then the anion is
introduced with the bias for the interionic distance using
$B(\vec{r}_{a},\vec{r}_{c})$. We found this sufficient to ensure that both
ions are sufficiently far from the walls.

Here, we discuss the case of confined electrolytes in a slit-pore whose
boundaries are at $z=\pm\frac{H}{2}$. The cation is inserted at a position
$\vec{r}_{c}$ with probability:

$C^{ins}(\vec{r}_{c})=c^{ins}_{x}(r_{x})c^{ins}_{y}(r_{y})c^{ins}_{z}(r_{z})\cdot
V,$ (22)

where we in fact consider no bias for the periodic directions, i.e.
$c^{ins}_{i}(r_{i})=L_{i}^{-1}$ for $i\in\\{x,y\\}$, and a one-dimensional
truncated-shifted Gaussian distribution function along the confinement:

$c^{ins}_{z}(r_{z})=\frac{1}{\mathcal{N}_{c}}\bigg{[}\exp\bigg{(}-\alpha_{c}r_{z}^{2}\bigg{)}-\exp\bigg{(}\alpha_{c}\frac{H}{4}^{2}\bigg{)}\bigg{]},$
(23)

where $\mathcal{N}_{c}$ is a normalization constant that ensures
$\int_{-H/2}^{H/2}c^{ins}_{z}(r_{z})dr_{z}=1/H$, and $\alpha_{c}$ determines
the sharpness of the distribution whose mean is at the center of the pore.
This choice of bias ensures that the insertion probability vanishes at the
confining wall boundaries, i.e. $c^{ins}_{z}(\pm\frac{H}{2})=0$. Without such
a bias, $C^{ins}(\vec{r}_{c})=1$ for randomly assigned 3D positions as in bulk
electrolytes.

For a trial ion-pair deletion, the flying cation is chosen among the
$N_{salt}+1$ cations with a probability:

$\begin{split}C^{del}(\vec{r}_{c})&=\frac{C^{ins}(\vec{r}_{c})}{\sum_{n=1}^{N_{salt}+1}C^{ins}(\vec{r}_{c,n})}\cdot(N_{salt}+1)\\\
&=\frac{c^{ins}_{z}(r_{z})}{\sum_{n=1}^{N_{salt}+1}c^{ins}_{z}(r_{z,n})}\cdot(N_{salt}+1).\end{split}$
(24)

The second equality is due to the random choice of the flying-cation position
along $x$ and $y$ directions. Again, $C^{del}(\vec{r}_{c})=1$ for randomly
assigned flying cation 3D positions. The flying anion is then chosen according
to Eq. 21. We note that such a bias $C(\vec{r}_{c})$ can also be applied to
the conventional instantaneous MC methods, since it is related to the
preparation of the positions of flying ions.

Once both the cation and anion have been chosen as described above, the
probability to accept the trial deletion ($N_{salt}+1\rightarrow N_{salt}$)
is:

$\begin{split}f_{des}(\vec{r}_{a},\vec{r}_{c})&=\min\bigg{[}1,\exp(-\beta\Delta
H_{N_{salt}+1\rightarrow
N_{salt}})\exp(-\beta\mu_{salt})\bigg{(}(N_{salt}+1)\frac{\Lambda_{s}^{3}}{L_{x}L_{y}H}\bigg{)}^{2}\frac{B^{ins}(\vec{r}_{a}|\vec{r}_{c})}{B^{del}(\vec{r}_{a}|\vec{r}_{c})}\frac{C^{ins}(\vec{r}_{c})}{C^{del}(\vec{r}_{c})}\bigg{]}\\\
&\equiv\min\bigg{[}1,\exp(+\beta\Delta M-\beta\mu_{salt})\bigg{]}.\end{split}$
(25)

Similarly, the trial insertion/deletion of solvent particles are generated
with the same biases Eqs. 22 and 24 and accepted with probabilities

$\begin{split}f_{ins}(\vec{r}_{s})&=\min\bigg{[}1,\exp(-\beta\Delta
H_{N_{solv}\rightarrow
N_{solv}+1})\exp(\beta\mu_{solv})\bigg{(}\frac{L_{x}L_{y}H}{\Lambda_{s}^{3}}\frac{1}{N_{solv}+1}\bigg{)}\frac{C^{del}(\vec{r}_{s})}{C^{ins}(\vec{r}_{s})}\bigg{]},\\\
&\equiv\min\bigg{[}1,\exp(-\beta\Delta M+\beta\mu_{solv})\bigg{]},\end{split}$
(26)

for insertion and

$\begin{split}f_{des}(\vec{r}_{s})&=\min\bigg{[}1,\exp(-\beta\Delta
H_{N_{solv}+1\rightarrow
N_{solv}})\exp(-\beta\mu_{solv})\bigg{(}(N_{solv}+1)\frac{\Lambda_{s}^{3}}{L_{x}L_{y}H}\bigg{)}\frac{C^{ins}(\vec{r}_{s})}{C^{del}(\vec{r}_{s})}\bigg{]}\\\
&\equiv\min\bigg{[}1,\exp(+\beta\Delta M-\beta\mu_{solv})\bigg{]}.\end{split}$
(27)

for deletion, respectively. Here, $N_{solv}$ is the number of solvent atoms
and $\mu_{solv}$ is their chemical potential. In this work, we set
$\Lambda_{s}=\sqrt{\Lambda_{+}\Lambda_{-}}=1$ in unit LJ length, and the sign
of $\Delta M$ in both trial moves follows the direction of trial insertions as
in Refs. 75 and 76.

Calculation of chemical potentials Grand-canonical simulations require the
knowledge of the chemical potential of the exchanged species (here the salt,
i.e. neutral ion pairs, and solvent particles), which is set by the
composition and the thermodynamic conditions in the bulk reservoir.
Configuration sampling to calculate the chemical potentials was done in the
$N_{solv}N_{salt}p^{*}T^{*}$ ensemble with $N_{solv}=5000$ and $N_{salt}=100$
at $p^{*}=p\sigma^{3}/\epsilon=1$ and $T^{*}=Tk_{B}/\epsilon=1$. The desired
pressure $p^{*}=1$ and temperature $T^{*}=1$ were maintained using the Nosé-
Hoover barostat and thermostat, with time constants of 1000 and 100 LJ units,
respectively. As described in Refs. 75 and 76, the chemical potentials were
computed using Crooks’ fluctuation theorem [82] and Bennett’s acceptance ratio
method [83, 84], resulting in: $\beta^{*}\mu^{*}_{solv}=-1.72\pm 0.06$, and
$\beta^{*}\mu^{*}_{salt}=-11.74\pm 0.09$ with $\beta^{*}=1/T^{*}$ (see Fig.S1
in SI). For the implicit solvent models, $\beta^{*}\mu^{*}_{salt,LJ}=-9.38\pm
0.01$ and $\beta^{*}\mu^{*}_{salt,WCA}=-8.59\pm 0.01$ for LJ and WCA
interactions between ions, respectively.

###  III.3 GCMD simulations

GCMD simulations were conducted in the
$\mu^{*}_{solv}\mu^{*}_{salt}A^{*}H^{*}T^{*}$ ensemble, i.e. the numbers of
both solvent and ion pairs fluctuate in a fixed confined space and their
statistics (in particular their average) are determined by their chemical
potentials (Fig. 2). Here, $A^{*}=L^{*}_{x}L^{*}_{y}$ and $H^{*}$ is the
distance between the wall atoms in each of the inner-most wall layers. GCMD
simulations are composed of two steps, including equilibrium MD to sample
configurations and non-equilibrium MD/MC steps to exchange solvent or ion pair
with reservoir. During equilibrium MD ($N_{eq}$ steps), the equations of
motion are integrated using the velocity Verlet algorithm with a timestep
$\delta t^{*}=\delta t\sqrt{\epsilon(m\sigma^{2})^{-1}}=0.005$. The desired
temperature $T^{*}=1$ is maintained using the Nosé-Hoover thermostat, with a
time constant of 100 LJ units.

During NEMD ($N_{f}$ steps), performed in the $NVE$ ensemble, we also use the
velocity Verlet algorithm but with $\delta t^{*}=0.01$, and the maximum
altitude $w^{*}_{max}$ (see Eq. 18) set to unity. While this may not seem the
most straightforward choice, we found it beneficial for the acceptance rate of
both solvent and ion-pair exchange in bulk electrolytes to use the 3D PPPM
estimate of the long-range part of electrostatic interactions (instead of the
4D expression in the presence of flying particles): In addition to its
simplicity of implementation, it results in higher acceptance rates for both
solvent and ion-pair exchanges in bulk electrolytes [76]. The short-range
Coulomb contribution was nevertheless calculated in 4D. Our implementation of
the H4D method for confined electrolytes in the LAMMPS simulation package [85]
is freely available as an update of our previous implementation of the H4D
method for bulk electrolytes [76] (see Data Availability statement). The
optimal choice of $N_{eq}$ and $N_{f}$ is system-specific, and we used
$N_{eq}=100$ and $N_{f}=800$ in most cases.

Figure 2: Illustration of grand canonical (GC) MD simulation of confined
electrolytes with explicit solvent particles in this study. GCMD was conducted
using the H4D method, extended to confined electrolytes, in order to enhance
the acceptance rate of the salt-pair and solvent exchange trial moves.
Chemical potentials used in this study were also computed using H4D:
$\beta^{*}\mu^{*}_{solv}=-1.72\pm 0.06$, and
$\beta^{*}\mu^{*}_{salt}=-11.74\pm 0.09$. See text for details.

NEMD details. We apply a bimodal biasing function for the flying-ion
distances, as described in Ref. 76. There are several choices for such a
biasing function (e.g., a Gaussian distribution with a 3D distance between
flying ions). Here, we use a bimodal distribution function for $b^{ins}$,
which helps to generate a pair of flying ions that are close to each other yet
without large overlap:

$b^{ins}(x_{a}|x_{c})=\sqrt{\frac{\alpha_{b}}{\pi}}\frac{\exp(-\alpha_{b}(x_{ac}-x_{b})^{2})+\exp(-\alpha_{b}(x_{ac}+x_{b})^{2})}{\text{erf}[\sqrt{\alpha_{b}}(x_{b}+L_{x}/2)]-\text{erf}[\sqrt{\alpha_{b}}(x_{b}-L_{x}/2)]},$
(28)

where $x_{ac}=x_{a}-x_{c}$. $x_{b}$ determines the mean separation of flying
ions, and $\alpha_{b}$ determines the sharpness of the biasing function.
Again, $b^{ins}$ is normalized such that
$\int_{-L_{x}/2}^{L_{x}/2}b^{ins}(x_{a}|x_{c})dx_{a}=1$. In this work, we used
$\alpha^{*}_{b}=0.5$ and $x^{*}_{b}=2$. The same functional form and
parameters are also used in the $y$ and $z$ directions. A small excluded
volume was also applied: $V^{*}_{ex}=0.125$ for early rejection (See details
in Ref. 76). In the case of implicit solvent models, trial insertions and
deletions of the ions are carried out instantaneously (this corresponds to an
infinite vertical velocity, with $N_{f}=0$) as in the conventional GCMC. For
the bias in the confined systems (see Eq. 23), $\alpha^{*}_{c}$ varies from
0.5 to 6 with decreasing separation distance. In order to avoid numerical
instabilities, one can increase $\alpha^{*}_{c}$ to obtain a narrow Gaussian
distribution near the center of the slit pore.

###  III.4 Computation of observables

We analyze the composition and structure of the confined liquid by computing
the ionic and solvent density profiles, $\rho^{*}_{\pm}(z^{*})$ and
$\rho^{*}_{\text{solv}}(z)$, using histograms with bin width $\Delta
z^{*}=0.1$, and symmetrized with respect to the pore center ($z^{*}=0$). Since
only neutral ion pairs or solvent particles are exchanged in GCMD, the
electroneutrality constraint

$\int_{-H^{*}/2}^{H^{*}/2}dz^{*}\rho^{*}_{q}(z^{*})=-2\Sigma^{*}_{s},$ (29)

with $\rho^{*}_{q}(z^{*})=\rho^{*}_{+}(z^{*})-\rho^{*}_{-}(z^{*})$ the charge
density, is always satisfied. We also compute the cumulative mean-charge
density from the charged surface

$Q^{*}(z^{*})=\Sigma^{*}_{s}+\int_{-H^{*}/2}^{z^{*}}\rho^{*}_{q}\big{(}z^{*}\big{)}dz^{*},$
(30)

which decays from $\Sigma^{*}_{s}$ at $z^{*}=-H^{*}/2$ to 0 at $z^{*}=0$ (the
centre of the slit-pore), if $H^{*}\gg\lambda^{*}_{D}$.

The decay of $Q^{*}(z^{*})$ also provides an estimate of the thickness of the
electrical double layer at the surface. Due to the short-range repulsion
between the ions and the wall, there is a region from which ions are excluded
where $\rho^{*}_{q}=0$ and $Q^{*}(z^{*})=\Sigma^{*}_{s}$. According to the
linearized PB equation, beyond this exclusion region, for dilute electrolytes
with weakly charged walls we expect an exponential decay as

$Q^{*}(z^{*}_{d})\approx\Sigma^{*}_{s}\exp[-(z^{*}_{d}-s^{*})/\lambda^{*}_{D}],$
(31)

where $z^{*}_{d}$ is the distance from the charged wall, and $s^{*}$ is an
effective width of the exclusion layer. We note again that this approximation
assumes $H^{*}\gg\lambda^{*}_{D}$, which is satisfied in all the cases
considered in the present work. As explained in Section II.3, for larger
$\Sigma^{*}_{s}$ and sufficiently large distances from the surface, one can
approximate the decay using an effective surface charge density
$\Sigma^{*}_{\text{eff}}=a_{\text{eff}}\Sigma^{*}_{s}$:

$Q^{*}(z^{*}_{d})\approx\Sigma^{*}_{\text{eff}}\exp[-(z^{*}_{d}-s^{*})/\lambda^{*}_{D}]$
(32)

with $\Sigma^{*}_{\text{eff}}=a_{\text{eff}}\Sigma^{*}_{s}$. Similar
expressions are obtained for the corresponding potential and ionic density
profiles. As a result of the exclusion of ions from both walls due to their
finite size, the effective width of the slit-pore is not the bare distance
$H^{*}$ between the inner-most layers of the confining walls, but an effective
size $H^{*}_{\text{eff}}=H^{*}-2s^{*}$. Based on the simulation data (see
Section IV), in the following we estimate $H^{*}_{\text{eff}}$ using
$s^{*}=0.9$ whenever comparing the simulation results to analytical
predictions involving the pore size.

##  IV Results and discussion

In this section, we discuss the results of GCMD simulations for the LJ
electrolytes with explicit solvent particles confined in a slit-like pore with
negatively charge surfaces. The system is in an equilibrium with a reservoir
setting the chemical potentials of the salt and of solvent particles, whose
values are computed as described in Section III.2. The salt concentration in
this reservoir ($\rho^{*}_{bulk}$=0.0144) is sufficiently small for Debye-
Hückel theory to apply, and the corresponding Debye screening length is
$\lambda^{*}_{D}=1.66$. As shown in Fig. S2 of the SI, we observe in the bulk
simulations that the charge-charge correlation function indeed decays
exponentially with the expected decay length.

Figure 3:  Local densities of cations $\rho^{*}_{+}$ (a-d), anions
$\rho^{*}_{-}$ (e-h), and solvent $\rho^{*}_{\text{solv}}$ (i-l) confined in
charged slit pores of different widths $H^{*}$ and surface charge densities
$\Sigma^{*}_{s}$, and the corresponding cumulative mean-charge density
$Q^{*}(z^{*})$ (m-p) defined by Eq. 30. Colors represent different $H^{*}$, as
indicated in the legend on the right. Horizontal dashed lines in
$\rho^{*}_{+}$ and $\rho^{*}_{-}$ represent the bulk salt concentration
$\rho^{*}_{\text{salt}}$ in the reservoir. In panels (m-p), the dotted and
dashed lines are the predictions from linearized PB theory with the bare and
effective surface charge densities, given by Eqs. 31 and 32, respectively.

Figure 3 displays the local densities of the cations (panels a-d), anions
(panels e-h), and solvent particles (panels i-l) confined in slit-pores with
different sizes ($H^{*}=12,~{}18,~{}24,~{}\text{and}~{}36$) and surface charge
densities ($\Sigma^{*}_{s}=-0.05,~{}-0.1,~{}-0.2,~{}\text{and}~{}-0.5$). Each
column displays the profiles for the same $\Sigma^{*}$, with the positions
$z^{*}$ shifted by $H^{*}/2$ to align the positions of the “left” walls (at
$z^{*}=-H^{*}/2$), and the colors in each panel correspond to the various pore
sizes. For all the considered $\Sigma^{*}_{s}$, the shifted profiles for the
various $H^{*}$ overlap. Deviations from the bulk behaviors emerge only at the
surfaces under these mild confinement conditions, for which
$H^{*}\gg\lambda^{*}_{D},l_{GC}$. The cumulative charge distribution, computed
from Eq. 30, is also reported in Figure 3 (panels m-p).

###  IV.1 Weakly charged pores

The leftmost column of Fig. 3 (a,e,i,m) displays the results of dilute LJ
electrolytes confined in a weakly charged pore ($\Sigma^{*}_{s}=-0.05$ and
$l^{*}_{GC}=3.2$). As expected, electrical double layers form at the charged
surfaces, with an excess of counterions (panel a) and a depletion of co-ions
(panel e) near the walls, with respect to the salt concentration far from
them. In addition, the ionic density profiles display oscillations near the
wall due to packing effects dominated by the solvent density (panel i), which
is much larger than that of the ions and displays at least three well defined
layers. This is in stark contrast to implicit solvent models where no such
layering is observed in the ion density profiles (see Section IV.4). The
cumulative mean-charge distribution $Q^{*}(z^{*})$ (panel m) follows an
exponential decay with characteristic length $\lambda^{*}_{D}$, well described
by the linearized PB result (Eq. 31), to zero at the center of the pore. In
this case the difference between the bare and renormalized surface charge
densities is small ($a_{\text{eff}}\approx 0.94$). Finally, we note that
$Q^{*}(z^{*})$ does not display any oscillations despite that of the ionic
density profiles. Such an observation is likely due to the simplicity of the
solvent model, which does not include an explicit charge distribution (dipolar
or more complex) which would contribute to $Q^{*}(z^{*})$ [86, 87, 88, 89,
42].

###  IV.2 Strongly charged pores

With more negative $\Sigma^{*}_{s}$ (and decreasing $l^{*}_{GC}$), the
structure of the EDLs evolves with further enrichment of cations (Figs. 3a-d)
and depletion of anions (Figs. 3e-h) near the wall. Due to the finite size of
the ions and solvent molecules, the increasing density of cations leads to the
progressive depletion of solvent particles at the interface (Figs. 3i-l). For
$\Sigma^{*}_{s}=-0.5$, the cation density is even larger than that of the
solvent in the first fluid layer in contact with the surface. Since the
solvent density near the wall is depleted (even more so that
$|\Sigma^{*}_{s}|$ is large) due to the increased density of counterions, so
is the average solvent density inside the pore (See Fig. S3 in the SI).

For large surface charge densities, the decay of $Q^{*}(z^{*})$ (Figure 3m-p)
departs significantly from the Debye-Hückel result (Eq. 31). Close to the
wall, it decays much faster due to the large concentration of ions in the
first layer, consistently with the Gouy-Chapman-Stern picture of the
interface. Beyond this layer, the decay is well described by the linearized PB
result with a renormalized surface charge density (Eq. 32 and the dashed lines
in Figure 3o-p). The failure from the bare DH prediction is not surprising,
since this is a regime where ion-wall interactions are large compared to ion-
ion interactions ($l^{*}_{GC}\lesssim\lambda^{*}_{D}$) [20, 21]. Importantly,
such an agreement suggests that (i) the effective charge describing the far-
field is well predicted by Eq. 12, and (ii) the decay length at long distance
is well described by $\lambda^{*}_{D}$.

###  IV.3 Donnan equilibrium

Figure 4:  Mean cation ($\bar{\rho}^{*}_{+}$, squares) and anion
($\bar{\rho}^{*}_{-}$, circles) densities in slit pores with different widths
$H^{*}$ and different surface charge densities: $\Sigma^{*}_{s}$ = -0.05 (a),
-0.1 (b), -0.2 (c), and -0.5 (d). Results are plotted as a function of the
inverse of the effective width of the pore, $H^{*}_{\text{eff}}=H^{*}-2s^{*}$,
where the parameter $s^{*}=0.9$ accounts for the excluded volume between the
liquid and the wall atoms (see Section III.4). The dotted and dashed blue
lines are the predictions from linearized PB theory with the bare (Eqs. 7 and
8) and renormalized (Eqs. 13 and 14) surface charge densities, respectively.
The horizontal black dashed line represents the bulk ion density
$\rho^{*}_{\text{salt}}$ in the reservoir.

The different composition of the confined solution with respect to the bulk
reservoir, which characterizes the Donnan equilibrium, can be quantified by
the mean cation and anion densities insides the pore. Since the linearized PB
theory predicts linear scalings as a function of the inverse of the distance
between the walls (see Eqs. 7 and 8), the results are reported in Fig. 4 as a
function of the inverse of the effective pore width,
$H^{*}_{\text{eff}}=H^{*}-2s^{*}$, where the parameter $s^{*}=0.9$ accounts
for the excluded volume between the liquid and the wall atoms (see Section
III.4). The predictions (Eqs. 7 and 8) are shown as dotted lines. They are
accurate for the smallest considered surface charge densities, but as
$|\Sigma^{*}_{s}|$ increases, they underestimate both the cation and anion
mean densities. Nevertheless, even for the larger values considered, the
scaling remains linear in $1/H^{*}_{\text{eff}}$ (which validates its choice
of the relevant definition of the pore width as well as the corresponding
value of $s^{*}$) and consistent with the limit
$\rho^{*}_{\pm}\to\rho^{*}_{salt}$ for very large pores (where interfacial
effects become negligible). In fact, all results are very well described by
the predictions of linearized PB with the effective charge
$\Sigma^{*}_{\text{eff}}=a_{\text{eff}}\Sigma^{*}_{s}$ (see Eqs. 13 and 14),
shown as dashed blue lines. This good agreement reflects the better accuracy
of the non-linear PB theory, whose result is implicitly taken into account via
the renormalized surface charge density, at least in the present regime where
$\lambda^{*}_{D}$ is small compared to $H^{*}_{\text{eff}}$.

###  IV.4 Comparison to implicit solvent models

Figure 5:  Local densities of cations $\rho^{*}_{+}$ and anions $\rho^{*}_{-}$
for the two implicit-solvent electrolyte models (implicit LJ, top, and WCA,
bottom; see Section III.1) confined in a slit of different width $H^{*}$ and
surface charge densities $\Sigma^{*}_{s}$, and the corresponding cumulative
mean-charge density $Q^{*}(z^{*})$ defined by Eq. 30. Colors represent
different $H^{*}$, as indicated in the legend on the right of Fig. 3.
Horizontal dashed lines in $\rho^{*}_{+}$ and $\rho^{*}_{-}$ represent the
bulk salt concentration $\rho^{*}_{\text{salt}}$ in the reservoir, while for
$Q^{*}$, the dotted and dashed lines are the predictions from linearized PB
theory with the bare and effective surface charge densities, given by Eqs. 31
and 32, respectively.

In order to investigate the solvent packing effect, we conducted additional
GCMD simulation of two dilute electrolytes with explicit ions (with short-
range LJ or WCA interactions, see Section III.1) and no explicit solvent
particles (only a permittivity screening the Coulomb interactions between
ions, as in the model with explicit solvent considered here). For bulk
electrolytes, all three models predict long-range screening behaviors well
described by Debye-Hückel theory, despite the differences in short-range
structures in ion-ion radial distribution functions (See Fig. S2 in SI).

Interfacial structure. In both implicit cases, local ion densities inside the
pores (Fig. 5) show no oscillations at all $\Sigma^{*}_{s}$’s, which further
supports the conclusion that the oscillations with an explicit solvent in Fig.
3 originate from the solvent packing effect. For the three smaller densities
considered ($|\Sigma^{*}_{s}|\leq 0.2$), the simulation results are very
similar to each other and to the explicit-solvent case, except for the above-
mentioned oscillations. In particular, the cumulative charge density is well
described by the linearized PB prediction with the bare charge for the lowest
$|\Sigma^{*}_{s}|$, and at long distance with the renormalized surface charge
for intermediate $|\Sigma^{*}_{s}|$ (but not at short distance).

The main difference between the two implicit models is observed for
$\Sigma^{*}_{s}=-0.5$. While the predictions of the implicit WCA model are
very close to the explicit solvent case, the implicit LJ model predicts a peak
in the anion concentration close to the cations adsorbed at the surface (at a
distance from the surface corresponding approximately to the shoulder of
$\rho^{*}_{+}$ next to the first peak). Such an excess of anions close to the
adsorbed cations results from the short-range attraction between ions, which
is absent in the implicit WCA case. Of course, since such ionic correlations
are not included in the PB theory, the latter is insufficient to predict the
cumulative charge, even at long distances since the presence of anions close
to the surface also modifies the electric field experienced by ions far from
it.

Excess ion density. The structure of EDL is closely related to the Donnan
equilibrium between the charged pore and the bulk reservoir. The predictions
of the three electrolyte models for the excess ion density
$\Delta\rho^{*}_{\text{ex}}$ with respect to the bulk reservoir (see Eq. 15)
are reported with different colors in Fig. 6, as a function of the inverse of
the effective width of the pore, $H^{*}_{\text{eff}}=H^{*}-2s^{*}$, for all
the considered surface charge densities, from $\Sigma^{*}_{s}=-0.05$ (circles)
to $\Sigma^{*}_{s}=-0.5$ (squares) (The mean cation and anion densities for
the implicit solvent models can be found in Fig. S4 in SI). For weakly charged
surfaces ($l^{*}_{GC}\gtrsim\lambda^{*}_{D}$), $a_{\text{eff}}\approx 1$ and
the excess ion density is small ($\Delta\rho^{*}_{\text{ex}}\approx 0$). For
larger $|\Sigma^{*}_{s}|$, both implicit models follow the linear trend
observed with the explicit solvent and predicted by the linearized PB theory
using a renormalized charge (see Eq. 15), also shown as dashed lines. However,
while the implicit WCA model closely follows the predictions with the explicit
solvent, the implicit LJ model consistently overestimates the salt
concentration, even more so that $|\Sigma^{*}_{s}|$ is large.

Figure 6: Mean excess ion density $\Delta\rho^{*}_{\text{ex}}$ of three
electrolytes models: Blue for the explicit solvent model, red for the implicit
WCA model, and green for the implicit LJ model. Four different markers
represent four different surface charge densities $\Sigma^{*}_{s}$. Results
are plotted as a function of the inverse of the effective width of the pore,
$H^{*}_{\text{eff}}=H^{*}-2s^{*}$, where the parameter $s^{*}=0.9$ accounts
for the excluded volume between the liquid and the wall atoms (see Section
III.4). The black dashed lines represent the prediction (Eq. 15) of the
linearized PB theory with a renormalized surface charge density
$a_{\text{eff}}\Sigma^{*}_{s}$ (see Eq. 12).

In order to better quantify the differences between the three electrolyte
models, we fit the simulation results of Fig. 6 by Eq. 15 treating the factor
$a_{\text{eff}}$ as a fitting parameter. The results are summarized in Table
1, together with the analytical predictions Eq. 12, or introducing the
nonlinear solution for the potential in the thin EDL approximation, Eq. 9, in
the Boltzmann distribution of the ions, and integrating the latter
numerically. Consistently with Fig. 6, these results indicate that Eq. 12 is
only accurate for sufficiently small surface charge density, and overestimates
$a_{\text{eff}}$ (hence underestimates $\Delta\rho^{*}_{\text{ex}}$) for
larger $|\Sigma^{*}_{s}|$. In contrast, the non-linear thin-EDL approximation
is accurate for all the considered cases (for which $\lambda^{*}_{D}$ is
smaller than $H^{*}$, as already noted). This suggests that the limitations of
Eq. 12 are not due to the breakdown of the mean-field description, but rather
to the charge renormalization approximation, which is not sufficient to
capture the overall excess ion density because it focuses on the behavior far
from the walls.

$\Sigma^{*}_{s}$ | $a_{\text{eff}}$ | $a_{\text{eff}}$ | $a_{\text{eff,fit}}$ | $a_{\text{eff,fit}}$ | $a_{\text{eff,fit}}$
---|---|---|---|---|---
| (Eq. 12) | (Eq. 9) | (explicit) | (implicit WCA) | (implicit LJ)
-0.05 | 0.94 | 0.75 | 0.73 | 0.99 | 0.95
-0.1 | 0.82 | 0.58 | 0.57 | 0.66 | 0.55
-0.2 | 0.61 | 0.38 | 0.37 | 0.44 | 0.30
-0.5 | 0.31 | 0.18 | 0.16 | 0.21 | 0.015
Table 1:  Charge renormalization factor $a_{\text{eff}}$ predicted by Eqs. 12
and 9 (see text for detail) and by numerically fitting the simulation results
for the excess ion density $\Delta\rho^{*}_{\text{ex}}$ for the three
electrolyte models (see Fig. 6) to Eq. 15.

Table 1 also confirms that the implicit WCA model overall provides a better
agreement with the explicit solvent than the implicit LJ model. The most
significant difference is observed for the largest surface charge density,
$\Sigma^{*}_{s}=-0.5$, where the implicit LJ model predicts
$a_{\text{eff}}\approx 0$, i.e. a perfect screening of the surface charge by
the nearby ions, which also include co-ions in that case (see Fig. 5). We
further note that when the field due to the surface charge is completely
screened, the salt concentration inside the pore is equal to that of the
reservoir (i.e. no Donnan exclusion). Finally, the good agreement between the
explicit solvent and implicit WCA model indicates that packing effects do not
play a significant role in the Donnan exclusion in the considered cases (even
though it does lead to oscillations in the ionic density profiles, as
discussed above).

##  V Conclusion

We investigated the Donnan equilibrium of coarse-grained dilute electrolytes
confined in charged slit-pores in equilibrium with a reservoir of ions and
solvent. This was achieved using an extension of a recently developed hybrid
Grand canonical / nonequilibrium molecular dynamics simulation method (J.
Chem. Phys. 151, 021101 (2019) and J. Chem. Phys. 159, 144802 (2023)) to
confined systems, which enhances the efficiency of solvent and ion-pair
exchange, using insertion/deletion via a fourth spatial dimension. We
restricted ourselves to the case of dilute electrolytes in the thin electric
double-layer limit and explored the influences of the pore size and of the
surface charge density, as well as packing effects due to an explicit solvent
by comparison with implicit solvent models, on the co- and counter-ion
distributions inside the pore, the cumulative charge density profiles at the
interface, and the mean ion concentrations of the confined liquid, which
differ from that in the electrolyte reservoir.

We showed that the validity range of linearized Poisson-Boltzmann theory to
predict the Donnan equilibrium of dilute electrolytes can be extended to
highly charged pores, by considering renormalized surface charge densities,
which can be computed analytically by comparing the full and linearized PB
equations, instead of the bare ones. By comparing with simulations of implicit
solvent models of electrolytes, we find that for the small salt concentration
considered here, an explicit solvent introduces oscillations in the ionic
density profiles, but has a limited effect on the excess salt concentration
inside the pore. In the low concentration and thin electric double-layer limit
considered here, the main limitations of the analytical predictions are not
due to the breakdown of the mean-field description, but rather to the charge
renormalization approximation, which is not sufficient to capture the overall
excess ion density because it focuses on the behavior far from the walls.

In the models of electrolytes considered in the present work, the solvent is
taken into account by its permittivity to screen the electrostatic
interactions between ions, with or without explicit solvent particles, is very
crude. While it has the advantage of allowing us to disentangle the effects of
long-range electrostatic interactions from short-range van der Waals
interactions, this is a very crude representation since the charge
distribution within solvent molecules results in more complex interactions
between themselves, with ions and with charged walls. Furthermore, the
computation of electrostatic interactions in the confined case assumes that
there is no dielectric contrast between the liquid and the wall, which is the
exception rather than the rule. The present Grand canonical simulation using
the H4D can also be applied with an explicit solvent such as the SPC/E water
model [75, 76], as well as more concentrated electrolytes than the one
considered here, or multivalent ions. This would provide insights into the
effects of ion-ion and ion-solvent correlations on the properties of the
confined electrolytes, and to explore the possible role of the Donnan
equilibrium in the observed scaling of Surface Force Balance measurements with
the concentration [30, 31, 32]. In addition, the present extension of the H4D
method to slit pores could be further extended to more complex geometries,
including disordered porous materials such as ion exchange membranes.

## Acknowledgements

We thank Luc Belloni for fruitful discussions. This project received funding
from the European Research Council under the European Union’s Horizon 2020
research and innovation program (grant agreement no. 863473).

## Data Availability Statement

Our implementation of the H4D method in the LAMMPS simulation package is
freely available at https://github.com/Jeongmin0658/h4d_lammps. All the data
presented in this work will be provided upon reasonable requests.

## References

  * Donnan [1911] F. G. Donnan, “The theory of membrane equilibrium in the presence of a non-dialyzable electrolyte,” Z. Electrochem 17, 572–581 (1911).
  * Nelson  _et al._ [2020] P. Nelson, D. Goodsell, K. Chen, and S. Bromberg,  _Biological Physics: Energy, Information, Life_ (Chiliagon Science, 2020).
  * Jardat  _et al._ [2009] M. Jardat, J.-F. Dufrêche, V. Marry, B. Rotenberg, and P. Turq, “Salt exclusion in charged porous media: a coarse-graining strategy in the case of montmorillonite clays,” Physical Chemistry Chemical Physics 11, 2023–2033 (2009).
  * Barry  _et al._ [2021] E. Barry, R. Burns, W. Chen, G. X. De Hoe, J. M. M. De Oca, J. J. De Pablo, J. Dombrowski, J. W. Elam, A. M. Felts, G. Galli, _et al._ , “Advanced materials for energy-water systems: the central role of water/solid interfaces in adsorption, reactivity, and transport,” Chemical Reviews 121, 9450–9501 (2021).
  * Hsiao and Hedström [2022] Y.-W. Hsiao and M. Hedström, “Salt effect on donnan equilibrium in montmorillonite demonstrated with molecular dynamics simulations,” The Journal of Physical Chemistry B 126, 8873–8881 (2022).
  * Ilgen  _et al._ [2023] A. G. Ilgen, K. Leung, L. J. Criscenti, and J. A. Greathouse, “Adsorption at nanoconfined solid–water interfaces,” Annual Review of Physical Chemistry 74, 169–191 (2023).
  * Briskot  _et al._ [2022] T. Briskot, N. Hillebrandt, S. Kluters, G. Wang, J. Studts, T. Hahn, T. Huuk, and J. Hubbuch, “Modeling the Gibbs–donnan effect during ultrafiltration and diafiltration processes using the poisson–boltzmann theory in combination with a basic Stern model,” Journal of Membrane Science 648, 120333 (2022).
  * Galama  _et al._ [2013] A. Galama, J. Post, M. C. Stuart, and P. Biesheuvel, “Validity of the boltzmann equation to describe donnan equilibrium at the membrane–solution interface,” Journal of membrane science 442, 131–139 (2013).
  * Siria, Bocquet, and Bocquet [2017] A. Siria, M.-L. Bocquet, and L. Bocquet, “New avenues for the large-scale harvesting of blue energy,” Nature Reviews Chemistry 1, 0091 (2017).
  * Gao  _et al._ [2022] K. W. Gao, X. Yu, R. M. Darling, J. Newman, and N. P. Balsara, “Increased donnan exclusion in charged polymer networks at high salt concentrations,” Soft Matter 18, 282–292 (2022).
  * Aydogan Gokturk  _et al._ [2022] P. Aydogan Gokturk, R. Sujanani, J. Qian, Y. Wang, L. E. Katz, B. D. Freeman, and E. J. Crumlin, “The donnan potential revealed,” Nature communications 13, 5880 (2022).
  * Kamcev and Freeman [2016] J. Kamcev and B. D. Freeman, “Charged polymer membranes for environmental/energy applications,” Annual review of chemical and biomolecular engineering 7, 111–133 (2016).
  * Shin  _et al._ [2022] S.-J. Shin, D. H. Kim, G. Bae, S. Ringe, H. Choi, H.-K. Lim, C. H. Choi, and H. Kim, “On the importance of the electric double layer structure in aqueous electrocatalysis,” Nature communications 13, 174 (2022).
  * Xu  _et al._ [2024] P. Xu, R. Wang, H. Zhang, V. Carnevale, E. Borguet, and J. Suntivich, “Cation modifies interfacial water structures on platinum during alkaline hydrogen electrocatalysis,” Journal of the American Chemical Society 146, 2426–2434 (2024).
  * Biesheuvel  _et al._ [2011] P. Biesheuvel, R. Zhao, S. Porada, and A. Van der Wal, “Theory of membrane capacitive deionization including the effect of the electrode pore space,” Journal of colloid and interface science 360, 239–248 (2011).
  * Misra  _et al._ [2019] R. P. Misra, J. P. de Souza, D. Blankschtein, and M. Z. Bazant, “Theory of surface forces in multivalent electrolytes,” Langmuir 35, 11550–11565 (2019).
  * Chen [2022] G. Chen, “Donnan equilibrium revisited: Coupling between ion concentrations, osmotic pressure, and donnan potential,” Journal of Micromechanics and Molecular Physics 7, 127–134 (2022).
  * Parsons [1990] R. Parsons, “The electrical double layer: recent experimental and theoretical developments,” Chemical Reviews 90, 813–826 (1990).
  * Bard, Faulkner  _et al._ [2001] A. J. Bard, L. R. Faulkner, _et al._ , “Fundamentals and applications,” Electrochemical methods 2, 580–632 (2001).
  * Markovich, Andelman, and Podgornik [2021] T. Markovich, D. Andelman, and R. Podgornik, “Charged membranes: Poisson–boltzmann theory, the dlvo paradigm, and beyond,” in  _Handbook of lipid membranes_ (CRC Press, 2021) pp. 99–128.
  * Herrero and Joly [2021] C. Herrero and L. Joly, “Poisson-boltzmann formulary,” arXiv preprint arXiv:2105.00720 (2021).
  * Park  _et al._ [2006] C. Park, P. A. Fenter, K. L. Nagy, and N. C. Sturchio, “Hydration and distribution of ions at the mica-water interface,” Physical review letters 97, 016101 (2006).
  * Ricci, Spijker, and Voïtchovsky [2014] M. Ricci, P. Spijker, and K. Voïtchovsky, “Water-induced correlation between single ions imaged at the solid–liquid interface,” Nature communications 5, 4400 (2014).
  * Brown, Goel, and Abbas [2016] M. A. Brown, A. Goel, and Z. Abbas, “Effect of electrolyte concentration on the Stern layer thickness at a charged interface,” Angewandte Chemie 128, 3854–3858 (2016).
  * Brown  _et al._ [2016] M. A. Brown, Z. Abbas, A. Kleibert, R. G. Green, A. Goel, S. May, and T. M. Squires, “Determination of surface potential and electrical double-layer structure at the aqueous electrolyte-nanoparticle interface,” Physical Review X 6, 011007 (2016).
  * Hartkamp  _et al._ [2018] R. Hartkamp, A.-L. Biance, L. Fu, J.-F. Dufrêche, O. Bonhomme, and L. Joly, “Measuring surface charge: Why experimental characterization and molecular modeling should be coupled,” Current opinion in colloid & interface science 37, 101–114 (2018).
  * Ma  _et al._ [2021] E. Ma, J. Kim, H. Chang, P. E. Ohno, R. J. Jodts, T. F. Miller III, and F. M. Geiger, “Stern and diffuse layer interactions during ionic strength cycling,” The Journal of Physical Chemistry C 125, 18002–18014 (2021).
  * Lee  _et al._ [2021] S. S. Lee, A. Koishi, I. C. Bourg, and P. Fenter, “Ion correlations drive charge overscreening and heterogeneous nucleation at solid–aqueous electrolyte interfaces,” Proceedings of the National Academy of Sciences 118, e2105154118 (2021).
  * Hunger  _et al._ [2022] J. Hunger, J. Schaefer, P. Ober, T. Seki, Y. Wang, L. Prädel, Y. Nagata, M. Bonn, D. J. Bonthuis, and E. H. Backus, “Nature of cations critically affects water at the negatively charged silica interface,” Journal of the American Chemical Society 144, 19726–19738 (2022).
  * Smith, Lee, and Perkin [2016] A. M. Smith, A. A. Lee, and S. Perkin, “The electrostatic screening length in concentrated electrolytes increases with concentration,” The journal of physical chemistry letters 7, 2157–2163 (2016).
  * Lee  _et al._ [2017] A. A. Lee, C. S. Perez-Martinez, A. M. Smith, and S. Perkin, “Underscreening in concentrated electrolytes,” Faraday Discussions 199, 239–259 (2017).
  * Elliott  _et al._ [2024] G. R. Elliott, K. P. Gregory, H. Robertson, V. S. J. Craig, G. B. Webber, E. J. Wanless, and A. J. Page, “The known-unknowns of anomalous underscreening in concentrated electrolytes,” Chemical Physics Letters 843, 141190 (2024).
  * Torrie and Valleau [1980] G. Torrie and J. Valleau, “Electrical double layers. i. Monte Carlo study of a uniformly charged surface,” The Journal of Chemical Physics 73, 5807–5816 (1980).
  * Torrie and Valleau [1982] G. Torrie and J. Valleau, “Electrical double layers. 4. limitations of the gouy-chapman theory,” The Journal of Physical Chemistry 86, 3251–3257 (1982).
  * Outhwaite, Bhuiyan, and Levine [1980] C. W. Outhwaite, L. B. Bhuiyan, and S. Levine, “Theory of the electric double layer using a modified poisson–boltzmann equation,” Journal of the Chemical Society, Faraday Transactions 2: Molecular and Chemical Physics 76, 1388–1408 (1980).
  * López-García, Horno, and Grosse [2011] J. J. López-García, J. Horno, and C. Grosse, “Poisson–boltzmann description of the electrical double layer including ion size effects,” Langmuir 27, 13970–13974 (2011).
  * Boström, Williams, and Ninham [2001] M. Boström, D. Williams, and B. Ninham, “Specific ion effects: why dlvo theory fails for biology and colloid systems,” Physical review letters 87, 168103 (2001).
  * Joly  _et al._ [2004] L. Joly, C. Ybert, E. Trizac, and L. Bocquet, “Hydrodynamics within the electric double layer on slipping surfaces,” Physical review letters 93, 257805 (2004).
  * Joly  _et al._ [2006] L. Joly, C. Ybert, E. Trizac, and L. Bocquet, “Liquid friction on charged surfaces: From hydrodynamic slippage to electrokinetics,” The Journal of chemical physics 125, 204716 (2006).
  * Netz and Orland [2000] R. R. Netz and H. Orland, “Beyond poisson-boltzmann: Fluctuation effects and correlation functions,” The European Physical Journal E 1, 203–214 (2000).
  * Hedley, Berthoumieux, and Kornyshev [2023] J. G. Hedley, H. Berthoumieux, and A. A. Kornyshev, “The dramatic effect of water structure on hydration forces and the electrical double layer,” The Journal of Physical Chemistry C 127, 8429–8447 (2023).
  * Schlaich  _et al._ [2024] A. Schlaich, J. O. Daldrop, B. Kowalik, M. Kanduč, E. Schneck, and R. R. Netz, “Water structuring induces nonuniversal hydration repulsion between polar surfaces: Quantitative comparison between molecular simulations, theory, and experiments,” Langmuir 40, 7896–7906 (2024).
  * Alexander  _et al._ [1984] S. Alexander, P. Chaikin, P. Grant, G. Morales, P. Pincus, and D. Hone, “Charge renormalization, osmotic pressure, and bulk modulus of colloidal crystals: Theory,” The Journal of chemical physics 80, 5776–5781 (1984).
  * Belloni [1998] L. Belloni, “Ionic condensation and charge renormalization in colloidal suspensions,” Colloids and Surfaces A: Physicochemical and Engineering Aspects 140, 227–243 (1998).
  * Aubouy, Trizac, and Bocquet [2003] M. Aubouy, E. Trizac, and L. Bocquet, “Effective charge versus bare charge: an analytical estimate for colloids in the infinite dilution limit,” Journal of Physics A: Mathematical and General 36, 5835 (2003).
  * Kreer, Horbach, and Chatterji [2006] T. Kreer, J. Horbach, and A. Chatterji, “Nonlinear effects in charge stabilized colloidal suspensions,” Physical Review E 74, 021401 (2006).
  * Trefalt, Behrens, and Borkovec [2016] G. Trefalt, S. H. Behrens, and M. Borkovec, “Charge regulation in the electrical double layer: ion adsorption and surface interactions,” Langmuir 32, 380–400 (2016).
  * Acharya and Lau [2020] P. Acharya and A. Lau, “Charge regulation of a surface immersed in an electrolyte solution,” The European Physical Journal E 43, 1–27 (2020).
  * Schlaich  _et al._ [2023] A. Schlaich, S. Tyagi, S. Kesselheim, M. Sega, and C. Holm, “Renormalized charge and dielectric effects in colloidal interactions: a numerical solution of the nonlinear poisson–boltzmann equation for unknown boundary conditions,” The European Physical Journal E 46, 80 (2023).
  * Brito, Nägele, and Denton [2023] M. E. Brito, G. Nägele, and A. R. Denton, “Effective interactions, structure, and pressure in charge-stabilized colloidal suspensions: Critical assessment of charge renormalization methods,” The Journal of Chemical Physics 159 (2023).
  * Valeau and Cohen [1980] J. Valeau and L. Cohen, “Primitive model electrolyte. i. grand canonical Monte Carlo computation,” The Journal of Chemical Physics 72, 5935–5941 (1980).
  * Smit [1995] B. Smit, “Grand canonical Monte Carlo simulations of chain molecules: adsorption isotherms of alkanes in zeolites,” Molecular Physics 85, 153–172 (1995).
  * Stern [2007] H. A. Stern, “Molecular simulation with variable protonation states at constant pH,” The Journal of chemical physics 126, 04B627 (2007).
  * Barr and Panagiotopoulos [2012] S. Barr and A. Panagiotopoulos, “Grand-canonical Monte Carlo method for donnan equilibria,” Physical Review E 86, 016703 (2012).
  * Michael  _et al._ [2020] E. Michael, S. Polydorides, T. Simonson, and G. Archontis, “Hybrid mc/md for protein design,” The Journal of Chemical Physics 153, 054113 (2020).
  * Mezei [1980] M. Mezei, “A cavity-biased (T, V, $\mu$) Monte Carlo method for the computer simulation of fluids,” Molecular Physics 40, 901–906 (1980).
  * Shelley and Patey [1994] J. Shelley and G. Patey, “A configuration bias Monte Carlo method for ionic solutions,” The Journal of Chemical Physics 100, 8265–8270 (1994).
  * Shelley and Patey [1995] J. Shelley and G. Patey, “A configuration bias Monte Carlo method for water,” The Journal of chemical physics 102, 7656–7663 (1995).
  * Garberoglio [2008] G. Garberoglio, “Boltzmann bias grand canonical Monte Carlo,” The Journal of chemical physics 128, 134109 (2008).
  * Shi and Maginn [2007] W. Shi and E. J. Maginn, “Continuous fractional component Monte Carlo: an adaptive biasing method for open system atomistic simulations,” Journal of chemical theory and computation 3, 1451–1463 (2007).
  * Moucka, Bratko, and Luzar [2015] F. Moucka, D. Bratko, and A. Luzar, “Electrolyte pore/solution partitioning by expanded grand canonical ensemble Monte Carlo simulation,” The Journal of Chemical Physics 142 (2015).
  * Panagiotopoulos [1989] A. Panagiotopoulos, “Exact calculations of fluid-phase equilibria by Monte Carlo simulation in a new statistical ensemble,” International Journal of Thermophysics 10, 447–457 (1989).
  * Soroush Barhaghi  _et al._ [2018] M. Soroush Barhaghi, K. Torabi, Y. Nejahi, L. Schwiebert, and J. J. Potoff, “Molecular exchange Monte Carlo: A generalized method for identity exchanges in grand canonical Monte Carlo simulations,” The Journal of chemical physics 149, 072318 (2018).
  * Fathizadeh and Elber [2018] A. Fathizadeh and R. Elber, “A mixed alchemical and equilibrium dynamics to simulate heterogeneous dense fluids: Illustrations for Lennard-Jones mixtures and phospholipid membranes,” The Journal of Chemical Physics 149, 072325 (2018).
  * Praprotnik, Delle Site, and Kremer [2005] M. Praprotnik, L. Delle Site, and K. Kremer, “Adaptive resolution molecular-dynamics simulation: Changing the degrees of freedom on the fly,” The Journal of chemical physics 123 (2005).
  * Baptista  _et al._ [2021] L. A. Baptista, R. C. Dutta, M. Sevilla, M. Heidari, R. Potestio, K. Kremer, and R. Cortes-Huerto, “Density-functional-theory approach to the hamiltonian adaptive resolution simulation method,” Journal of Physics: Condensed Matter 33, 184003 (2021).
  * Duane  _et al._ [1987] S. Duane, A. D. Kennedy, B. J. Pendleton, and D. Roweth, “Hybrid Monte Carlo,” Physics letters B 195, 216–222 (1987).
  * Mehlig, Heermann, and Forrest [1992] B. Mehlig, D. Heermann, and B. Forrest, “Hybrid Monte Carlo method for condensed-matter systems,” Physical Review B 45, 679 (1992).
  * Boinepalli and Attard [2003] S. Boinepalli and P. Attard, “Grand canonical molecular dynamics,” The Journal of chemical physics 119, 12769–12775 (2003).
  * Chen and Roux [2014] Y. Chen and B. Roux, “Efficient hybrid non-equilibrium molecular dynamics-Monte Carlo simulations with symmetric momentum reversal,” The Journal of Chemical Physics 141, 09B612_1 (2014).
  * Radak  _et al._ [2017] B. K. Radak, C. Chipot, D. Suh, S. Jo, W. Jiang, J. C. Phillips, K. Schulten, and B. Roux, “Constant-pH molecular dynamics simulations for large biomolecular systems,” Journal of chemical Theory and Computation 13, 5933–5944 (2017).
  * Ross  _et al._ [2018] G. A. Ross, A. S. Rustenburg, P. B. Grinaway, J. Fass, and J. D. Chodera, “Biomolecular simulations under realistic macroscopic salt conditions,” The Journal of Physical Chemistry B 122, 5466–5486 (2018).
  * Prokhorenko  _et al._ [2018] S. Prokhorenko, K. Kalke, Y. Nahas, and L. Bellaiche, “Large scale hybrid Monte Carlo simulations for structure and property prediction,” npj Computational Materials 4, 1–7 (2018).
  * Nilmeier  _et al._ [2011] J. P. Nilmeier, G. E. Crooks, D. D. Minh, and J. D. Chodera, “Nonequilibrium candidate Monte Carlo is an efficient tool for equilibrium simulation,” Proceedings of the National Academy of Sciences 108, E1009–E1018 (2011).
  * Belloni [2019] L. Belloni, “Non-equilibrium hybrid insertion/extraction through the 4th dimension in grand-canonical simulation,” The Journal of Chemical Physics 151, 021101 (2019).
  * Kim, Belloni, and Rotenberg [2023] J. Kim, L. Belloni, and B. Rotenberg, “Grand-canonical molecular dynamics simulations powered by a hybrid 4D nonequilibrium MD/MC method: Implementation in LAMMPS and applications to electrolyte solutions,” The Journal of Chemical Physics 159, 144802 (2023).
  * Attard [1995] P. Attard, “Ion condensation in the electric double layer and the corresponding poisson-boltzmann effective surface charge,” The Journal of Physical Chemistry 99, 14174–14181 (1995).
  * Barrat and Hansen [2003] J.-L. Barrat and J.-P. Hansen,  _Basic concepts for simple and complex liquids_ (Cambridge University Press, 2003).
  * Buyukdagli, Achim, and Ala-Nissila [2012] S. Buyukdagli, C. Achim, and T. Ala-Nissila, “Electrostatic correlations in inhomogeneous charged fluids beyond loop expansion,” The Journal of Chemical Physics 137 (2012).
  * Buyukdagli and Blossey [2016] S. Buyukdagli and R. Blossey, “Beyond poisson–boltzmann: fluctuations and fluid structure in a self-consistent theory,” Journal of Physics: Condensed Matter 28, 343001 (2016).
  * Guo, Haji-Akbari, and Palmer [2018] J. Guo, A. Haji-Akbari, and J. C. Palmer, “Hybrid Monte Carlo with lammps,” Journal of Theoretical and Computational Chemistry 17, 1840002 (2018).
  * Crooks [1999] G. E. Crooks, “Entropy production fluctuation theorem and the nonequilibrium work relation for free energy differences,” Physical Review E 60, 2721 (1999).
  * Bennett [1976] C. H. Bennett, “Efficient estimation of free energy differences from Monte Carlo data,” Journal of Computational Physics 22, 245–268 (1976).
  * Shirts  _et al._ [2003] M. R. Shirts, E. Bair, G. Hooker, and V. S. Pande, “Equilibrium free energies from nonequilibrium measurements using maximum-likelihood methods,” Physical review letters 91, 140601 (2003).
  * Plimpton [1995] S. Plimpton, “Fast parallel algorithms for short-range molecular dynamics,” Journal of computational physics 117, 1–19 (1995).
  * Bonthuis, Gekle, and Netz [2011] D. J. Bonthuis, S. Gekle, and R. R. Netz, “Dielectric profile of interfacial water and its effect on double-layer capacitance,” Phys. Rev. Lett. 107, 166102 (2011).
  * Schlaich, Knapp, and Netz [2016] A. Schlaich, E. W. Knapp, and R. R. Netz, “Water dielectric effects in planar confinement,” Phys. Rev. Lett. 117, 048001 (2016).
  * Gonella  _et al._ [2021] G. Gonella, E. H. Backus, Y. Nagata, D. J. Bonthuis, P. Loche, A. Schlaich, R. R. Netz, A. Kühnle, I. T. McCrum, M. T. Koper, _et al._ , “Water at charged interfaces,” Nature Reviews Chemistry 5, 466–485 (2021).
  * Nickel  _et al._ [2024] O. Nickel, L. J. Ahrens-Iwers, R. H. Meißner, and M. Janssen, “Water, not salt, causes most of the seebeck effect of nonisothermal aqueous electrolytes,” Physical Review Letters 132, 186201 (2024).