Transferable Ion Force Fields in Water from a Simultaneous Optimization of Ion Solvation and Ion–Ion Interaction

The poor performance of many existing nonpolarizable ion force fields is typically blamed on either the lack of explicit polarizability, the absence of charge transfer, or the use of unreduced Coulomb interactions. However, this analysis disregards the large and mostly unexplored parameter range offered by the Lennard-Jones potential. We use a global optimization procedure to develop water-model-transferable force fields for the ions K+, Na+, Cl–, and Br– in the complete parameter space of all Lennard-Jones interactions using standard mixing rules. No extra-thermodynamic assumption is necessary for the simultaneous optimization of the four ion pairs. After an optimization with respect to the experimental solvation free energy and activity, the force fields reproduce the concentration-dependent density, ionic conductivity, and dielectric constant with high accuracy. The force field is fully transferable between simple point charge/extended and transferable intermolecular potential water models. Our results show that a thermodynamically consistent force field for these ions needs only Lennard-Jones and standard Coulomb interactions.


Contents
S5 Fitting the MSD S11 S6 Details on calculating the dielectric decrement S11 S7 Observables for different water models S12 S8 Comparison to other force fields S12 References S14 S1 Optimization strategy The goal of the optimization strategy is to find the chloride parameters in a way that reproduces the experimental solvation free energy of every cation/anion pair, and provides the best fit to the logarithmic activity derivative a cc as a function of concentration. Therefore, we start by selecting one set of parameters for the chloride ion and calculating the free energy isolines of sodium and potassium. This gives one isoline for sodium and one isoline for potassium for every choice of chloride. The details of the isoline calculation are discussed in Section S3. For a number of points along these isolines, we calculate the logarithmic activity derivative of NaCl and KCl as a function of concentration. To avoid having to calculate a cc curves of NaBr and KBr for every possible choice of sodium and potassium parameters, we first perform a partial optimization of the sodium and potassium parameters for every choice of chloride parameters. To that end, we calculate the partial square deviation k (defined identically to Eq. (11) from the main text) for each of the a cc curves, by summing the square deviation over all concentrations used, but only over the salts NaCl and KCl. This gives us one set of best-performing sodium and potassium parameters for every initial choice of chloride parameters. For these best-performing sodium and potassium parameters, we calculate the logarithmic activity derivative of NaBr and KBr as a function of concentration for a number of points along the NaBr and KBr free energy isolines. We now optimize the bromide parameters by calculating the square deviation k (Eq. (11) from the main text), summing over the concentrations used and over all four salt combinations. This gives us one value of k for every choice of chloride parameters, corresponding to the total square deviation from the experimental activity derivatives when using the best possible sets of sodium, potassium, and bromide parameters. We find the final optimized force field by repeating the scheme above for different choices of the initial chloride parameters.

S2 Experimental solvation free energy
The solvation free energy of ion pairs is measured experimentally at a temperature of 298.15 K. A summary of experimental values of the enthalpy ∆H, entropy ∆S and Gibbs free energy ∆G is given in Table S1, in chronological order.   Some data sets listed in Table S1 report two independent values and calculate the third, making them thermodynamically consistent automatically, whereas some report three independent values. These data sets are typically not exactly thermodynamically consistent, but the differences are rarely more than 1 k B T . Exceptions include early estimates of the Br free energy and, notably, the book by Marcus 10 . One possible reason for the inconsistency in Marcus' data is that at least one of his papers reports the 'solvation energy proper', referring to the transfer from the 1 atm gas phase to a single ion in solution, using the wrong sign for the pressure correction 15 . The affected data set also appears to be thermodynamically inconsistent, even though according to the text ∆G is calculated from ∆H and ∆S. This inconsistency might well be due to the sign error, because the deviations are of the order of the value of the pressure correction.

S4
In his book, however, Marcus uses standard states, so the pressure correction is not to blame 10 . It is clear that the enthalpy data is taken from the NBS tables compiled by Wagman et al. 7 , shifted by z(433.2) kJ/mol, where z is the ion valency. Taking the entropy change for the transition of the uncharged species in gas to the charged species in solution also from the NBS tables yields values for ∆G that are very close to the listed ∆G values 10  Interestingly, using ∆H and ∆S from Marcus 10 , which corresponds to the transition of the charged species in the gas phase to the charged species in the hypothetical ideal 1 mol/l liquid phase, gives almost exactly the same values as the ones found independently by Tissandier et al. 11 . A subsequent estimate by Fawcett 12 using the same approach for the electronic contribution as Marcus and Loewenschuss 8 , as well as estimates by Schmid et al. 13 and Kelly et al. 14 give very similar results as well. Therefore, we believe that the Tissandier/Marcus (∆H, ∆S) data is probably the most accurate data set available.
Since our simulations -like most molecular dynamics simulations -are performed at a temperature of 300 K, we calculate the solvation free energy ∆G at 300 K from the entropy change ∆S and enthalpy change ∆H. Because the volume change during solvation of an ion pair is negligible 16 , we can compare the simulated solvation free energy F directly to the experimental solvation free energy ∆H − T ∆S. Averaging the values from Marcus 10 and Tissandier et al. 11 gives the final values listed in Table S2.  [16][17][18] . These values were also used by Mamatkulov and Schwierz 19 . Given the relatively small deviations, however, these different reference values are generally not expected to have had a major impact on the final force field parameters in those papers. The calculation of the derivatives ∂H/∂λ of single ions is described in the main text. The single-ion free energy is obtained using Simpson's integration rule. The numerical integration is verified and found to be accurate by comparing to the energy obtained using a 12-point open Gaussian integration scheme. From a thermodynamic point of view, the only relevant quantity is the solvation free energy of ion pairs. The free energy isoline is defined as the line in parameter space where the simulated solvation free energy of the pair reproduces the experimental value. To obtain the free energy isolines, we calculate the free energy for several single anions and cations. Note that when solvating ion pairs from the vacuum into the bulk water, the electrostatic energy originating in bringing a single ion across the water interface vanishes, and therefore we do not take it into account. Figure S1 shows the resulting solvation free energy for the SPC/E water model. Black dots show the grid on which the simulations where performed and the values between the black dots are obtained using cubic interpolation. Since water is charge-asymmetric it favors anionic versus cationic solvation [20][21][22][23] . It is also clear from Fig. S1 that solvation of ions with a smaller σ are favored as predicted by the Born model 24 .  Now the solvation free energies of ion pairs are compared to the experiment. We fix the parameters of one ion of the pair, and calculate the ion pair solvation free energy as a function of the parameters of the other ion. As an example, Figure S2a shows the solvation S7 free energy F of NaBr, reduced by the experimental solvation free energy (∆H − T ∆S) exp . Here the sodium parameters are fixed and the bromide parameters are varied. The solid lines are fits according to

S3 Free energy isolines
where a, b and c are fit parameters. Different colors denote different values of σ Br . From the roots of the fits (shown as blue crosses), we obtain the isoline, i.e. pairs of σ Br and Br which reproduce the experimental solvation free energy for fixed values of σ Na and Na . We now fit the isoline according to where d, e, and f are again fit parameters. The isoline and the fit are shown in Figure S2b and the best fit parameters are shown in table S3. Along these isolines, we calculate the activity coefficients. In Figure S3 we show the activity derivatives a cc for several parameter sets using the SPC/E water model. Note that all parameters shown reproduce the experimental solvation free energy as they are all drawn from isolines. However, only a few also agree with the experimental activities. The best parameter sets are shown as solid dots. As mentioned in the main text we find that especially the sodium chloride activity is sensitive to changes in the Lennard Jones parameters.
In Fig. S4, we show the estimates of the error of the individual values of a cc for the optimal force field parameters.

S4 Radial Distribution functions
We show the radial distribution functions of the optimized force field at two different concentrations in Fig. S5.

S6 Details on calculating the dielectric decrement
We follow the procedure detailed in Ref. 25 to obtain the dielectric spectra from each salt simulation. This procedure begins with the linear response relationship between the electric susceptibility χ(ω) and the equilibrium auto-correlation function for fluctuations of the polarization M (t), For our system the total instantaneous polarization is the sum of two contributions M (t) = M W (t) + M I (t) for the water molecules and ions, respectively. However, in a periodic simulation cell, it is convenient to express the ionic polarization contribution in terms of the ionic current J = dM I /dt due to toroidal shifts that occur when charged particles traverse the simulation boundaries 26 . Manipulating the above expression, we can express χ(ω) in terms of time correlations between these two contributions to the total polarization, and arrive at The first term χ W (ω) is due to the water dipole moment correlations and is given by, where we have defined the water dipole time correlation function as The second term χ IW (ω) arises from cross-correlations between the water dipoles and salt ion currents. This term can ultimately be expressed as with the ion-water cross-correlation function defined as The final term χ I (ω) contributing to the total susceptibility stems from the ion current S11 auto-correlation function. However, this term diverges at zero frequency due to the finite DC-conductivity of ionic solutions 25 . Thus, we calculate a regularized susceptibility given by, with the ion current auto-correlation function defined as Plots of the real (χ ) and imaginary (χ ) parts of the total frequency-dependent susceptibility are shown for each ion pair for different concentrations in Fig. S7. As noted in the main text, the dielectric constant used in calculating the dielectric decrement, denoted ∆ε, is inferred from the lowest frequency point of the real contribution to the total susceptibility for each salt type at every concentration.

S7 Observables for different water models
The mass density, ionic conductivity, and dielectric decrement for water models other than SPC/E are shown in Fig. S8. The values of the water diffusion constant and the dielectric constant at zero ion concentration of the different water models are listed in Table S5.

S8 Comparison to other force fields
The optimal chloride parameters obtained from our optimization are close to the Dang chloride used earlier 17,43 . However, the solvation free energy of the Dang NaCl pair is off by several k B T , so the present optimization provides a significantly better agreement. The   [30][31][32][33][34][35][36][37] . The zero concentration dielectric constants for the different water models are ε TIP3P = 94, ε TIP4P = 50 and ε TIP4P/ε = 80 38,39 . Water self-diffusion constant, normalized by the diffusion constant of pure water, together with the experimental data 40 . The pure water diffusion constant D 0 of the different water models is listed in Table S5. 3.77 ± 0.09 (3.31) 53.9 (50) TIP4P/ε 2.23 ± 0.01 (2.1) 79.7 (80) solvation free energy of the CHARMM and Mamatkulov and Schwierz 19 ion pairs is closer to the experimental value, but still off by several k B T . Because the force field by Mamatkulov and Schwierz 19 is optimized at 300 K with respect to the solvation free energy from Marcus 10 at 298.15 K, a large part of this deviation is caused by the slight differences between the experimental reference values and the temperature difference. This deviation is probably not the primary problem, however. Instead, the strategy employed by Mamatkulov and Schwierz 19 is to pick parameters for the chloride ion in TIP3P which reproduce the singleion solvation free energy of the Dang chloride in SPC/E water. In focussing on singleion solvation free energies, their optimization strategy is fundamentally different from our approach. In particular, in contrast to the solvation free energy of ion pairs, single-ion solvation free energies do sensitively depend on the water force field. Relatedly, interface potential differences, which are important for single-ion solvation free energies, also depend on the water force field. Overall, these issues lead to a poor choice of the chloride parameters and the spurious conclusion that ion force fields are not transferable between SPC/E and TIP3P.