Osmotic Force Balance Evaluation of Aqueous Electrolyte Osmotic Pressures and Chemical Potentials

Concentrated aqueous salt solutions are ubiquitous in problems of biological and environmental relevance. The development of accurate force fields that capture the interactions between dissolved species in solution is crucial to simulating these systems to gain molecular insights into the underlying processes under saline conditions. The osmotic pressure is a relatively simple thermodynamic property connecting the experimental and simulation measurements of the associative properties of the ions in solution. Milner [C. Gillespie and S. T. Milner, Soft Matter, 16, 9816 (2020)] proposed a simulation approach to evaluate the osmotic pressures of salts in solution by applying a restraint potential to the ions alone in solution and determining the resulting pressure required to balance that potential, referred to here as the osmotic force balance. Here, we expand Milner’s approach, demonstrating that the chemical potentials of the salts in solution as a function of concentration can be fitted to the concentration profiles determined from simulation, additionally providing an analytical expression for the osmotic pressure. This approach is used to determine the osmotic pressures of 15 alkali halide salts in water from simulations. The cross interactions between cations and anions in solution are subsequently optimized to capture their experimental osmotic pressures.


■ INTRODUCTION
Aqueous electrolytes play a vital role in biological processes and in environmental waste remediation.Potassium 1 and sodium 2 channel proteins, for instance, regulate muscle and neuron activity, utilizing selectivity filters to promote the flow of specific cations across cell membranes. 3In bacterial colonies, on the other hand, ion channels conduct long-range electrical signals across biofilms to coordinate metabolic states among the cells. 4Multivalent metal ions like Ca 2+ and Zn 2+ act to moderate processes including cell fusion, 5 apoptosis, 6 enzymatic activity, 7 and protein binding to nucleic acids and other proteins. 8On the macro-scale, osmoreceptor nerve cells help regulate the osmotic balance in the body by moderating the intake and discharge of water and salt via thirst and renal excretion. 9−14 The osmotic pressure of ions provides the impetus for flow across semipermeable membranes.In the reverse osmotic purification of water from salt solutions, the pressure of a waste effluent stream must be increased so that the pressure difference between the effluent and purified water is greater than the osmotic pressure.If the pressure drops below this threshold, water can flow back to the effluent.The growth in the treatment of produced water, water extracted from underground formations as a byproduct of oil and gas production, 15 for either reuse in oil and gas extraction or for beneficial use (e.g., agriculture) is expected to outpace reinjection for the foreseeable future, accounting for more than half of the $10 billion water management market by 2025. 16Membrane separations have received significant attention for produced water treatment due to the fact they can produce higher-quality water via a less chemically intensive process while being cost competitive with alternate strategies. 17,18Salts comprise one of the dominant constituents of produced water. 15As such, a better understanding of the osmotic properties of salts in solution can beneficially impact the design of water purification systems.For example, simulations have shown that the pressure required to affect water flow through membranes against the osmotic pressure of aqueous NaCl solutions is significantly underestimated if the salt is assumed to obey the ideal gas law. 19he development of molecular simulation protocols and models to accurately capture the osmotic properties of salt solutions is an active area of ongoing research. 20Perhaps the earliest simulation approach for the evaluation of the osmotic pressure from molecular simulations was proposed by Murad and Powles. 21In their method, two semipermeable membranes are used to partition the simulation box into pure solvent and mixture cells.The osmotic pressure is determined by the pressure difference between the mixture and solvent cells.Luo and Roux 22 later adapted this technique to examine the osmotic equation-of-state of aqueous electrolyte solutions.Rather than model an explicit membrane, which can introduce solution structuring that can muddle interpretation of finitesized simulations, 23 they introduced half-harmonic potentials to act as virtual walls that act only on the ions but not the water.The osmotic pressure is determined by the force of the ions on the walls.Kohns et al. 24 used Luo and Roux's approach to evaluate the activity of water in NaCl solution.Saxena and Garci ́a, 25 on the other hand, used this approach to reparametrize divalent ion interactions to minimize clustering in the simulation to bring the osmotic pressure into agreement with experiment.Alternately, the osmotic pressure was evaluated from the free energies of various components in aqueous electrolyte solutions.The osmotic pressure is evaluated from these calculations as the pressure required to raise the chemical potential of water in the mixture to that of pure water at atmospheric pressure.Benavides et al., 26 for example, used thermodynamic integration to evaluate the free energies of NaCl in water, from which the chemical potential of water and osmotic pressure can be inferred.Nezbeda and co-workers, 27,28 on the other hand, used a combination of open ensemble and isothermal−isobaric ensemble simulation techniques to evaluate a wide range of solution properties of aqueous NaCl solutions including its osmotic pressure.
A difficulty with the techniques described above is that multiple simulations must be performed to examine the concentration dependence of the osmotic equation-of-state.Moreover, following the free energy simulation approaches, the osmotic pressure frequently is only indirectly inferred after processing of results for the chemical potentials of the ions and water in the mixture.Milner 23 recently proposed a simulation technique that circumvents these difficulties, referred to here as the osmotic force balance.In the osmotic force balance, an external potential acting along a single axis of the simulation box is applied to the solutes but not the solvent.This external potential is balanced by the osmotic pressure gradient that results from the nonuniform distribution of solutes in the simulation box.Milner and co-workers used this approach to evaluate the osmotic pressure of NaCl, 23 poly(ethylene oxide) 29 in water, and mixtures of benzene and pyridine. 30he osmotic force balance is an extension of the sedimentation equilibrium approach proposed by Hansen and co-workers 31,32 to evaluate the equation-of-state of colloidal particles.In these simulations, an effective gravitational field is applied to the particles in solution, giving rise to a concentration gradient in the simulation.Hansen's approach has been used in simulations and experiments to determine the equation-ofstate of hard spheres, 31,32 charged colloids, 32,33 Janus colloids, 34,35 granular suspensions, 36 and supercooled water. 37he main difference between the Milner and Hansen approaches is that Hansen considered one-component systems, while Milner selectively applied the external potential to individual components of a mixture to establish osmotic equilibrium.
Here, we expand the osmotic force balance method to determine the free energies of aqueous electrolytes over a range of concentrations from a single molecular simulation.We also explore the connection between calculations performed with the chemical potential of water fixed and calculations performed with the bulk pressure fixed to obtain an improved comparison between simulations and experiment.The methods developed here are then applied to simulations of 15 alkali halide salts in aqueous solution (excluding the fluoride salts) to evaluate their osmotic pressure, activity coefficients, and solution density.The Lennard-Jones interactions between cation−anion pairs are then optimized so that the simulations reproduce the experimental osmotic equationsof-state for each salt.
■ THEORY Osmotic Force Balance.Consider an electroneutral system of solutes (e.g., ions) in solution with an external potential applied to the solutes along the z-axis, U(z) but not to the solvent.Here, we assume that this potential is zero at the origin (i.e., U(0) = 0 at z = 0) and monotonically increases with distance from the origin.As such, the solutes will congregate near the minimum in the applied potential.The variation in the solute concentration subsequently establishes a gradient in the osmotic pressure, Π(z).Millner 23 demonstrated that the osmotic pressure at equilibrium satisfies the force balance where the left-hand side corresponds to the gradient of the osmotic pressure and the right-hand side is the opposing force acting on the solutes by the applied field.In this expression, C s (z) is the solute concentration at z, while ν is one for a neutral solute and is equal to the sum of the stoichiometric coefficients for a salt, e.g., ν = ν + + ν − for salt X ν+ Y ν− .As an aside, C s is determined as N s /V, where N s is the number of solute molecules and V is the volume.The number of cations and anions in solution are N + = ν + N s and N − = ν − N s , respectively.This expression reduces to that utilized by Hansen and co-workers 31,32 when the solute concentration is replaced by that of the single component in the system and the potential gradient term is replaced by the gravitational acceleration constant.The osmotic pressure at z is evaluated by integration of eq 1 where we have assumed Π(∞) → 0 as C s (∞) → 0 since the applied potential grows unbounded with distance from the origin.Given the salt concentration as a function of the position, we can determine the osmotic equation-of-state for the solute in solution, Π(C s ).From a practical point of view, the solute concentration profile would be evaluated via molecular simulation.In these simulations, the chemical potentials of both water and the solutes are expected to be uniform throughout the simulation box at equilibrium.While water's chemical potential, μ w , arises from its interactions with other water molecules and the solutes, the solute's chemical potential is where μ s (z|μ w ) is the solute's chemical potential arising from interactions of the solutes with water and the other solutes at constant μ w and μ s tot is the total solute chemical potential that includes the effect of the external potential.So while μ s tot is independent of position within the simulation at equilibrium, μ s is not.The osmotic pressure is subsequently determined by assuming that μ s (z|μ w ) is equal to the solute's chemical potential in a bulk solution at concentration C s (z).Specifically, at fixed temperature (dT = 0) and water chemical potential (dμ w = 0), we obtain from the Gibbs−Duhem equation where Π is the osmotic pressure, which is substituted for P when μ w is constant.When the total solute chemical potential is fixed (dμ s tot = 0), it follows from eq 3 that Substituting eq 5 into eq 4 yields When the derivative is taken with respect to z, this expression reduces to eq 1, demonstrating that the osmotic pressure obtained from the force balance conforms to that obtained following a thermodynamic path.To be thermodynamically consistent, however, we expect the results obtained from the osmotic force balance to be independent of the applied external potential so that C s (z) is representative of the bulk solute concentration at a given osmotic pressure.
Noninteracting Solutes.In the case of a noninteracting solute, the chemical potential is i k j j j j j y where μ s∞ ex is the excess chemical potential of the solute at infinite dilution, R is the ideal gas constant, and C 0 is a reference solute concentration.Substituting this expression into eq 5 and integrating, the solute concentration as a function of z is where C max is the maximum solute concentration at the external potential's minimum.Substituting this equation into eq 2, the osmotic equation-of-state is found to be Thus, as expected, the osmotic pressure for noninteracting solutes is governed by the ideal gas law.
Molar Description of Aqueous Electrolyte Solutions.The excess chemical potential of an aqueous electrolyte in solution at constant μ w can be described by the inclusion of an activity coefficient 28,38,39 i k j j j j j y where The activity coefficient expression (eq 10b) is termed a modified Debye−Huckel equation (modDH).The first term on the right-hand side of eq 10b corresponds to the extended Debye−Huckel expression for the free energy of salts in solution.While B is assumed here to be a fitting parameter, the constant A is where q j indicates the charge on the cation (j = +) or anion (j = −), ε is the solvent's dielectric constant, and ε 0 is the permittivity of free space.The summation on the right-hand side of eq 10b is an empirical polynomial correction to the extended Debye−Huckel expression up to the power i max with fitted constants α i .Integrating eq 5 from the potential minimum at the origin [where U(0) = 0] to z and substituting in eq 10a yields The B and α i parameters in this expression can be fitted to the simulation results for C s (z).On the other hand, the maximum salt concentration, C max , is determined by ensuring the fit gives the correct number of ions in solution where L j denotes the length of the simulation box in the j (=x, y, or z) direction, whose coordinates run from −L j /2 to L j /2 with the origin at the center of the box.Inversion of eq 12 to find a simple analytical expression for C s (z) is not possible, in general.Nevertheless, eq 13 can be numerically integrated by substituting values of C s less than C max into eq 12 and determining the value of z that satisfies this expression.In the case of a harmonic potential where k is the spring constant, z is determined as Once eq 12 has been fitted to the simulation results, the osmotic pressure is determined by substituting eq 10 into eq 4 and integrating From the point of view of the force balance (eq 1), it is most natural to assume that the osmotic pressure is a function of C s .This is not necessarily the case following alternate thermodynamic approaches for describing osmotic equilibrium.

Molal Description of Aqueous Electrolyte Solutions.
−41 The osmotic pressure is then determined as the pressure increment required to raise water's chemical potential in the mixture to that of pure water at atmospheric pressure.These approaches effectively fix the system at atmospheric pressure, which is distinct from the force balance.In this case, water's chemical potential is determined from the solute's chemical potential at fixed temperature (dT = 0) and pressure (dP = 0) using the associated Gibbs−Duhem relationship where m s (=N s /(M w kg N w )) is the solute molality and M w kg is the molecular weight of water in kg/mol.More generally, M i j is the molecular weight of solution component i in units of j = g or kg per mol.Thus, when pressure is fixed, the solute molality is the natural variable to describe the solution composition following eq 17.As such, the majority of experimental results for the activity of water in salt solutions and the osmotic pressure are reported as a function of molality.
The salt's chemical potential in solution can be determined by adapting the modDH expression (eq 10) to a molal dependence at fixed pressure where Considering the relationship between the solute concentration and molality in the low concentration limit, A and B from eq 10 can be related to where X is A or B and ρ w0 is the density of pure liquid water in g/cm 3 . 39As with eq 10, the parameters B ̃and α̃i are fitting parameters.
For a constant pressure description of electrolyte solutions, the osmotic pressure is determined as the pressure required to raise the free energy so that water's activity is 1, so that the chemical potential of water in the mixture at pressure P + Π is equal to that of pure water at pressure P. The activity of water at fixed pressure, a w (m s |P), is determined by substitution of eq 18 into eq 17 and integration along a line of constant pressure.
where μ w0 is the chemical potential of pure water at pressure P.
Assuming that the solution is incompressible, the osmotic pressure is determined as where V̅ w is the partial molar volume of water in the mixture, which can be estimated from the solution density as described below.Alternately, water's partial molar volume in salt solutions does not vary significantly from that of pure water, less than ∼3% difference over the salt concentration ranges considered here.The osmotic pressure is then well approximated using the molar volume of pure water in the absence of any information about water's volumetric properties in mixtures.Solution Densities.An additional thermodynamic variable that can be derived from osmotic force balance simulations is the dependence of the density on the solution composition.Specifically, the water concentration as a function of position, C w (z), is obtained from these simulations along with the solute concentration.The solution density profile within the simulation box can subsequently be evaluated, albeit at a constant μ w .The density at a constant P can then be determined from knowledge of the osmotic pressure and an estimate for the solution compressibility.While the solution concentration changes with the pressure, its molality does not.Here, then, we assume that the density at constant μ w or P is described by the expression where ρ w0 , the density of pure water, and the coefficients θ i j are fitting parameters determined along lines of constant pressure (j = P) or water chemical potential (j = μ, where the subscript w has been dropped for simplicity).The solution density at constant P can be approximated as where κ eff is the effective compressibility of the solution.We estimate the compressibility of the salt solutions in the Results and Discussion.The partial molar volume of water as a function of the salt molality can be estimated from eq 22 as which we use to evaluate the water's activity from the osmotic pressure using eq 21.

Journal of Chemical Theory and Computation
Fitting the Alternate the Osmotic Pressure.Here, we describe our procedure for fitting constant μ w or P descriptions of the osmotic pressure to our simulation results.In this fitting, A is determined from Debye−Huckel theory (eq 11), while the A and A ̃terms and B and B ̃terms are assumed to be related to one another through eq 19.This approximation is reasonable given that the terms appear in the extended Debye−Huckel contribution for the electrostatic free energy best apply at a low concentration where the difference between the constant μ w and P is minimal.An initial value of B = 1.5 M −1/2 was used for each fitting and then varied to minimize the difference between the osmotic pressures evaluated following eqs 16 and 21.The α i terms appearing in eq 10b are determined by a least-squares fit of eq 12 to the salt concentration profile determined from simulation under the mass balance constraint (eq 13).The α̃i terms appearing in eq 18b are then determined by a leastsquares fit of eq 21 to the osmotic pressures evaluated from eq 2 by numerical integration of the simulation concentration profile.We found that a value of i max = 2 provided an excellent description of our simulation results over the simulated range of salt concentrations, up to ∼3.5 M (corresponding to ∼4 molal).

■ MOLECULAR SIMULATION METHODS
Osmotic force balance molecular dynamics simulations of aqueous electrolyte solutions were conducted in the canonical ensemble 42 using GROMACS 2020.7. 43We considered salts composed of all combination of the alkali metal cations Li + , Na + , K + , Rb + , and Cs + with the halide anions Cl − , Br − , and I − , for a total of 15 monovalent 1:1 electrolytes in water (LiCl, LiBr, LiI, NaCl, NaBr, NaI, KCl, KBr, KI, CsCl, CsBr, CsI, RbCl, RbBr, and RbI).The fluoride anion proved to be problematic as a result of solubility issues that arose during the simulations, resulting in the fluoride salts being neglected.The temperature was set to 25 °C and was controlled using the Nose−Hoover thermostat. 44,45Water was modeled using the TIP4P/2005 potential. 46The dielectric constant of the model was taken to be 59.1 for fitting to the free energy models above. 47The ions were modeled using the potential parameters reported by Joung and Cheatham 48 parametrized for TIP4P/Ew water, 49 which have also been shown to accurately capture the single-ion properties of the monovalent salts in TIP4P/2005 water. 50Cross interactions of the ions with water were evaluated using Lorentz−Berthelot combining rules (Table 1). 42While initial cross interactions between the cations and anions were evaluated using Lorentz−Berthelot combining rules, 42 those interactions were subsequently optimized to reproduce the experimental osmotic pressures as described below.Lennard-Jones interactions were evaluated out to 9 Å with a long-range mean field correction for potential truncation effects.Electrostatic interactions were evaluated using Particle Mesh Ewald Summation. 51The internal holonomic constraints for water were fixed using the LINCS algorithm. 52The equations-of-motion were integrated with a 2 fs time step.Simulations were equilibrated for 25 ns, followed by a production run of 200 ns to determine equilibrium averages.
The primary simulation box used was a rectangular box of 30 × 30 × 100 Å along the x, y, and z directions.The simulations included 60 cations and 60 anions (N s = 60) and approximately 3000 waters, except in the case of CsI where N s = 40 to ensure the maximum concentration fell below the solubility limit.A harmonic potential (eq 14) was applied to the ions to confine the ions near the center of the z-axis while forming a reservoir of pure water on either end of the simulation box along z.The spring constant k was adjusted so that the salt concentration at the z origin was ∼3.5 M.An equilibrated snapshot from a simulation of NaBr in water is shown in Figure 1 illustrating the impact of the harmonic restraint on the distribution of cations and anions in the simulation.Since the simulations were conducted in the canonical ensemble, to ensure the chemical potential of water was effectively equal to that of pure water at 25 °C and atmospheric pressure, the number of waters in the final simulation for each ion pair was adjusted in a series of shorter preliminary simulations so that the density of water at either end of the box was 0.997 ± 0.0005 g/cm 3 , the density of pure TIP4P/2005 water at ambient conditions (the resulting bulk pressures lie within ±10 atm of atmospheric pressure).The density of pure water, ρ w0 , was determined by fitting eq 22 to our results.In addition to simulations with a 30 × 30 × 100 Å simulation box, we also simulated NaBr and RbCl in boxes a The optimal value of χ is also reported.Notably for CsI, the unmodified potential (χ = 0) produced the optimal result.The units for the Lennard-Jones diameter and well-depth are Å and K, respectively.S1 in the Supporting Information.Following Milner, 23 we fit the simulation results for the osmotic pressure to the experimental equations-of-state by modifying the repulsive portion of the cross interaction between unlike ions.This is achieved multiplying the repulsive r −12 (r represents the distance between ions) portion of the cross interaction by 1 + χ, where χ (>−1) is an adjustable optimization parameter, while not modifying the attractive r −6 portion of the potential = + The effect of this fitting scheme is to modify the Lennard-Jones well-depth and diameter as and (1 ) The best value of χ was determined by minimizing the mean square difference between the experimental and simulation results for the osmotic pressure as a function of the salt concentration.Experimental results for the osmotic pressures of the monovalent salts were taken from the extensive tables compiled by Hamer and Wu. 39Initial simulations are conducted with the cross interaction well-depth and diameter obtained from Lorentz−Berthelot combining rules (χ = 0).A series of additional simulations with differing values of χ were then run to find the optimal interaction.We typically started with large variations in χ to bound the optimal result, followed by additional simulations to zero in on a best fit.This would typically take ∼4 additional simulations to find a reasonable estimate of the optimum value following our initial simulation with χ = 0.
To compare the salt solution densities we obtained from the osmotic force balance simulations against those at atmospheric pressure, we conducted simulations of bulk 1, 2, and 3 molal aqueous solutions of LiI, KCl, and CsBr in the isothermal− isobaric ensemble.The pressure was maintained with a Parinello−Raman barostat. 53The ions were modeled by using the optimized potentials from the force balance simulations.The numbers of waters and ions in each of these simulations is reported in the Supporting Information (Table S2).These simulations were equilibrated for 25 ns, followed by a 200 ns run to evaluate the thermodynamic averages.
■ RESULTS AND DISCUSSION Thermodynamic Consistency.A necessary condition for the application of the osmotic force balance to accurately capture the thermodynamic properties of aqueous electrolyte solutions is that the results derived are independent of the specifics of the simulation boxes used.One potential consistency check is to fit the chemical potentials from one simulation of a given electrolyte and predict the concentration profile in a different simulation box by using the fitted chemical potential.Milner 23  while the harmonic constants of the external potential in the two simulations are related to one another as When the concentration profiles in the two simulations are plotted as a function of the dimensionally reduced zcoordinate, 2z/L z , they should collapse on top of one another if they are consistent.
To demonstrate thermodynamic consistency of the osmotic force balance, the concentration profiles for two representative salts, NaBr and RbCl, for three different simulations spanning a factor of 4 in volume are reported in Figure 2.These profiles exhibit a maximum at the center of the box (z = 0 Å) that tails off with increasing distance from the minimum in the applied potential.While at first glance these distributions may appear to be Gaussian, they are not normally distributed as a result of intermolecular interactions, as demonstrated below.More importantly, the concentration profiles of the different simulations collapse onto one another for each salt, indicative of thermodynamic consistency.The simulation of NaBr with L z = 28.23 Å (Figure 2a) exhibits the most significant differences between the concentration profiles, although the difference is small.This is the smallest simulation cell considered; however, we attribute this discrepancy to the fact that consistency is expected to breakdown as the external potential forces get larger as the box gets smaller.Nevertheless, the differences of this simulation with those performed using larger boxes are small, giving us confidence that we can obtain reliable results for the chemical potential as long as we use large enough simulation cells.Resultantly, we fit the chemical potentials below to simulations performed using the largest box considered, L z = 100 Å.
A second consistency condition for the application of the osmotic force balance is that the simulations are locally electroneutral.Specifically, the charge densities of cations, q + C + , and anions, q − C − , should cancel each another as a function of z (i.e., q + C + (z) + q − C − (z) = 0).Milner 23 demonstrated that the bulk concentrations of cations and anions can differ when the osmotic pressure is determined by the forces acting on a membrane-like restraint wall that induces ion structuring at its interface.It can be difficult in this case to even ascribe a specific salt concentration to the observed pressure, complicating the determination of the osmotic equation-of-state.Milner, however, did not demonstrate local electroneutrality for the osmotic force balance technique itself.The concentration profiles of the cations and anions comprising the NaBr and RbCl systems with L z = 100 Å are reported in Figure 3.The anion and cation concentration profiles of both salts match one another over the entire box length, which is expected for monovalent salts if they are locally electroneutral.This holds even down to submolar concentrations (Figure 3 inset).Indeed, statistical differences between the anion and cation concentrations are observed only below 1 mM (0.001 M) in these simulations, which can be attributed to deficiencies in sampling such low concentrations.Indeed, 1 mM corresponds to one ion pair per 55,000 waters, which is well above the number of waters in our simulations (∼3000 waters).Similarly good results are found for simulations of the other ions in water (not reported here).
When fitting the modDH chemical potential expression (eq 10) to our simulation concentration profiles below, we only consider data for concentrations above 10 mM, which is an order of magnitude above the concentration for which electroneutrality fails as a result of poor sampling.Taken together, the demonstrated box size independence and electroneutrality give us confidence that thermodynamically consistent results for the osmotic pressure and chemical potential of ions in water can be obtained from osmotic force balance simulations.
Osmotic Pressure Evaluation.The concentration profiles for NaBr and RbCl using the Lennard-Jones potentials (χ = 0) determined from simulation are compared to the fits of the modDH equation for the chemical potential in Figure 4a,b, respectively.These concentration profiles are described with quantitative accuracy by the modDH equation down into the sub-Molar concentration regime (Figure 4a,b inset).Similarly good agreement is found for all the other simulated salts as reported in the Supporting Information (Figures S1−S13).The fit parameters of the modDH equation (eq 10) to the concentration profiles for the salts with χ = 0 are reported in the Supporting Information (Table S3).
The osmotic pressures associated with the fits of the concentration profiles for NaBr and RbCl with χ = 0 are reported in Figure 5, while results for all other simulated salts are reported in the Supporting Information (Figures S14− S26).At low salt concentrations, both salts effectively behave ideally, although slight negative deviations are observed as a result of screening, as described by the Debye−Huckel limiting law.The osmotic pressures determined for both NaBr and RbCl using the unmodified Lennard-Jones potentials are greater than those in the corresponding experiments at higher concentrations.Nevertheless, both the simulation and experimental pressures are greater than would be obtained assuming the salts behaved like an ideal gas at an elevated concentration (C s ≳ 2 M), resulting from multibody interactions of the salts in solution.Interestingly, NaBr exhibits a markedly greater pressure at elevated concentrations than RbCl from both simulation and experiment, indicative of distinct differences between how these salts behave in solution.
To bring the osmotic pressures from simulation into agreement with experiment, the cross interactions between the anions and cations were modified following eq 25.As reported in Figure 5, the simulated osmotic pressures for NaBr and RbCl can be brought into excellent agreement with experiments following this approach.Similarly good agreement is obtained for all of the other simulated salts (Figures S14− S26 in the Supporting Information).The optimized Lennard-Jones cross-interaction parameters obtained following this scheme for all the simulated salts are reported in Table 1.(The unmodified interactions for CsI were notably found to be optimal.)As with the original Lennard-Jones interactions, the modDH expression provides an excellent quantitative description of the concentration profiles of NaBr and RbCl obtained by using the optimized Lennard-Jones interactions.
(Figure 4c,d.Results for the other simulated salts reported in Figures S1−S13 in the Supporting Information.)The fitted parameters of the concentration-dependent modDH expression for the chemical potential to the simulation results for all of the salts using the optimized Lennard-Jones interactions are reported in Table 2. Interestingly, while the predictions for osmotic pressures above using the unmodified Lennard-Jones cross interactions are too repulsive, the optimized values of χ for NaBr (χ = −0.41)and RbCl (χ = 0.50) differ in sign.Intuitively, for an overly repulsive pressure, we might expect the optimized χ would be negative to reduce the repulsion between the anions and cations, in line with what is observed for NaBr.RbCl, however, breaks this expectation.The osmotic pressure is the result of not only the cross interaction of anions with cations, however, but also from the interactions of the anions and cations with themselves and with water.We suspect that to achieve agreement between simulation and experiment at higher concentrations, the change in the cross interactions plays a nontrivial secondary role in the interactions of the anions and cations with their similarly charged brethren.Moreover, our simulations neglect any effects of polarization on intermolecular interactions, which could lead to additional sources of error in the simulations.Here, we have optimized the interactions to obtain agreement with experiment, but not the molecular origin of this behavior, which would require deeper investigation.Indeed, a range of negative and positive results for χ were found for all of the salts with no discernible correlation between the anion and cation pairs, suggesting potential new avenues of investigation.
Electrolyte Solution Thermodynamics.Before the free energies of the salts at atmospheric pressure can be determined, we must determine the volumetric properties of the solution to evaluate the activity of water from the osmotic pressure (eq 21).This, in turn, requires an estimate of the compressibility of the salt solutions (eq 23).The aqueous densities of LiI, KCl, and CsBr determined from the osmotic force simulations using the optimized Lennard-Jones inter- Table 2. Parameters of the Modified Debye−Huckel Expressions for the Free Energy of the Simulated Salts in Aqueous Solution; Table (a) Reports Fits Determined Directly from the Osmotic Force Balance Simulations at Constant μ w ; The Units of A, B, α 1 , and α 2 Are M −1/2 , M −1/2 , M −1 , and M −2 , Respectively; Table (b) Reports Fits at Constant P to the Osmotic Pressures; The Units of A ̃, B ̃, α1, and α2 Are molal −1/2 , molal −1/2 , molal −1 , and molal −2 , Respectively  actions, ρ(m s |μ w ), are compared against experiment in Figure 6.The reported simulation densities are lower than those of their experimental counterparts in this figure.Perhaps, this is not surprising given that the unmodified potentials were assessed by their ability to reproduce single ion properties, but not other solution properties. 50Indeed, the nonpolarizable Madrid model for aqueous electrolytes required optimization of not only the cross anion−cation Lennard-Jones interaction but also the anion−water oxygen and cation−water oxygen interactions and scaled-charges as well to reproduce a widerange of solution properties. 26,54That said, the predicted solution densities for KCl (Figure 6b) and CsBr (Figure 6c) only slightly underpredict the experimental densities, although the predictions for LiI (Figure 6a) are significantly off.Notably, classical simulations of Li + ions have reported a range of first-shell coordination numbers in water, indicative of uncertainty in the models. 55Given that the experiments were conducted at atmospheric pressure, while the simulations are conducted under osmotic stress, we should expect some error.
Removing the effect of the osmotic pressure on the solution densities, however, would only reduce the predicted densities and increase the difference between simulation and experiment.
A more appropriate comparison to assess the accuracy of the predicted densities is to compare the results from the force balance simulations against bulk salt solution simulations performed at ambient pressure.To this end, atmospheric pressure simulation results for LiI, KCl, and CsBr solutions at concentrations of 1, 2, and 3 molal using the optimized Lennard-Jones cross interactions are reported in Figure 6.While qualitatively similar to the results obtained at constant μ w , the densities at atmospheric pressure are slightly lower, with the difference between the two simulation series increasing with increasing concentration.This difference reflects the growing compression of the solution as a result of the osmotic pressure with an increasing concentration.While the solution compressibility depends on the concentration and the specific salts in solution, it is not expected to differ significantly from that of water.We subsequently assume we can use an effective value of the compressibility, κ eff , across all salts and solution concentrations to obtain an estimate of the solution density at atmospheric pressure.Excellent agreement is observed between the predictions for ρ(m s |P) using eq 23 and the bulk atmospheric pressure simulations using a value of κ eff = 3.95 × 10 −5 bar −1 .Interestingly, κ eff is only 20% lower than κ w for TIP4P/2005 water (=4.65 × 10 −5 bar −1 ).In the absence of information about the compressibility of these solutions, reasonable results can be obtained assuming κ eff = κ w .Fits of the densities at constant μ w and P to eq 22 for all of the simulated salts using the optimized Lennard-Jones potentials are reported in Table 3. Results for simulations using the unmodified Lennard-Jones interactions (χ = 0) are reported in the Supporting Information (Table S4).
Since the fits of the constant P expression of the modDH theory (eq 18) were performed to reproduce the simulation osmotic pressures, which themselves had been fitted to the experiment, the theory is expected to accurately capture the activity of water at ambient conditions given its simple proportionality with the osmotic pressure (results not shown).It is not assured, however, that the activity coefficients of the salts in solution will match experiment.The constant P activity coefficient fits of NaBr and RbCl are subsequently reported in Figure 7 to assess the accuracy of the simulations in capturing the experimental values.Activity coefficients for all the other simulated salts are reported in Figures S27−S39 in the Supporting Information.For NaBr and RbCl, the simulation results are in qualitative agreement with experiment.In the case of NaBr (Figure 7a), both the simulation and experimental results exhibit an initial steep drop, followed by a minimum in the neighborhood of a 0.5 to 1 molal solution.For RbCl (Figure 7b), the simulations and experiments exhibit the same initial drop, but afterward they both decrease more slowly up to 4 molal.The initial drop in the activity coefficients for both salts is a direct result of the leading m s 1/2 dependence dictated by the Debye−Huckel limiting law.This drop is more significant for the simulations as a result of the dielectric constant for TIP4P/2005 water (ε = 59.1) being lower than that for real water (ε = 78.3 56), resulting in a greater value of A ãs dictated by eq 11 for TIP4P/2005 water.More interestingly, the simulation results for both NaBr and RbCl are shifted downward from the experiment by 0.1 to 0.2 at higher concentrations.The deeper minimum in the activity coefficient  CsBr.We compare results from the osmotic force balance simulations (μ w constant), osmotic force balance simulations corrected to 1 bar pressure (P constant), experiment (experiment), and bulk simulations at 1 bar pressure (simulation).The lines were fit to eq 22.The figure symbols are defined in the legend of (a).could be remedied by either using a water model that exhibits the correct dielectric constant for water or by using scaled-ion charges to mimic dielectric screening effects as has been proposed by a number of investigators. 26,54,57,58The fits of the modDH model at ambient pressure (eq 18) to the osmotic pressures determined using the optimized interactions above are reported in Table 2. Fits for simulations using the unmodified Lennard-Jones interactions (χ = 0) are reported in the Supporting Information (Table S3).

■ CONCLUSIONS
Here, we utilized osmotic force balance simulations for measuring the osmotic pressures of ions in solution to examine the osmotic pressures of a wide range of alkali halide salts up to concentrations of ∼4 M. The utility of the osmotic force balance technique is that the osmotic pressure over a wide range of concentrations can be evaluated from a single simulation.Moreover, the application of an external potential applied to the ions to perform these simulations does not require specialized coding of the available molecular simulation packages.The cation and anion Lennard-Jones crossinteraction parameters were optimized to reproduce the experimental osmotic pressures of all of the simulated salts.We were able to fit the concentration profiles obtained from the osmotic force balance simulations utilizing an assumed form for the free energies of ions in solution, the modified Debye−Huckel equation.This theory quantitively described the salt concentration profiles from low to high concentrations, in addition to providing an analytical expression for the salt osmotic pressures.We also provided a framework for mapping osmotic force balance simulation results to measurements of the activity of water that is typically determined experimentally at ambient pressures, thereby facilitating an improved comparison between simulations and experiment.
Accurate force field parameters are essential for meaningful simulation studies of molecular association in solution.Force fields, however, are frequently only parametrized to reproduce solute properties at infinite dilution.The framework laid out here is not limited to monovalent electrolyte solutions but can be applied to a broader range of systems.For example, it has been noted that interactions for both monovalent 59,60 and multivalent 25 salts can have erroneous aggregation propensities, confounding predictions of charge-mediated properties in solution.Fitting ion potentials to properties such as the osmotic pressure could thereby potentially improve the description of salts at concentrations close to their solubility limits.Indeed, using an off the shelf force field for the fluoride ion here proved problematic, resulting from exaggerated solute  The units of ρ w0 , θ 1 j , θ 3/2 j , and θ 2 j are g/cm 3 , g/(cm 3 molal), g/(cm 3 molal 3/2 ), and g/(cm 3 molal 2 ), respectively.aggregation.Following the work of Vega and co-workers, 26 it would potentially be useful to also fit additional solution properties like the salt solution density, which could also necessitate adjustments of ion−water interactions as well.
−63 ■ ASSOCIATED CONTENT * sı Supporting Information The Supporting Information is available free of charge at https://pubs.acs.org/doi/10.1021/acs.jctc.3c00982.Details for the osmotic force balance simulations and isothermal−isobaric ensemble simulations; fits of the salt chemical potentials for the simulations using the unmodified (χ = 0) Lennard-Jones cross interactions to the modDH expressions at constant μ w and P; fits of the salt solution densities at constant μ w and P for the simulations using the unmodified (χ = 0) Lennard-Jones cross interactions to eq 22; simulation and fitted concentration profiles from the osmotic force balance simulations using the unmodified (χ = 0) and optimized Lennard-Jones cross interactions for all the simulated salts except NaBr and RbCl; osmotic pressures determined from the osmotic force balance simulations using the unmodified (χ = 0) and optimized Lennard-Jones cross interactions for all the simulated salts except NaBr and RbCl; and activity coefficients of the salts at atmospheric pressure using the optimized Lennard-Jones cross interactions for all the simulated salts except NaBr and RbCl (PDF) ■

Figure 1 .
Figure1.Harmonic potential, φ, applied to the Na + and Br − ions in aqueous solution in a 30 × 30 × 100 Å cell overlaid on top of the final snapshot from the simulation.The Na + and Br − ions as shown as the blue and red van der Waals spheres, respectively, while water is illustrated by the background stick representation.The simulation snapshot was visualized using VMD.64
proposed an alternate check that avoids having to fit the chemical potentials by comparing the concentration profiles from two simulations (1 and 2) with different sizes where the concentrations of water and salt in the two simulations are equal N

Figure 3 .
Figure 3. Cation and anion concentration as a function of z for (a) NaBr and (b) RbCl in a simulation box with dimensions 30 × 30 × 100 Å with χ = 0.The concentration profiles have been averaged about the line of symmetry (z = 0) to improve statistical accuracy.The inset figures show the concentration on a logarithmic scale to more clearly see the ion concentrations at down to sub mM concentrations.The simulation lines are defined in the legend of each figure.

Figure 4 .
Figure 4. Comparison of the simulation concentration profiles of NaBr and RbCl against the fits of the modDH theoretical expression (eq 10a).(a,c) Results for NaBr using the original (χ = 0) and optimized (χ = −0.41)Lennard-Jones cross interactions, respectively.(b,d) Results for RbCl using the original (χ = 0) and optimized (χ = 0.50) Lennard-Jones cross interactions, respectively.The concentration profiles have been averaged about the line of symmetry (z = 0) to improve statistical accuracy.The inset figures in (a,b) show the concentration on a logarithmic scale to more clearly see the ion concentrations at down to sub mM concentrations.The figure symbols are defined in the legend of (a).

Figure 5 .
Figure 5. Osmotic pressure of (a) NaBr and (b) RbCl in water as a function of the salt concentration.Results obtained from simulation using original (χ = 0) and optimized (χ = −0.41 and 0.50 for NaBr and RbCl, respectively) are compared against experiment and the ideal gas law.The experimental results are taken from Hamer and Wu.39The figure symbols are defined in the legend of each figure.
(a) Coefficients for μ s (C s |μ w ) along a Line of Constant μ w (eq 10 with A = 1.7964)

Figure 6 .
Figure 6.Density of the aqueous salt solutions as a function of the solution molality.Figures report results for(a) LiI, (b) KCl, and (c)CsBr.We compare results from the osmotic force balance simulations (μ w constant), osmotic force balance simulations corrected to 1 bar pressure (P constant), experiment (experiment), and bulk simulations at 1 bar pressure (simulation).The lines were fit to eq 22.The figure symbols are defined in the legend of (a).
Figure 6.Density of the aqueous salt solutions as a function of the solution molality.Figures report results for(a) LiI, (b) KCl, and (c)CsBr.We compare results from the osmotic force balance simulations (μ w constant), osmotic force balance simulations corrected to 1 bar pressure (P constant), experiment (experiment), and bulk simulations at 1 bar pressure (simulation).The lines were fit to eq 22.The figure symbols are defined in the legend of (a).

Figure 7 .
Figure 7. Activity coefficients of (a) NaBr and (b) RbCl in water at ambient pressure by fitting eq 20 to the osmotic pressures of the salts at 25 °C.The simulations were conducted using the optimized Lennard-Jones cross interactions.The experimental results are taken from Hamer and Wu.39The figure symbols are defined in the legend of (a).

Table 1 .
Cation−Anion Lennard-Jones Cross-Interaction Parameters for the Joung−Chetham 48 Potential Determined by Lorentz−Berthelot Combining Rules and Optimized to Reproduce the Experimental Osmotic Pressure Following eq 25 a