Ionic force field optimization based on single-ion and ion-pair solvation properties

Aqueous electrolyte solutions are of fundamental importance not only in physical chemistry but also for biological function and technological applications. In biology, the presence of ions, specifically of the monovalent ions $K+$⁠, $Na+$⁠, and $Cl−$⁠, has significant effects on the stability, structure, and function of nucleic acids and proteins and the regulation of biomolecular processes.^1–3 In technological applications, ions can play an important role in chemical reactions by influencing their rates, as well as in controlling the solubility of various cosolutes.^4,5 For salt concentrations larger than $∼10mM$ salt effects are typically ion specific even for simple bulk properties such as the osmotic pressure, which, in turn, can be highly relevant for transport and function in biomolecular systems.^1,3 The molecular understanding and prediction of the complex and often context-dependent effects of aqueous electrolytes pose a challenging task to the scientific and, in particular, theoretical community.^6,7

The successful molecular modeling of ionic effects typically involves computer simulations in which the ionic and water degrees of freedom are explicitly resolved and evolved by a set of effective, classical interactions: the simulation force field. Quite commonly used force fields are pairwise additive and nonpolarizable to keep the parameter space small. On that level, the atoms are characterized by (partial) Coulombic point charges $qi$⁠, excluded-volume radii, and dispersion attraction strengths. In standard protocols, the nonelectrostatic interaction for atoms $i$ and $j$ at a distance $rij$ is modeled by a pairwise Lennard-Jones (LJ) interaction of the form

with two free parameters, the interaction length $σij$ and the energy scale $εij$⁠, per pair of atoms. The whole set ${σij,εij,qi}$ with $i,j=1,…,M$⁠, defines the total force field behind the nonbonded inter- and intramolecular molecular dynamics (MD) interactions for $M$ atomic species. Typically, the vast number of parameters is reduced by using heuristic mixing rules for the cross interactions $(i≠j)$ so that the only remaining parameters are the diagonal coefficients $σii$ and $εii$⁠. The common mixing rules are $εij=εiiεjj$⁠, and either $σij=(σii+σjj)/2$⁠, constituting the Lorentz–Berthelot (LB) mixing rules, or $σij=σiiσjj$⁠, defining the geometric mixing rules.⁸ It is assumed that the force fields take implicitly into account the polarizability, as well as many-body effects. In fact, potentials which account for electric polarizability do not strictly seem to be required for modeling ion pairing. It has been shown that for mono- and divalent ions even the first hydration shell is not significantly polarized compared to water molecules in the bulk,⁹ though it should be kept in mind that a force field with more parameters has the principal possibility to be more accurate. For water usually simple point charge models (e.g., TIP3P or SPC/E) are used in which oxygen and hydrogen atoms are resolved.^10,11 The latter are connected by rigid intramolecular bonds and carry partial charges optimized in such a way that a few important water properties (density, structure, surface tension, and dielectric constant) are well reproduced.

Throughout the years, there have been numerous attempts to systematically optimize ionic force fields. In principle, to properly describe the interactions between ions and between ions and water, a high level of quantum theory is required, which however, turns out to be computationally too demanding for many-ion systems. Because of this, but also due to the weak computer power in the early days of force field development, it has become a common habit to parametrize the empirical force fields of ions based on single ion properties, such as ion solvation free energies or hydration structure in small water clusters, see for example Ref. 12. Force fields based on this procedure, though, often fail to reproduce realistically the ion-ion fluid structure and thermodynamics of the electrolytes at nonvanishing concentrations, even for simple ionic solutions such as NaCl.^13,14 As has been recognized in the literature, there is a strong sensitivity of solution thermodynamics^15–17 and contact ion pairing^18,19 to small changes in the effective pair potential between the interacting ions. Ion force field development is, thus, a difficult task and remains an active field of research.^20,21 Recently, two studies have revisited and systematically explored single ion hydration free energies by scanning through a wide region of the LJ parameter space ${σ,ε}$⁠.^22,23 This was triggered by the observation that while traditionally force field parameters have been adjusted in order to correctly reproduce single-ion free energies of solvation, one experimental observable is clearly not sufficient to unequivocally fix the two parameters $ε$ and $σ$⁠. In Ref. 22, crystal lattice energies have been used as a second independent optimization target. In Ref. 23, the single-ion solvation entropy and the effective solution ion size [as constructed from the ion-water radial distribution function (rdf)] have been used, though the simultaneous optimization of two parameters, especially for the cations, was problematic. It was seen that the three observables considered, namely, the free energy of solvation, the entropy of solvation, and the effective ion size, roughly matched the experimental values on a whole curve in the ${σ,ε}$ parameter plane, not allowing to single out one of the cationic parameter combinations as truly optimal.²³ It would be desirable to nail down the final cationic parameter set by benchmarking to collective, thermodynamic solution properties such as the electrolyte activity or osmotic pressure. Weerasinghe and Smith²⁴ introduced and carried out this idea for NaCl solutions by reproducing Kirkwood–Buff (KB) integrals as determined by experiments, ensuring that a good representation of solution activity is obtained. The same approach has been used recently to investigate the cation specific binding with protein surface charges.²⁵ However, the parametrization, which involves a free fit parameter in the mixing rule, does not conserve the free energy of solvation of single ions. It is therefore an interesting question, whether one can simultaneously describe single ion properties (such as the free energy of solvation) and ion-pairing properties by choosing optimized parameters for ${σ,ε}$ alone or whether an additional parameter has to be introduced, e.g., in the form of a generalized mixing rule as was recently suggested.²⁵ 

An alternative method for benchmarking MD force fields has been introduced by Hess et al.²⁶ and Kalcher and Dzubiella.²⁷ Here, effective, MD-derived ion-ion pair potentials are used in a many-body corrected virial route to obtain osmotic pressures. The electrolyte structure at a given concentration, which forms an input to the virial equation, can be obtained directly from a MD simulation or approximately, but with much less computational effort, from simulations with implicit solvent^14,26 or hypernetted chain (HNC) integral equation theory.²⁸ It was shown that the KB derived NaCl force field²⁴ and a few alkali-Cl force fields proposed by Dang²⁹ could reproduce experimental osmotic coefficients in SPC/E water quite well, while others badly failed.^26,27 The reason for the failure of some of the force fields when considering ion-pairing properties did not become clear.

In this work, we explore the optimization of ionic force fields based on single-ion and ion-pair thermodynamic properties, using the infinite-dilution solvation free energy and the finite-concentration electrolyte osmotic pressure as benchmarks. As for NaCl reasonable force fields exist, we use for $Na+$ and $Cl−$ the parameters given by Dang,^29,30 and focus on the salts KCl, NaI, KF, CsCl, and CsI in SPC/E water and systematically vary the force fields of $K+$⁠, $Cs+$⁠, $F−$⁠, and $I−$⁠. To satisfy the experimental ion solvation free energies, we confine our search in LJ parameter space to the experimental equisolvation free energy lines in $ε−σ$ space as calculated previously.²³ In this paper we do not vary systematically the mixing rule and in most simulations use the LB mixing rules. We apply the procedure proposed by Kalcher and Dzubiella²⁷ to calculate the electrolyte osmotic coefficients at a finite concentration of 1 M for a wide range of LJ parameters and compare to experiments. Our calculations are accompanied by HNC integral equation calculations that allow to efficiently cover and investigate a wide range of electrolyte concentrations. We systematically explore (a) whether ionic force field optimization is possible consistently for both single ion solvation energies and collective electrolyte thermodynamics without loosening the parameter space constraint given by the mixing rule, and (b) how the thermodynamics of electrolyte solutions react to a change of the LJ parameters along the equisolvation free energy path. As a main result, this simultaneous optimization of the force field seems possible for the cations $Cs+$ and $K+$⁠, while for the anions $I−$ and $F−$ the results are less promising. This suggests that in the long run, and in order to consistently describe finite-concentration electrolyte thermodynamics, the mixing rules have to be modified and systematically optimized. By calculating osmotic coefficients for a solution of CsI we also investigate the transferability of ion parameters for $Cs+$ and $I−$ that have been separately optimized by matching osmotic coefficients for CsCl and NaI. Modulo the previously mentioned restricting comments on the optimization of $I−$ parameters, ion parameters seem transferable in the sense that trends in the osmotic coefficients of certain ions also are found when those ions are assembled into different ion pairings.

We perform atomistic simulations using the MD package GROMACS^31,32 in the $(N,p,T)$ canonical ensemble, for which the particle number $N$⁠, as well as the pressure $p=1$ bar and temperature $T=300K$ are held constant using a Berendsen barostat and thermostat.³³ The simulation box is cubic, with an edge length of $L≃4nm$⁠, and periodically repeated in all three dimension and includes explicit ions and a total number of about 2000 SPC/E water molecules.¹⁰ Finite size effects are not significant for these sizes, as shown by a previous study on similar systems and sizes.²⁷ The three dimensional particle-mesh Ewald sum is used for the electrostatics³⁴ with a grid spacing in Fourier space of 0.12 nm in all three directions. We use an interpolation order of four, a distance cutoff of 0.9 nm for the real-space interactions, and a relative strength of the electrostatic interaction at the cutoff of $10−5$⁠. Typical times for the simulations for gathering statistics are 150 ns for low concentrations and 40–50 ns for moderate concentrations.

Five different salt solutions were simulated, CsCl, KCl, NaI, KF, and CsI at densities of 0.3 and $1M$ with 12 and 39 pairs of ions in the solution, respectively. Here, the concentration $M$ denotes molarity (mol/l). For the ions, charged and nonpolarizable spheres were used interacting with the LJ potential [Eq. (1)]. For the SPC/E water the LJ parameters are $εOO=0.6500kJ/mol$ and $σOO=0.3166nm$ and are assigned only to the oxygen molecule of water (no parameters are related to the two hydrogen atoms). Point charges of $qO=−0.8476e$ and $qH=+0.4238e$ are assigned to the oxygen and hydrogen atoms, respectively. For the ions we varied both parameters $εiO$ and $σiO$ and present the analysis in Sec. III. Only for sodium $(Na+)$ and chloride $(Cl−)$ we used fixed parameters, those given by Dang,^29,30 as they have been proven to be consistent in determining thermodynamic properties^25,27 and give reasonable hydration energies.^22,23 The LJ parameters are for $Na+$⁠, $εNa,O=0.5216kJ/mol$ and $σNa,O=0.2876nm$ and for $Cl−$⁠, $εCl,O=0.5216kJ/mol$ and $σCl,O=0.3785nm$⁠. In this notation, the subscripts $iO$ denote the parameters between ion $i$ and the oxygen atom of the SPC/E water model. For the cross interactions between two ions we use the LB mixing rules (except where noted otherwise).

For NaCl we also tried a different force field that was optimized based on KB integrals.²⁴ The parameters of this force field are for $Na+$⁠, $εNa,O=0.342kJ/mol$ and $σNa,O=0.279nm$ and for $Cl−$⁠, $εCl,O=0.547kJ/mol$ and $σCl,O=0.373nm$⁠. That force field involves a modified mixing rule for the relation between ion-ion and ion-water interaction, as it was not possible to fit the experimental data without breaking the geometric mixing rule.²⁴ 

We begin with a brief description of the derivation of the effective ionic (infinite dilution) pair potentials for the salts studied here, for more details see Ref. 27. The pair potentials are derived from the rdfs obtained within finite-concentration MD simulations. The rdf between a pair of atoms or ions $i$ and $j$ at distance $r$ is defined as $gij(r;ρ)$ at a given salt concentration $ρ$⁠. The potential of mean force (pmf) $wij(r;ρ)$ at concentration $ρ$ results from the rdf through a Boltzmann inversion,^35,36

where $β=1/kBT$ is the inverse thermal energy. The pmf can be decomposed into a short-ranged and a long-ranged contribution, $wijsr(r;ρ)$ and $wijlr(r;ρ)$⁠, respectively.^37–40 The long-ranged part of the pmf is a nonspecific Debye–Hückel potential and can be subtracted from $wij(r;ρ)$ leading in this way to the short-ranged part of the pair potential as detailed previously,²⁷ 

In the low concentration limit, the pmf between two ions reduces to their effective pair potential and the decomposition described in Eq. (2) can be written as

where the potential is split into the short-ranged part of the pair potential, $Vijsr(r)$⁠, and the usual Coulombic part. In this equation, $zi,zj$ are the valencies for the two ion types, respectively, and $λB(ρ)$ is the concentration dependent Bjerrum length with an infinite-dilution (pure water) value of $λB(0)=0.78nm$ for SPC/E water (about 10% larger than the real water value).⁴¹ The key assumption of our derivation is now that the short-ranged part of the pair potential, $Vijsr$⁠, can be extracted from the finite concentration pmf, $Vijsr(r)≃wijsr(r;ρ)$⁠, as calculated in Eq. (3), at not too high concentration. This is a good approximation, as long as the density is smaller than the density where the hydration layers of ions begin to overlap, which has been found to be between 0.5 and $1M$⁠.²⁷ On empirical grounds, the above procedure for the derivation of the ionic pair potentials works well for rdfs generated at a finite concentration of $ρ≃0.3M$⁠. Here, $Vijsr(r)≃wijsr(r;ρ)$ is well fulfilled and accurate rdfs can be sampled with good statistics.²⁷ In the following, the total effective pair potential $Vijeff(r)$ is used as an input to the pressure calculations by the virial route and as an input to the HNC method.

The optimization of the ionic force fields for the salts studied in this work is done by comparing the derived osmotic coefficients to their experimental values. We use the virial route to calculate osmotic coefficients as was done previously.²⁷ The osmotic pressure $Π=2ρϕ(ρ)/β$ of the ionic solution is defined through the osmotic coefficient $ϕ(ρ)$ at a concentration $ρ$⁠. The osmotic coefficient is given through the virial equation,³⁹ 

where the indices $i$ and $j$ represent the two salt components and $gij(r;ρ)$ needs to be evaluated at the respective concentration.

The virial route, as implemented in this work, is not exact as it employs the infinite dilution pair potential $Vijeff$⁠. Accordingly, many-body contributions to the ion-ion interactions for higher densities as induced by the water are not included. It has been suggested, though, that many-body contributions to the pair-potential can be qualitatively included by taking into account the concentration dependence of the water dielectric constant $ε(ρ)$⁠.^14,26 Thus, the long-range part in the pair potential $Vijeff(r)$ has to be altered by using $ε(ρ)$ instead of the infinite dilution limit $ε(0)$⁠. The following correction has been shown to lead to agreement of the virial route with the exact compressibility route up to a concentration of roughly $≃2M$⁠,²⁷ 

where again $λB(ρ)=βe2/4πε0ε(ρ)$ is the Bjerrum length in the aqueous electrolyte solution of concentration $ρ$⁠, and $ε(0)∼72±1$ for SPC/E water,²⁷ consistent with previous studies.⁴¹ The input parameter $ε(ρ)$ is directly calculated from the MD simulations and fitted through the function

where the values of the constant $A$ for each salt are shown in Table I.

A more efficient yet more approximate evaluation of the variation of the osmotic coefficient over a wide concentration range is possible with integral equation theory based on the HNC closure.¹⁵ The latter relates the pair correlation function $gij(r)$ between two particles $i$ and $j$ to the pair potential in an approximative way,³⁹ through

where $cij(r)$ is the direct correlation function, and $hij(r)=gij(r)−1$ is the total correlation function. The HNC equation is closed by the Ornstein–Zernicke (OZ) equation of liquid state theory, which relates $hij(r)$ and $cij(r)$⁠.⁴² For an $N$-component mixture with particle number densities $ρn$ the OZ equation is given as $hij(r)=cij(r)+∑k=1Nρk∫dr⃗′cik(r⃗−r⃗′)hjk(r′)$⁠. The HNC approach uses the $ε(ρ)$-dependent effective pair potential [Eq. (6)] from the MD simulations as input. Output is the liquid structure in the form of the rdfs $gij(r)$⁠. The osmotic coefficient of the salt solutions under consideration can then be calculated by the virial equation [Eq. (5)]. This approach has been used before and details can be found elsewhere.²⁸ The derived osmotic coefficient-concentration curves are in good agreement with the MD-derived ones but start to show significant deviations above $1.5–2M$⁠.

We note that experiments as well as atomistic MD simulations treat the system at the so-called Lewis–Randall (LR) level while the implicit HNC theory uses the McMillan–Mayer (MM) level.³⁶ At the MM level, the thermodynamic functions are calculated at constant chemical potential of the solvent, thus the MM and LR approaches correspond to different statistical ensembles. The pressure is given by the total pressure of the solution within LR and by the osmotic pressure of the solutes in equilibrium with the solution within MM. In order to compare the results between MD and the HNC approach used here, we follow the conversion,

where $m$ and $ρ$ are the molality and molarity of the solution, $Ms$ is the molar mass of the solute, and $ρ0$ and $ρ′$ are the mass densities of the pure solvent and the solution, respectively.²⁸ We thus use Eq. (9) to convert the HNC results to the LR level. Throughout the paper the osmotic coefficients shown are those corresponding to the LR level, and we use the notation $ϕ$ without a subindex.

As a crucial ingredient to our strategy, all LJ parameters investigated by us lie on the curve that reproduces the experimental free energy of solvation for a given single ion, which has been calculated previously.²³ This way, we can check whether parameter combinations exist that simultaneously reproduce single ion solvation as well as collective solution properties. All LJ ionic parameters employed in this work are depicted in Fig. 1.

The procedure we have followed here for the optimization of ionic force fields is the following: for each salt solution and each parameter set used, (a) we start with MD simulations of about 150 ns at a concentration of $0.3M$ and (b) we simulate the same system for about 30–50 ns at a concentration of $1M$⁠. We determine the rdfs between the ions and extract the effective pair potentials from the low-concentration simulation according to the methodology outlined in Sec. II B. The effective pair potential and the rdf lead to the osmotic coefficient for the specific solution at $1M$ according to the virial route and Eq. (5). Note that for step (a) the simulation time is longer as more statistics need to be gathered for the computation of the pmfs at $0.3M$⁠, compared to the rdfs at $1M$ concentration. In order to obtain the variation of the osmotic coefficient $ϕ(ρ)$ for a whole concentration range, we follow the HNC approach outlined above.

Optimization of the LJ parameter for each ion is attempted by comparing the osmotic coefficient at $1M$ as calculated from the MD simulations to the corresponding experimental values. The $ϕ(ρ)$ curves from the HNC calculations only serve as an indication whether the overall behavior of $ϕ$ is reasonably compared to the experimental curves and is not used for parameter optimization.

The salts that were modeled in this work are CsCl, KCl, NaI, and KF, with the goal to obtain optimized force fields for $Cs+$⁠, $K+$⁠, $I−$⁠, and $F−$⁠. As a check on the transferability of the obtained parameters, we also considered a solution of CsI with parameters optimized for CsCl and NaI. NaCl from Dang in SPC/E has been found to describe the osmotic properties well when compared to experiments,²⁷ while the individual ions also yield reasonable values for the solvation free energies.^22,23 For this reason, we do not attempt to optimize the force fields of $Na+$ and $Cl−$ and rather consider them as a given reference. We begin with MD simulations of CsCl, KCl, (for which the $Cl−$ force field is that from Ref. 30), and NaI (for which the $Na+$ force field is that from Ref. 29). MD simulations are performed for those three salt solutions for all ionic parameters of $Cs+$⁠, $K+$⁠, and $I−$ summarized in Fig. 1.The optimal ionic parameters for $Cs+$⁠, $K+$⁠, and $I−$ are then estimated from the comparison of the resulting osmotic coefficients to the experimental values. At the next step, we use the optimal LJ parameters for $K+$ in simulations of KF solutions and pursue a similar parameter space survey for $F−$ with the goal of finding an optimal parameter set for $F−$⁠. The choice of the salts NaI and KF is mainly motivated by the fact that the standard force fields for those salts gave very poor description of the osmotic coefficients in previous investigations.²⁷ As a consistency check, we finally take the optimized LJ parameters for $Cs+$ and $I−$ and consider CsI solutions and compare the calculated osmotic coefficient $ϕ$ for CsI to experiments.

A well-known problem of ionic force fields⁴³ has to be mentioned. Our MD simulations show that very low $εiO$ parameters for the cations (⁠$Cs+$ and $K+$⁠) lead to unphysical clustering of the ions for CsCl and KCl even at low concentrations of $0.3M$⁠. Accordingly, no pmfs could be extracted and the corresponding parameters (Cs1, Cs2, K1, K2, K7, and K9) are neglected for the calculation of the osmotic coefficient and for the optimization procedure of the $Cs+$ and $K+$ force fields. On the other hand, very low $εiO$ for the anions (⁠$I−$ and $F−$⁠) do not lead to similar clustering for NaI and KF, respectively, and can still be considered as candidates for an optimized force field. Interestingly, ion clustering was observed for all salts at a concentration $1M$ in some regions of the LJ parameter space. Specifically, the ions in KF clustered for the parameter sets F12, F13, for which the K11 parameter was used. For F1, F2, we could not obtain reasonable results, as the ion-water system could not be energetically optimized within MD, thus those parameters had to be eliminated as well. Apart from these restrictions, we have used all other parameters in Fig. 1 for the MD simulations. Parameter sets for which results are not shown for the osmotic coefficients in Fig. 5 are the ones that led to clustering of the ions in the aqueous environment.

A robust measure of electrolyte structural properties is the ion-ion rdf, which shows distinct structural signatures that differ among different salt solutions. Here, we have calculated the rdfs for all systems at two concentrations, 0.3 and $1M$⁠, respectively. The rdfs for the two concentrations and the same LJ parameters do not differ significantly, apart from the height of the rdf peaks. The heights of the contact and the solvent-separated peak indicate different hydration properties of the ions. All curves exhibit strong electrostatic screening and reach the asymptotic value of 1 below a distance of 2 nm for the $1M$ concentration. This is consistent with the small screening length of about 0.25 nm at $1M$⁠. As shown in Fig. 2, for CsCl and KCl the first peak in the cation-anion rdfs is much higher than the second one, indicating close contact of the anion and cation and predominant direct ion pairing. For KF and NaI, the first peak in the cation-anion rdfs is of similar height or lower than the second one, indicating that water enters between the anion and cation indicative of indirect ion pairing. Due to the variety of LJ parameters used in this study, the relative strength of direct and indirect ion pairing for the different salts shows rich behavior.^18,19 The height of the first rdf peak at $0.3M$ for CsCl ranges from about 6–30, for KCl from 5 to 25, brought about by variations of the cationic force field parameters, as illustrated in Fig. 2. For the ion pairs that show pronounced solvent-separated pairing, NaI and KF, the influence is much less; here the height of the first rdf peak varies for NaI from 2.5 to 3.5 and for KF from 9 to 12. Note that these variations are induced by changing the anionic force fields. This demonstrates that by changing force field parameters that keep the single-ion properties invariant (as judged by the ion-water rdf or the solvation free energy), the ion-pairing properties can be significantly affected.¹⁷ The position of the contact peak also shows interesting behavior. As an example, we present the trends for some representative LJ parameters for the KCl cation-anion rdfs: for the order of the bare LJ radius $σKO$⁠, $σKO(K11)<σKO(K8)<σKO(K6)<σKO(K5)<σKO(K3)$⁠, the ordering of the contact peak position between $K+$ and $Cl−$⁠, $rKClcp$ is $rKClcp(K3)<rKClcp(K5)<rKClcp(K6)<rKClcp(K8)<rKClcp(K11)$⁠, as is visible from Fig. 2. We find similar trends for CsCl, NaI, and KF, as well. Ions with smaller bare LJ radius thus show larger ion-ion separation, which clearly has to do with the fact that the interaction strengths $(εiO)$ of the single ions are varied along with $σiO$⁠. In addition to the peaks of the rdfs we also study the minima in the rdfs. The positions of the first and second minima in the rdfs, $r1$ and $r2$⁠, and the distance at which the rdfs vanish, $r0$⁠, are given for representative salt parameters at a concentration of $0.3M$ in Table II.

We have also compared the $g(r)s$ from our MD simulations with results from HNC for a few different salt parameters and observed only small differences; see Fig. 3. The reason for the deviations is the approximate treatment of statistical mechanics in HNC. The differences are mostly visible in the height of the first and second peaks in the $g(r)$⁠.

We derive short-ranged pair potentials for all salt parameters that do not lead to ion clustering from Eq. (3). Examples are shown in Fig. 4. In accordance with the rdfs, the pair potentials for NaI and KF reveal deeper second minima, while the first two minima for CsCl and KCl show roughly similar depth. The second minimum of KF is broader than for the other salts, indicating an unusual water configuration in the solvent-separated ion pair. The $Cs+–Cs+$ and $K+–K+$ pair potentials show smaller oscillations compared to the $Na+–Na+$ potential. This has been observed before, and was rationalized by the fact that $Na+$ has tightly bound hydration shells giving rise to energy barriers when two sodium cations approach.²⁷ In the case of the anion-anion pmfs, the $F−–F−$ pmf shows deeper minima and stronger oscillations than the $Cl−–Cl−$ or the $I−–I−$ potential. Going from the small fluoride to the large iodide, one observes a trend toward a soft-spherelike potential. Similar to the rdfs, we do not observe a simple dependence of the position and depth of the minima in the pmfs with the variation of the LJ parameters, again due to the simultaneous change of both the LJ radius and interaction strength along the lines of constant solvation free energies. Note again that there is substantial variation in the cation-anion potentials for the different force-field parameters, which gives hope to be able to fit the osmotic coefficients. The scattering in the $Cl−–Cl−$ potentials (note that the $Cl−$ force field is not changed) is due to bad sampling statistics as the ions typically do not get very close.

Having determined the rdfs and short-ranged pair potentials for different LJ parameters, we next calculate the osmotic coefficient $ϕ$ using the virial route based on the rdfs and pair potentials from MD simulations (not HNC) as explained in Sec. II The results for CsCl, KCl, NaI, and KF at a concentration of $1M$ are summarized in Fig. 5 where we have added lines as guides to the eye in order to bring out the main features of the results more clearly. The error bars for the calculated $ϕ$ are in the range ±0.01–0.05. The experimental values of the osmotic coefficients^44,45 for each salt are also shown, on which we base our optimization of the ionic force fields. Inspection of this figure reveals that for the salts CsCl, KCl, and KF $ϕ$ shows a maximum for intermediate values of $εiO$⁠. For NaI, no significant variation of $ϕ$ with $εiO$ is observed. We first focus on the cations $Cs+$ and $K+$ in CsCl and KCl, respectively, and panel (a) in Fig. 5. As more clearly shown by the lines added as guides to the eye for the CsCl and KCl data, there are two crossings between the simulated and the experimental values for $ϕ$⁠, so in principle there are two optimal force field sets for these cations. Note that we also included $ϕ$ for the Dang parameters for KCl and CsCl^27,30 as open symbols. Interestingly, though they are somewhat off from the curve corresponding to the experimental solvation free energy (see Fig. 1), they show quite good agreement with the experimental data for $ϕ$ within the error.

Inspection of panel (b) reveals that the situation for the anions $I−$ and $F−$ is distinctively different. For all iodide parameters, the corresponding $ϕ$ values are close to the experimental value (within the error), but show very little variation. This suggests that varying the force field for $I−$ has no considerable effect on the osmotic coefficients. A very similar result was obtained for chloride in NaCl in previous simulations.²⁴ However, this does not seem to be generally true for anions, as KF in the same panel shows a much larger variation in $ϕ$⁠, ranging from about 0.4 up to 0.7. For most of the results for KF shown in Fig. 5(b), the optimized K11 parameter set from Table III was used. However, as seen in the figure, the use of the equally optimal force field K5 does not lead to qualitative changes. For all $F−$ parameters, the osmotic coefficients of KF are too low compared to the experimental osmotic coefficient and even for the best $F−$ parameter combinations are too low by about 0.25, which is larger than the error bars. We conclude that while force fields for $Cs+$ and $K+$ can be derived with relative ease, the situation is more complicated for the anions considered by us; while $I−$ works quite well (which seems to be coincidental since the variation of $ϕ$ with $εIO$ is very small), the optimization for $F−$ fails though here the variation of $ϕ$ with $εIO$ is pronounced.

Using the HNC approach, we have also calculated the variation of the osmotic coefficient with the concentration for all salts and LJ parameters investigated here. Curves for some of the parameter sets used are shown in Fig. 6 together with experimental data. For some of the LJ parameters for fluoride, $ϕ$ for KF diverges above a certain concentration; for those cases the corresponding curves are cut above the concentration of $2M$⁠. Inspection of the overall shape of the curves reveals that there is good agreement with the experimental curves for some LJ parameters for CsCl, KCl and NaI, in contrast to KF. This is not unexpected, since for KF simulation results for $ϕ$ did not match the experimental value at $1M$ for any parameter combination. Note that there are small differences between the HNC results (lines) and MD results at $0.3M$ and $1M$ (symbols) for $ϕ$⁠, which are of the order of the error inflicted by the virial approximation and the simulation numerical scatter of about ±0.05, see our previous discussion.²⁸ 

Based on the MD osmotic coefficient results, we suggest optimal LJ parameters for $Cs+$⁠, $K+$ for use with $Cl−$ parameters from Dang,³⁰ and optimal parameters for $I−$ for use with $Na+$ from Dang;²⁹ we also suggest a force field for $F−$ derived from simulations of KF using our optimized force field for $K+$⁠. For most ions, two parameter sets are suggested and summarized in Table III. For the cations, $Cs+$ and $K+$⁠, we take two values closest to the crossing points of the fitting curves for $ϕ$ with the experimental data, see Fig. 5, namely, Cs6, Cs9, K5, and K11. For $K+$ a slight ambiguity exists, as for example also the K8 force field matches the experimental data. We rather chose K11 because that parameter is further away from the parameter sets K7 and K9 for which the ions were found to cluster. For iodide, our choice of parameters is less obvious; basically all parameters within the error bars are equally acceptable.

For KF, we performed two distinct sets of MD simulations, the first with the K11 force field and the second with the K5 parameters. Interestingly, both $K+$ give comparable results for $ϕ$⁠, see Fig. 5(b), meaning that the redundancy found with KCl seems to be also present for KF. However, none of the KF parameters reproduces experimental values for $ϕ$⁠, so we chose a single force field for $F−$ that shows the least deviation.

In Fig. 4, pmfs for the optimal parameters are compared with a pmf of a nonoptimal LJ parameter set. As expected, for $Cs+$ and $K+$ the cation-anion pmf of the nonoptimal force field shows larger deviations from the optimal force field results; for $I−$ all anion-cation pmfs are quite similar.

We briefly return to the discussion on peak heights and ion pairing.^18,19 For the optimal force fields in Table III, we find for the height of the first contact peak in the rdfs at $0.3M$ concentration, values of about 9.5 and 5.5 for K5Cl and K11Cl, respectively, 14.8 and 8.6 for Cs6Cl and Cs9Cl, respectively, about 2.3 for both NaI1 and NaI4, respectively, and 3.3 for K11F5. The order of the peak height is $KCl∼CsCl>KF>NaI$⁠, consistent with previous theories on ion pairing,^18,19 according to which the tendency to form direct ion pairs goes down as the ion sizes become more dissimilar. We note that this ordering of contact pair formation probability is only realized for the optimized force fields, nonoptimal force field combinations can easily lead to partial or complete reversal of this ordering.

We so far were occupied with finding force field parameters that match experimental osmotic coefficients best. We now turn to a separate question and check whether the optimized force fields presented in the previous section are transferable. To that end, we perform a set of MD simulations for CsI in water, for which the parameters for $Cs+$ and $I−$ are taken to be the optimized force fields from Table III. The hope can be not to match experimental osmotic coefficients for CsI perfectly, as the iodide parameters are not perfect by themselves (when compared with experimental osmotic coefficient data for NaI). Rather, we intend to check whether the trends of the simulated osmotic coefficients of CsI reflect the properties of the $Cs+$ and $I−$ force fields. The specific salts modeled are Cs6I4 and Cs9I4 (we did not check the I1 parameter set, as variations of iodide parameters do not seem to considerably affect the osmotic coefficients, as was seen from Fig. 5 for NaI). The MD and HNC results for the osmotic coefficient are shown in Fig. 7(a). Note that for Cs6I4 we only show the MD data for $0.3M$⁠, since the simulation data at $1M$ show bad convergence behavior, probably due to the vicinity to the crystallization transition for this particular force field (note that experimentally CsI has the smallest maximal solubility of all considered salts in this paper, of about $3M$⁠, so such problems might be anticipated). We find overall good agreement between the MD and HNC results with deviations in $ϕ$ smaller than 0.05. There are quite sizeable deviations between the theoretical predictions and the experimental data. However, note that we have used the optimized parameters for both $Cs+$ and $I−$⁠, based on osmotic coefficients for CsCl and NaI, and that the deviations from experimental data for CsI in Fig. 7(a) are of the same order as the deviations observed for CsCl and NaI in Fig. 5. Furthermore, the deviations from experimental data go in the expected direction, namely, the simulation prediction for Cs6I4 lies below the experimental value, while Cs9I4 lies above [similar to the data for Cs6Cl and Cs9Cl in Fig. 5(a)]. We tentatively conclude from these data that the force fields obtained by our optimization strategy are transferable, meaning that if one fixes the force fields of ions $A+$ and $D−$ by optimizing the osmotic coefficients of the salts AB and CD, the osmotic coefficient of the salt AD should come out approximately right without further adjustment. Such reasoning of course assumes that one can optimize the salts AB and CD and therefore does not apply to $F−$⁠, where the optimization of KF fails in the first place. So one sees that while optimization and transferability are logically distinct operations, transferability presumes successful optimization of force field parameters. We will come back to this point in the conclusions.

As an additional check of our methodology, we also used a different force field for NaCl that was designed to reproduce experimentally determined KB integrals (and thus experimental activity coefficients).²⁴ We performed simulations with those NaCl parameters using the mixing rules proposed, namely, geometric mixing including a freely adjusted mixing coefficient (and not the LB rule which we were using up to this point). The MD and HNC results for the osmotic coefficient are shown in Fig. 7(b), together with the experimental data and those derived previously from the NaCl Dang parameters.^27,29,30 It is evident from this figure that the KB NaCl force fields²⁴ give good agreement with experimental values for the osmotic coefficient, indicating that activity coefficient and osmotic coefficient data probe the same ionic features, namely, the pairing characteristics in aqueous solution.

We have used a combination of MD simulations and statistical mechanics analysis to systematically explore and optimize force field parameters for ions in aqueous solution, with particular emphasis on the interplay of single-ion and ion-pair thermodynamics.

The ion parameters considered by us, namely, the LJ radius and strength, were confined to a curve on which the experimental single ion solvation free energy is reproduced.²³ For a whole number of specific force fields on those curves, we constructed effective ionic pair potentials which led, through the virial route, to electrolyte osmotic coefficients. These were then compared to experimental values at $1M$ concentration. We have used the MD-virial route derived osmotic coefficients to optimize the LJ parameters for the single ions $Cs+$⁠, $K+$⁠, $I−$⁠, and $F−$ based on data for CsCl, KCl, NaI, and KF, respectively, treating $Na+$ and $Cl−$ as reference ions. This optimization strategy works well for the cations $Cs+$ and $K+$⁠, for which we obtain, due to the peculiar shape of the MD-based osmotic coefficient curves, two optimized force fields for each ion. This proves that at least for the cations, the simultaneous description of single-ion and ion-pairing thermodynamics is possible if one systematically explores the full LJ parameter space without the need to modify the combination rules. For anions, on the other hand, the optimization is more problematic. We could not get reasonable match of experimental osmotic coefficient data varying the LJ parameters of $F−$ in KF, simply because the maximum in the simulated osmotic coefficient is significantly below the experimental value. Iodide can be more or less satisfactorily optimized based on NaI data, but this seems to be mere coincidence, since NaI osmotic coefficients show very little dependence on the LJ parameter, which indicates a basic limitation of our approach. With these restrictions in mind, we checked for the transferability of our optimized force fields by considering the osmotic coefficient of a CsI solution, i.e., an ion combination that was not targeted in the optimization process. Based on the far-from-perfect performance of the $I−$ force field in NaI, we judge the transferability properties of the force fields as satisfactory since in particular the trends of the CsCl osmotic coefficients upon variations of the $Cs+$ force field are fully recovered when looking at the trends of the CsI osmotic coefficients.

These results suggest that for truly accurate nonpolarizable ion parameters, one apparently needs to lift the constraint of the simple mixing rules and introduce an additional scaling parameter in the mixing rule, in agreement with previous approaches.^24,25 However, it seems desirable and possible to introduce such a mixing parameter only in the cation-anion interaction, so that the anion-water and the cation-water interactions stay the same as used in the single-ion solvation free energy optimization (the cation-cation and anion-anion interactions are much less important). That way, one would correctly describe the single-ion properties (as embodied in the infinite-dilution solvation free energy) as well as ion-pairing properties (as important for osmotic and activity coefficients). As an unfortunate byproduct, one would have to optimize this mixing parameter for each ion pair, requiring substantial simulation efforts.

The authors wish to thank D. Horinek for discussions. We acknowledge support by the Ministry for Economy and Technology (BMWi) in the framework of the AiF project “Simulation and prediction of salt influence on biological systems,” by the NIM Cluster of Excellence in Munich and the Emmy-Noether Program of the Deutsche Forschungsgemeinschaft (DFG). M.F. acknowledges support from the “Gender Issue Incentive Funds” of the Cluster of Excellence in Munich. The Leibniz Rechenzentrum Munich is acknowledged for supercomputing access.