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

solvation properties

Maria Fyta1 , Immanuel Kalcher1 , Joachim Dzubiella1 , Luboš Vrbka2 , and Roland R. Netz1
1
arXiv:0912.4526v1 [physics.chem-ph] 22 Dec 2009

Physics Department (T37), Technical University of Munich, 85748 Garching, Germany

2
Institute for Physical and Theoretical Chemistry,
University of Regensburg, 93040 Regensburg, Germany

Abstract
Molecular dynamics simulations of ionic solutions depend sensitively on the force fields employed
for the ions. To resolve the fine differences between ions of the same valence and roughly similar
size and in particular to correctly describe ion-specific effects, it is clear that accurate force fields
are necessary. In the past, optimization strategies for ionic force fields either considered single-ion
properties (such as the solvation free energy at infinite dilution or the ion-water structure) or ion-
pair properties (in the form of ion-ion distribution functions). In this paper we investigate strategies
to optimize ionic force fields based on single-ion and ion-pair properties simultaneously. To that
end, we simulate five different salt solutions, namely CsCl, KCl, NaI, KF, and CsI, at finite ion
concentration. The force fields of these ions are systematically varied under the constraint that the
single-ion solvation free energy matches the experimental value, which reduces the two-dimensional
{σ, ε} parameter space of the Lennard Jones interaction to a one dimensional line for each ion. From
the finite-concentration simulations, the pair-potential is extracted and the osmotic coefficient is
calculated, which is compared to experimental data. We find a strong dependence of the osmotic
coefficient on the force field, which is remarkable as the single-ion solvation free energy and the
ion-water structure remain invariant under the parameter variation. Optimization of the force field
is achieved for the cations Cs+ and K+ , while for the anions I− and F− the experimental osmotic
coefficient cannot be reached. This suggests that in the long run, additional parameters might
have to be introduced into the modeling, for example by modified mixing rules.

PACS numbers:

1
I. INTRODUCTION

Aqueous electrolyte solutions are of fundamental importance not only in physical chem-
istry but also for biological function and technological applications. In biology, the presence
of ions, specifically of the monovalent ions K+ , Na+ , Cl− , has significant effects on the stabil-
ity, structure, and function of nucleic acids and proteins and the regulation of biomolecular
processes [1, 2, 3]. In technological applications, ions can play an important role in chem-
ical reactions by influencing their rates, as well as in controlling the solubility of various
co-solutes [4, 5]. For salt concentrations larger than ∼ 10 mM 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 molec-
ular understanding and prediction of the complex and often context-dependent effects of
aqueous electrolytes poses 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 pair-wise additive and non-polarizable 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 non-electrostatic interaction
for atoms i and j at a distance rij is modeled by a pairwise Lennard-Jones (LJ) interaction
of the form
" 12  6 #
σij σij
VLJ (rij ) = 4εij − , (1)
rij rij

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 intra-molecular 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 6= 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 [8]. It is assumed that the force fields take implicitly into account

2
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 [9], 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, 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 be-
tween 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 proper-
ties, 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 real-
istically the ion-ion fluid structure and thermodynamics of the electrolytes at non vanishing
concentrations, even for simple ionic solutions such as NaCl [13, 14]. As has been recog-
nized in the literature, there is a strong sensitivity of solution thermodynamics [15, 16, 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 [22], crystal lattice energies have been used as a second independent
optimization target. In [23], the single-ion solvation entropy and the effective solution ion
size (as constructed from the ion-water radial distribution function) 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

3
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 [23]. 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 [24]. The same approach has been used recently to investigate the cation-specific
binding with protein surface charges [25]. 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 [25].
An alternative method for benchmarking MD force fields has been introduced by Hess et
al. [26] and Kalcher and Dzubiella [27]. 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 ob-
tained directly from an MD simulation or approximately, but with much less computational
effort, from simulations with implicit solvent [14, 26] or hypernetted-chain (HNC) integral
equation theory [28]. It was shown that the KB derived NaCl force field [24] and a few alkali-
Cl force fields proposed by Dang [29] could reproduce experimental osmotic coeffcients 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 equi-solvation free energy lines in ε−σ
space as calculated previously [23]. In this paper we do not vary systematically the mixing

4
rule and in most simulations use the LB mixing rules. We apply the procedure proposed
by Kalcher and Dzubiella [27] 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 equi-solvation 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.

II. METHODS

A. Simulation details

We perform atomistic simulations using the Molecular Dynamics (MD) package GRO-
MACS [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 [33]. The simulation box is cubic, with an edge length of L ≃ 4
nm, and periodically repeated in all three dimension and includes explicit ions and a total
number of about 2000 SPC/E water molecules [10]. Finite size effects are not significant
for these sizes, as shown by a previous study on similar systems and sizes [27]. The three
dimensional particle-mesh Ewald sum is used for the electrostatics [34] with a grid spacing

5
in Fourier space of 0.12 nm in all three directions. We use an interpolation order of 4, a
distance cutoff of 0.9 nm for the real-space interactions, and a relative strength of the elec-
trostatic 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 M and 1 M with 12 and 39 pairs of ions in the solution, respectively. Here, the
concentration M denotes molarity (mol/l). For the ions, charged and non-polarizable spheres
were used interacting with the LJ potential (Eq. 1). For the SPC/E water the LJ parameters
are εOO = 0.6500 (kJ/mol) and σOO = 0.3166 (nm) 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
Section 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+ , εN a,O = 0.5216 (kJ/mol) and σN a,O = 0.2876 (nm) and for Cl− , εCl,O = 0.5216
(kJ/mol) and σCl,O = 0.3785 (nm). 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 Lorentz-Berthelot mixing rules (except where noted otherwise).
For NaCl we also tried a different force field that was optimized based on Kirkwood-
Buff integrals [24]. The parameters of this force field are for Na+ , εN a,O = 0.342 (kJ/mol)
and σN a,O = 0.279 (nm) and for Cl− , εCl,O = 0.547 (kJ/mol) and σCl,O = 0.373 (nm).
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 [24].

B. Effective ionic pair potentials

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 [27]. The pair potentials are
derived from the radial distribution functions (rdfs) obtained within finite-concentration
MD simulations. The rdf between a pair of atoms or ions i and j at distance r is defined

6
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]:

βwij (r; ρ) = − ln[gij (r; ρ)] = β wijsr (r; ρ) + wijlr (r; ρ) ,

 
(2)

where β = 1/kB T 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, 38, 39,
40]. The long-ranged part of the pmf is a non-specific 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 [27],

wijsr (r; ρ) = wij (r; ρ) − wijDH (r; ρ). (3)

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:

βVijeff (r) = βVijsr (r) + zi zj λB (0)/r (4)

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.78 nm for SPC/E water (about 10% larger than
the real water value) [41]. 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 (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 betwen 0.5 and 1 M [27]. 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.3 M. Here, Vijsr (r) ≃ wijsr (r; ρ) is well fulfilled and
accurate rdfs can be sampled with good statistics [27]. 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.

C. Virial route to the osmotic coefficient φ(ρ)

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

7
route to calculate osmotic coefficients as was done previously [27]. 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 [39]:

π X ∞ dβVijeff (r) 3
Z
φ(ρ) = 1 − ρ gij (r; ρ) r dr, (5)
3 i,j 0 dr

where the indices i and j represent the two salt components and gij (r; ρ) needs to be evalu-
ated 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 ≃ 2 M [27]:

zi zj
Ṽijeff (r) = Vijeff (r) − [λB (0) − λB (ρ)] (6)
r

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 [27], consistent with previous studies
[41]. The input parameter ε(ρ) is directly calculated from the MD simulations and fitted
through the function
ε(0)
ε(ρ) = (7)
(1 + Aρ)
where the values of the constant A for each salt are shown in Table I.

CsCl KCl NaI KF CsI

A [1/M] 0.23 0.24 0.34 0.19 0.26

TABLE I: Values for the A parameter in Eq. 7 as calculated from this and previous work [27].

8
 D. The hypernetted chain (HNC) approach

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 [15]. The latter relates the pair correlation function gij (r) between two particles i
and j to the pair potential in an approximative way [39], through:
h i
gij (r) = exp − β Ṽijef f + hij (r) − cij (r) (8)

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) [42]. For an N-component mixture with particle
number densities ρn the OZ equation is given as hij (r) = cij (r) + N
P R ′
k=1 ρk d~r 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 form of the radial distribution
functions 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 [28]. 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 [36]. 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:
ρ0 mρ0
φM M = φLR (1 + mMs ) = φLR , (9)
ρ′ ρ
where m and ρ are the molality and molarity of the solution, Ms the molar mass of the
solute, ρ0 and ρ′ the mass densities of the pure solvent and the solution, respectively [28].
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

9
φ without a subindex.

E. Choice of parameters and optimization procedure

As a crucial ingredient to our strategy, all Lennard-Jones 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 [23]. This way, we can check whether parameter
combinations exist that simultaneously reproduce single ion solvation as well as collective
solution properties. All Lennard-Jones 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 fol-
lowing: for each salt solution and each parameter set used, (a) we start with MD simulations
of about 150 ns at a concentration of 0.3 M and (b) we simulate the same system for about
30-50 ns at a concentration of 1 M. We determine the rdfs between the ions and extract the
effective pair potentials from the low-concentration simulation according to the methodology
outlined in II B. The effective pair potential and the rdf leads to the osmotic coefficient for
the specific solution at 1 M 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.3 M, compared to the rdfs at 1 M 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 1 M 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, 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
and activity well when compared to experiments [27], while the individual ions also yield
reasonable values for the solvation free energies [22, 23]. For this reason, we do not attempt

10
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 [27]. 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 [43] 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.3 M. 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.

11
III. RESULTS AND DISCUSSION

A. Structural properties

1. Radial distribution functions

A robust measure of electrolyte structural properties is the ion-ion radial distribution

function (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 M
and 1 M, 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 1 M concentration. This is consistent with the small screening length of about
0.25 nm at 1 M. 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.3 M for CsCl ranges
from about 6 to 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 [17]. 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− , rcp cp cp
KCl is rKCl (K3) < rKCl (K5) <
cp cp cp
rKCl (K6) < rKCl (K8) < rKCl (K11), as is visible from Fig. 2. We find similar trends for

12
CsCl, NaI, and KF, as well. Ions with smaller bare Lennard Jones 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.3 M in Table II.

Cs6Cl Cs9Cl K5Cl K11Cl NaI1 NaI4 K11F5

r0 (nm) 0.29 0.30 0.27 0.28 0.28 0.28 0.29

r1 (nm) 0.43 0.43 0.40 0.40 0.41 0.40 0.43

r2 (nm) 0.67 0.68 0.63 0.63 0.64 0.65 0.66

TABLE II: The rdf characteristics at 0.3 M for the optimized ion parameters; see Tab. III. The
distance (r0 ) at which the rdf vanishes to zero and that of the first (r1 ) and second (r2 ) minimum
in the cation-anion rdfs for representative parameters of Cs+ , K+ , I− , and F− in CsCl, KCl, NaI,
and KF aqueous environments, respectively.

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).

2. Short-ranged potentials of mean force

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 rdf’s, 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 [27]. In the case of the anion-anion pmfs, the F− -F− pmf shows deeper minima and

13
stronger oscillations than the Cl− -Cl− or the I− -I− potential. Going from the small fluoride
to the large iodide, one observes a trend towards a soft-sphere like 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.

B. Osmotic coefficients

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 Section II. The results
for CsCl, KCl, NaI, and KF at a concentration of 1 M 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

14
field for I− has no considerable effect on the osmotic coefficients. A very similar result was
obtained for chloride in NaCl in previous simulations [24]. 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 combination 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 1 M 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 [28].

C. Optimum ionic force field parameters

Based on the MD osmotic coefficient results, we suggest optimal LJ parameters for Cs+ ,
K+ for use with Cl− parameters from Dang [30], and optimal parameters for I− for use with
Na+ from Dang [29]; 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

15
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.

ion/label σiO (nm) εiO (kJ/mol)

Cs6 0.333 0.5

Cs9 0.325 1.0

K5 0.31 0.41

K11 0.293 1.26

I1 0.45 0.1

I4 0.425 0.32

F5 0.3665 0.1

TABLE III: The optimized Lennard-Jones parameters for Cs+ , K+ , I− , and F− as extracted by
the MD simulations and the comparison of the osmotic coefficients with experimental data.

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. 5b, 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 non-optimal LJ
parameter set. As expected, for Cs+ and K+ the cation-anion pmf of the non-optimal 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.3 M 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,

16
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, non optimal
force field combinations can easily lead to partial or complete reversal of this ordering.

D. Transferability of the optimized ionic force fields

We so far were occupied with finding force field parameters that match experimental os-
motic 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 be 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 crystalization transition for this particular force field (note that ex-
perimentally 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 coeffi-
cients 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. 5a). We tentatively conclude from these data that

17
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 Kirkwood-Buff integrals (and
thus experimental activity coefficients) [24]. We performed simulations with those NaCl
parameters using the mixing rules proposed, namely geometric mixing including a freely
adjusted mixing coefficient (and not the Lorentz-Berthelot 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 pa-
rameters [27, 29, 30]. It is evident from this figure that the Kirkwood-Buff NaCl force
fields [24] give good agreement with experimental values for the osmotic coefficient, indi-
cating that activity coefficient and osmotic coefficient data probe the same ionic features,
namely the pairing characteristics in aqueous solution.

IV. SUMMARY AND CONCLUSIONS

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 Lennard-Jones radius and strength,
were confined to a curve on which the experimental single ion solvation free energy is re-
produced [23]. 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 co-
efficients. These were then compared to experimental values at 1 M 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, respec-

18
tively, 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
satisfactorly optimized based on NaI data, but this seems to be mere coincidence, since
NaI osmotic coefficients shows 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 non-polarizable ion parameters, one appar-
ently 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]. How-
ever, 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 un-
fortunate byproduct, one would have to optimize this mixing parameter for each ion pair,
requiring substantial simulation efforts.

19
 Acknowledgments

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 Forschungsgemein-
schaft (DFG). MF acknowledges support from the ’Gender Issue Incentive Funds’ of the
Cluster of Excellence in Munich. The Leibniz Rechenzentrum Munich is acknowledged for
supercomputing access.

[1] C. Anderson and M. Record, Annu. Rev. Biophys. Biophys. Chem. 46, 657 (1995).
[2] R. L. Baldwin, Biophys. J. 71, 2056 (1996).
[3] B. Alberts, A. Johnson, J. Lewis, M. Raff, K. Roberts, and P. Walter, Molecular Biology of
the Cell (Garland Science, New York, 2002).
[4] W. Kunz, Pure Appl. Chem. 78, 1611 (2006).
[5] A. Kumar, Chem. Rev. 101, 1 (2001).
[6] K. D. Collins, Methods 34, 300 (2004).
[7] P. Jungwirth and D. J. Tobias, Chemical Reviews 106, 1259 (2006).
[8] D. Frenkel and B. Smit, Understanding Molecular Simulation: From Algorithms to Applica-
tions (Academic Press, 1996).
[9] C. Krekeler and L. Delle Site, J. Phys.: Cond. Matt. 19, 192101 (2007).
[10] H. J. C. Berendsen, J. R. Grigera, and T. P. Straatsma, J. Phys. Chem. 91, 6269 (1987).
[11] W. L. Jorgensen, J. Chandrasekhar, J. D. Madura, R. W. Impey, and M. L. Klein, J. Chem.
Phys. 79, 926 (1983).
[12] D. E. Smith and L. X. Dang, J. Chem. Phys. 100, 3757 (1994).
[13] A. P. Lyubartsev and A. Laaksonen, Phys. Rev. E 55, 5689 (1997).
[14] B. Hess, C. Holm, and N. van der Vegt, J. Chem. Phys. 124, 164509 (2006).
[15] P. S. Ramanathan and H. L. Friedman, J. Chem. Phys. 54, 1086 (1971).
[16] B. M. Pettitt and P. J. Rossky, J. Chem. Phys. 84, 5836 (1986).
[17] L. X. Dang, B. M. Pettitt, and P. J. Rossky, J. Chem. Phys. 96, 4046 (1992).

20
[18] K. D. Collins, Biophys. J. 65, 65 (1997).
[19] C. J. Fennell, A. Bizjak, V. Vlachy, and K. A. Dill, J. Phys. Chem B 113, 6782 (2009).
[20] J. Ponder and D. Case, Adv. Protein Chem. 66, 27 (2003).
[21] M. Patra and M. Karttunen, J. Comp. Chemistry 25, 678 (2004).
[22] I. S. Joung and T. E. Cheatham III, J. Phys. Chem. B 112, 9020 (2008).
[23] D. Horinek, S. Mamatkulov, and R. Netz, J. Chem. Phys. 130, 124507 (2009).
[24] S. Weerasinghe and P. E. Smith, J. Chem. Phys. 119, 11342 (2003).
[25] B. Hess and N. van der Vegt, Proc. Nat. Acad. Sci. 106, 13296 (2009).
[26] B. Hess, C. Holm, and N. van der Vegt, Phys. Rev. Lett 96, 147801 (2006).
[27] I. Kalcher and J. Dzubiella, J. Chem. Phys. 130, 134507 (2009).
[28] L. Vrbka, W. Kunz, I. Kalcher, R. Netz, J. Dzubiella, and M. Lund, J. Chem. Phys. 131,
154109 (2006).
[29] L. X. Dang, J. Am. Chem. Soc. 117, 6954 (1995).
[30] L. X. Dang, J. Chem. Phys. 96, 6970 (1992).
[31] H. J. C. Berendsen, D. van der Spoel, and R. van Drunen, Comp. Phys. Comm. 91, 43 (1995).
[32] E. Lindahl, B. Hess, and D. van der Spoel, J. Mol. Mod. 7, 306 (2001).
[33] H. J. C. Berendsen, J. P. M. Postma, W. F. van Gunsteren, A. DiNola, and J. R. Haak, J.
Chem. Phys. 81, 3684 (1984).
[34] U. Essmann, L. Perera, M. L. Berkowitz, T. Darden, H. Lee, and L. G. Pedersen, J. Chem.
Phys. 103, 8577 (1995).
[35] J. E. Mayer, Equilibrium Statistical Mechanics (Pergamon, Oxford, 1968).
[36] W. G. McMillan and J. E. Mayer, J. Chem. Phys. 13, 276 (1945).
[37] S. Gavryushov, J. Phys. Chem. B 110, 10888 (2006).
[38] S. Gavryushov and P. Linse, J. Phys. Chem. B 110, 10878 (2006).
[39] J. P. Hansen and I. R. McDonald, Theory of Simple Liquids (Academic Press, London, 1986),
2nd ed.
[40] J. Perkyns and B. Pettitt, J. Chem. Phys 97, 7656 (1992).
[41] P. Kusalik and I. Svishchev, Science 265, 1219 (1994).
[42] L. Ornstein and F. Zernike, Proc. Acad. Sci. (Amsterdam) 17, 793 (1914).
[43] P. Auffinger, T.E. Chetham, and A.C. Vaiana, J. Chem. Theory Comput. 3, 1851 (2007).
[44] A. Heydweiller, Annalen der Physik 335, 873 (1910).

21
[45] W. Hamer and Y.-C. Wu, J. Chem. Phys. Refer. Data 1, 1047 (1972).

22
 1.5
1.5 Cs10 (a) (b)
K11
K10
Cs9 K9 1
εiO (kJ/mol)

εiO (kJ/mol)
1
Cs8 K8
Cs7 K7
Dang Dang
0.5 Cs6 K6 0.5
Cs5 K5
Cs4 Cs3 K4 K3
Cs2 Cs1 K2
K1
0 0
0.28 0.3 0.32 0.34 0.36 0.38 0.4 0.42 0.44
(c) 0.46 0.28F130.3 0.32 0.34 0.36 0.38 0.4 0.42 0.44
(d)0.46
3
σiO (nm) σiO (nm)
I8
1 2.5

F12 2
εiO (kJ/mol)

εiO (kJ/mol)
I7
I6 1.5
F11
0.5 Dang I5 F10
I4 1
F9
I3 I2 Dang F8 0.5
F6 F5
I1 F7 F4 F3 F2
0 F1 0
0.28 0.3 0.32 0.34 0.36 0.38 0.4 0.42 0.44 0.46 0.28 0.3 0.32 0.34 0.36 0.38 0.4 0.42 0.44 0.46
σiO (nm) σiO (nm)

FIG. 1: All Lennard-Jones εiO and σiO parameters for the ions (a) Cs+ , (b) K+ , (c) I− , and (d)
F− used in the optimization procedure in this work. All parameters lie on the curves on which
the experimental single ion solvation free energies are reproduced [23]. The open symbols are the
respective LJ parameters from Dang [29, 30].

23
 25
30 (a) ++
Cs -Cl
-- (b) +
K -Cl
-

20

20 Cs3 K3 15
gij (r)

gij (r)
Cs6 K5
Cs7 K6
Cs9 K8 10
Cs10 K11
10
5

0 0
0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1
(c) r (nm) Na -I
+ - (d) r (nm) +
K -F
-

3 10
I1
I4 K11F4
gij (r)

I5

gij (r)
2 K11F5
I7 K11F7
I8 K11F9
K11F11 5
1

0 0
0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1
r (nm) r (nm)

FIG. 2: Cation-anion radial distribution functions (rdfs) at a concentration of 0.3 M for various
parameter sets for (a) CsCl, (b) KCl, (c) NaI, and (d) KF. The Cl− and Na+ force fields are taken
from Dang [29, 30], the K+ force field in (d) is K11.

24
 8
+ -
4 Cs -Cl
0
6
Cs6
βVij(r)

+ + Cs9
3 Cs -Cs
0
6
- -
3 Cl -Cl

0
0.2 0.4 0.6 0.8 1 1.2
r (nm)

FIG. 3: Radial distribution functions (rdfs) at a concentration of 0.3 M for the Cs6 parameter set.
Red (solid) lines correspond to the MD results and blue (dashed) lines to the rdfs as generated
through the HNC approach (based on the pmfs from MD). From top to the bottom, the rdfs for
the cation-anion, cation-cation, and anion-anion are shown.

25
 8 (a) 8 (b)
+ - + -
4 Cs -Cl 4 K -Cl
0 0
6 6
Cs6 K5

βVij (r)
βVij (r)

+ + + + K11
3 Cs -Cs Cs9 3 K -K

sr
sr

Cs3 K3
0 0
6 6
- - - -
4 Cl -Cl 4 Cl -Cl
2 2
0 0

8 (c) 8 (d)
+ - + -
4 Na -I 4 K -F
0 0
6 6
I1 K11F5
βVij (r)

βVij (r)

+ + I4 + +
3 Na -Na 3 K -K K11F3
sr

sr

I7
0 0
6 6
- - - -
4 I -I 4 F -F
2 2
0 0
0.2 0.4 0.6 0.8 1 1.2 0.2 0.4 0.6 0.8 1 1.2
r (nm) r (nm)

FIG. 4: Short-ranged potentials of mean force (pmfs) at a concentration of 0.3 M for the optimal
LJ parameters (see Section III C) for (a) Cs (sets Cs6, Cs9), (b) K (sets K5, K11), (c) I (sets I1,
I4) and (d) KF (set K11F5). In each panel, from top to the bottom, the pmfs for the cation-anion,
cation-cation, and anion-anion are shown. For comparison, in all panels the pmfs for a non-optimal
parameter set is included: (a) set Cs3, (b) K3, (c) I7, and (d) K11F3.

26
 1.1 CsCl
KCl
Cs7
1 K6 Cs8 K10 Cs10
K5
φ 0.9 KCl K11
CsCl K4 K8 Cs9

K3 Cs6
0.8
Cs3
Cs4 Cs5
0.7 (a)
0 0.5 1 1.5
1.2 I5 εI7iO(kJ/mol)
I1 I3 I8
1.1 I6
I2 I4 NaI
1
NaI KF
0.9
K11F
φ 0.8 F8 K5F
F9
0.7 F5 F6
F7
0.6
F11
0.5 F3
F10
0.4 F4
(b)
0.3
0 0.5 1 1.5
εiO(kJ/mol)

FIG. 5: The osmotic coefficient φ as calculated from the MD simulations by means of the virial route
for CsCl, KCl, NaI and KF in a 1 M solution. In (a) the magenta(dashed) and turquoise(dotted)
lines represent the experimental results for CsCl and KCl, respectively. In (b) the green (dotted)
and orange(dashed) lines are the experimental values for NaI and KF, respectively. Representative
results for KF for the K5 parameter set are also shown. The open symbols are the osmotic
coefficients for the Dang parameters obtained in a previous work [27, 29, 30]. The solid lines are
guides to the eye.

27
 1.1 1.1
(a) (b)

1 1

Cs4 K3 0.9
0.9 Cs5 K4
φ HNC

φ HNC
Cs6 K5 0.8
Cs7 K6
0.8 Cs8 K8
Cs9 K10 0.7
Cs10 K11
0.7 exp exp
K5 - MD 0.6
Cs6 - MD
Cs9 - MD K11 - MD
0.6 0.5
0 0.5 1 1.5 2 0 0.5 1 1.5 2
1.3 Concentration (M) Concentration (M) 1.1
(c) (d)
1
1.2 F3
I1 F4 0.9
I2 F5
F6
φ HNC

φ HNC
I3
1.1 I4 F7 0.8
I5 F8
I6 F10
I7 exp 0.7
1 I8 F5 - MD
exp 0.6
I1 - MD
I4 - MD
0.9 0.5
0 0.5 1 1.5 2 0 0.5 1 1.5 2
Concentration (M) Concentration (M)

FIG. 6: The osmotic coefficient φ as calculated through the HNC approach (lines, see text) of
representative LJ parameter sets for (a) CsCl, (b) KCl, (c) NaI, and (d) KF as a function of
concentration (M). The experimental curves [44, 45] are also shown. The symbols are the values
directly from MD for the parameters sets proposed in Table III. In panel (d) all KF results are
obtained by using the K11 parameter set for the potassium ion.

28
 (a) (b)
Dang - MD
1.1 kbff - MD 1.1
Dang - HNC
1 kbff - HNC
exp
φ 0.9 φ
1
0.8
Cs6I4 - MD
0.7 Cs9I4 - MD
Cs6I4 - HNC
Cs9I4 - HNC
0.6 exp
0.9
0 0.5 1 1.5 2 0 0.5 1 1.5 2
Concentration (M) Concentration (M)

FIG. 7: MD results (symbols) for the osmotic coefficient, plotted together with the HNC results
(broken lines) as a function of the salt concentration (in M). Panel (a) shows results for CsI
and serves as a check on the transferability of the proposed force fields for Cs+ and I− . Panel
(b) depicts the results for NaCl from the Kirkwood-Buff derived force field (kbff) [24] and the
Dang parameters [27, 29, 30]. In both panels, the black (solid) lines correspond to the respective
experimental curves [44, 45].

29