Using interpolation for fast and accurate calculation of ion-ion interactions.

We perform extensive molecular dynamics (MD) simulations between pairs of ions of various diameters (2-5.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

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.[1,2] 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 simple spherical solutes, extensive sampling is required to get converged results. Various sampling techniques have been introduced to mitigate the computational expense:[3,4] constraints,[5−8] umbrella sampling,[9] and weighted histogram methods,[10] 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[11−15] but not quantitative details.[13,15] An alternative is to use implicit-solvent models (such as Poisson–Boltzmann (PB) or generalized Born (GB)),[16−19] in which the solvent is approximated as a dielectric continuum.[19−26] 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.

Table 1

Dang Lennard-Jones Parameters for Selected Alkali Metal and Halogenide Ions with Charges of Cations at +1, Charges of Anions at −1, and ϵLJ = 0.1 kcal mol–1

cation	σLJ (Å)	anion	σLJ (Å)
Na+	2.5865	Cl–	4.4045
K+	3.3345	Br–	4.6265
Rb+	3.5305	I–	5.1705
Cs+	3.8865

However, implicit-solvent models do not capture the subtle features of PMFs, which require some representation of particulate water[29,30] (see Figure 1). To correct for the energy errors in implicit modeling, radii can be adjusted empirically, but this can sacrifice transferability. Other approaches have been applied to improving implicit-solvent models. Implicit-solvent models can allow interstitial high dielectrics (e.g., water), but they give incorrect PMFs between nonbonded atoms.[31] 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.[32] 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.,[33] and it was shown that this significantly improves the predictions. Geney et al.[34] 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.[35] 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.[36] 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. Mongan et al.[37] 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[38] 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.

Figure 1

PMFs of a model NaCl from explicit-solvent MD simulation (symbols) and implicit-solvent theories: PB (continuous blue line) and GB (dashed green line). The TIP3P water model[27] was used in MD simulations. PB calculations were done with APBS 1.3 software,[28] and GB results were obtained according to the GB model of Still et al.[23] 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).

Our goal here is an approach to computing ion–ion PMFs that is both fast and accurate. 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. Thus, i-PMF gets its physical accuracy from this precomputation stage. 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.[39−42]

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.

Theoretical Methods

Details of the Ion–Ion PMF Simulations by MD in Explicit Solvent

We performed MD simulations, mostly on the TITAN supercomputer at Oak Ridge National Laboratory,[43] using version 4.6.2 of GROMACS.[44] Each simulation consisted of two charged solute particles and 635 water molecules in a rhombic dodecahedron box. Periodic boundary conditions were applied. All simulations were performed in the isothermal–isobaric ensemble at 298.15 K and 1 atm by applying the Nosé–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)[45] 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 Å–1), 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 kcal·mol–1 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.[46−49] 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 performed simulations using the TIP3P[27] water model. Lorentz–Berthelot mixing rules were used for σLJ and ϵLJ.

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.[5−7,50] The SHAKE algorithm was used to apply the interparticle constraint. 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.[6] Uncertainties in the average mean force values were estimated from the limiting value of the block averages,[51] 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 call these the “presimulated PMFs”. 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. All results apply to systems at 298.15 K and 1 atm. The LJ energy parameter was the same for all ions (ϵLJ,ion = 0.1 kcal mol–1).

Interpolation Procedure for Capturing the Presimulated MD PMFs

Our procedure handles three independent variables. By definition, 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(x, y), that correspond to points on a two-dimensional (n × m) grid of points (x, y). 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 x ≤ xa ≤ x and x ≤ yb ≤ x; let us denote these four functional points in the following wayIn the simplest case of interpolation on the grid square—bilinear interpolation—the interpolation formula is[52]where t ∈ [0, 1] and u ∈ [0, 1] areIn 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 f′ = ∂f/∂x, f′ = ∂f/∂y, f′ = ∂f2/∂x∂y 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 form[52]and we need to find 16 coefficients c by first solving the equations for the square grid points. The function f was in our case the PMF, and x and y 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 (Q = +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, we used the function “interp2” implemented in the GNU Octave[54,53] software (version 3.6.2). 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. Figure 2 shows a cartoon representation of the proposed algorithm for PMF prediction.

Figure 2

Illustration of the i-PMF process from the bottom upward. 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.

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 2–3 s on a PC. This is a significant speed-up compared to MD simulations. 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.

Comparison of Different Interpolation Schemes

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

Figure 3

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). Symbols denote the test MD results obtained independently of the interpolation scheme. The cubic spline interpolation algorithm best captures the magnitude and shape of PMF features.

Table 2

Comparison of the Predicted Positions and Values of the First and Second Minimum and the First Maximum of the PMF for Potassium Iodide Obtained with MD and the i-PMF Method Using Nearest-Neighbor (NN), Bilinear (BL), Piecewise Cubic Hermite Polynomial (PC), and Cubic Spline Interpolation with Smooth First and Second Derivatives (CS)

	first minimum		first maximum		second minimum
method	r (Å)	PMF (kcal mol–1)	r (Å)	PMF (kcal mol–1)	r (Å)	PMF (kcal mol–1)
MD	3.6	–1.22(3)	4.6	0.37(3)	6.0	–0.36(2)
NN	3.6	–1.38(3)	4.7	0.28(3)	6.0	–0.33(3)
BL	3.6	–0.92(4)	4.6	0.35(3)	6.0	–0.31(3)
PC	3.6	–1.16(3)	4.6	0.40(3)	6.0	–0.33(3)
CS	3.6	–1.21(4)	4.6	0.39(3)	6.0	–0.34(3)

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 χ2 parameter (see the Supporting Information file for the definition). Smaller χred2 values indicate better agreement between the predicted and simulated PMF. The values of χred2 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 χred2 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 103 kcal·mol–1 (and larger), but the error in the case of +1:–1 ion pairs is small (∼10–2 kcal·mol–1). Even though the differences between the predicted and simulated PMFs are only around 1 kcal·mol–1, the χred2 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 χred2 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 χred2, then the value of χred2 for bilinear interpolation lowers to 6.1. This smaller value compared to the χred2 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 χred2 = 1 suggests that the extent of the agreement between observations and estimates is in accord with the error variance. χred2 > 1 indicates that the fit has not fully captured the data (or that the error estimate has been underestimated), while χred2 < 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 S1–S4 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 our studies are limited. 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. The present approach could readily be generalized but at substantial computational expense. 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.

Results and Discussion

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.

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. We selected this example to demonstrate that the interpolation method works in “extreme” 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 χred2 values for the cases displayed in Figure 4 are given in Table S2 of the Supporting Information file. In the cases where χred2 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, the cubic spline interpolation gives the best results, as already discussed.

Figure 4

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

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 χred2 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.

Figure 5

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.

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.[55−64] Again, we find excellent agreement between the interpolation formulas and the MD simulations; see Figure 6. We see the expected charge–charge repulsion. Interestingly, because the solvent is water, the repulsions are very weak. Examples of PMFs for other like-charged pairs of ions from Table 1 are given in Figure S4 of the Supporting Information file.

Figure 6

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.

Predicting Association Constants for Ion Pairing by Using Interpolated PMFs

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. This quantity is related to the equilibrium constant for the process.

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⊖)−1·([B–]/c⊖)−1, where [X] denotes the molar concentration of species X and c⊖ = 1 mol dm–3 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, one can calculate Ka. This, however, requires running long explicit-solvent MD simulations for a given ion pair in order to obtain reasonable statistics.[65−70] An alternative is to calculate the association constants from the radial distribution function. 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, we treat ions as free. Concentrations of free ions and of the ion pair can be estimated by integrating the gAB(r).[50] The ratio of the concentration of the ion pair to the total concentration of the electrolyte (ctot) is equal toand 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/(NA4πRmax3), 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[50,71] 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.

Table 3

Association Constants of Selected Salts Obtained from the Predicted (Kai-PMF) and from MD Simulation (KaMD) PMFs with the Estimated Error of the Last Digit Reported in Parentheses

	Cl–		Br–		I–
	Kai-PMF	KaMD	Kai-PMF	KaMD	Kai-PMF	KaMD
Na+	0.148(7)	0.141(4)	0.118(8)	0.112(3)	0.061(3)	0.064(2)
K+	0.40(1)	0.435(9)	0.40(1)	0.446(9)	0.42(1)	0.418(8)
Rb+	0.47(1)	0.49(1)	0.52(1)	0.55(1)	0.59(2)	0.55(1)
Cs+	0.65(2)	0.64(2)	0.78(2)	0.77(1)	0.99(2)	0.91(2)

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.[72] 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, we make a simple test using our interpolated PMFs. 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[73] 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–1 is small), whereas CsI and NaCl (two big and two small ions) have the smallest solubilities (S–1 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–1 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 7

Inverse values of the experimentally determined solubilities,[73]S–1 (panel a), and ion-association constants from predicted PMFs, Ka (panel b), for selected alkali metal halides.

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.

Conclusions

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.