salts of monovalent ions are highly soluble in aqueous solution , and the solvated ions have important functional roles in screening charge , moving charge , and also influencing the structure and dynamics of biomolecules such as proteins and nucleic acids .
given this , accurate modeling of biomolecular structure , dynamics , and function requires some representative model of ionic interactions and mobile ions . in this work ,
we focus on the most common molecular mechanical models , specifically those based on a pairwise additive potential that includes a coulombic treatment of the electrostatic interactions and a lennard - jones representation of dispersionattraction and core repulsion . in this formulation ,
the potential energy ( uij ) between any pair of nonbonded atoms ( i and j ) in a system composed of the ions and water molecules is usually expressed as 

 here , qi and qj are the point charges of the atoms , rij is the distance between atoms , and rmin , ij and ij are the van der waals radius and well depth of the lennard - jones potential . to represent the water , simplified point charge models such as tip3p,(8 ) spc / e,(9 ) and tip4pew(10 ) are utilized .
although the models are simple , determination of the lennard - jones parameters for the ions is fairly subtle , and therefore , a wide variety of different parameter sets have been proposed .
many of these are summarized in our previous work that introduced a new set of force field parameters for the alkali ( li , na , k , rb , cs ) and halide ions ( f , cl , br , i ) where the specific choice of lennard - jones parameters depends on the specific choice of water model.(11 ) here we further investigate the structural and dynamic properties of these ions in comparison with other existing parameter sets or force fields .
new properties calculated include activity coefficients , diffusion coefficients , populations of ion clusters , residence times of ionion and ionwater dimers , association constants , and solubility .
although many of these properties have been determined experimentally , some properties ( such as ionwater residence times and ion clusters populations ) are elusive or disputed and can only be estimated through simulation .
we compare our results with existing experiment and previous simulations , noting that much of the prior simulation work is fragmental as it tends to focus on calculating the properties themselves rather than assessing the influence of the force field .
the intent here is to provide insight into the validity and range of applicability of simple pairwise atomic molecular mechanical models of alkali and halide ions .
this information can help guide the development of more accurate models of ions in aqueous solution .
assuming there are only solute molecules a and solvent molecules b in a molecular system , the chemical potential a , is as follows 

 where xa , a , k , and t refer to the mole fraction , activity coefficient of molecule a , boltzmann constant , and temperature , respectively.(20 ) the chemical potential in the standard state composed of pure a molecules is a0 .
equation 2 can be rewritten in more convenient molality units ( ma [ mol / kgsolvent ] ) instead of mole fraction(21 ) 
 where m0 = 1 mol / kg and the new activity coefficient in terms of molality ( a , m ) is related to the original activity coefficient through the relation , a , m = axb .
contributions to the activity term ( ktln a , m ) can be estimated by subtracting chemical potentials calculated at different concentrations .
this is not trivial since , as the concentration rises ( up to 1 molality ( m ) ) , the experimentally measured activity coefficients for nacl and kcl drop down to 0.6 , which is equivalent to a very small energetic change in the chemical potential ( 0.30 kcal / mol at 300 k ) . at low concentrations , due to the sensitivity of the chemical potential to changes in the activity coefficient , only very small errors in the calculated chemical potential are tolerable .
as typically errors in calculating free energies ( or the chemical potential ) in conventional simulations are around 1 kcal / mol in most practical applications,(22 ) reducing the error to the order of 0.01 kcal / mol essentially requires an enormous number of ensembles or , equivalently , much longer simulations . to overcome this problem , previous investigations estimated activities using the potentials of mean force acting on the counterions .
the benefit of this approach is that the explicit water is no longer necessary once the potentials of mean force are constructed .
to more fully simulate the relevant ensemble , we chose instead to use fully explicit simulations following kokubo and co - workers.(21 ) chemical potentials can be calculated by applying the widom method(23 ) as expanded into ntp ensembles by shing and chung.(24 ) when a test particle of a is inserted in the system with na  1 particles of a(solute ) and nb particles of b(solvent ) , the chemical potential of a can be expressed as 

 where a and qa are thermal de broglie wavelength and the molecular partition function of particle a ; v is the volume of the system ;  = 1/kt ; and  is the potential energy of the test particle interacting with the rest of the pre - existing atoms .
equation 4 includes two terms as shown below , and the excess term can be simplified assuming the correlation between energy and volume is weak . 


the excess term is equivalent to the solvation free energy of the test particle , which can be determined by several methods . in this research
, we chose to use thermodynamic integration ( ti ) . in order to calculate activity coefficients , equate eq 3 and eq 4 and rearrange to obtain 

 activity coefficients go to 1 as the molality goes to zero
thus , the reference state of the most diluted solution (  ) can be introduced 

 subtraction of eq 7 from eq 8 gives 

 the calculation of each term in eq 9 is straightforward , except for the last term , which is the difference of the solvation free energy of the inserted particle at a specified concentration from that of the most diluted solution . as mentioned ,
previous attempts to calculate mean activity coefficients of alkalihalide salts were limited to implicit solvent methods using the potential mean force between molecules .
lyubartsev and laaksonen derived the potential mean force using smith and dang s nacl(27 ) in flexible spc water,(28 ) while gavryushov and linse used spc / e water and aqvist(29 ) and dang s parameters .
lenart(14 ) used koneshan and rasaiah parameters(30 ) for van der waals interactions of ions , and the potentials in water were adjusted for a specifically designed implicit solvent model .
experimental methods to measure mean activity coefficients are various , but the results are all very similar.(31 ) we used the values summarized the crc handbook(32 ) for comparison .
the more intuitive method is direct calculation of the average mean square displacement of ions for a unit time period .
calculation of the velocity autocorrelation function also provides a link to the diffusion coefficient.(15 ) however , given long simulations , the two methods give equivalent results for the diffusion coefficients . in this research
information regarding the dynamics of ions in water is very limited and relies mostly on simulation .
for example , the correlation time ( or the residence time ) of ionwater and ionion pairs can easily be calculated by simulation .
previous estimates of the residence times applied the autocorrelation function designed by impey et al.(15 ) more recently , laage and hynes raised the issue that residence times calculated by that method are very sensitive to the tolerance time.(33 ) the tolerance time
was originally introduced to the calculation in order to represent the transition state between the paired ( reactant ) and unpaired ( product ) states as the reactant state moves to become the product state . in other words , this catches the case when a pair only briefly becomes unpaired before pairing again , with a time period less than the tolerance time .
the updated method by laage and hynes suggests calculating the new time correlation function ( r(t ) ) based on the stable states picture model.(34 ) to apply this method , the states of stable reactant , r , ( or paired ) and stable product , p , ( or unpaired ) need to be defined . in the r state ,
molecule pairs are separated by almost the radius of the first shell ( rmax1 ) , which is the first peak of their radial distribution function ( rdf , g ) .
the two molecules in the pair reaches the p state when they are separated by almost the radius of the second shell ( rmax2 ) , which is the second peak of the rdf .
the probability , p(t ) , is the probability that a molecule pair in r state at time 0 transforms into the p state after time t : 

 the residence time (  ) of the molecule pair is defined by the equation above .
although experimental measurements of the solubility of alkalihalide crystals in water are easily found , it is not simple to calculate solubility in simulation . to compare with experiment ,
to do this , the chemical potential of the crystal and solutions at various concentrations is calculated and the concentration at which the two potentials are equal indicates saturation concentration . instead of this more tedious approach , we used a more direct method to derive the equilibrium state .
sufficiently large crystals in contact with an almost saturated ion solution will reach an equilibrium state given sufficient sampling . to avoid finite size surface effects of a small crystal ,
we note that this approach , although direct , will likely require longer simulations to fully equilibrate the heterogeneous mixture .
however , the direct method is advantageous since there are no underlying assumptions ( as are required in the calculation of chemical potentials ) . as a result
molecular dynamics simulations were performed with the sander module of amber9 using the amber nonpolarizable potential energy function shown in eq 1 .
general conditions applied in all of the molecular dynamics simulations , unless otherwise stated , are discussed below .
a 2  buffer was built into the pairlist , and rebuilding of the list was triggered if any particle moved more than 1.0  .
long - ranged electrostatic interactions were estimated with the particle mesh ewald method , and a homogeneous density approximation was applied to correct for the cutoff of the lennard - jones interactions .
particle grids for the ewald summation were approximately one grid point per angstrom in each dimension .
the atomic charges for the grid points were smoothed onto the grid using a fourth - order b - spline .
temperature and pressure were regulated by using weak coupling algorithm.(41 ) the coupling time for temperature and pressure control was 1 ps . the target temperature and pressure were 298 k and 1 bar , respectively .
all of the water models used here are rigid , and the shake algorithm(42 ) was applied to constrain all bonds to hydrogen with a tolerance of 10  .
the initial coordinates of water and ions were randomly positioned in a rectangular periodic box , which was prepared to have density of 1 g / cm .
the initial box size was about 3840  , and the number of waters and ions in each case is displayed in table 1 .
the systems were minimized by steepest descent for 1000 steps , followed by a two - step equilibration . in the first step ,
md simulation was performed for 40 ps in the nvt ensemble . in the second step ,
a cationanion pair was inserted in the system slowly to calculate the hydration free energy of the ion pair .
lennard - jones interactions of the inserted ion pair were turned on . in the second step ,
the charge of the anion was activated . the perturbation steps adopted the following mixing rule of the potential energy : u( ) = f()u0 + ( 1  f())u1 , where u0 and u1 are the potential energy of the initial state and the final state .
intermediate states between perturbations were determined adjusting the parameter ,  ( 0    1 ) .
for the first perturbation , f( ) = ( 1  )i = 06(6 + i)!/6!/i ! , and for the second and third perturbations , f( ) = 1   was applied , respectively .
numerous intermediate states between perturbations were defined by substituting  with these numbers : 0.01 , 0.05 , 0.1 , 0.15 , 0.2 , 0.25 , 0.3 , 0.35 , 0.4 , 0.45 , 0.5 , 0.55 , 0.6 , 0.65 , 0.7 , 0.75 , 0.8 , 0.85 , 0.9 , and 0.95 for the first perturbation and 0 , 0.05 , 0.1 , 0.15 , 0.2 , 0.25 , 0.3 , 0.35 , 0.4 , 0.45 , 0.5 , 0.55 , 0.6 , 0.65 , 0.7 , 0.75 , 0.8 , 0.85 , 0.9 , 0.95 , and 1 for the second and the third perturbations . using the previously equilibrated coordinates ,
all of the states were simulated in the ntp ensemble to calculate < du / d>. the initial 100 ps for each of the simulations was discarded .
the < du / d > values were fit into a cubic spline as a function of  , and this was used to calculate 01du / dd numerically , which is equivalent to eq 6 . shown
are the molalities ( m ) of kcl and nacl solutions simulated when calculating the activity coefficients , including the number of explicit ion pairs ( nion ) and water molecules ( nw ) in the periodic unit cell .
the initial coordinates of the water and ions were randomly placed in a rectangular periodic box keeping the density of the system to be 1 g / cm .
fifteen hundred explicit water molecules and 27 ion pairs were added to make the concentration 1 m. the box sizes were approximately 3640  .
note that lif was omitted since its solubility is lower than 1 m.(32 ) each system was minimized and equilibrated in the same way as described in the previous section , and then md simulations were performed for 10 ns .
each curve was smoothed using a bzier curve to better obtain the characteristic radii of the rdf curves at dg / dr = 0 .
the three shortest radii were determined , and they were assigned sequentially as the radius of the first shell ( rmax1 ; radius of the first peak ) , the radius of first coordination shell ( rmin1 ; radius of the first minimum after rmax1 ) , and the radius of the second shell ( rmax2 ; radius of the second peak ) .
the first shell was refined by repeatedly fitting the curve into quadratic equations.(11 ) the radius at the vertex of the quadratic fit was substituted for rmax1 .
the relationship , 6d = limtd < r>/dt was utilized to calculate the diffusion coefficient ( d).(44 ) average mean - square displacement ( < r > ) of each type of molecule up to 10 ps was plotted versus time every 1 ps .
the curves were fit to a straight line that passes through the origin by least - squares .
the diffusion coefficients were then corrected for finite periodic box size effects as per yeh and hummer.(45 ) this was done by measuring the diffusion coefficients of ions with various sizes of the system , specifically 1500 water molecules and 27 ion pairs , 1000 water molecules and 18 ion pairs , and 500 water molecules and 9 ion pairs .
these simulations were extended to 10 ns each . for the analysis of ion clusters ,
cluster populations were analyzed from the simulation of 1500 water molecules and 27 ion pairs .
calculation of the residence time for a paired molecule requires the precise determination of the r and p states , which are based on the distance of the two molecules .
the upper boundary of the r state ( rr ) was the radius satisfying g(rr ) = (g(rmax1)g(rmin1 ) ) .
a pair of molecules at the same or shorter distance than rr was considered to be in the r state . in a similar way
, the lower boundary of the p state ( rp ) was determined by the radius satisfying g(rp ) = (g(rmin1)g(rmax2 ) ) .
a pair of molecules separated by a distance equal or larger than rp were defined to be in the p state .
originally , the lower boundary of the p state suggested by laage and hynes was the radius corresponding to the half - energy barrier in mean force potential curve . however , our definition of the boundary is equivalent to their definition because radial distribution is exponentially proportional to the potential of mean force ( i.e. , g(r )  e , where w(r ) is the potential of mean force ) .
the residence time ,  , was calculated by integrating the time correlation function eq 10 . 

 because r(t ) can not be evaluated infinitely , the equation was integrated only from 0 to 1 ns and the rest of the integration was carried out under the assumption of exponential decay : 

 integration of r(t ) over t was performed using the trapezoidal rule , and eq 12 was numerically solved to obtain . rectangular crystals of alkalihalide salts were built using their lattice constants as described previously.(11 ) when building the nacl - type crystals , the unit crystal was aligned exactly parallel to the rectangular axes ( figure 1a ) . on either side of the yz planes of the crystals ,
an ion and water mixture with the density of 1.52.0 g / cm was prepared .
presaturated ion solutions were composed of roughly 20% more ion pairs than expected to be soluble in terms of molal concentration ; 862 ( licl ) , 266 ( nacl ) , 206 ( kcl ) , 336 ( rbcl ) , 901 ( libr ) , 397 ( nabr ) , 246 ( kbr ) , 304 ( rbbr ) , 534 ( lii ) , 531 ( nai ) , 386 ( ki ) , 336 ( rbi ) , 1626 ( csf ) ion pairs and 2000 water molecules . for the cscl - type crystals ( cscl , csbr , and csi )
, the unit crystal was initially aligned to the rectangular axes then rotated 45 along the z axis in order to reduce the net dipole along the x axis ( figure 1b ) .
this prevented the strong net dipole from artificially breaking the crystal immediately after the start of the simulation .
on either one of the yz planes of the crystals , an ion and water mixture with the density of 1.52.0 g / cm was apposed ; 120% saturated solutions and various other solutions at different concentrations were tried .
however , we encountered several problems with cscl - type crystals as explained in the next section . coexisting salt crystal and ion solution in a periodic box .
the horizontal axis is the x axis , and ion solution was placed on one of the yz planes of the crystals .
( a ) nacl - type crystal and ( b ) cscl - type crystal .
the equilibration and subsequent dynamics were normally not very sensitive to the precise setup of the crystal solution and its initial density .
however , in a few cases , the crystal instantaneously disintegrated and melted into the aqueous solution due to instability in the built crystal ; in these cases , the results were discarded . as the simulation cell is periodic ,
the initial setup constituted infinite crystal sheets with aqueous ion solution between the crystal sheets .
the system was minimized by using the steepest descent method for 1000 steps and equilibrated at a constant temperature and pressure for 40 ps while harmonically restraining the coordinates of crystal ions to their initial positions with a force constant of 50 kcal / mol /  .
pressure was controlled by adjusting the periodic box size anisotropically so as not to break the connectivity of the crystals between neighboring boxes .
that is , pressure applied to x , y , and z directions were calculated separately , and the three edges of the box were rescaled according to the three normal pressures .
anisotropic rescaling intrinsically requires longer equilibration time , but it is indispensible in nonisotropic systems like this .
the coupling time for the pressure regulation was 1 ps . after removing the restraints , md simulation was continuously performed and the crystals gradually grew or dissolved as they moved toward the equilibrium state .
the distribution of the water and ions along the x axis changes as the crystal grows or is dissolved . to obtain the distribution of molecules at different times in a consistent manner , a static reference needed to be defined . on the time scale of the md simulations performed , the inner crystal layers remained intact .
therefore , we used the geometric center of the inner crystal layers as the reference point for the x axis distribution . from the reference point , the histogram of the molecules along the x axis was collected up to half of the x axis edge of the periodic box ( xmax ) with a step of 0.1  .
this final histogram was converted to number density per unit volume by dividing by the volume of each bin of the histogram .
the bulk molar concentration of the molecules in the solution was assumed to be the average molar concentration of the molecules located between xmax  10  and xmax from the reference point . the volume and total potential energy of the system were traced and averaged .
the system was considered to reach the equilibrium state when both the volume and potential energy leveled off . to avoid introducing biases
when determining the equilibrium state , running averages over 60 ns of the two properties were recorded every 100 ps along with their standard deviation .
the ranges that these averages cover with the tolerance of their standard deviations were compared each other .
when the average overlap of the ranges from the last 40 ns onto the very last range is above 95% , we considered the system to be in equilibrium state ( see figure s1 in the supporting information ) .
the final concentration of the solution phase was considered to be the average concentration of the last 30 ns .
activity coefficients of nacl and kcl solutions at varying concentrations ( 5 m ) were calculated with tip3p water(8 ) and the tip3p - compatible ions .
only this limited set was investigated as these simulations were very costly and effectively required the equivalent computer time of 11  700 days  4 processors based on the performance of available clusters . yet
, these results still provide general insight into how much the activity coefficients from simulations vary with respect to experiment . to calculate activity coefficients ,
the number of explicit atoms in the periodic box was adjusted to simulate a variety of concentrations , as is shown in table 1 . in eq 9 ,
all of the other parameters were determined ( see table 1 ) except for the excess chemical potential .
as the activities are very sensitive to small changes in the excess chemical potential terms , it is necessary to discuss how the error was estimated . in general , the errors of < du / d>i from the ith ti window can be estimated if the correlation length ( c ) of < du / d>i is known.(46 ) given this , we denoted the error as ew , i , and it was estimated using the relationship 

 where ( rmsd)i is the rms deviation of < du / d>i and ti is the length of the simulation required to obtain the rmsd . although we did not precisely estimate the correlation length ( c ) of < du / d>i for all of the ti windows , we expect that the correlation lengths are shorter than 1 ps ( for example , please refer to the decay of the autocorrelation curve shown in figure s2 in the supporting information ) .
when we assume that the correlation lengths are less than 1 ps , the error estimates should be valid .
as the simulations proceeded , it was noted that some windows had inherently low errors , dropping below 0.00001 kcal / mol within a few hundreds of picoseconds of simulation , whereas other windows required tens of nanoseconds before the error drops below 0.1 kcal / mol .
the contribution of each window to the free energy change of a perturbation is roughly inversely proportional to the total number of windows ( nw ) of the perturbation because the ti windows are almost evenly distributed .
therefore , an error in free energy calculation of each perturbation ( ep , j ) is 

 where the subindex j denotes jth perturbation .
finally , the error of the total free energy calculation ( etotal ) can be estimated as 

 we differentiated the time of simulations for the higher concentrations ( > 1 m ) and the lower concentrations ( < 1 m ) because we wanted higher precision at the lower concentrations .
the finally estimated total errors fell around 0.082 kcal / mol for the higher concentrations and 0.034 kcal / mol for the lower concentrations .
however , real errors might be even smaller because we assumed that c was 1 ps , which is likely an overestimate .
the calculated excess chemical potentials are listed in table s1 and displayed in figure 2 . in order to calculate activity coefficients of ion pairs applying eq 9 , the difference in chemical potential from the most diluted solution
due to the practical difficulty in carrying out simulations of high volume systems , with explicit solvent , we did not simulate solutions at less than 0.037 m. in fact , this is not sufficiently dilute , and we could not assume that the activity coefficient of the solution is unity at this concentration . when they are compared to experimental activity coefficients,(32 ) the chemical potential at this concentration is almost 0.1 kcal / mol more negative than that from the most diluted solution ( i.e. ,  0.1 kcal / mol  ktln m ) .
, we assumed that there is 0.1 kcal / mol difference in chemical potential at 0.037 m from the most diluted solution .
in actuality , the values of the activity coefficients at the other concentrations are greatly affected by this number . however , choosing this number is reasonable , and regardless of the precise value , the overall shape of the curve of activity coefficients at various concentrations is preserved .
the activity coefficients calculated with this assumption are listed in tables 2 and s1 . calculated ( calc ) and experimental(32 ) ( ref ) mean activity coefficients of nacl and kcl versus square root of molality .
the x axis denotes the square root of molality of an ion pair , and the y axis indicates its mean activity coefficients .
closed circles and closed squares denote the activity coefficients of tip3p - compatible ion pairs , nacl and kcl , in tip3p water .
open circle and open squares are the experimental mean activity coefficients of nacl and kcl in water .
the contributions to the chemical potential from the three perturbations are indicated as vdw , elec1 , and elec2 ( table s1 in the supporting information ) .
the vdw generally increased as the concentration of the solution increased for both nacl and kcl .
for elec1 , the insertion of sodium ions was less favored as concentration increased ; however , it became more favored for potassium ions beyond 1 m. overall , the contributions to the differences in the chemical potentials observed were relatively small ( columns 2 , 3 , and 4 in table s1 in the supporting information ) with the exception of the hydration free energies .
the calculated activity coefficients deviate from the experimental values as the concentration increases ( table 2 ) .
however , the curve of nacl was placed above that of kcl , which qualitatively corresponds to the experimental curves .
although we will discuss solubility further below , the solubility of nacl and kcl determined by simulations in tip3p water was actually lower than their experimental solubility .
however , we did not observe crystal formation on the time scale of these simulations , suggesting the solutions were supersaturated above certain concentrations .
lyubartsev and laaksonen calculated the activity coefficients of smith and dang s nacl(27 ) in flexible spc water.(28 ) the results are very impressive as their numbers are very close to the experimental mean activity coefficients , even up to 5 m. both observations imply that potential energies among smith and dang s na , cl , and spc / e water are reasonably balanced and can properly estimate the activity coefficient and the solubility .
however , this does not sufficiently support the idea that all of the smithdanggarrett ions are appropriately balanced .
we note that it is also possible that the flexible spc water model provides better effective interaction potentials than spc / e .
gavryushov and linse applied a similar methodology to lyubartsev and laaksonen , except that they used spherical boundary conditions and the spc / e fixed - point water model .
they determined mean activity coefficients of aqvist s ions(29 ) ( na and k ) in addition to smith and dang s ions .
note that the activity coefficients of nacl generally correspond to the true activity , regardless the source of the ion parameters , while those of kcl tend to poorly reproduce the experimental activity coefficients .
this likely results from inaccurate estimates of the solubility of kcl in simulation with available parameters compared to experiment .
lenart et al.(14 ) simulated ion solutions with a specially designed implicit solvent model that included a distance - dependent dielectric permittivity , and this led to improved results .
however , it should be noted that they did not use a standard combining rule for the lennard - jones potential and instead defined separate lennard - jones interactions between all of the atom pair types .
this suggests that current pair potentials still have room for improvement in the context of doing away with the simple combining rules .
experimental and computed diffusion coefficients of alkali and halide ions at infinite dilution are listed in table 3 and compared with those from pre - existing ion models .
the diffusion coefficient increases as the size of ion increases with characteristic peaks at rb and br .
although previous work claims to reproduce the qualitative trend of the diffusion coefficients with the smithdanggarrett ion parameters , we were not able to reproduce those results .
we only found the characteristic peaks with the water - model - specific cation parameters in spc / e .
the differences likely relate to differences in how the long - range electrostatic interactions were determined and the correction for the finite size effects .
moreover , in the current work , larger periodic unit cells were applied , which may also partially explain the difference . until our recent work ,
yet , it is well - known that tip3p water diffuses too rapidly and this in turn leads to ions that are too swift in tip3p water .
consequently , experimental diffusion coefficients are better reproduced in tip4pew and spc / e water , although there still are significant discrepancies .
these discrepancies partially result from concentration effects as the simulations were performed in 1 m solution while the experimental numbers are from dilute solution . when the concentration of the ions in solution increases , the diffusion coefficients of the ions decrease.(18 ) diffusion coefficients of the ions also appear to depend on the choice of counterion as ions tend to move swifter when they are solvated with fast - moving counterions ( see table s2 in supporting information ) .
generally , if the counterion is larger , the increment in the diffusion coefficient is greater . considering the relative diffusion rates of the ions and water molecules ,
swift molecules tend to pull or push the ions and thus enhance the diffusion coefficients of the ions .
also , spc / e - compatible br diffuses most rapidly with k and rb .
looking at the residence time between cation and anion ( shown in the next section ) , the results suggest that high residence times impede the movements of both the cations and anions ( because the paired ions are effectively heavier ) . generally ,
the larger ions have higher diffusion coefficients ; however , the relative values are influenced by ion pairing and the diffusion rates of other molecules in the solution .
diffusion coefficients were calculated at 1 m concentration with the counter ions cl and na . note that the experimental diffusion coefficients(32 ) represent dilute solution .
the diffusion coefficients before correcting for finite size effects are also provided in the supporting information .
the population of these clusters ( displayed in table s4 in the supporting information ) is closely related to the other dynamic properties of the solution.(51 ) from the data , the association constant for each single cationanion pair was calculated .
the true association constant is the ratio of the activity of ion pair to the activities of each of the single ions ( k = apair/(acation  aanion ) ) . due to the overwhelming computational demands of estimating activities for each of the ion models
, we assumed the activities were unity and substituted mole fraction , molar concentration , and molal concentration .
the association constants in terms of molar concentration are similar and are reported in the supporting information .
chen and pappu(52 ) performed a similar analysis using qvist ions(29 ) where the association constants were calculated in terms of molal concentration .
when these numbers are compared with fuoss conductimetry results,(53 ) the results suggest that the combination of qvist cations with the anions from chandresekhar et al.(54 ) and lybrand et al.(55 ) is excellent . however , this analysis neglected the activity coefficients .
as shown in chen and pappu s table 3 , the association constant is clearly greater when concentration is low , which proves activity coefficients can not be ignored .
the association at the low concentration is indeed close to the actual association constant because activity coefficients are close to unity at low concentrations
. if ion - dipole effects are neglected , the activity coefficients of ion pairs are almost unity .
if we take experimental mean activity coefficients for single ions from nacl and kcl , they are 0.779 and 0.768 at 0.1 m. therefore , 1.65 ( = 1/0.779 ) for nacl and 1.70 ( = 1/0.768 ) for kcl should be multiplied to their calculations . in other words ,
association constants of water - model - specific ions , and smithdanggarrett ions are shown in table 4 , and generally the values were far less than fuoss association constants even though these numbers are probably underestimated .
fuoss estimated the numbers based on the conductometric data assuming that the paired ions play no effective role in conductance of the ionic solution .
however , consistent with chen and pappu , we assume that the distances between ion pairs are within rmin1 .
however , gurney radii from fuoss data ( see fuoss table 2 ) , which are arguably the criterion determining free ion and ion pair , are generally much greater than our rmin1 values , especially for li- , na- , and k - associated salts , although the distances for some rb- and cs - associated salts are slightly shorter than our rmin1 values . for this reason , his nonconductive
ion pairs possibly include some solvent - shared pairs as well as contact ion pairs .
the association constants were recalculated using the criteria provided by fuoss to see how much the association constants change depending on the criteria ( table s5 in the supporting information ) . although significant differences from fuoss estimation were still observed , the new association constants agreed much better with fuoss numbers . with the same assumptions ,
most of the association constants of qvist s ions calculated by chen and pappu would exceed fuoss estimation .
it is not clear why the qvist ions , compared to the other ion models , show such significant differences .
ion concentrations listed in table s4 of supporting information were used to calculate the constants in terms of both molal concentration .
the unit of the constants is m. the deviations were calculated from the deviations of the individual populations in table s4 of supporting information .
the residence time between two molecules is an indicator of the lifetime of pair association at a particular temperature .
we calculated residence times of cationanion , watercation , and wateranion pairs in 1 m solutions using water - model - specific and smithdanggarrett ion models .
the characteristic radii from the rdf were used in the determination of the residence time and are listed in the supporting information ( table s3 ) .
the residence times between li and f are considered to be very long ( compared to the md simulation time scale ) , and thus we could not estimate their values .
also , we could not calculate the residence time between tip4pew - compatible cs and f and smithdanggarrett cs and f in spc / e water because of an insufficient size for the ensemble .
consistent trends in the cationanion residence times across a particular column or row of the table ( as a function of the ion sizes ) are not observed . when the cation is fixed ( along the column in table s6a in the supporting information )
, the curves of the residence times generally monotonically increase , monotonically decrease , or display a convex shape with a single maximum . for the tip3p - compatible cations
li and na , the residence times monotonically drop as the size of the paired anions increases . k and rb have maximal residence time with br , whereas cs does with i. the trend is similar with the tip4pew - compatible cations . however , here na has a maximal residency time with br , while k , rb , and cs have their peaks with i. spc / e - compatible cations behave slightly differently from the ions of the other models .
na ( and possibly li ) ions have both local minimum and maximum with cl and br , respectively .
however , smithdanggarrett cations do not have such maxima in the middle of the series .
when the cation is smaller than na , cationanion residence times continuously decrease as the size of the anion increases .
when the cation is larger than k , residence times tend to increase as a function of the anion size .
these behaviors consistently can be observed in the three different types of water models . on the other hand , in the series of anions ( along the row in table s6a in the supporting information ) , the curves of the residence time generally show two patterns : monotonous decrease or concave shape having one minimum .
the minimums are located at slightly different positions depending on the anion and ion models .
all the tip3p - compatible anions have the longest residence time when they are paired with li .
the residence times in the series gradually decreased as the size of the pairing cations increases excluding i , which has a slight increase compared with cs .
however , in the series of br and i , the shortest residence times were reached at rb
. minimums of residence times also can be observed with all the spc / e - compatible anions ; f with k , and the other anions with na . for the smithdanggarrett anions , such minimums
unexpectedly , the residence time of i with these parameters in spc / e shows an increasing pattern as a function of cation size .
overall , the residence times are complicated and depend on both the ion parameters and the applied water model .
residence times between cations and water molecules ( table s6b in the supporting information ) show little dependence on the size of the counterions . unlike cationanion residence times
, the average watercation residence times generally are not strongly sensitive to the choice of ion force field , although this appears to be less true with li and na .
however , the differences among anionwater residence times with different ion models are small , even with the smallest anion , f ( table s6c in the supporting information ) .
these observations are confusing because , while significant influence of the water and ion parameter model is observed for anioncation interactions , little is observed with ionwater residence times .
this likely results from the heavy focus on ionwater interactions , as compared to ionion , in force field development .
although our work in the development of water - model - specific ion models did include ionion interactions , specifically , the crystal lattice energy and lattice constant , the ionion residence times lack consistency . on the other hand ,
residence times of ionwater pairs ( wi ) are closely related with the diffusion coefficients of ions ( di ) and water molecules ( dw ) .
this can be seen when di and dw/wi1/2 are plotted in an ad hoc manner to show relationships among the data ; when this is done , a roughly linear correlation is observed ( figure s3 in the supporting information ) .
interestingly , the slopes of the linear relationship are almost identical no matter what ion or water models are used . however , the diffusion rates clearly depend on the concentration of the ions and also the nature of the cationanion and waterwater interactions
. unexpected crystallization of alkalihalide salts at concentrations below their saturation limit motivated our development of new ion parameters.(11 ) a necessary test of these parameters is alkalihalide pair solubility .
transitions in the crystal structures , detected by sudden change of volume or potential energy , were observed .
also , the structures of the newly growing crystals on the surface of the pre - existing crystal were different from their natural crystal structure .
another problem was exposed in salts paired with fluoride ions . the volume and the potential energy of the systems fluctuated at low frequency .
although the systems may satisfy the condition of equilibrium as described in the methods section , in this case , the equilibrium states were suspicious because of the low - frequency fluctuations of the two properties .
for these reasons , we only show the valid solubility results in table 6 . depending on the choice of ion pair and ion parameters , simulations required about 110530 ns for the solution crystal coexisting systems to reach their equilibrium states .
this simulation time could , in principle , be reduced if the initial concentration and density of the solution opposed with the crystal are very close to its saturation condition . notwithstanding , it is presumed that at least 100 ns of simulation are required .
the equilibrium concentrations of the cation and the anion in the same system were almost identical as expected .
thus , the concentrations of the paired salt shown in table 6 are the average of the cation and the anion .
the chemical potentials of crystals can be assumed to be constants because the dependency of chemical potential of crystal on its size is usually negligible .
the low solubility of the tip3p - compatible ions implies that the chemical potential of the crystals are too low in general because at low concentration the deviation of the activity coefficients from the experimental results is small , as shown in figure 2 , and the deviation of chemical potential of the solution from the real solution might be small .
although the absolute solubility values of the tip3p - compatible ions deviate significantly from experiment , the qualitative solubility trends are most accurately reproduced with these parameters compared to the other models ( table 6 ) .
smithdanggarrett ions also showed low solubility in tip3p water , except for nai , which had a higher solubility than true nai solution .
however , the qualitative trends were not observed with smithdanggarrett ions in tip3p water : nacl and nai were relatively good , but the other salts showed significant deviation .
tip4pew - compatible ions well reproduced the saturation concentrations with na - associated or cl - associated salts .
spc / e - compatible ions showed similar behaviors with na - associated or cl - associated salts , but the solubilities were not as accurate as tip4pew ions .
smithdanggarrett ions did not behave much differently in spc / e water than they did in tip3p water .
they showed relatively good solubilities with nacl and nai . since the ion models are identical for the smithdanggarrett ions , the chemical potential of the crystals are presumably the same in both water models .
in addition to this , as the solubility is not strongly affected by the water models , the chemical potential of the solution as a function of concentration is likely also almost identical regardless of the choice of water model . as a result , we can presume that the activity coefficients do not change greatly as a function of concentration in different water models as long as the ion models are identical . however , absolute hydration free energies can be different.(11 ) the various simulation sets except tip3p
their solubilities were generally well reproduced with nacl . however , as the relative distance from nacl becomes farther , the deviation of the solubilities seems to increase . to better understand the deviations in solubility compared to experiment
lattice energy is the potential energy difference generated when ions in a crystal are separated to infinite distance . a plot of the difference in experimental and calculated lattice energies versus the natural log of the ratio of calculated and experimental molal concentrations exposes trends for the different ion model parameter sets ( figure 3 ) .
if good agreement between the experimental and calculated values in both properties is obtained , the points will tend to cluster about the origin .
the data suggest that the errors in the lattice energies for the various salts calculated using smithdanggarrett ions are generally positive ( calculated < experiment ) in contrast to water - model - specific ions.(11 ) interestingly , regardless of the sign of the deviation in the lattice energy , most of the data points are located above x axis .
the data in the figure also suggest that greater absolute deviations between experiment and calculation of the lattice energy tend to lead to greater deviations in the solubility
. however , there are trends in the deviations of the lattice energy that depend on the choice of cation . with the water - specific ion parameters , na - associated salts were less sensitive to the deviation of lattice energy ( solid line ) while rb - associated salts were much more sensitive ( dashed line ) .
although there is a slight discrepancy regarding the extent of the effect , the solubility of smithdanggarrett salts was similarly affected by the deviations of lattice energy , and the effects were also distinctive depending on the cations . in conclusion
we tried to find correlations of solubility with lattice constants , but we could not find obvious correlations .
these results suggest that the parametrization of the tip3p - compatible ions can be revisited to improve salt solubilities by focusing more heavily on reproducing the lattice energy .
x axis indicates the difference of the lattice energy calculated ( lecalc ) from the experimental lattice energy ( leref).(11)y axis indicates the logarithmic ratio of experimental saturation concentration ( mref ) to the calculated saturation concentration ( mcalc ) in terms of molality .
each point denotes the alkalihalide salts , as shown in the legend on the right - hand side .
different colors specify the ion and water models : blue for tip3p - compatible ions in tip3p water , magenta for tip4pew - compatible ions / tip4pew water , yellow for spc / e - compatible ions in spc / e water , cyan for smithdanggarrett ions in tip3p water , and violet for smithdanggarrett ions in spc / e water . among them , water - model - specific ions were fit to y = ax curves , where a is a fitting parameter .
the solid line , the dotted line , and the dashed line are the fit curves for na - associated salts , k - associated salts , and rb - associated salts .
figure 4 shows the correlations between the association constants of alkalihalide pairs ( listed in table 4 ) and their solubility . generally , the two properties were inversely proportional , such that ion pairs with higher association constants have a greater tendency to form crystals .
fuoss association constants do not show this correlation ( figure 4b ) ; however , as discussed , their gurney radii might not be able to distinguish contact ion pairs from solvent - shared pairs , which means the reported association constants also include non - crystal - forming associations . in any case , the correlation between the calculated cluster populations and the saturation concentrations is useful since the more easily calculated association constant can be utilized to estimate solubility .
correlation of the association constant ( k ) of cationanion pair to the saturation concentration ( m ) .
calculated molal saturation concentrations of salts are inversely proportional to the association constant of the ion pair ( a ) .
the solid line indicates the fitted curve , mcalc = 0.593k ( r = 0.7618 ) .
however , the experimental saturation concentrations(32 ) did not show correlation with fuoss association constants(53 ) ( b ) .
the results suggest that high cationanion residence times drive the ions in solution toward the crystal form , whereas high watercation and wateranion residence times tend to dissolve the crystal ( see the legend of figure 5 ) . although the results appear reasonable ,
these observations can not be experimentally validated as we do not have an experimental measure of the residence times .
correlation of the residence times to the saturation concentration . under the assumption that the inverse of residence time is proportional to the dissociation rates , and
that the pairing between ionwater and ionion is a first - order reaction implying that the association constants are the ratio of the rates of the forward and backward reactions , the calculated molal saturation concentrations of alkalihalide salts ( mcalc ) were plotted against c - a/w - a  w - c . here c - a , w - c , and w - a are residence times of cationanion , watercation , and wateranion , respectively .
the solid line indicates the fitted curve of all the data points , y = 0.460x ( r = 0.7656 ) .
the color scheme is same as that in figure 3 . by investigating these correlations , the links among the properties  including diffusion coefficients , residence times , association constants , and solubility 
figure s3 in the supporting information shows that the diffusion coefficients of the ions are largely dependent on the choice of water model .
also , looking at the smithdanggarrett ions in table s6b , c in the supporting information , the ionwater residence times also depend on the choice of water model .
the dependency of the residence times on the choice of water model was generally in the order , wateranion < watercation < cationanion , as shown in table s6 in the supporting information .
when we assume that the wateranion and watercation residence times can be correctly calculated , this suggests a critical need to balance the ionion interactions .
this balance between the residence times is also closely related to the solubility , as shown in figure 5 , and even to the association constants based in figure 4 .
originally , water - model - specific ion models were optimized using multiple properties reflecting ionion interactions . in spite of the effort ,
cationanion residence times , association constants , and solubility show large variability depending on the choice of water model .
this , along with the sensitivity of the solubility to the lattice energy ( figure 3 ) , again emphasizes the importance of ionion interactions in the parametrization of ion force fields .
the relevance of each of these properties ( diffusion coefficients , residence times , lattice energy , association constants , and solubility ) to the activity coefficients was not assessed in detail here as the limited set of calculations performed was extremely computationally demanding .
the activity coefficients of the smithdanggarrett ions calculated by gavryushov and linse(12 ) suggest that the calculated nacl activity coefficients are relatively good compared to kcl .
these results suggest that , if solubility is accurately simulated , the activity coefficient should also be accurately reproduced .
therefore , the tip3p - compatible ions should not be used for simulations at high salt concentrations .
although smithdanggarrett ions are fine with nacl , their kcl is not appropriate for simulating high salt environments .
for this reason , we recommend tip4pew and spc / e water models together with their water - model - specific ions for simulating high salts concentration .
the dynamic and energetic properties of alkali and halide ions were calculated in simulation using various different ion and water force fields including our previously developed water - model - specific ion parametrizations.(11 ) mean activity coefficients were determined for nacl and kcl in tip3p water .
qualitatively , the mean activity coefficients correspond to experiment ; however , the deviations tend to increase as the salt concentration increases .
although extremely computationally demanding , to better determine which model is most appropriate for reproducing the experimental activity coefficients , all of the pairs of alkalihalide salts should be investigated using explicit solvent with the various parameter sets . however , these and previous results do clearly suggest that accurate estimation of solubility is critical to accurately calculate activity .
diffusion coefficients appear to be largely dependent on the choice of water model . generally , faster moving water models resulted in higher diffusion coefficients for the solvated ions , but at the same time , the diffusion constants were almost inversely proportional to the square root of the residence time between the ion and the water molecule .
we also calculated the association constants of various alkalihalide salts in terms of molar and molal concentration .
however , the lower association constants could be explained by two reasons , specifically ( 1 ) since fuoss experimental data use a greater cutoff length defining ion pairs which may include solvent - separated pairs , and ( 2 ) if activity coefficients are properly considered , our calculated association constants in terms of molar and molal concentration would increase . the residence times between ion and water and ion and ion were also calculated .
ionion residence times varied significantly according to the choice of water model , while the ionwater residence times did not .
the tip3p - compatible ions showed qualitatively correct differences among the various salts , but the solubilities were generally significantly underestimated .
the tip4pew and spc / e ion parameter models well reproduced the solubilities of na - associated salts and cl - associated salts .
smithdanggarrett ions also qualitatively reproduced the solubilities of nacl ; however , the deviations of the other ions in the parameter sets from experimental were generally greater than observed with the tip4pew- and spc / e - compatible ion parameters .
some of the dynamic and energetic properties of alkali and halide ions show a strong dependence on the choice of water model .
ion force fields can then be developed on top of that foundation with care to balance the ionion interactions as many pairwise properties such as the association constant , the cationanion residence time , and the solubility are very sensitive to these . in summary , the results suggest that our spc / e and tip4pew water - model - specific ion parameters fairly accurately model the behavior of the monovalent ions in aqueous solution across a wide range of concentrations , whereas the tip3p - compatible ion parameters should be avoided at high salt concentrations .
 the dynamic and energetic properties of the alkali and halide ions were calculated using molecular dynamics ( md ) and free energy simulations with various different water and ion force fields including our recently developed water - model - specific ion parameters . 
 the properties calculated were activity coefficients , diffusion coefficients , residence times of atomic pairs , association constants , and solubility . through calculation of these properties 
 , we can assess the validity and range of applicability of the simple pair potential models and better understand their limitations . 
 due to extreme computational demands , the activity coefficients were only calculated for a subset of the models . 
 the results qualitatively agree with experiment . 
 calculated diffusion coefficients and residence times between cationanion , watercation , and wateranion showed differences depending on the choice of water and ion force field used . 
 the calculated solubilities of the alkalihalide salts were generally lower than the true solubility of the salts . 
 however , for both the tip4pew and spc / e water - model - specific ion parameters , solubility was reasonably well - reproduced . finally , the correlations among the various properties led to the following conclusions : ( 1 ) the reliability of the ion force fields is significantly affected by the specific choice of water model . 
 ( 2 ) ionion interactions are very important to accurately simulate the properties , especially solubility . 
 ( 3 ) the spc / e and tip4pew water - model - specific ion force fields are preferred for simulation in high salt environments compared to the other ion force fields . 
Introduction
Methods
Results and Discussion
Conclusion