# Alchemical Osmostat for Monte Carlo Simulation: Sampling Aqueous Electrolyte Solution in Open Systems

Ambroise De Izarra, François-Xavier Coudert, Alain Fuchs, Anne Boutin

► To cite this version:

Ambroise De Izarra, François-Xavier Coudert, Alain Fuchs, Anne Boutin. Alchemical Osmostat for Monte Carlo Simulation: Sampling Aqueous Electrolyte Solution in Open Systems. Journal of Physical Chemistry B, 2023, 127 (3), pp.766-776. 10.1021/acs.jpcb.2c07902. hal-03938290

HAL Id: hal-03938290

https://hal.science/hal-03938290

Submitted on 13 Jan 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Alchemical Osmostat for MC Simulation:
Sampling Aqueous Electrolyte Solution in Open
Systems

Ambroise de Izarra,†,‡ François-Xavier Coudert,*,‡ Alain H. Fuchs,‡ and Anne
Boutin*,†

†PASTEUR, Département de Chimie, École Normale Supérieure, PSL University,
Sorbonne Université, CNRS, 75005, Paris, France

‡Chimie ParisTech, PSL University, CNRS, Institut de Recherche de Chimie Paris, Paris,
France

E-mail: fx.coudert@chimieparistech.psl.eu; anne.boutin@ens.psl.eu
Abstract

Molecular simulations involving electrolytes are usually performed at a fixed amount of salt ions in the simulation box, reproducing macroscopic concentration. Although this statement is valid in the bulk, the concentration of an electrolyte confined in a nanoporous materials such as MOFs or zeolites is greatly affected and remains a *priori* unknown. The nanoporous material in equilibrium with the bulk electrolyte exchange water and ions at a given chemical potential $\Delta\mu$ in the semigrand canonical ensemble that must be calibrated in order to determine the concentration in the nanoporous material. In this work, we propose an algorithm based on non-equilibrium candidate Monte Carlo (NCMC) moves to ultimately perform MC simulations in contact with a saline reservoir. First, we adapt the Widom insertion technique to calibrate the chemical potential by alchemically transmuting water molecules into ions by using NCMC moves. The chemical potential defines a Monte Carlo osmostat in the semigrand constant volume and temperature ensemble ($\Delta\mu, N, V, T$) to be added in a Monte Carlo simulation where the number of ions fluctuates. In order to validate the method, we adapted the NCMC move to determine the free energy of water solvation and subsequently explore thermodynamics of electrolyte solvation at infinite dilution in water. Finally, we implemented the osmostat in MC simulations initialized with bulk water that are driven towards electrolytes of similar concen concentration as the saline reservoir. Our results demonstrate that alchemical osmostat MC simulation is a promising tool for use to sample electrolyte insertion in nanoporous materials.
Introduction

In computational chemistry, atomistic Monte Carlo (MC) simulations are routinely used
to study the structure and equilibrium properties of condensed phases, solid or liquid, and
their interfaces. Water is an important fluid in nature, and aqueous electrolytes at various
concentrations are a key component of many chemical and biological systems, as well as used
in applications such as batteries, desalination, and pollution remediation. Consequently, the
insertion of ionic species in a dense fluid such as water is required for molecular simulation in
many situations involving open systems, such as the osmotic equilibrium through membranes
or the uptake of electrolytes in nanoporous materials. The Grand Canonical Monte Carlo
(GCMC) simulation method is well-suited for simulating the molecular insertion or exchange
of neutral species in fluid phases, and is used ubiquituously in the simulation of gas and liquid
adsorption in porous materials.¹ However, its use for ionic species is particularly inefficient: in
addition to the fact that finding suitable cavities large enough in the fluid to accomodate ion
insertion would appear infrequently, as for any large molecule, the strong organization of the
solvent induced by the electric charge gives the insertion move a low probability. Secondly,
electrical neutrality of the system imposes the insertion of pairs of charge-balancing ions (or
more), which further decreases the probability of a successful Monte Carl insertion move.

Yet, the convergence and accuracy of a MC-based molecular simulation depends greatly
on the ability to carry out insertion or deletion of ions, in order to describe phenomenan
like the adsorption of electrolytes at interfaces or in porous materials, without an explicit
description of the electrolyte reservoir (which is computational expensive, and induces finite-
size effects). The design of adapted MC moves that have higher probabilities of particle
insertion or deletion in dense fluid is therefore required,² for which a chemical potential is
imposed and dictate composition evolution in response to this chemical potential. A large
number of biasing techniques have been developed to this end, and without attempting an
exhaustive review of this vast literature, we summarize below some of the most relevant
works to provide the reader with context for our own work.
Bias insertion techniques$^{3-6}$ are MC moves that have been designed to search the favorable energetic location,$^{5,7}$ suitable cavities$^{8}$ or distance-biased insertion holes$^{3}$ in the system to enhance sampling of ion insertion in water. An important example is the method of thermodynamic integration originally proposed by Kirkwood$^{9}$ which is suited for calculations of free energies between two given states. In past work, it has been applied to the calculation of the chemical potential of ions in aqueous solutions.$^{10-12}$ Its basic principle is the definition of a thermodynamic path where an ion is gradually added to the electrolyte solution by coupling its interaction with a scaling parameter $\lambda$, varying in the interval $[0,1]$, during a MC or molecular dynamics (MD) simulation. Other advanced techniques such as multicanonical methods,$^{13-15}$ metadynamics for chemical potential calculation$^{6}$ or advanced configurational bias moves$^{16,17}$ have been developed.

Another important class of methods are the fractional component Monte Carlo moves$^{18,19}$ where gradual insertion and deletion of particles are performed through use of a coupling parameter $\lambda$ with an adaptative bias potential — bearing some similarity in essence to the thermodynamic integration methods. These methods have been used in single and multi-component systems in dense phases, as a way to improve phase equilibria simulations and to allow the direct calculation of chemical potentials and partial molar properties. In fractional component MC simulations, the number of particles is a continuous variable through $\lambda$ in a way that particle insertion or deletion is done with a so-called *fractional particle*. More precisely, a set of Metropolis-like acceptance rules can be derived for $\lambda$: if a change of the coupling parameter occurs so that $\lambda > 1$, the current fractional component particle is converted into a real particle and a new fractional particle is inserted with coupling parameter of $\lambda - 1$. Likewise, if the coupling parameter becomes $\lambda < 0$, the fractional particle is removed from the system and a remaining particle is designed as the new fractional particle with a coupling parameter of $\lambda + 1$. To further improve the convergence of these simulations, an appealing strategy is to improve the acceptance probability of particle insertions or deletions by using some local relaxation process, in order to minimize energy penalty and remove non-
physical forces by collective motion in the surrounding molecules. This relaxation is typically
performed through the introduction of a number of MD steps. Therein, Çağin and Pettitt
developed an hybrid method²⁰ where an MC move is used to change the composition of the
system, associated with deterministic MD trajectories for relaxation. Another example of
scheme has been proposed by Boinepalli and Attard²¹ where a MD step in the microcanoni-
cal ensemble is performed at a given state associated to a certain value of coupling parameter
λ of the current fractional molecule. For further details, we refer the interested reader to the
recent and comprehensive review by Rabhari and co-workers on fractional component Monte
carlo methods.²²

Recently, Nilmeier and co-workers²³ introduced a new type of hybrid MC move denoted
*nonequilibrium candidate Monte Carlo (NCMC)* based on non equilibrium molecular dynam-
ics. In this method, a candidate move is proposed through a finite-time process where the
system is driven out of equilibrium, and accepted with a probability that preserves micro-
scopic reversibility. In this process, the acceptance rule follows a Metropolis-type probability,
where the total work performed during the out-of-equilibrium process is taken into account
(instead of the difference of free energy). The generation of out-of-equilibrium states induces
an additional computational cost, compensated by the reduced structural correlation that
can occur in dense systems such as electrolytes and the versatile nature of the algorithm.
Along this line, Ross and co-workers²⁴ adapted the NCMC scheme through implementation
of a Monte Carlo osmostat added to MC simulations that can sample from a semigrandcanon-
ical isobaric isothermal ensemble in which the number of NaCl salt pairs (that we will called
later an electrolyte unit for generalization) varies to study ion distribution in water around
biomolecules. The NCMC moves enable high acceptance rates for which water molecules
are alchemically transformed into a pair of salt ion (i.e. an electrolyte unit) or vice-versa.
MC simulation in the semigrand canonical ensemble requires the calibration of the chemical
potential Δμ used as a thermodynamic parameter governing the osmostat in equilibrium
with an infinitely sized reservoir of electrolyte at a specified salt concentration. We invite
the readers to refer to the original publication of Ross and Chodera for more details. In their work, the self-adjusted mixture sampling<sup>25</sup> sample uniformly salt concentration during a single simulation, follow by application of the Bennett acceptance ratio to extract relevant free energy difference from NCMC move. The chemical potential can be thus extracted from relation between macroscopic electrolyte concentration and free energies.

In this work, we proposed a variant of the osmostat implementation that we implemented in RASPA,<sup>26</sup> a software dedicated to simulation of adsorption in porous materials. Our ultimate goal would be to employ in the future our osmostat version for the study of various electrolyte insertion in porous materials in equilibrium with a saline reservoir. The paper is organized as follows. First, we proposed another version of chemical potential calibration in the semigrand canonical ensemble ($\Delta\mu, N, V, T$) by using the Widom insertion technique with NCMC move to transform alchemically water into an electrolyte unit. Then, in order to assess the accuracy of our chemical potential calibration, NCMC move is used to study thermodynamics of electrolyte solvation at infinite dilution by calculating the free energy of water solvation with our modified Widom method. Finally, the chemical potential is used as the thermodynamic parameter governing the osmostat to gradually perturb bulk water that transition into an electrolyte that reaches concentration dictated by the previously calibrated chemical potential.

## General simulation methodology

We ran molecular dynamics (MD) and Monte Carlo (MC) simulations using the RASPA software.<sup>26</sup> We used a combination of Lennard-Jones and electrostatic potentials. Water molecules are represented with the rigid TIP4P model<sup>27</sup> and the electrolytes are modeled as charged particles that differs by their Lennard-Jones parameters and charges. We considered the following monovalent electrolytes which are: LiCl, NaCl, KCl, RbCl, CsCl, KF, KBr and KI. We also explored the following divalent electrolytes to probe influence of charge ion:
$MgCl_2$, $CaCl_2$, $SrCl_2$ and $BaCl_2$. The details of the force field are given in the SI.²⁸,²⁹

We used Lorentz-Berthelot mixing rules for cross-terms in Lennard-Jones potential, and Ewald summation was used to account efficiently for electrostatic interactions. We used a cutoff of 8.0 Å for both the Lennard-Jones potential and the separation between real space and Fourier space in the Ewald summation. MD simulations were performed with a timestep of 1 fs and all simulations used a Nosé-Hoover thermostat, setting the temperature at 300 K. In order to explore the phase space by using MC simulations, translations were performed both on ions and water while rotations were applied only to water molecules, in a ratio of 1:1:1 respectively. NCMC moves are performed less frequently than these “standard” MC moves so that the relative probability to perform NCMC moves was set to $10^{-3}$. The reason is that ion insertion/deletion achieved during a NCMC move would appear less frequently than simple MC moves to rotate or translate a chosen particle. Also, short MD steps performed during alchemical transformation of water molecules into ions (or vice versa) accounts also for collective motion of particles to accommodate insertion of water (or ions). Consequently, the execution time of a NCMC move is much higher than the execution time of a “standard” MC move, as discussed below. In our validation of the method, we considered systems containing at most 150 molecules to keep the computational cost reasonable, due to the sequential nature of the RASPA code.

## Calibration of the chemical potential in the semigrand canonical ensemble

Let us consider a system made of $N$ molecules with $N_{H_2O}$ water molecules, $N_{C^{m+}}$ salt cations and $N_{A^{n-}}$ salt anions where $m$ and $n$ design respectively the absolute charge of the cation and of the anion. The number of electrolyte units (i.e. the number of salt pair of NaCl where $m$ and $n$ are equal to one for instance) is determined as follows:
$$N_{C_n A_m} = \frac{N_{C^{m+}} + N_{A^{n-}}}{m + n} \qquad (1)$$

The partition function for such a system having $N_{C_n A_m}$ electrolyte units reads:

$$Z(N_{C_n A_m}, N, V, T) = \frac{N!}{N_{C^{m+}}! N_{A^{n-}}! N_{H_2O}!} \int \mathrm{d}\mathbf{r} e^{-\beta U[\mathbf{r}, N_{C_n A_m}]} \qquad (2)$$

The chemical potential in the semigrand canonical ensemble $\Delta\mu_{(n+m)H_2O \to C_n A_m}$ (denoted $\Delta\mu_{C_n A_m}$ for clarity) is the difference of chemical potential between withdrawing $(n+m)$ water molecules and inserting an electrolyte unit $C_n A_m$.²⁴ The number of water molecules and electrolyte units become respectively $N'_{H_2O} = N_{H_2O} - (n+m)$ and $N'_{C_n A_m} = N_{C_n A_m} + 1$ which implies that the number of anions and cations increases respectively to $N'_{A^{n-}} = N_{A^{n-}} + m$ and $N'_{C^{m+}} = N_{C^{m+}} + n$. Finally, the chemical potential $\Delta\mu_{C_n A_m}$ is expressed as follows:

$$\begin{aligned} \Delta\mu_{C_n A_m} &= -\frac{1}{\beta} \ln \left( \frac{Z(N'_{C_n A_m}, N, V, T)}{Z(N_{C_n A_m}, N, V, T)} \right) \\ \Delta\mu_{C_n A_m} &= -\frac{1}{\beta} \ln \left( \frac{\prod_{i=0}^{m+n-1} (N_{H_2O} - i)}{\prod_{j=1}^{n} (N_{C^{m+}} + j) \prod_{k=1}^{m} (N_{A^{n-}} + k)} \right) \left\langle e^{-\beta W_{(n+m)H_2O \to C_n A_m}} \right\rangle_{N_{C_n A_m}} \end{aligned} \qquad (3)$$

where $\langle e^{-\beta W_{(n+m)H_2O \to C_n A_m}} \rangle_{N_{C_n A_m}}$ is the Boltzmann factor of the work performed for transforming alchemically water molecules into an electrolyte unit (denoted alchemical work from now on), averaging in the canonical ensemble of the system containing $N_{C_n A_m}$ electrolyte units. This adapted version of the Widom technique is slightly different from the original one applied to GCMC simulations as it requires only the difference of potential energy before and after insertion of a particle.³⁰ In our case, we employ NCMC moves²³,²⁴ to compute the alchemical work $W_{(n+m)H_2O \to C_n A_m}$ (denoted $W_{C_n A_m}$ for clarity) performed during the alchemical transformation: $(n+m)$ water molecules were chosen randomly in the box and alchemically transformed into an electrolyte unit $C_n A_m$ over an alchemical path divided in
T = 200 segments, each one consisting of a NCMC step where the system is perturbated and then relaxed. For a given segment t, the alchemical path is decribed by two parameters: the state of the system r_t and the values of the interpolated non-bonded parameters λ_t of the molecules undergoing alchemical transformation. In the same manner as in the work of Ross and co-workers,²⁴ the alchemical path is a linear interpolation of the nonbonded parameters of water and salt ions, applied to randomly chosen (n+m) water molecules to be transformed into ions:

$$ \lambda_t = (1 - f_t)\lambda_{\text{water}} + f_t\lambda_{\text{ion}} \quad \text{and} \quad f_t = \frac{t}{T} \qquad (4) $$

A segment t within the alchemical path consists in a perturbation of the system by updating non-bonded parameters (r_t, λ_t)→(r_t, λ_{t+1}), followed its relaxation from (r_t, λ_{t+1})→(r_{t+1}, λ_{t+1}) by using 200 steps of MD in the NVE ensemble. A NCMC move requires thus 40 000 MD steps which represents a duration of 40 ps. At each segment, the change of potential energy after perturbation of the system is recorded and summed, yielding to the alchemical work W_{C_nA_m} performed on the system to alchemically transform (n + m) water molecules into an electrolyte unit C_nA_m along the alchemical path:

$$ W_{C_nA_m} = \sum_{t=1}^{T} (U[r_t, \lambda_{t+1}] - U[r_t, \lambda_t]) \qquad (5) $$

Hence, a single NCMC move yields to one value of W_{C_nA_m} and we performed 3000 NCMC moves to estimate W_{C_nA_m} and obtain the chemical potential Δμ_{C_nA_m} during a MC simulation by using our modified Widom technique. We present typical execution time of a NCMC move compared to standard moves implemented in RASPA in the Supplementary Information (Table S8). A NCMC move is computationally expensive as its execution time is seven order of magnitude higher than those for translation and rotational moves, at the
benefit of a proper equilibration of the system during the alchemical transformation. The calibration process of the chemical potential is performed for different electrolyte over a wide range of concentration (Table 1).

<table><caption>**Table 1:** Composition of systems of 150 particles after transformation of water into electrolyte units. For each concentration, $\Delta\mu_{C_nA_m}$ is calculated according to equation 3. The upper table presents concentration of monovalent electrolytes ($n = 1$ and $m = 1$) while the lower table presents concentration of divalent cation electrolytes ($n = 1$ and $m = 2$).</caption><thead><tr><th>Number of electrolyte units (N<sub>CA</sub>)</th><th>N<sub>H<sub>2</sub>O</sub></th><th>ratio N<sub>H<sub>2</sub>O</sub>/N<sub>CA</sub></th><th>Concentration (M)</th></tr></thead><tbody><tr><td>1</td><td>148</td><td>148</td><td>0.37</td></tr><tr><td>6</td><td>138</td><td>23</td><td>2.41</td></tr><tr><td>11</td><td>128</td><td>11.63</td><td>4.77</td></tr><tr><td>20</td><td>110</td><td>5.5</td><td>10.09</td></tr><tr><td>30</td><td>90</td><td>3</td><td>18.50</td></tr><tr><td>40</td><td>70</td><td>1.75</td><td>31.72</td></tr></tbody><thead><tr><th>Number of electrolyte unit (N<sub>CA<sub>2</sub></sub>)</th><th>N<sub>H<sub>2</sub>O</sub></th><th>ratio N<sub>H<sub>2</sub>O</sub>/N<sub>CA<sub>2</sub></sub></th><th>Concentration (M)</th></tr></thead><tbody><tr><td>1</td><td>147</td><td>147</td><td>0.38</td></tr><tr><td>6</td><td>132</td><td>22</td><td>2.52</td></tr><tr><td>10</td><td>120</td><td>12</td><td>4.62</td></tr><tr><td>18</td><td>96</td><td>5.33</td><td>10.40</td></tr><tr><td>25</td><td>75</td><td>3</td><td>18.50</td></tr><tr><td>32</td><td>54</td><td>1.68</td><td>32.89</td></tr></tbody></table>

This concentration are chosen according to experimental studies where adsorption of electrolyte in zeosils are studying at concentration up to 20 mol/L.³¹⁻³³ Each system contains 150 particles and the volume of the cubic simulation of side $L = 16.646$ Å is fixed to match the density of bulk water at atmospheric pressure in absence of ions. The number of electrolyte units $N_{C_nA_m}$ and the number of water $N_{H_2O}$ shown in Table 1 are the targeted number of species after removal of ($n + m$) water and insertion of an electrolyte unit $C_nA_m$.

The Figure 1 shows the evolution of the alchemical work and chemical potential for various electrolytes as a function of the electrolyte concentration. We could not explore KI electrolyte for concentration greater than 10 M (Figure 1-B, dark blue) because the largest anions $I^-$ ($\sigma_{I^-} = 5.167$ Å) could not be inserted in the box with other species. For all electrolyte, the alchemical work and the free energy are the highest (i.e., more negative)
<figure><img src="image_2.png" /><figcaption>Figure 1: A-B) For monovalent electrolytes CA: Evolution of the alchemical work performed during alchemical transformation $W_{CA}$ (left) and associate chemical potential $\Delta\mu_{CA}$ (right). For A), the cation is varied and the anion is fixed to $Cl^{-}$. For B), the cation is fixed at $K^{+}$ and the anion is varied. C) For divalent electrolytes $CA_{2}$: Evolution of the alchemical work $W_{CA_{2}}$ (left) and associate chemical potential $\Delta\mu_{CA_{2}}$ (right).</figcaption></figure>
after insertion of a single electrolyte unit ($c = 0.37$ M), which is consistent with similar calculations performed on NaCl electrolyte.²⁴ Upon increase of electrolyte units initially present in solution, the alchemical work required to transform water into ions decreases in absolute value as well as the chemical potential. Previous simulations of potassium halide aqueous solutions³⁴ showed that the hydration number (i.e. the number of water molecules in the first hydration shell) of cations and anions drops while the number of ion pair increases slightly. This trend is consistent with the radial distribution functions (Figures S4 to S15 in the Supporting Information) and the number of water molecules in the first solvation shell around ions (Tables S9 to S20). The hydration shells around ions are disturbed when concentration of electrolyte increases, which lower solvation energy of salt ions. Another common trend highlighted is that the alchemical work and chemical potential increases in absolute value when the radius of the ion diminishes, because smaller ions have a higher solvation free energy.³⁵ The absolute difference of chemical potential after insertion of an electrolyte unit is about ~ 215 kJ/mol between LiCl and CsCl, and about ~ 494 kJ/mol between MgCl₂ and BaCl₂, which is more than twice the value for monovalent electrolytes.

We recall that Li⁺ and Mg²⁺ have similar sizes of respectively $σ_{Li⁺} = 1.505$ Å and $σ_{Mg²⁺} = 1.630$ Å, as well as Cs⁺ and Ba²⁺ with a respective size of $σ_{Cs⁺} = 3.883$ Å and $σ_{Ba²⁺} = 3.820$ Å. A higher charge of the cation leads thus to a stronger energy of solvation. This observation is in qualitative agreement with previous experiments³⁶⁻³⁸ and simulations³⁵ showing that only the first solvation shell is perturbed by monovalent cation, which is not the case for divalent cations where perturbations occur in the second solvation shell,³⁵ that could account for the higher difference of chemical potential for divalent cation electrolytes.
<figure><img src="image_3.png" /><figcaption>Figure 2: Evolution of the alchemical work as a function of the NCMC moves for alkali-chlorides electrolyte solution at two different concentration. Left) $N_{CA} = 1$ is equivalent to 0.37 M. Right) $N_{CA} = 30$ is equivalent to 18.45 M.</figcaption></figure>

The evolution of the alchemical work for all alkali-chloride electrolyte solutions as a function of the NCMC moves is given in Figure 2 and in the SI for every electrolytes (Figures S1 to S3). The numerical data and an estimation of the error for both the alchemical work and the chemical potential are also provided in SI (Tables S4 to S7). The evolution of the alchemical work as a function of NCMC moves remains constant but displays more fluctuations when the number of electrolyte units increases in the system, leading to higher error calculated for these cases. Overall, errors of the alchemical work calculated for monovalent electrolytes are below of ~1 kJ/mol, except for divalent cations where error reach at most ~4 kJ/mol. Along this line, errors associated with the calculation of the chemical potential are reasonable and reach at most ~5.5 kJ/mol (~1.3 kcal/mol) at high electrolyte concentration. However, the error associated with the calibration process using our method is higher than those of Ross and Chodera for which error of $\Delta\mu_{NaCl}$ when $N_{NaCl}$ is four times lower than ours. However, our calibration process is rather straightforward to implement and errors are small enough so that it represents another way for calculating chemical potential in solution with low electrolyte concentration.
# Thermodynamics of electrolyte solvation at infinite dilution

We consider one of the previous system for which we calculated the chemical potential $\Delta\mu_{C_nA_m}$ where only one electrolyte unit ($N_{C_nA_m} = 1$) is inserted in the box (first line of upper and lower Table 1). Electrolyte solvation at infinite dilution is thus studied using a moderate box size. These small systems contain enough water molecules to capture the thermodynamics of electrolyte solvation because experiments$^{36-38}$ and MD simulations$^{35}$ demonstrated that ions disturb hydrogen bonds between water molecules only in the first and second solvation shell. At the end of each NCMC move, we capture the distance separating the ion pair in order to construct the radial distribution function and deduce the potential of mean force of ion pair at infinite dilution, as shown in Figure S16 in the Supporting Information. This approach is more suitable than analyzing $(N, V, T)$ simulations where only one ion pair is present because ion pair formation (i.e., short distance between ions) can be more easily sampled.

We first calculated the free energy of electrolyte solvation in water at infinite dilution. This quantity is the variation of the Helmholtz free energy $\Delta F_{\text{vacuum} \to C_nA_m}$ (abbreviated later by $\Delta F_{C_nA_m}$ for clarity) for bringing an electrolyte unit $C_nA_m$ from the vacuum to the TIP4P liquid water. We divided the solvation process in two steps. First, water molecules are removed from vacuum and solvated in bulk water. Then, water molecules are alchemically transformed into ions as described in the previous section, so that the free energy of solvation reads:

$$ \Delta F_{C_nA_m} = \Delta\mu_{C_nA_m} + \Delta\mu_{(n+m)H_2O} - (n+m)k_B T \quad (6) $$

where $\Delta\mu_{(n+m)H_2O}$ is the work required to solvate $(n+m)$ TIP4P water molecules (called
subsequently water solvation work) while the term $(n+m)k_B T$ is the energy cost for removing $(n+m)$ water molecules from vacuum. We previously computed $\Delta\mu_{C_n A_m}$ where a single electrolyte unit $C_n A_m$ is present in the system after water transmutation. However, the solvation free energy of $(n+m)$ water molecules $\Delta\mu_{(n+m)H_2O}$ needs to be calculated. Therein, we prepared a box having the same side $L$ as previously, containing $150-(n+m)$ water molecules and used our modified Widom insertion technique using NCMC moves. We randomly selected $(n+m)$ positions in the simulation box, satisfying the condition that the chosen positions are located at least from a distance of $\sigma_{Ow}/2$ from the center of water oxygens in order to avoid initial interpenetration of beads resulting in divergent Lennard-Jones potential and thus blow up of the system during relaxation process by $(N, V, E)$ molecular dynamics. The work required to solvate $(n+m)$ TIP4P water molecules is the excess Widom potential defined as:

$$ \Delta\mu_{(n+m)H_2O} = -\frac{1}{\beta} \ln \left( \langle e^{-\beta W_{(n+m)H_2O}} \rangle_{(150-(n+m))H_2O} \right) \quad (7) $$

where $\langle e^{-\beta W_{(n+m)H_2O}} \rangle_{(150-(n+m))H_2O}$ is the Boltzmann factor of the water solvation work averaging in the canonical ensemble of the system containing $150-(n+m)$ water molecules. However, the NCMC move is slightly modified as follows: the calculation of the solvation work $W_{(n+m)H_2O}$ was performed in two stages,³⁹ growth and charging. The molecules are grown according to a path of $T = 200$ segments where each segment $t$ consists in a perturbation of the system by updating the Lennard-Jones parameters followed by relaxation through 200 NVE MD steps:

$$ \epsilon_t^{ij} = f_t \epsilon_{\text{final}}^{ij} \qquad \sigma_t^{ij} = f_t \sigma_{\text{final}}^{ij} \qquad f_t = \frac{t}{T} \qquad (8) $$
Secondly, another path of $T = 200$ segments is repeated where the electrostatic interaction is turned on, in the same way as for the growth of the molecules. The water solvation work is summed during the process of growth and charging in the same way as equation 5.

<figure><img src="image_4.png" /><figcaption>Figure 3: Evolution of alchemical work $W_{2H_2O}$ (in black) for solvating 2 water molecules in bulk water. The contribution of vdW and electrostatic interactions are displayed in blue and red, respectively.</figcaption></figure>

For simplicity, we choose to determine the free energy solvation of $n + m = 2$ water molecules and performed 3000 NCMC moves to calculate $W_{2H_2O}$ for subsequently deducing the free energy of water solvation in bulk $\Delta\mu_{2H_2O}$. The average value of the water solvation work is found to be $W_{2H_2O} = -48.767 \pm 0.232$ kJ/mol, resulting in a solvation free energy for 2 water molecules of $\Delta\mu_{2H_2O} = -53.010 \pm 0.491$ kJ/mol. The solvation free energy of a single water molecule obtained from NCMC move is thus $\Delta\mu_{2H_2O}/2 = -26.505$ kJ/mol (or -6.334 kcal/mol), which is consistent with previous studies$^{24,40}$ and confirms that the NCMC method is versatile as it can be used to perform both alchemical transformation and particle insertion.

The evolution of $W_{2H_2O}$ and its decomposition in the work performed during growth and charging as a function of NCMC moves is presented in Figure 3. It is clear that the variation of the non-bonded parameters between two alchemical segment is small enough to prevent blow-up of the system during the relaxation process by NVE MD as shown by the fluctuating but constant value of the work performed during growth ($19.542 \pm 0.184$
kJ/mol) and charging (−68.309 ± 0.112 kJ/mol). Finally, one can deduce the free energy of solvation for electrolytes at infinite dilution, summarized in left column of Table 2. Overall, our calculations are in qualitative agreement with experimental results of Schmidt *et al*<sup>41</sup> in that the free energy are all negatives, and become less negative for larger ions.<sup>35</sup> Also, Ross and Chodera calculated solvation free energy of NaCl at infinite dilution in TIP4P water and found −713.79 kJ/mol<sup>24</sup> which are within 5% of our estimates (−733.4 KJ/mol).

<table><caption>**Table 2:** Solvation free energy of electrolytes $\Delta F_{C_nA_m}$, internal energy $\Delta U_{C_nA_m}$ and entropy $\Delta S_{C_nA_m}$ versus experiments at infinite dilution.</caption><thead><tr><th rowspan="2">C<sub>n</sub>A<sub>m</sub></th><th colspan="2">ΔF<sub>C<sub>n</sub>A<sub>m</sub></sub> (kJ/mol)</th><th colspan="2">ΔU<sub>C<sub>n</sub>A<sub>m</sub></sub> (kJ/mol)</th><th colspan="2">ΔS<sub>C<sub>n</sub>A<sub>m</sub></sub> (J/mol/K)</th></tr><tr><th>This work</th><th>Ref<sup>41</sup></th><th>This work</th><th>Ref<sup>41</sup></th><th>This work</th><th>Ref<sup>41</sup></th></tr></thead><tbody><tr><td>LiCl</td><td>-853 ± 3</td><td>-849</td><td>-886 ± 2</td><td>-898</td><td>-110 ± 16</td><td>-164.1</td></tr><tr><td>NaCl</td><td>-738 ± 2</td><td>-744</td><td>-768 ± 3</td><td>-783</td><td>-100 ± 14</td><td>-132.8</td></tr><tr><td>KCl</td><td>-671 ± 2</td><td>-671</td><td>-695 ± 2</td><td>-700</td><td>-80 ± 10</td><td>-96.7</td></tr><tr><td>RbCl</td><td>-658 ± 2</td><td>-649</td><td>-680 ± 2</td><td>-675</td><td>-73 ± 10</td><td>-85.9</td></tr><tr><td>CsCl</td><td>-635 ± 2</td><td>-626</td><td>-656 ± 2</td><td>-650</td><td>-70 ± 14</td><td>-81.0</td></tr><tr><td>KF</td><td>-838 ± 2</td><td>-799</td><td>-881 ± 2</td><td>-847</td><td>-143 ± 13</td><td>-159.4</td></tr><tr><td>KBr</td><td>-658 ± 1</td><td>-644</td><td>-678 ± 3</td><td>-669</td><td>-67 ± 12</td><td>-80.9</td></tr><tr><td>KI</td><td>-603 ± 2</td><td>-609</td><td>-622 ± 2</td><td>-629</td><td>-63 ± 12</td><td>-62.8</td></tr><tr><td>BaCl<sub>2</sub></td><td>-1999 ± 4</td><td>-1984</td><td>-2061 ± 3</td><td>-2065</td><td>-206 ± 20</td><td>-273.9</td></tr><tr><td>SrCl<sub>2</sub></td><td>-2099 ± 3</td><td>-2112</td><td>-2175 ± 5</td><td>-2203</td><td>-253 ± 23</td><td>-306.3</td></tr><tr><td>CaCl<sub>2</sub></td><td>-2217 ± 4</td><td>-2239</td><td>-2299 ± 3</td><td>-2335</td><td>-273 ± 23</td><td>-324.4</td></tr><tr><td>MgCl<sub>2</sub></td><td>-2495 ± 4</td><td>-2566</td><td>-2582 ± 3</td><td>-2682</td><td>-290 ± 24</td><td>-391.1</td></tr></tbody></table>

In addition, we estimated separately the variation of internal energy $\Delta U_{C_nA_m}$ and deduced the variation of entropy $\Delta S_{C_nA_m}$ from the variation of free solvation energy $\Delta F_{C_nA_m}$. The variation of free energy is the sum of the variation of internal energy and entropy

$$ \Delta F_{C_nA_m} = \Delta U_{C_nA_m} - T\Delta S_{C_nA_m} \qquad (9) $$

By identifying equation 9 with equation 6, the expression for the variation of internal energy writes<sup>35</sup>
$$ \Delta U_{C_n A_m} = U_{C_n A_m + (150 - (n+m)\text{H}_2\text{O})} - U_{(150 - (n+m)\text{H}_2\text{O})} - \frac{1}{2}(n+m)k_B T \quad (10) $$

We calculated separately the internal energies of the system containing only 150 – (n+m) water molecules and with an electrolyte unit C_nA_m in addition to 150 – (n + m) water molecules by performing MD simulations in the NVT ensemble for 20 millions steps to gather statistics after equilibration during 1 millions steps. By taking the difference between the variation of Helmholtz free energy and variation of internal energy, we obtain an estimation of the variation of entropy ΔS_{C_nA_m} for electrolyte solvation at infinite dilution. The results for the variation of internal energy and of entropy are displayed respectively in the middle and left columns of Table 2. Again, the values for ΔU_{C_nA_m} and ΔS_{C_nA_m} are all negative and are more negative for smaller ions, except for CsCl that displays a more negative entropy than RbCl, likely due to accumulated uncertainty from ΔU_{C_nA_m} and ΔF_{C_nA_m}. Not only ΔS_{C_nA_m} are underestimated for all electrolytes (less negative) compared to experimental data, but also values of entropy for divalent cation electrolytes show a larger difference with experiments than monovalent electrolytes which is likely due to the limitation of the empirical force field used in this work. It is well known that empirical force fields that model ions with their charge and a Lennard-Jones center do not describe properly divalent cations in solution. For instance, Mg²⁺ and Zn²⁺ have almost the same radius (0.72 vs 0.74 Å), but it has been shown that their hydration number differs greatly.⁴² This is because divalent cations present strong polarization and charge transfer effects that can not be modeled with simple empirical force fields. Better force fields has been developed to model properly divalent cations in solution. Li and Merz proposed to add another term scaling as r⁻⁴ to the Lennard-Jones potential⁴³ and obtained a better estimation of the hydration free energy, ion oxygen distance as well as coordination number in system involving multivalent cations than using empirical force field with regular 6-12 Lennard-Jones potential. Polarizable force field were successfully applied
to divalent cations for which the charge and the van der Waals centers are distributed to multiple site represented by dummy atoms.⁴⁴ Also, the Electronic Continuum Correction is another approach that takes into account of the polarization of the divalent cation in a mean field approach: the charge of the cation are rescaled by the inverse square root of the electronic part of solvent dielectric constant.⁴⁵ In this work, we have made the choice of using the refined empirical force field proposed by Mamatkulov and coworkers,²⁹ not only because it reproduces empirically experimental energy solvation at infinite dilution for divalent cation electrolytes, but also for the transferability of the force field for monovalent salt, water and later, porous materials.

## MC simulation in the semigrand canonical ensemble: from bulk water to electrolyte solution

We implemented the osmostat in RASPA so that electrolyte units can be added or removed during a simulation, in the same way as in the work of Ross and co-workers.²⁴ By setting up the osmostat, a system is driven to chemical equilibrium with a saline reservoir at a given electrolyte concentration that required the calibration of $\Delta\mu_{C_nA_m}$ previously done in the first section. The maximal number of electrolyte units $N_{C_nA_m}^{\max}$ that can be present in the system is the floor function of $N_{H_2O}/(n+m)$ where $N_{H_2O}$ is the number of water molecules in the bulk system chosen as an initial configuration (150 water molecules). We briefly present a description of the osmostat implementation in RASPA.

* First, we choose whether to insert or delete an electrolyte unit

$$P_{\text{insert}} = \begin{cases} 1 & \text{if } N_{C_n A_m} = 0 \\ \frac{1}{2} & \text{if } 0 < N_{C_n A_m} < N_{C_n A_m}^{\max} \\ 0 & \text{if } N_{C_n A_m} = N_{C_n A_m}^{\max} \end{cases} \qquad (11)$$
$$P_{\text{delete}} = \begin{cases} 0 & \text{if } N_{C_n A_m} = 0 \\ \frac{1}{2} & \text{if } 0 < N_{C_n A_m} < N_{C_n A_m}^{\max} \\ 1 & \text{if } N_{C_n A_m} = N_{C_n A_m}^{\max} \end{cases} \qquad (12)$$

The choice of $P_{\text{insert}}$ and $P_{\text{delete}}$ ensures that electrolyte insertion is attempted when no electrolytes remain in the system while electrolyte deletion is attempted when the maximum number of electrolyte units in the system is reached.

* Then, an NCMC move is proposed for alchemically transmute chosen particles during which the alchemical work is recorded in the same fashion as in the section where calculation of chemical potential was performed. Again, the alchemical path involves 200 segments where each segment consists in perturbation through interpolation of non-bonded parameters of the chosen particles, followed by relaxation with 200 NVE MD steps. The alchemical work for transmuting water into an electrolyte unit is the same as defined previously $W_{C_n A_m}$ while the alchemical work for removing an electrolyte is $W_{-C_n A_m}$.

* Let us derive again the probability of accepting and deleting an electrolyte unit $C_n A_m$.²⁴ We recall that the density probability of a system of water and electrolytes in equilibrium with an osmostat at chemical concentration $\Delta\mu_{C_n A_m}$ in the semi-grand canonical ensemble writes

$$\pi(r, N_{C_n A_m}, N, \Delta\mu_{C_n A_m}, V, T) \propto \frac{N!}{N^{C_m+!} N_{A^n-!} N_{H_2O}!} e^{-\beta(U[\mathbf{r}, N_{C_n A_m}] - \Delta\mu_{C_n A_m} N_{C_n A_m})} \quad (13)$$

The probability of acceptance of a trial NCMC move in which (n+m) water molecules are transmuted into an electrolyte unit, denoted $\alpha(N_{C_n A_m} \rightarrow N_{C_n A_m} + 1)$, is the ratio
of the corresponding probability densities (equation 13):

$$ 
\begin{aligned}
\alpha(N_{C_n A_m} \rightarrow N_{C_n A_m} + 1) &= \frac{\pi(r, N_{C_n A_m} + 1, N, \Delta\mu_{C_n A_m}, V, T)}{\pi(r, N_{C_n A_m}, N, \Delta\mu_{C_n A_m}, V, T)} \\
\alpha(N_{C_n A_m} \rightarrow N_{C_n A_m} + 1) &= \frac{\prod_{i=0}^{m+n-1} (N_{H_2O} - i)}{\prod_{j=1}^{n} (N_{C^{m+}} + j) \prod_{k=1}^{m} (N_{A^{n-}} + k)} e^{-\beta(W_{C_n A_m} - \Delta\mu_{C_n A_m})}
\end{aligned}
\quad (14) $$

In the same manner, the probability of acceptance of a trial NCMC move in which an electrolyte unit is transmuted into (n + m) water molecules denoted $\alpha(N_{C_n A_m} \rightarrow N_{C_n A_m} - 1)$ writes

$$ 
\begin{aligned}
\alpha(N_{C_n A_m} \rightarrow N_{C_n A_m} - 1) &= \frac{\pi(r, N_{C_n A_m} - 1, N, \Delta\mu_{C_n A_m}, V, T)}{\pi(r, N_{C_n A_m}, N, \Delta\mu_{C_n A_m}, V, T)} \\
\alpha(N_{C_n A_m} \rightarrow N_{C_n A_m} - 1) &= \frac{\prod_{i=0}^{n-1} (N_{C^{m+}} - i) \prod_{j=0}^{m-1} (N_{A^{n-}} - j)}{\prod_{k=1}^{m+n} (N_{H_2O} + k)} e^{-\beta(W_{C_n A_m} + \Delta\mu_{C_n A_m})}
\end{aligned}
\quad (15) $$

Finally, the equations 14 and 15 lead to the MC moves implemented for insertion or deletion of an electrolyte unit that define the osmostat
$$ acc(N_{C_n A_m} \rightarrow N_{C_n A_m} + 1) = \min \left( 1, \frac{\prod_{i=0}^{m+n-1} (N_{H_2O} - i)}{\prod_{j=1}^{n} (N_{C^{m+}} + j) \prod_{k=1}^{m} (N_{A^{n-}} + k)} e^{-\beta(W_{C_n A_m} - \Delta\mu_{C_n A_m})} \right) \\ acc(N_{C_n A_m} \rightarrow N_{C_n A_m} - 1) = \min \left( 1, \frac{\prod_{i=0}^{n-1} (N_{C^{m+}} - i) \prod_{j=0}^{m-1} (N_{A^{n-}} - j)}{\prod_{k=1}^{m+n} (N_{H_2O} + k)} e^{-\beta(W_{-C_n A_m} + \Delta\mu_{C_n A_m})} \right) \quad (16) $$

We test our implementation of the osmostat in the semi-grand canonical ensemble by imposing a series of chemical potential for CsCl electrolyte $\Delta\mu_{CsCl}$ to MC simulations starting from bulk water. In SI, supplementary MC simulations were performed by giving extreme values for $\Delta\mu_{CsCl}$ to probe behaviour of osmostat in asymptotic cases: a highly negative $\Delta\mu_{CsCl}$ yields a system full of water while setting a highly positive $\Delta\mu_{CsCl}$ results in a box full of ions. In Figure 4, we present the evolution of the concentration of CsCl as a function of MC cycles for values of chemical potential of -567.8, -561.2, -551.1 and -525.6 kJ/mol (Table S6) that defines a saline reservoir at increasing concentration of respectively 2.41, 4.77, 10.09 and 18.50 M. Starting from bulk water, the systems in contact with the saline reservoir exchange ions and water molecules until equilibrium concentration is reached. A higher concentration of the reservoir required more MC steps to drive the system towards equilibrium. Only ~50 MC cycles are required to reach equilibrium in contact with reservoirs of 2.41 and 4.77 M while ~120 and 200 cycles are needed for equilibration with reservoirs of 10.06 and 18.46 M respectively.
<figure><img src="image_5.png" /><figcaption>Figure 4: Evolution of CsCl concentration for different chemical potential $\Delta\mu_{CsCl}$ as a function of MC cycles, for which the concentration of the saline reservoir is indicated in black for each case. The distribution of salt equilibrium is given on the right of the time-series plots.</figcaption></figure>

Once equilibrium is reached, the concentration in the systems fluctuates around the con-
centration of the saline reservoir with average values of $c = 2.07$, 5.01, 9.98 and 19.31 M,
which illustrates that NCMC moves are very sensitive to the calibration precision of the
chemical potential: a variation of $\sim 6.5$ kJ/mol of chemical potential induces the concentra-
tion to double from 2.41 to 4.77 M. Importantly, salt fluctuations completely drops when
concentration increases above saturation concentration in water which is associated to higher
fluctuations of the alchemical work during the transformation (Table S4 and Figure S1) and
more rejection of NCMC moves. Nonetheless, our osmostat implementation in RASPA is
able to efficiently sample the electrolyte concentration within a small amount of MC cycles.

Conclusion

A technique for insertion of electrolytes in water has been implemented in this work, based on
non-equilibrium candidate Monte Carlo moves as proposed in the previous work of Ross and
co-workers. At constant number of particles in the system, water molecules and electrolytes
are exchanged at a given chemical potential, so that MC steps that insert or delete ions
define an osmostat in contact with the system.

First, the calibration of the osmostat has been handled using the Widom insertion
technique requiring alchemical work extracted from proposed NCMC moves where water
molecules are alchemically transformed to an electrolyte unit. Our method is straightforward
to implement, at the expense of a slightly higher error on the computed alchemical potential
compared to the algorithm of Ross and co-workers. Then, we calculated solvation free energy
of water by decoupling Lennard-Jones and electrostatic interactions during insertion of water
molecules in the bulk using NCMC moves. Previous chemical potential for inserting a single
electrolyte unit in bulk water by alchemical transformation is used to grasp thermodynamics
of electrolyte solvation at infinite dilution. Our value of water solvation free energy as well as
electrolyte solvation free energies are in good agreement with the literature, demonstrating
that NCMC move is a versatile tool that can be used for both alchemical transormation and
particle insertion in dense fluid. Finally, the MC moves for inserting/removing an electrolyte
unit are similar to those of Ross and co-workers and show that our implementation of the
osmostat yields the correct electrolyte concentration in system initialized with pure bulk
water.

We chose to develop and test our own version of the osmostat in RASPA to have access
of functionalities already implemented in the software for studying adsorption of electrolytes solutions in nanoporous materials. We plan to use our osmostat to simulate adsorption of various electrolyte solutions in zeolithes or MOFs that requires calibration of osmostat of different electrolytes solution over a wide range of concentration, which was the main goal of the current work. This opens the route for probing thermodynamics of electrolyte intrusion, pression intrusion and structure of intruded electrolyte solution in nanoporous material which will be the object of the next study. However, because of the choice of our empirical force field that is not optimal for divalent cations, we do not focus on reproducing precisely the structure of divalent cation electrolytes although we expect to grasp the thermodynamics and thus concentration of divalent cation electrolyte in porous materials. In addition, a version of RASPA with implementation of the osmostat will be provided therein. As a last remark, the osmostat was tested on relatively small boxes of monovalent or divalent-cation aqueous electrolytes as a proof-of-concept because RASPA is a sequential program and it would be time consuming to consider larger systems. In this regards, zeolites/MOFs are in most cases periodic so that only a supercell of several repeating units is required for future MC simulations where only few hundreds of free particles are adsorbed which is fortunately feasible with RASPA.

## Acknowledgement

We acknowledge funding from Agence Nationale de la Recherche under project MESAMM (ANR-19-CE05-0031).

## Supporting Information Available

Force field parameters, evolution of alchemical work as a function of NCMC moves, error and raw data for alchemical work and chemical potential as well as radial distribution functions, water solvation and tests of the osmostat for extreme cases are provided in supplementary
References

(1) Ren, E.; Coudert, F.-X. Thermodynamic exploration of xenon/krypton separation based on a high-throughput screening. *Faraday Discuss.* **2021**, *231*, 201–223.

(2) Delgado-Buscalioni, R.; Coveney, P. USHER: an algorithm for particle insertion in dense fluids. *J. Chem. Phys.* **2003**, *119*, 978–987.

(3) Orkoulas, G.; Panagiotopoulos, A. Z. Chemical potentials in ionic systems from Monte Carlo simulations with distance-biased test particle insertions. *Fluid Phase Equilib.* **1993**, *83*, 223–231.

(4) Cui, S.; Harris, J. Solubility of sodium chloride in supercritical water: A molecular dynamics study. *J. Phys. Chem.* **1995**, *99*, 2900–2906.

(5) Shelley, J.; Patey, G. A configuration bias Monte Carlo method for ionic solutions. *J. Chem. Phys.* **1994**, *100*, 8265–8270.

(6) Perego, C.; Giberti, F.; Parrinello, M. Chemical potential calculations in dense liquids using metadynamics. *Eur. Phys. J. Spec. Top.* **2016**, *225*, 1621–1628.

(7) Snurr, R. Q.; Bell, A. T.; Theodorou, D. N. Prediction of adsorption of aromatic hydrocarbons in silicalite from grand canonical Monte Carlo simulations with biased insertions. *J. Phys. Chem.* **1993**, *97*, 13742–13752.

(8) Mezei, M. A cavity-biased (T, V, µ) Monte Carlo method for the computer simulation of fluids. *Mol. Phys.* **1980**, *40*, 901–906.

(9) Kirkwood, J. G. Statistical mechanics of fluid mixtures. *J. Chem. Phys.* **1935**, *3*, 300–313.
(10) Hummer, G.; Pratt, L. R.; Garcia, A. E. Free energy of ionic hydration. *J. Phys. Chem.* **1996**, *100*, 1206–1215.

(11) Ferrario, M.; Ciccotti, G.; Spohr, E.; Cartailler, T.; Turq, P. Solubility of KF in water by molecular dynamics using the Kirkwood integration method. *J. Chem. Phys.* **2002**, *117*, 4947–4953.

(12) Lynden-Bell, R.; Rasaiah, J.; Noworyta, J. Using simulation to study solvation in water. *Pure Appl. Chem.* **2001**, *73*, 1721–1731.

(13) Berg, B. A.; Neuhaus, T. Multicanonical ensemble: A new approach to simulate first-order phase transitions. *Phys. Rev. Lett.* **1992**, *68*, 9.

(14) Orkoulas, G.; Panagiotopoulos, A. Z. Phase behavior of the restricted primitive model and square-well fluids from Monte Carlo simulations in the grand canonical ensemble. *J. Chem. Phys.* **1999**, *110*, 1581–1590.

(15) Fenwick, M. K.; Escobedo, F. A. Expanded ensemble and replica exchange methods for simulation of protein-like systems. *J. Chem. Phys.* **2003**, *119*, 11998–12010.

(16) Siepmann, J. I.; Frenkel, D. Configurational bias Monte Carlo: a new sampling scheme for flexible chains. *Mol. Phys.* **1992**, *75*, 59–70.

(17) Laso, M.; de Pablo, J. J.; Suter, U. W. Simulation of phase equilibria for chain molecules. *J. Chem. Phys.* **1992**, *97*, 2817–2819.

(18) Nezbeda, I.; Moučka, F.; Smith, W. R. Recent progress in molecular simulation of aqueous electrolytes: Force fields, chemical potentials and solubility. *Mol. Phys.* **2016**, *114*, 1665–1690.

(19) Shi, W.; Maginn, E. J. Continuous fractional component Monte Carlo: an adaptive biasing method for open system atomistic simulations. *J. Chem. Theo. Comput.* **2007**, *3*, 1451–1463.
(20) Çağin, T.; Pettitt, B. M. Molecular dynamics with a variable number of molecules. *Mol. Phys.* **1991**, *72*, 169–175.

(21) Boinepalli, S.; Attard, P. Grand canonical molecular dynamics. *J. Chem. Phys.* **2003**, *119*, 12769–12775.

(22) Rahbari, A.; Hens, R.; Ramdin, M.; Moultos, O.; Dubbeldam, D.; Vlugt, T. Recent advances in the continuous fractional component Monte Carlo methodology. *Mol. Simul.* **2021**, *47*, 804–823.

(23) Nilmeier, J. P.; Crooks, G. E.; Minh, D. D.; Chodera, J. D. Nonequilibrium candidate Monte Carlo is an efficient tool for equilibrium simulation. *Proc. Nat. Acad. Sci.* **2011**, *108*, E1009–E1018.

(24) Ross, G. A.; Rustenburg, A. S.; Grinaway, P. B.; Fass, J.; Chodera, J. D. Biomolecular simulations under realistic macroscopic salt conditions. *J. Phys. Chem. B* **2018**, *122*, 5466–5486.

(25) Tan, Z. Optimally adjusted mixture sampling and locally weighted histogram analysis. *J. Comput. Graph. Stat.* **2017**, *26*, 54–65.

(26) Dubbeldam, D.; Calero, S.; Ellis, D. E.; Snurr, R. Q. RASPA: molecular simulation software for adsorption and diffusion in flexible nanoporous materials. *Mol. Simul.* **2016**, *42*, 81–101.

(27) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. Comparison of simple potential functions for simulating liquid water. *J. Chem. Phys.* **1983**, *79*, 926–935.

(28) Chowdhuri, S.; Chandra, A. Hydration structure and diffusion of ions in supercooled water: Ion size effects. *J. Chem. Phys.* **2003**, *118*, 9719–9725.
(29) Mamatkulov, S.; Fyta, M.; Netz, R. R. Force fields for divalent cations based on single-ion and ion-pair properties. *J. Chem. Phys.* **2013**, *138*, 024505.

(30) Widom, B. Some topics in the theory of fluids. *J. Chem. Phys.* **1963**, *39*, 2808-2812.

(31) Confalonieri, G.; Ryzhikov, A.; Arletti, R.; Quartieri, S.; Vezzalini, G.; Isaac, C.; Paillaud, J.-L.; Nouali, H.; Daou, T. J. Structural interpretation of the energetic performances of a pure silica LTA-type zeolite. *Phys. Chem. Chem. Phys.* **2020**, *22*, 5178-5187.

(32) Huve, J.; Daou, T. J.; Nouali, H.; Patarin, J.; Ryzhikov, A. The effect of nanostructures on high pressure intrusion-extrusion of water and electrolyte solutions in hierarchical nanoboxes of silicalite-1. *New J. Chem.* **2020**, *44*, 273-281.

(33) Ryzhikov, A.; Ronchi, L.; Nouali, H.; Daou, T. J.; Paillaud, J.-L.; Patarin, J. High-pressure intrusion-extrusion of water and electrolyte solutions in pure-silica LTA zeolite. *J. Phys. Chem. C* **2015**, *119*, 28319-28325.

(34) Soper, A. K.; Weckström, K. Ion solvation and water structure in potassium halide aqueous solutions. *Biophys. Chem.* **2006**, *124*, 180-191.

(35) Irudayam, S. J.; Henchman, R. H. Prediction and interpretation of the hydration entropies of monovalent cations and anions. *Mol. Phys.* **2011**, *109*, 37-48.

(36) Näslund, L.-Å.; Edwards, D. C.; Wernet, P.; Bergmann, U.; Ogasawara, H.; Pettersson, L. G.; Myneni, S.; Nilsson, A. X-ray absorption spectroscopy study of the hydrogen bond network in the bulk water of aqueous solutions. *J. Phys. Chem. A* **2005**, *109*, 5995-6002.

(37) Omta, A. W.; Kropman, M. F.; Woutersen, S.; Bakker, H. J. Negligible effect of ions on the hydrogen-bond structure in liquid water. *Science* **2003**, *301*, 347-349.
(38) Collins, K. D.; Neilson, G. W.; Enderby, J. E. Ions in water: Characterizing the forces that control chemical processes and biological structure. *Biophys. Chem.* **2007**, *128*, 95–104.

(39) Grossfield, A.; Ren, P.; Ponder, J. W. Ion solvation thermodynamics from simulation with a polarizable force field. *J. Am. Chem. Soc.* **2003**, *125*, 15671–15682.

(40) Jorgensen, W. L.; Blake, J. F.; Buckner, J. K. Free energy of TIP4P water and the free energies of hydration of CH₄ and Cl⁻ from statistical perturbation theory. *Chem. Phys.* **1989**, *129*, 193–200.

(41) Schmid, R.; Miah, A. M.; Sapunov, V. N. A new table of the thermodynamic quantities of ionic hydration: values and some applications (enthalpy-entropy compensation and Born radii). *Phys. Chem. Chem. Phys.* **2000**, *2*, 97–102.

(42) Duboué-Dijon, E.; Mason, P. E.; Fischer, H. E.; Jungwirth, P. Hydration and ion pairing in aqueous Mg2+ and Zn2+ solutions: force-field description aided by neutron scattering experiments and ab initio molecular dynamics simulations. *The Journal of Physical Chemistry B* **2017**, *122*, 3296–3306.

(43) Li, P.; Song, L. F.; Merz Jr, K. M. Parameterization of highly charged metal ions using the 12-6-4 LJ-type nonbonded model in explicit water. *The Journal of Physical Chemistry B* **2015**, *119*, 883–895.

(44) Saxena, A.; Sept, D. Multisite ion models that improve coordination and free energy calculations in molecular dynamics simulations. *Journal of chemical theory and computation* **2013**, *9*, 3538–3542.

(45) Leontyev, I.; Stuchebrukhov, A. Accounting for electronic polarization in non-polarizable force fields. *Physical Chemistry Chemical Physics* **2011**, *13*, 2613–2626.
TOC Graphic

<figure><img src="image_6.png" /></figure>