we 
 describe here a method that can rapidly and accurately compute 
 the potentials of mean force ( pmfs ) between spherical univalent ions 
 in water .
the pmf between an ion a and ion b represents the reversible 
 work ( as a function of distance between the two ions , r ) , averaged over all the configurations of the surrounding solvent , 
 needed to bring a and b from infinite separation to a distance r. pmfs are central to many solvation 
 processes .
for example , some pairs of ions form a stable contact 
 pair or salt bridge ; others form a stable solvent - separated pair and have intervening bridging 
 water molecules .
the pmf is relevant to understanding the 
 solubilities of salts in water , the equilibria and dynamics when a 
 drug molecule approaches  and binds to  a protein binding 
 site , and how some parts of a protein molecule approach other parts 
 as the protein molecule folds .
currently , the gold standard 
 for computing the pmfs between large 
 molecules  such as are typical in biology  is to use explicit - solvent 
 simulations with semiempirical force fields .
however , these simulations 
 are often computationally too expensive for many of the calculations 
 that would be of practical interest . even in the case of
various sampling techniques have been introduced to mitigate the computational 
 expense : constraints , umbrella sampling , and weighted histogram methods , for 
 example . on the other hand ,
treating water explicitly in analytical 
 theories is challenging , because of the need to properly account for 
 the orientation ( angular ) effects of water .
for simple salt solutions , 
 integral - equation theories often predict qualitatively correct trends but not quantitative details .
an alternative is to 
 use implicit - solvent models ( such as poisson  boltzmann ( pb ) 
 or generalized born ( gb ) ) , in which the solvent is approximated as a dielectric 
 continuum .
implicit - solvent models have the advantage of computational speed : 
 a pmf costs around 100 h of md simulation time in explicit water but 
 only minutes on a personal computer for pb pmf or fractions of a second 
 for gb .
pmfs of a model nacl from explicit - solvent md simulation ( symbols ) 
 and implicit - solvent theories : pb ( continuous blue line ) and gb ( dashed 
 green line ) .
pb calculations were done with apbs 1.3 software , and gb results were obtained according to the 
 gb model of still et al .
the permittivity 
 of water was set to 78.54 in the pb and gb calculations .
parameters 
 from dang et al . were used for na and cl ions ( see table 1 ) . however
, implicit - solvent models do not capture the subtle 
 features 
 of pmfs , which require some representation of particulate water ( see figure 1 ) . to correct for the energy 
 errors in implicit modeling
implicit - solvent models can 
 allow interstitial high dielectrics ( e.g. , water ) , but they give incorrect 
 pmfs between nonbonded atoms .
salt bridges , 
 for example , treated via gb models are typically too stable compared 
 to explicit - solvent results , which leads to oversampling such conformations 
 in md simulations .
much work has been 
 done to overcome these difficulties , such as the following .
a method 
 for including a single explicit bridging water between the pairs of 
 the salt bridge charges was proposed by yu et al . , and it was shown that this significantly improves the predictions . 

geney et al . introduced an empirical 
 correction to the dielectric radii of hydrogen atoms of charged protein 
 groups , bringing the contact depth of gb pmfs in better agreement 
 with explicit solvent for the tested salt bridges .
recently , nguyen 
 et al . did more systematic fits to improve 
 the gb solvent model parameters for protein simulations . to account 
 for the proper coupling between nonpolar and polar solvation energies 
 in implicit - solvent models , dzubiella et al .
proposed a variational formalism where the gibbs free energy of 
 the system is expressed as a functional of the solvent volume exclusion 
 function .
the method captures the sensitivity of solvent expulsion 
 in cases of simple spherical solutes .
proposed an extension of the gb model to correctly describe 
 the solvent - excluded volume of each pair of atoms , and that improved 
 the nonbonded pmfs .
an empirical model of chen and brooks captures the context dependence of the effective 
 surface tension , and the authors showed that it can resolve deficiencies 
 in pmfs of nonpolar peptide side chain analogues .
we first perform md simulations in explicit 
 solvent as our sort of gold standard for what we aim to capture in 
 a fast computational method .
we describe here a way  which we 
 call interpolation pmf ( i - pmf)to determine 
 any ion  ion pmf from this set of presimulated md explicit - water 
 pmf simulations .
i - pmf gets its speed because we simply capture 
 the explicit results by simple interpolations of the precomputed results 
 that we store in look - up tables .
we show here that this simple approach 
 results in fast and accurate pmf determinations .
i - pmf is in the spirit 
 of a recent solvation model called sea ( semi - explicit assembly ) in 
 which solvation free energies of single spheres are 
 computed through an explicit - simulation precompute step , and then 
 assembled into the solute of interest at run - time by summing regional 
 free energies .
we describe below the two components of the i - pmf method : 
 first , 
 we performed extensive md simulations of small spherical ions of different 
 charges and radii , and harvested their pmfs .
second , we describe our 
 interpolation algorithm , the i - pmf , that allows us to capture the 
 presimulated md results with interpolation fits .
next , we show that 
 the i - pmf method gives good agreement with md and experimental results .
we performed md simulations , mostly on the titan 
 supercomputer at oak ridge national laboratory , using version 4.6.2 of gromacs .
each simulation consisted of two charged solute particles and 635 
 water molecules in a rhombic dodecahedron box .
hoover 
 thermostat ( coupling constant of 2 ps ) and parrinello  rahman 
 barostat ( coupling constant of 10 ps ) .
the equations of motion were 
 integrated using the leapfrog algorithm with a time step of 0.2 ps . 

the smooth version of the particle mesh ewald ( pme ) was used to treat electrostatics ( grid spacing of 0.12 , 
 pme order of 4 , real - space cutoff of 9  , and ewald screening 
 parameter of 0.347  ) , and a 9  cutoff 
 for lennard - jones ( lj ) interactions was used along with long - ranged 
 dispersion corrections for energy and pressure .
the lj size parameter 
 ( lj ) of the solutes ranged from 2 to 5.5  
 in increments of 0.5  , while the lj energy parameter ( lj ) was set to 0.1 kcalmol for all 
 solutes .
the choice for the lj range and for the 
 lj value was based on the values for the alkali 
 metal ( na , k , rb , cs ) and halogenide ( cl , br , i ) ion lj parameters of the dang force field .
the solute particles were each given a formal charge of 1 
 or + 1 .
we simulated all possible combinations of solute pairs ( of 
 different sizes and charges ) , a total of 136 combinations ( + 1:1 , 
 + 1:+1 , and 1:1 ion pairs ) for a given water model . 

we obtained the 
 pmfs from md simulations in which constraints were applied to hold 
 the solute pairs at a series of fixed interparticle separations during 
 the course of the simulation .
the 
 closest distance between the solute particles was 1.6  , and 
 the largest interparticle distance was 12  . in the interval 
 from 1.6 to 5  ,
the separation step was 0.1  , from 5 to 
 8  it was 0.2  , and 0.4  increments were used for 
 the rest ( up to 12  ) .
therefore , to construct the pmf for a 
 given solute pair , 60 independent simulations were performed where , 
 after energy minimization and 0.1 ns equilibration , the trajectory 
 over 3 ns was collected ( in total , this amounted to 136  60 
 = 8160 independent md runs ) .
the individual pmf was calculated by 
 integrating the average mean force over the solute separation distance , 
 while the entropic force due to the increase in phase space with solute 
 separation was added to the average mean force prior to integration .
uncertainties in the average mean force values 
 were estimated from the limiting value of the block averages , and the errors in pmf were calculated by integrating 
 the variances of the mean force . in short , in the presimulation 
 step of the i - pmf method , we performed 
 explicit - solvent ( tip3p water ) md simulations to obtain the pmfs for 
 different combinations of solute charges ( q1 , q2 ) and lj size parameters ( lj,1 and lj,2 ) .
we tabulated the values of the 
 presimulated pmfs for given ( q1 , q2 ) and ( lj,1 , lj,2 ) at the end of the supporting information file . to validate that the i - pmf method is predictive , we 
 performed additional 
 test md runs for selected charge and size combinations of ions with 
 the same protocol .
the lj energy parameter was the same for all ions ( lj , ion = 0.1 kcal mol ) .
a pmf is the free energy as a function of distance r between the two ions .
in addition , here we aim to account for two 
 other variables at the same time , lj,1 and lj,2 , the diameters of the two ions , each of which has a given 
 fixed charge .
below is the procedure for performing this interpolation 
 from our md simulation data , and we use this process to determine 
 pmfs for any given pair of charges on the ions .
here is a short summary 
 of how interpolations are performed on two variables at the same time . 

suppose that we are given a matrix of functional values , f(xi , yj ) , that correspond to points on a two - dimensional 
 ( n  m ) grid of points ( xi , yj ) .
here , i goes from 1 to n and j goes from 1 to m. in a two - dimensional interpolation , we want to estimate , by interpolation , 
 the function f at some point not tabulated ( xa , yb)we 
 want to find f(xa , yb ) .
the point ( xa , yb ) falls into a grid square ; i.e. , 
 four tabulated points surround the desired interior point .
say that xi  xa  xi+1 and xj  yb  xj+1 ; let us denote these four functional points in the following 
 way1 in the simplest case of interpolation on the 
 grid square 
t  [ 0 , 1 ] and u  [ 0 , 1 ] are3 in order to obtain a smoother interpolated 
 surface , higher - order interpolation techniques are applied .
for example , 
 in bi - cubic spline interpolation , we enforce the 
 smoothness of some of the derivatives as the interpolating point cross 
 grid square boundaries .
to do a bicubic interpolation within a grid 
 square , one needs to know the function f and the 
 derivatives fx = f/x , fy = f/y , fxy = f/xy at each of the four corners of the square . 

the values of the derivatives at the grid points are determined globally 
 by one - dimensional splines .
the interpolating function has the following 
 form4and we need to find 16 coefficients cmn by first solving the equations 
 for the square grid points .
the function f was in 
 our case the pmf , and xi and yj were the lj . besides bilinear and cubic spline interpolation , we tested 
 the performance of two additional interpolation techniques ( nearest - neighbor 
 and piecewise - cubic - hermite - polynomial interpolations ) . for
a more 
 detailed discussion regarding the two - dimensional interpolation algorithms , 
 see refs ( 52 ) and ( 53 ) .
now , we describe 
 the i - pmf interpolation methodology for capturing the pmf for any 
 particular pair of ion sizes and charges and water model . for a given 
 pair of solutes with charges q1 and q2 ( qi = + 1 or 1 ) and for a given interparticle separation 
 distance r , a matrix of pmf values corresponding 
 to different combinations of the solutes lj
size parameters 
 ( lj,1 and lj,2 ) was built from 
 the presimulated pmfs .
the span of lj,1 and lj,2 was from 2 to 5.5  in increments of 0.5  , 
 so the matrix dimension was 8  8 . in total ,
60 such matrices 
 were constructed , each belonging to a given interparticle separation 
 distance r  [ 1.6,12 ]  ( see section 
 above ) . for every matrix , a two - dimensional interpolation fit through 
 the grid points was performed . for that purpose
a matrix of interpolated values was constructed 
 and used as a look - up table to determine the value of the pmf at a 
 given r for an arbitrary pair of lj,1 , lj,2 
 [ 2 , 5.5 ]  values . lj,1 and lj,2 are different from the size parameters 
 lj,1 and lj,2 of the ions used 
 in the grid .
we stored the 
 interpolated values of the pmf for every lj,1 in the 
 range from 2 to 5.5  in increments of lj,1 = lj,2 = 0.01 
  ( this means that the matrix of interpolated values was of dimension 
 351  351 at each r ) .
although , for example , 
 lj values of ions in dang force field are given 
 to five significant digits ( see table 1 ) , no 
 significant differences in predicted pmfs were observed if smaller 
 lj were used to store the interpolated data .
we tested different 
 interpolation algorithms . by comparing the 
 pmf obtained via an interpolation scheme and a full explicit - solvent 
 simulation
, we concluded that the cubic spline interpolation method 
 worked best in the majority of cases .
details of this comparison are 
 shown below . for our initial i - pmf implementation
, it takes 
 approximately 8 
 s on a personal computer ( pc ) to build tables of presimulated values , 
 interpolating through the grid points , making the look - up tables of 
 interpolated values , and extracting from them the pmf for a desired 
 combination of charges and ion sizes .
after the interpolation is performed 
 a single time and the look - up tables are available , the pmf can be 
 extracted in 23 s on a pc .
our md pmfs required a cluster of computers 
 with many processors ( the value of the pmf at each r is its own simulation ) , and in total , it takes more than 100 h of 
 cpu time to get a single pmf .
presimulating 
 a systematic series of ion size and charge combinations on a supercomputer , 
 oak ridge national lab s titan in this case , provides us with 
 a data set of accurately calculated pmfs .
we perform independent two - dimensional 
 interpolations over these solute size combination grids to determine 
 the interaction strength for two arbitrarily sized solutes across 
 the full range of particle separation distances .
the estimated pmf 
 is finally assembled from this interpolation data for this specific 
 solute pair .
figure 3 shows 
 the results of the i - pmf method utilizing 
 four different interpolation algorithms .
a size asymmetric + 1:1 
 system mimicking potassium iodide ( lj,+ = 3.3345 
  , lj, = 5.1705  ; see table 1 ) at infinite dilution was selected as a representative 
 case .
a similar analysis of i - pmf predictions in other cases shown 
 in this paper is provided in the supporting information .
the simulated pmf shows a deep first minimum , corresponding to 
 the contact ion pair ( rcp  3.6 
  ) , a desolvation barrier at r  4.6 
  , and further peaks involving the solvent shared / separated pairs 
 ( cf . table 2 ) .
we see that all of the interpolation 
 algorithms give a pmf that qualitatively agrees with the simulated 
 one . nevertheless , the performance of the nearest - neighbor interpolation 
 is the least accurate of the four interpolation algorithms over a 
 broad range of interparticle separation distances .
it overestimates 
 the contact pair ( cp ) peak and underestimates the desolvation barrier 
 ( numerical data are collected in table 2 ) . 

the bilinear interpolation underestimates the depth of the contact 
 pair minimum , and the values at short r ( repulsive 
 wall ) are shifted to somewhat larger r. the value 
 of the cp which follows from piecewise cubic hermite polynomial interpolation 
 lies on the border of the experimental error , but the rest of the 
 pmf agrees within the error bars with the md results .
the pmf obtained 
 via the cubic spline interpolation through the grid points gives the 
 best agreement with the simulated pmf .
the positions of the characteristic 
 peaks are within the error bars ( table 2 ) . here , we want 
 to find the best interpolation method .
the goodness 
 of the interpolation schemes for the i - pmf method can be expressed 
 in terms of the reduced  parameter ( see the supporting information file for the definition ) . 

the values of red are 
 22.0 for bilinear interpolation , 8.2 for nearest - neighbor interpolation , 
 1.8 for piecewise cubic hermite polynomial interpolation , and 0.4 
 for cubic spline interpolation algorithm .
the large value of red in the case of the bilinear interpolation algorithm 
 originates in the differences between the predicted and simulated 
 pmf at small interparticle separation distances ( the repulsive wall 
 of the pmf ) .
the values of the pmf on the steep repulsive wall are 
 on the order of 10 kcalmol ( and 
 larger ) , but the error in the case of + 1:1 ion pairs is small 
 ( 10 kcalmol ) . 
 even though the differences between the predicted and simulated pmfs 
 are only around 1 kcalmol , the red becomes extremely large due to the small error 
 in the simulated pmf .
this region of the pmf is not very important 
 when one discusses the structural aspects of ion  ion and ion 
water 
 interactions ( the red was indeed defined 
 in such a way as to avoid points on the repulsive wall ; see the supporting information file ) . if two points on 
 the steep repulsive wall , lying just before the cp minimum ( r < 3.6  ) , are not taken into consideration when 
 calculating red , then the value of red for bilinear interpolation lowers to 6.1 .
this 
 smaller value compared to the red for 
 the nearest - neighbor case agrees with the previous visual judgment 
 of a better overall prediction of the bilinear interpolation compared 
 to a nearest - neighbor fit . as a rule , a value of red = 1 suggests that the extent of the agreement between 
 observations and estimates is in accord with the error variance .
red > 1 indicates that the fit has not fully 
 captured 
 the data ( or that the error estimate has been underestimated ) , while 
 red < 1 speaks in favor of an unexpectedly 
 good fit .
we have concluded that ( in the majority of cases ) the most 
 satisfactory interpolation algorithm for predicting the pmf of opposite - charged 
 as well as like - charged pmfs involves the cubic spline interpolation 
 ( see tables s1s4 of the supporting information file ) .
the method gives quantitative agreement within the error 
 bars in most of the relevant cases displayed in the rest of the figures 
 and in data presented in the supporting information . in the following subsections ,
we show a few characteristic cases 
 that further demonstrate the quality of the proposed pmf prediction 
 scheme . here
we use only a single value 
 of the 
 lj energy parameter ( lj ) , regardless of the charge 
 or size of the ions .
this is true for the selected ions in the dang 
 ion force field ( cf .
table 1 ) , but in principle , 
 each ion could require its own lj value .
in that case , the presimulated pmfs would need to be collected 
 for different combinations of lj . instead of performing 
 a two - dimensional interpolation , a four - dimensional interpolation 
 through the pmf values corresponding to ( lj,1 , lj,1 ; lj,2 , lj,2 ) would 
 need to be conducted .
comparison of i - pmf predicted pmfs using different interpolation 
 algorithms for pmf prediction in the case of a model ki salt in tip3p 
 water : nearest - neighbor interpolation ( blue line ) , bilinear interpolation 
 ( orange line ) , piecewise cubic hermite polynomial interpolation ( green 
 line ) , and cubic spline interpolation with smooth first and second 
 derivatives ( red line ) .
the cubic spline interpolation 
 algorithm best captures the magnitude and shape of pmf features .
here , we 
 compare the results of the i - pmf method to the results 
 of the precomputed md simulations .
we find that i - pmf works quite 
 well across the whole range of solutes , from small ions with high 
 charge density to large ones having small charge density .
it captures 
 the position and the depth of the characteristic wells corresponding 
 to the contact and solvent shared / separated pairs and the height of 
 the peaks of the desolvation barriers .
the method is even capable 
 of capturing the complex behavior of the pmf in the case of small 
 ions having a large charge density ( <3  ) , where the path 
 from the solvent shared to solvent separated pair involves multiple 
 transition barriers .
i - pmf prediction ( lines ) and md simulation ( symbols ) pmfs 
 for model 
 size symmetric cations and anions .
solute sizes are lj,+ = lj, = 2.75 ( blue ) , 3.25 ( violet ) , 3.75 
 ( purple ) , 4.25 ( brown ) , 4.75 ( green ) , and 5.25  ( yellow ) ; inset , 
 2.25  . in all cases , we find that the predictions are nearly 
 indistinguishable from md simulations , and i - pmf even captures the 
 structural features of  extreme  ion pairing cases where 
 the ion sizes are unphysically small ( inset ) . up to a distance of the cp minimum ( i.e. , for the values 
 on the 
 repulsive wall ) , the cubic hermite polynomial interpolation performs 
 slightly better than cubic spline interpolation ( not shown ) . in an 
 extreme case of lj = 2.25  ,
the cubic hermite 
 polynomial interpolation works better up to the distance of the desolvation 
 barrier .
however , we need to stress that in these cases the selected 
 size parameters of ions are small ( 2.25  ) and do not correspond 
 to practical model situations .
cases 
 ( inset in figure 4 ) . improved performance of 
 the interpolation could be achieved by providing a finer grid at the 
 size edges , i.e. , between 2 and 2.5  , and equivalently between 
 5 and 5.5  for larger solutes .
this would , however , require 
 an increase in the number of pmfs calculated in the presimulation 
 step .
figure s1 of the supporting information shows how the interpolated pmf values change with increasing size 
 of the ions ( size symmetric case ) for few selected interparticle separation 
 distances .
comparison to the simulated values shows that the deviations 
 between predicted and simulated values are the largest for solutes 
 between 2 and 2.5  .
the red values 
 for the cases displayed in figure 4 are given 
 in table s2 of the supporting information file .
in the cases where red is smaller 
 for the cubic hermite polynomial interpolation compared to cubic spline 
 interpolation , this originates from the data on the repulsive wall . 
 for the relevant domains of the pmf
we also find that , 
 in addition to working for ions having identical 
 sizes , i - pmf works for ions with different sizes .
figure 5 shows pmfs for combinations of small and large 
 cations ( na and cs , respectively ) with small 
 and large anions ( cl and i , 
 respectively ) , using ion parameters listed in table 1 .
results of predicted pmfs for nacl and nai ion pairs in 
 tip3p water are given in panel a , and pmfs for cscl and csi , in panel 
 b. displayed md simulation data ( symbols ) were obtained independently 
 of the prediction scheme . in all cases , the predicted pmf ( cubic spline 
 interpolation through a presimulated grid of pmf values ) agrees quantitatively 
 with the reference md simulations .
results for other combinations 
 of ions given in table 1 are shown in figure 
 s3 of the supporting information file . 

the method performs somewhat worse in the case of kcl , kbr , rbbr , 
 and csbr where some parts of the pmf are either slightly over- or 
 underestimated compared to corresponding md results .
the differences 
 are , however , small , and the agreement can still be considered very 
 good .
the quality of the prediction was estimated for all the shown 
 cases in terms of the red parameter 
 ( table s1 of the supporting information file ) .
out of all four interpolation algorithms tested here , the 
 cubic spline interpolation performs the best .
we find that i - pmf gives 
 more accurate pmfs ( figure 5a ) than the simple 
 pb or gb models in figure 1 with similar computational 
 performance .
we also computed pmfs for like - sign pairs , such 
 as na  na and cl  cl ion pairs .
previously , such quantities have been studied 
 by integral - equation 
 theories and computer simulations .
again , we find excellent agreement between the interpolation formulas 
 and the md simulations ; see figure 6 .
. examples of pmfs 
 for other like - charged pairs of ions from table 1 are given in figure s4 of the supporting information file .
i - pmf prediction ( lines ) and md simulation ( symbols ) pmfs for salts 
 with dissimilarly sized cations and anions . displayed are sodium ( panel 
 a ) and cesium ( panel b ) chlorides ( green ) and iodides ( purple ) .
we 
 find i - pmf gives quantitative agreement in all cases , indicating that 
 more complex mixed - size solute pairings are no more challenging than 
 the previous size - symmetric solute pairings . in this section
, we explore the application of our pmf 
 modeling to ion solubilities in water . a quantitative measure of the 
 ion pair formation is given by the corresponding change in standard 
 free energy .
i - pmf predicted ( lines ) and md simulation ( symbols ) pmfs 
 for na 
na ( yellow down - triangles and 
 line ) , cl  cl ( green up - triangles 
 and 
 line ) , and na  cl ( brown circles 
 and line ) .
these results indicate that i - pmf can successfully predict 
 all the potential pair interactions ( including like - charge interactions ) 
 present in a given salt solution simulation .
the i - pmf nacl curve 
 also can be compared with the implicit results in figure 1 to highlight the method s ability to capture 
 the detailed pairing features seen in explicit md simulations .
the apparent association constant 
 for the process of forming an 
 ion pair is in the case of univalent ions ( a + b  ab ) defined as ka = ( [ ab]/c)([a]/c)([b]/c ) , where 
 [ x ] denotes the molar concentration of species x and c = 1 mol dm is the standard 
 molar concentration .
the change in standard free energy of ion pair 
 formation is related to the equilibrium constant via a thermodynamic 
 relation g = rt ln ka . 
 by conducting equilibrium md simulations where the free energy change 
 is obtained by counting the time that the ions spend in the associated 
 and dissociated states
this , however , requires running long explicit - solvent md simulations 
 for a given ion pair in order to obtain reasonable statistics .
both approaches imply that a definition 
 of associated / dissociated state is given . here
, we define an ion pair 
 as a configuration of a cation and anion where the two are not separated 
 more than a certain distance rcp . we take 
 this distance to be the first minimum in the cation  anion radial 
 distribution function gab(r ) = exp[pmfab(r)/rt ] ( r denoting the gas constant ) , corresponding to 
 the so - called contact pair distance . for all r > rcp
concentrations of 
 free ions and of the ion pair can be estimated by integrating the gab(r ) .
the ratio of the concentration of the ion pair to the total concentration 
 of the electrolyte ( ctot ) is equal to5and the concentration of 
 the free ions is 
 [ a ] = [ b ] = ctot  [ ab ] .
rmax is the 
 largest interparticle separation distance ( 12  in our case ) . 

the total concentration of the electrolyte is ctot = 3/(na4rmax3 ) , 
 where na is the avogadro constant . in table 3 , we compare predicted association 
 constants for ion pair formation , obtained from integrating the appropriate 
 parts of the predicted pmfs ( eq 5 ) , with the ka from reference md simulation pmfs .
we report 
 the values for all combinations of cation  anion pairs given 
 in table 1 .
prior to integration , a cubic spline 
 interpolation was used to smooth the original pmfs .
it was demonstrated 
 in previous studies that this does not significantly 
 affect the value of ka . the error in the 
 association
constant was determined as ka = |ka  ka|/2 , where ka and ka correspond to constants obtained from a pmf increased 
 and decreased by its variance , respectively . from table 3 , we see that the agreement 
 of the predicted association constants with the ones obtained by integration 
 of the simulated reference pmfs falls in most cases within the limits 
 of the estimated errors .
the predicted ka is systematically smaller than the reference one in the cases of 
 kcl , kbr , and rbbr , and larger in the case of rbi and csi .
small differences 
 in predicted and simulated pmfs accumulate in the integration process , 
 resulting in larger ka differences than 
 one might expect through visually comparing both pmfs .
also , the rmax in eq 5 should in 
 principle be  , but we were restricted by the dimensions of 
 the simulation box .
differences in pmfs at larger interparticle separation 
 distances amplify by integration and are reflected in estimated ka . in figure s6 of the supporting 
 information ,
we show predicted association constants for size 
 symmetric + 1:1 ion pairs as a function of ion sizes .
we see 
 that the predicted ka agree well with 
 constants determined from reference md simulation pmfs .
finally , 
 we explore the  law of matching water affinities  
 of k. d. collins .
in short , this is the 
 idea that small  small ion pairs are insoluble in water ( because 
 of strong electrostatic attractions ) , that large 
large ion 
 pairs are poorly soluble in water ( because they act largely like hydrophobic 
 spheres ) , but that large  small ion pairs are very soluble in 
 water . here
since 
 we can consider the cp state as a precipitated form of the salt , the ka correlates with the solubility of the salt 
 in question .
figure 7 therefore shows inverse 
 values of the experimentally determined solubilities in panel a and our predicted ka values in panel b. experiments show that nai and cscl ( one ion big 
 and the other small ) have large solubilities ( s is small ) , whereas csi and nacl ( two big and two 
 small ions ) have the smallest solubilities ( s is large ) . with respect to a given cation ,
the predicted 
 association constants show a qualitatively similar pattern as experiments : 
 for cs salts ka diminishes 
 from csi to cscl , for rb and k there is not 
 much difference with respect to the anion , while for na the constants increase going from nai to nacl . for potassium salts , 

experiments show a more significant increase in s for kcl than is reflected in the predicted ka values . for salts with a common anion ,
the 
 predicted ka values increase with increasing 
 size of the cation ( from na to cs ) .
such a 
 trend is expressed also in experimental data for i and br salts but not for cl .
it seems that the predicted ka value 
 for kcl is too small and the one for csi too large . in that respect , 
 the predicted ka for cscl does not follow 
 the law of matching water affinities .
it should be noted that these ka results are only for the set of ion parameters 
 listed in table 1 .
many other sets of ion parameters 
 have been proposed , and alternate sets may show different behavior .
figure s7 of the supporting information shows a density map of predicted ka values 
 for a systematic variation in the sizes of anion  cation pairs 
 ( from 3 to 5.5  ) .
we see that ka is the largest in the case of two big ions and the smallest for 
 the combination of one big and one small ion , in agreement with collins 
 proposed behavior .
inverse values of the experimentally determined solubilities , s ( panel 
 a ) , and ion - association constants from predicted pmfs , ka ( panel b ) , for selected alkali metal halides .
we 
 performed extensive md simulations of the pmfs between ions 
 of different sizes and charge in tip3p water . by applying an interpolation 
 scheme ( i - pmf ) on this presimulated pmf data set , we can calculate 
 pmfs accurately and quickly for an arbitrary combination of size and 
 charge of the two ions .
we tested various interpolation formulas and 
 found that the cubic spline interpolation algorithm gives the most 
 accurate representation of the majority of md results .
the accuracy 
 of the prediction is scalable based on the density of the presimulated 
 pmf grid .
improved agreement could come from a denser grid of pmfs 
 or an expanded set of sizes .
the advantage of the interpolations is 
 about a 5 order of magnitude gain in speed , without loss of predictive 
 accuracy ( within the error bars ) .
our interpolated pmfs are consistent 
 with experimentally observed ion - pair solubilities and with collins 
 law of matching affinities .
the i - pmf method may be useful where there 
 is a need to calculate accurate pmfs rapidly .
 we 
 perform extensive molecular dynamics ( md ) simulations between 
 pairs of ions of various diameters ( 25.5  in increments 
 of 0.5  ) and charge ( + 1 or 1 ) interacting in explicit 
 water ( tip3p ) under ambient conditions . 
 we extract their potentials 
 of mean force ( pmfs ) . 
 we develop an interpolation scheme , called i - pmf , 
 that is capable of capturing the full set of pmfs for arbitrary combinations 
 of ion sizes ranging from 2 to 5.5  . 
 the advantage of the interpolation 
 process is computational cost . whereas it can take 100 h to simulate 
 each pmf by md 
 , we can compute an equivalently accurate i - pmf in seconds . 
 
 this process may be useful for rapid and accurate calculation of the 
 strengths of salt bridges and the effects of bridging waters in biomolecular 
 simulations . 
 we also find that our data is consistent with collins 
  law of matching affinities  of ion solubilities : small  small 
 or large  large ion pairs are poorly soluble in water , whereas 
 small  large are highly soluble . 
Introduction
Theoretical
Methods
Results and Discussion
Conclusions