Theory of Ionic Hydration: Insights from Molecular Dynamics

Subsequently, Latimer et al.15 showed that the experimental hydration free energies of alkali and halide ions could be fitted by eq 1 with R = Rion + ...
0 downloads 0 Views 144KB Size
Download PDF
7958

J. Phys. Chem. B 1999, 103, 7958-7968

Theory of Ionic Hydration: Insights from Molecular Dynamics Simulations and Experiment C. Satheesan Babu† and Carmay Lim*,†,‡ Institute of Biomedical Sciences, Academia Sinica, Taipei 11529, Taiwan R.O.C, and Department of Chemistry, National Tsing Hua UniVersity, Hsinchu 300, Taiwan R.O.C. ReceiVed: June 28, 1999

The Born model of ionic solvation prescribes the solvation free energy of ions in a dielectric continuum and is widely used. Despite the fact that the solvent molecules around the ions are highly structured giving rise to molecular phenomena such as dielectric saturation and electrostriction, the Born model yields accurate results compared to experimental data if a proper radius for the ion, the effective Born radius, Reff, is chosen. On the basis of molecular dynamics simulations of ions of varying charge and size in TIP3P water, Reff is identified as the mean of the ionic radius (Rion) and the first peak position in the ion-water pair distribution function (Rgmax). This relationship was verified using experimentally available Rgmax for ions of varying size (0.4-3.0 Å) and charge (-3e to +4e). The new interpretation of Reff implies that the correct expression for the solvation free energy is a combination of two Born free energies, which incorporate the effects of dielectric saturation and electrostriction implicitly. The new solvation free energy formula was used to derive unknown Rgmax for a number of ions whose free energy data are available as well as new expressions for the solvation entropies and enthalpies. An empirical relationship between Rgmax and Rion is established. Potential applications of the findings of this work to free energy simulations, solvation in different solvents, solvation of molecular systems, and biomolecules are discussed.

Introduction Continuum dielectric theory has been applied to many problems in solution chemistry and biophysics.1,2 The Born model,3 proposed in 1920, represents the simplest continuum theory of ionic solvation. It gives the excess solvation free energy of a spherical ion of radius R and charge qe (where e is the magnitude of the electronic charge) in a continuum solvent medium of dielectric constant ; i.e.,

∆GBorn )

-q2e2 1 12R 

[

]

(1)

A key assumption in deriving the Born formula (eq 1) is that the response of the solvent molecules to the charge distribution of the solute is linear; this gives rise to the quadratic dependence of the Born free energy on the solute charge. However, it is well-known that solvent dipoles in the first solvation shell of the ion are highly structured due to the electric field of the ion, and thus the dipoles around it need not respond linearly.4-9 This gives rise to phenomena such as dielectric saturation and electrostriction.10-12 Dielectric saturation arises from the orientational immobilization of solvent dipoles around the ion, thereby decreasing the dielectric permittivity of the solvent in the immediate vicinity of the ion, which in turn decreases the solvation free energy.4 On the other hand, electrostriction stems from an increase in the solute-solvent electrostatic interaction energy near the ion, thereby decreasing the excess volume of the solute, and thus increasing the solvation free energy.5,9,11,12 Both these molecular effects, which appear to be manifestations of the same microscopic phenomenon, act in opposite directions and contribute to the nonlinear response of the solvent.5,9 These † ‡

Academia Sinica. National Tsing Hua University.

solvent molecular effects become particularly important when estimating free energies of transfer of ions from one solvent to another (e.g., transfer from a hydrogen-bonding solvent like water to a non-hydrogen-bonding solvent such as propylene carbonate, and transfer from one solvent isotope to another, e.g., from H2O to D2O).11 Despite the seemingly crude assumptions inherent in the Born model, it can reproduce experimental solvation free energies with a proper choice of the radius R in eq 1. This radius, termed the “effective Born radius (Reff)”, is larger than the crystal ionic radius (see below). There have been many attempts aimed at understanding how a continuum representation of the solvent in the Born model can treat solvation accurately simply by adjusting the radius R. An early attempt in understanding the molecular origin of the effective Born radius originated with the work of Voet,13 who showed that the experimental hydration enthalpies of alkali and alkaline earth metal ions as well as closed-shell trivalent ions were consistent with eq 1 if R ) Rion+ ∆, where Rion is the Goldschmidt ionic radius14 and the empirical constant ∆ is 0.85 Å for the alkali cations but 0.1 Å for the cations studied. The distance Rion+ ∆ was interpreted as the distance from the center of a cation to the center of the electric dipole in a neighboring water molecule. Subsequently, Latimer et al.15 showed that the experimental hydration free energies of alkali and halide ions could be fitted by eq 1 with R ) Rion + ∆, where Rion here is the Pauling ionic radius and the empirical constant ∆ is 0.85 Å for alkali cations but 0.1 Å for halide anions. They interpreted Rion + ∆ as the radius of the cavity formed by the water dipoles around the ion; i.e., for cations, it is the ion-oxygen distance while for anions, it is the ionhydrogen distance of a neighboring water molecule (Figure 1a,b). The origin of the Rion + ∆ radius obtained by Voet13 and Latimer et al.15 can be understood within the framework of the

10.1021/jp9921912 CCC: $18.00 © 1999 American Chemical Society Published on Web 08/27/1999

Theory of Ionic Hydration

J. Phys. Chem. B, Vol. 103, No. 37, 1999 7959

Figure 1. (a) Definition of the ionic and solvation regions around a cation. Rgmax is the distance from the center of the ion to the peakposition in the ion-oxygen rdf. (b) Same for an anion with Rgmax defined as the distance from the center of the ion to the peak-position in the ion-hydrogen rdf. (c) The dielectric constant as a function of the distance from the ion, represented by a step function (solid line) or a continuous (dotted curve) function (see text for details). The solvent molecules in the first shell are represented by dipoles, R1gmax is the radius corresponding to the peak position in the first solvation shell, and R2gmax corresponds to that for the second solvation shell; 0, 1, and 2 represent the dielectric constants in various solvation regions.

Debye-Pauling model, which was originally introduced to improve the Debye-Hu¨ckel theory of electrolytes.11,16 In the Debye-Pauling model, the solvation free energy of an ion of charge qe is given by11

∆GDP )

[

(

) (

)]

1 -q2e2 1 1 1 1 1 2 Rion 0 Rion Rion + ∆  Rion + ∆ (2)

Here, the first term arises from the free energy about the ion in the gas phase. The second term is for the free energy in the annular space between Rion and Rion+ ∆ which has a lower dielectric constant 0 compared to the bulk value of . The last term corresponds to the free energy outside a sphere of radius Rion+ ∆, in a solvent medium with dielectric constant . When 0 ) 1, the dielectric permittivity of free space, eq 2, becomes

∆GDP )

-q2e2 1 1 2(Rion + ∆)

[

]

(3)

Equation 3 is in accord with the observations of Voet and Latimer et al. It implies that the region around the spherical ion with a thickness of ∆ has such extreme dielectric saturation that 0 attains the dielectric permittivity of free space. Further

extensions of this model involved refinement of the dielectric response in the solvation shells and assumed a distancedependent dielectric function instead of a step function as in eq 3.17-21 Figure 1c illustrates the various functional forms of the dielectric constant  for the Born model, the Debye-Pauling model, and refinement of these models. Rashin and Honig22 gave an alternative interpretation of the effective Born radius. They defined Reff as the radius of a spherical cavity that contains negligible solvent charge density. On the basis of electron density distributions in ionic crystals, they concluded that covalent radii for cations and ionic radii for anions provided a reasonable measure of the cavity radius that was introduced by Latimer et al.15 When these measures of cavity size were used in eq 1, calculated and experimental solvation enthalpies for 31 ions ranging in charge from -1 to +4 agree within 7.5% if all radii were empirically increased by 7%.22 On the other hand, from liquid-state statistical mechanical theories23,24 and computer simulation studies of model systems for monovalent ions,4,25 Reff was identified with Rgmax, the first peak position of the ion-solvent radial distribution function (rdf) (Figure 1). The “cavity” defined by Rgmax contains some solvent charge density and thus differs from that defined by Rashin and Honig, except for a hard-sphere ion in a hard-sphere monatomic solvent, where the cavity, defined by the sum of the hard-sphere radii of the ion and solvent, contains zero solvent density. The interpretation that Reff ) Rgmax can, to some extent, explain the different solvation properties of various solvents. It also accounts for the asymmetry in the solvation free energies of cations and anions of equal charge magnitude and ionic size,11 since the orientation of the solvent dipoles for negative ions is opposite to that for positive ions (Figure 1). In contrast, the Born expression (eq 1) with Reff ) Rion predicts equal solvation free energies of positive and negative ions of equal charge and size. However, it is known from neutron diffraction experiments and simulation studies of ion-solvent rdfs that Rgmax is considerably longer than the Rion + ∆ obtained by Voet and Latimer et al.21,25,26 Furthermore, from the Debye-Pauling model and the empirical analyses, it is possible that two distances (as opposed to only one distance in Reff ) Rgmax) are involved in the free energy formula. On the other hand, the analyses of Voet and Latimer et al. failed to account for the solvation free energies in non-hydrogen bonding solvents and solvent isotopes.11 Although analyses incorporating Rion + ∆ or a distance-dependent dielectric constant take into account the effects of dielectric saturation, it is not clear how they account for electrostriction. In this work, Reff is identified as the mean of two radii; viz., Rion and Rgmax; i.e.,

Reff ) (Rion+ Rgmax)/2

(4)

This new interpretation was based on simulation studies of ions whose potential energy parameters were derived from experimental solvation free energies using free energy perturbation simulations. The ions studied include alkali metal cations, alkaline earth metal dications, and hypothetical mono- and dinegative ions. The derivations of the ion-water parameters and the simulation methodology are outlined in the Methods. The structural features and energetics of the solvation shells around the ions and the new interpretation of Reff (eq 4) and the solvation free energy are presented in the Results. Equation 4 is verified for “realistic” ions using experimental Rgmax for ions of varying size (0.4-3.0 Å) and charge ranging from -3e to +4e. Furthermore, an empirical relationship between Rion and

7960 J. Phys. Chem. B, Vol. 103, No. 37, 1999

Babu and Lim

TABLE 1: Potential Energy Parameters, Computed Hydration Free Energies, and First-Shell Coordination Numbers for Alkali and Alkaline Earth Metal Ions in TIP3P Water Based on Free Energy Perturbation Simulations ion Li+

Na+ K+ Rb+ Cs+ Mg2+ Mn2+ Ca2+ Sr2+ Ba2+ CsMg2-

-10-2min,a kcal/mol 1.828 0.277 0.033 0.017 0.008 87.51 61.03 44.97 11.82 4.710 0.008 87.51

Rmin,a Å

-∆Gcalc,b kcal/mol

-∆Gexpt,c kcal/mol

1.137 1.868 2.658 2.956 3.395 0.787 0.920 1.326 1.741 2.124 3.395 0.787

122.2 ( 98.8 ( 0.8f 80.9 ( 1.0e 75.5 ( 0.9e 67.7 ( 0.7e 455.2 ( 0.3f 437.9 ( 0.3f 380.8 343.3 ( 0.2f 313.0 ( 0.2f 87.7 ( 0.2f 776.9 ( 0.8f

122.1 98.2 80.6 75.5 67.8 455.5 437.8 380.8 345.9 315.1

0.6e

nd 4.5 5.9 6.6 7.0 8.0 6.0 6.0 8.0 8.5 9.1 7.0 6.0

a min and Rmin are single-ion Lennard-Jones parameters; the corresponding parameters for pair interactions are derived using the arithmetic combining rules, as in the CHARMM program.31 b The computed free energies for alkali metal ions are from Åqvist,27 whereas those for alkaline earth metals and hypothetical ions are from this work. c Experimental values from Burgess.29 d Coordination numbers computed from eq 5 up to the first minimum in the ion-oxygen rdf. e For SPC water.27 f ForTIP3P water.

Rgmax is established. The new solvation free energy formula is used to derive experimentally unavailable Rgmax and new expressions for the solvation enthalpies and entropies. A qualitative explanation of why eq 4 works is presented in the Discussion, which also highlights potential applications of eq 4 to the solvation of polyatomic systems. The key results of this work are summarized in the Conclusion. Methodology Ion-Water Parameters. The ions studied include alkali metal ions (Li+, Na+, K+, Rb+, Cs+), alkaline earth metal ions (Mg2+, Ca2+, Sr2+, Ba2+), Mn2+, as well as two hypothetical negative ions, Cs- and Mg2-. To make a connection between the size of the ion and its associated first solvation shell with the experimental solvation free energy, it is important that the ion-water potential model reflects the bulk solvation free energy. Åqvist27 has pioneered the application of free energy perturbation simulations to derive ion-water potential parameters that reproduce the experimental solvation free energies and ion-solvent structures of alkali and alkaline earth metal ions in SPC28 water. The experimental solvation free energies in Åqvist’s parametrization were taken from Burgess,29 who employed a proton solvation free energy of -260.5 kcal/mol. We have transferred Åqvist’s parameters to the TIP3P30 water model and computed the free energies of divalent cations relative to Ca2+ in TIP3P water. The excellent agreement between the computed and experimental free energies of divalent cations (this work) and monovalent ions (Åqvist’s work27) supports the transferability of Åqvist’s parameters to other solvent models. For the hypothetical negative ions, the van der Waals parameters for Mg2+ and Cs+ were used for Mg2- and Cs-, respectively. Table 1 lists the ion-water potential energy parameters, the computed hydration free energies, and the corresponding experimental values.29 Simulation Protocol. Constant volume MD simulations of the aforementioned ions with charge ranging from -2 to +2 in TIP3P water were performed using the CHARMM program31 at a mean temperature of 300 K. Each ion was fixed at the center of a 20 Å radius water sphere, constructed from an equilibrated

configuration of TIP3P water molecules at a pressure of 1 atm and an experimental density of 0.0334 molecules/Å3. The boundary water molecules were subjected to a modified mean field potential, which accounts for interactions with infinite bulk water not treated explicitly in the simulation.32 All atoms were propagated according to Newton’s equations using a Verlet integrator and a time step of 2 fs. Nonbonded interactions were truncated by applying an atom-based force-switching function in the region between 17 and 19 Å, and a nonbond cutoff at 20 Å (the radius of the simulation system). The nonbond pair list was updated every 10 steps. Each ion-water system was equilibrated for 20 ps and subjected to 100 ps of production dynamics, during which coordinates were stored every 0.1 ps for analysis. Relative hydration free energies of Cs- and Mg2were computed by free energy simulations using Cs+ and Mg2+, respectively, as the starting point (Table 1). The free energy simulations employed a thermodynamic integration protocol31 with 11 windows for the sampling. At each window the perturbation energies were collected for a 20 ps period after allowing 5 ps for equilibration of each window point. The average free energies were obtained from forward and reverse perturbation runs (Table 1). Analysis. We define the radius of the first solvation shell as the distance from the ion to the first minimum in the ion-oxygen distribution functions. This differs from Rgmax, which corresponds to the first peak in the rdf. The running coordination numbers were computed from the corresponding ion-solvent rdfs and are defined as

nX(R) ) 4πFX

∫0RgX(r)r2 dr

(5)

where FX is the bulk density of species X (i.e., O or H). The electrostatic solvation energies were computed from the ionoxygen and ion-hydrogen rdfs using the following equation:33

∫0RuelX(r) gX(r)r2 dr

ESolv X (R) ) 4πFX

(6)

In eqs 5 and 6, X denotes either O or H and uelX is the electrostatic energy of the ion with X, i.e., uelX ) e2qIqX/rIX where rIX is the distance between the ion and X, qI is the ionic charge, and qX is the charge on O or H; for TIP3P water,30 qH ) 0.417 and qO ) -0.834. The ion-solvent energies (solvation energies) were computed as the sum of ion-hydrogen and ionoxygen energies. Results Ion-Solvent Structure. Figure 2 shows the ion-oxygen rdfs for some selected monovalent ions (Li+, Na+, and Cs+), divalent ions (Mg2+ and Ca2+), and two hypothetical negative ions, Mg2and Cs-. As expected, the water molecules are strongly oriented around divalent ions with a peak height of 26 for Mg2+ and 16 for Ca2+. The monocations have lower peak heights, broader solvation shells, and less prominent second solvation shells compared to the divalent ions. The structural features for the mono- and divalent ions using the new potential energy parameters are similar to those found in earlier simulation studies of ionic solutions.34-37 In particular, the peak positions in the rdfs are close to those observed in the SPC water simulations.27 A feature that distinguishes the Mg2--O and Cs--O rdfs (Figure 2c) from the corresponding Mg2+-O (Figure 2b) and Cs+-O rdfs (Figure 2a) is the sharper and higher solvation peaks for the anions compared to the respective cations. Moreover, the peak positions of the Mg2--O and Cs--O rdfs are shifted toward shorter distances compared to the peak positions of the

Theory of Ionic Hydration

Figure 2. Ion-oxygen radial distribution functions for (a) Li+, Na+, and Cs+, (b) Mg2+ and Ca2+, and (c) hypothetical negative ions Mg2and Cs-. The numbers in (b) and (c) correspond to the peak height of the first peak.

Figure 3. First solvation shell of dipositive and dinegative magnesium ions.

Mg2+-O and Cs+-O rdfs. The observed differences reflect the difference in packing of the first-shell water molecules around the ion. To illustrate this, snapshots of representative configurations of Mg2+ and Mg2- are shown in Figure 3. For Mg2-, one hydrogen from each water molecule in the first hydration shell points to the ion, and hydrogen bonds with the dianion, whereas the other water hydrogen forms hydrogen bonds with the second-shell water molecules. Thus, the negative ion seems to promote more structure compared to its positively charged counterpart. This observation has also been found in previous works.6,25,46-48 First-Shell Structure and Energetics. Since the first solvation shell dictates the value of Rgmax, the structural and

J. Phys. Chem. B, Vol. 103, No. 37, 1999 7961

Figure 4. First-shell solvation properties of alkali metal ions: (a) Radial distribution functions for Li+, Na+, K+, Rb+, and Cs+; (b) running coordination numbers n(r) computed from the ion-oxygen rdfs; (c) solvation energies computed from the sum of ion-hydrogen and ion-oxygen energies according to eq 6.

energetic features of the first shell were analyzed. The rdfs, running coordination numbers, and solvation energies of the first solvation shells for alkali monocations, alkaline earth dications, and hypothetical monoanions are shown in Figures 4-6, respectively. The coordination numbers for the first hydration shell of the ions are given in Table 1. The first hydration shell of Li+, which extends from 1.8 to 2.5 Å, has a coordination number of 4.5, suggestive of a tetrahedral-like arrangement, and the minimum in the solvation energy for Li+ occurs at 2.35 Å, the outer-edge of the first solvation layer (Figure 4). In contrast, the first hydration shell of Na+ is octahedral, and its larger coordination number (6) relative to Li+ may account for the similar ion-solvent energy minima computed for Li+ and Na+. As the ionic size increases from Na+ to Cs+, the width of the rdfs increases, the first solvation shell becomes more diffuse, and the magnitude of ion-solvent energy minimum decreases, reflecting the weaker ion-solvent interactions. The rdfs of divalent cations (Figure 5) are sharper than those of monocations (Figure 4) owing to the higher electric field of the dications. Mg2+ and Mn2+ exhibit similar solvation properties, which are distinct from the other dications. Both ions form a regular octahedral structure, and the energetic features of the first solvation shells of Mg2+ and Mn2+ show a flat valley corresponding to the region between the first and second solvation shells. The unique first-shell solvation structural

7962 J. Phys. Chem. B, Vol. 103, No. 37, 1999

Figure 5. Same as in Figure 4 for Mg2+, Mn2+, Ca2+, and Ba2+.

features of Mg2+ and Mn2+ may play an important role in controlling the hydrolytic process. For example, most enzymes that require Mg2+ as a cofactor are also active with Mn2+. Ca2+, rather than the smaller Mg2+ or Mn2+, has the most favorable first-shell solvation energy among all the dications, due partly to its larger first-shell coordination number compared to Mg2+ and Mn2+ (see Table 1). This correlates with the role of Ca2+ as a stabilizing ion in calcium-binding enzymes such as calmodulin. For the hypothetical negative ions the distribution functions plotted in Figure 6a are the ion-hydrogen rather than ionoxygen rdfs. The first-shell solvation structures of Mg2- and Cs- differ from those of Mg2+ and Cs+, respectively, as noted above. Although both Mg2+ and Mg2- ions have a first-shell coordination number of 6, the running coordination number rises sharply at 1.7 Å for Mg2- (Figure 6b) whereas it rises more gradually at a longer distance of 1.9 Å for Mg2+ (Figure 5b). This indicates that water molecules are more strongly attracted to Mg2- than to the Mg2+ ion (see Figure 3 and above). Consequently, the first-shell solvation energy of Mg2- (Figure 6c) is more negative than that of Mg2+ (Figure 5c) even though both ions have the same ionic size. The solvation free energy for Mg2- is also significantly more negative (by 321 kcal/mol) than that for Mg2+ (Table 1). This is in sharp contrast to the Born formula (eq 1), which predicts equal solvation free energies for Mg2- and Mg2+. For Cs-, the ion-solvent distribution is more liquidlike with a diffuse first shell, indicating weak solvation forces. As observed for Mg2-, the first-shell solvation energy of Cs- is more negative than its oppositely charged

Babu and Lim

Figure 6. Same as in Figure 4 for Mg2- and Cs- except that the rdfs in (a) are ion-hydrogen rdfs.

counterpart, Cs+, and the solvation free energy for Cs- is more negative (by 20 kcal/mol) than that for Cs+ (Table 1). Effective Born Radius and Its Relationship to Rion and Rgmax. Experimental effective Born radii, ReffE, were derived from experimental hydration free energies (Table 1) using eq 1 with  ) 78.5 for water at 298 K. Rgmax for cations and anions were computed from the first peak positions of the ion-oxygen and ion-hydrogen rdfs, respectively (see Figures 4-6). These two sets of radii are listed in Table 2 along with the Goldschmidt ionic radii,14 RionG. Also shown in column 4 of Table 2 are ReffC, which are calculated by taking the mean of RionG and the computed Rgmax radii (see eq 4). The ReffC and ReffE values are in excellent agreement for the ions studied here. These results suggest that R in eq 1 can be replaced by (Rion + Rgmax)/2, so that the corrected solvation free energy is given by

∆GCorr )

-q2e2 1 1Rion + Rgmax 

[

]

(7)

Solvation free energies computed using eq 7 (column 6, Table 2) agree with the experimental free energies (column 5, Table 1) to within 3.4%, with Li+ showing the largest percentage deviation. To further verify the internal consistency of eqs 6 and 7, the absolute solvation free energy of the proton, ∆Gabs(H+), was calculated for all the cations studied using the corrected absolute solvation free energies (column 6, Table 2) and experimental

Theory of Ionic Hydration

J. Phys. Chem. B, Vol. 103, No. 37, 1999 7963

TABLE 2: Computed Solvation Thermodynamics of Ions ion

RionG,a Å

Rgmaxb Å

ReffC,c Å

ReffE,d Å

∆Gabs,e kcal/mol

∆Grel,f kcal/mol

∆Gabs(H+),g kcal/mol

Mg2+

0.78 0.91 1.06 1.27 1.43 0.78 0.98 1.33 1.49 1.65 1.65 0.78

2.06 (2.09 ( 0.04) 2.08 (2.19 ( 0.01) 2.39 (2.42 ( 0.05) 2.57 (2.64 ( 0.04) 2.73 (2.90) 2.00 (2.08 ( 0.07) 2.39 (2.36 ( 0.06) 2.67 (2.80 ( 0.08) 2.81 (2.89 ( 0.10) 3.03 (3.14 ( 0.08) 2.19 0.86

1.42 1.50 1.73 1.92 2.08 1.39 1.69 2.00 2.15 2.34 1.92 0.82

1.44 1.50 1.72 1.90 2.08 1.34 1.67 2.03 2.17 2.42 1.86 0.84

-461.6(1.3) -438.5(0.2) -380.0(0.2) -341.4(1.3) -315.2(0.0) -117.9(3.4) -97.3(1.0) -81.9(1.7) -76.2(1.0) -70.0(3.1) -85.4(2.7) -799.4(2.9)

65.5 83.2 140.2 175.1 205.9 138.4 162.3 179.9 185.0 192.7

-263.6 -260.8 -260.1 -258.3 -260.5 -256.3 -259.6 -261.8 -261.2 -262.7

Mn2+ Ca2+ Sr2+ Ba2+ Li+ Na+ K+ Rb+ Cs+ CsMg2a

Goldschmidt ionic radius from Cotton and Wilkinson.14 b Computed from the rdfs in this work; values in parentheses are experimental values from Marcus.26 c ReffC ) (RionG + Rgmax)/2, using computed Rgmax in column 3. d Computed from ∆Gexpt in Table 1 using eq 1 with  )78.5. e Computed from eq 7; values in parentheses are percentage deviations of the computed free energies from experimental values. f Relative experimental solvation free energies from Rosseinsky.39 g Computed from eqs 8 and 9.

Figure 7. Plot of q2/2X versus q2/2ReffE, where X is ReffC (open circles), Rion (filled diamonds), or Rgmax (filled circles). The solid line corresponds to a plot of q2/2ReffE versus itself. The data points from left to right are Cs+, Rb+, K+, Cs-, Na+, Li+, Ba2+, Sr2+, Ca2+, Mn2+, Mg2+, and Mg2-.

relatiVe (or conventional) solvation free energies ∆Grel39 (column 7, Table 2) using the following equations:

∆Grel(Mn+) ) ∆Gabs(Mn+) - n∆Gabs(H+)

(8)

∆Grel(Xn-) ) ∆Gabs(Xn-) + n∆Gabs(H+)

(9)

The computed proton solvation free energy ∆Gabs(H+) (column 8, Table 2) is constant to within 3% with a mean value of -260.5 kcal/mol. Note that this is the same value employed to obtain the experimental absolute free energies29 in Table 1 (see Methods). This observation further establishes the accuracy of the relationship in eq 7. Equation 7 shows that the solvation free energy is linear with respect to 1/(Rion+ Rgmax), and not 1/Rion or 1/Rgmax. This relationship and the accuracy of eq 7 in predicting solvation free energies are illustrated in Figure 7, which shows a plot of q2/2X (where X ) ReffC, Rion, or Rgmax) vs q2/2ReffE for all the ions studied. The solid line in Figure 7 represents the “experimental” reference line, i.e., a plot of q2/2ReffE vs itself. It is remarkable that the q2/2ReffC data points (open circles) lie exactly on the “experimental” line. In contrast both the q2/2Rion (filled diamonds) and q2/2Rgmax (filled circles) cation data points deviate significantly from the experimental line and the deviations magnify with increasing q2/2ReffE. The latter implies that using Rion or Rgmax in the Born expression (eq 1) would give erroneous free energies for small, multivalent cations. Nevertheless, the q2/2Rgmax values are closer to the experimental line than the respective q2/2Rion points, indicating that Rgmax accounts

for experimental solvation free energies better than Rion. Indeed, Rgmax was identified as the effective Born radius for monovalent ions in earlier works.23-25 For negative ions, represented by hypothetical Mg2- and Cs-, Rion is close to Rgmax since for anions Rgmax is determined by the first peak of the ion-hydrogen rdf. This is consistent with the observation that adding a relatively small empirical constant of 0.1 Å to the ionic radii15 reproduced the experimental hydration free energies of halide ions. But even in the case of negative ions, the average of Rion and Rgmax is closer to ReffE than either Rion or Rgmax (see Table 2). Thus, the actual dielectric boundary between the ion and the solvent starts not at Rion or Rgmax as originally thought but at (Rion + Rgmax)/2. Extension to Realistic Systems. The above analyses were made for ion-water models that reproduced experimental hydration free energies. Since Rgmax can be experimentally derived from X-ray or neutron diffraction and EXAF measurements,26,38 the applicability of eq 4 for realistic systems can be verified. Table 3 lists Rgmax as well as Pauling (RionP) and Goldschmidt (RionG) ionic radii14,26,39 for ions (including polyatomic ones) of varying charge (+4e to -3e) and size (0.43.0 Å). The Rgmax radius in Table 3 is an average of several ion-water distances determined from X-ray and neutron diffraction experiments as well as computer simulation studies of model systems.26 The Pauling and Goldschmidt ionic radii are identical for halide anions but differ by as much as 0.14 Å for cations (the largest deviation occurring for La3+). Two sets of experimental free energies are given in Table 3. One set, ∆GexptM, is taken from Marcus,21 based on relative ∆Grel values in the NBS compilations40 and ∆Gabs(H+) ) -252.4 kcal/mol. The other set, ∆GexptFK, is taken from Friedman and Krishnan11 based on relative ∆Grel values taken mainly from the compilations of Rosseinsky39 and Latimer41 and ∆Gabs(H+) ) -260.5 kcal/mol. Apart from ∆Grel and ∆Gabs(H+) differences between the two sets of experimental free energies, another difference is that the ∆GexptM values have been corrected by 1.9 kcal/mol to account for the compression of the space available to the ion on its transfer from the gaseous to aqueous state.21,42 Using ∆GexptM and ∆GexptFK in eq 1 with  ) 78.5, effective Born radii ReffE,M and ReffE,FK were obtained, respectively. The experimental Reff values are in close agreement with the ReffC,P and/or ReffC,G radii computed using eq 4 with Rgmax and RionP or RionG, respectively (see Table 3). Cr3+ shows the largest percentage difference between experimental and computed Reff values (16%) due to its relatively short Rgmax (1.97 Å). The accuracy of eq 7 in predicting solvation free energies of (realistic) ions is illustrated in Figure 8. The solid line in Figure

7964 J. Phys. Chem. B, Vol. 103, No. 37, 1999

Babu and Lim

TABLE 3: Comparison of Effective Born Radii Computed Using Experimental Rgmax with Two Sets of Experimental Effective Born Radiia ion

Rgmax, Å

Rion,P Å

RionG, Å

Li+ Na+ K+ Rb+ Cs+ Ag+ Be2+ Mg2+ Ca2+ Sr2+ Ba2+ Mn2+ Fe2+ Co2+ Ni2+ Cu2+ Zn2+ Cd2+ Hg2+ Sn2+ Al3+ Y3+ La3+ Ce3+ Pr3+ Nd3+ Sm3+ Eu3+ Gd3+ Tb3+ Dy3+ Er3+ Tm3+ Lu3+ Cr3+ Fe3+ In3+ Tl3+ Th4+ FClBrINO3ClO4SO42SeO42SeO42PO43-

2.08 2.36 2.80 2.89 3.14 2.42 1.75 2.09 2.42 2.64 2.90 2.19 2.11 2.11 2.06 2.11 2.10 2.30 2.42 2.33 1.89 2.37 2.53 2.55 2.54 2.47 2.45 2.45 2.39 2.40 2.37 2.36 2.36 2.34 1.97 2.03 2.16 2.23 2.53 1.68 2.24 2.42 2.70 2.21 2.75 2.87 2.87 3.00 2.11

0.74 1.02 1.38 1.49 1.70 1.15 0.40 0.72 1.00 1.25 1.36 0.83 0.78 0.75 0.69 0.73 0.75 0.95 1.02 0.93 0.53 1.01 1.18 1.14 1.14 1.12 1.09 1.07 1.06 1.04 1.03 1.00 0.99 0.97 0.62 0.65 0.79 0.88 1.06 1.33 1.81 1.96 2.20 1.79 2.40 2.30 2.30 2.43 2.38

0.78 0.98 1.33 1.49 1.65 1.13 0.34 0.78 1.06 1.27 1.43 0.91 0.83 0.82 0.68 0.72 0.69 1.03 0.93 0.45 0.90 1.04

0.53 0.53 0.81 0.91 1.33 1.81 1.96 2.20

-GexptM, kcal/mol

-GexptFK, kcal/mol

113.5 87.2 70.5 65.7 59.8 102.8 572.4 437.4 359.7 329.8 298.8 420.7 439.8 457.7 473.2 480.4 467.3 419.5 420.7 356.1 1081.5 824.6 751.7 764.8 775.6 783.9 794.7 803.1 806.6 812.6 818.6 835.3 840.1 840.1 958.4 1019.4 951.2 948.9 1389.8 111.1 81.3 75.3 65.7 71.7 102.8 258.1 258.1 215.1 660.9

123.5 98.3 80.8 76.6 71.0 114.5 582.3 455.5 380.8 345.9 315.1 437.8 456.4 479.5 494.2 498.7 484.6 430.5 436.3 371.4 1103.3

1037.0 1035.5 973.2 975.9 103.8 75.8 72.5 61.4

ReffC,P, Å

ReffC,G, Å

ReffE,M, Å

ReffE,FK, Å

1.41 1.69 2.09 2.19 2.42 1.78 1.08 1.40 1.71 1.95 2.13 1.51 1.45 1.43 1.38 1.42 1.42 1.63 1.72 1.63 1.21 1.69 1.85 1.85 1.84 1.80 1.77 1.76 1.73 1.72 1.70 1.68 1.68 1.66 1.29 1.34 1.47 1.56 1.80 1.51 2.02 2.19 2.45 2.00 2.58 2.58 2.58 2.72 2.25

1.43 1.67 2.06 2.19 2.39 1.77 1.05 1.44 1.74 1.96 2.17 1.55 1.47 1.46 1.37 1.42 1.39 1.67 1.68

1.44 1.88 2.32 2.49 2.74 1.59 1.15 1.50 1.82 1.99 2.19 1.56 1.49 1.43 1.39 1.36 1.40 1.56 1.56 1.84 1.36 1.79 1.96 1.93 1.90 1.88 1.86 1.84 1.83 1.82 1.80 1.77 1.76 1.76 1.54 1.45 1.55 1.55 1.89 1.47 2.02 2.18 2.49 2.29 1.59 2.54 2.54 3.05 2.23

1.33 1.67 2.03 2.14 2.31 1.43 1.13 1.44 1.72 1.90 2.08 1.50 1.44 1.37 1.33 1.31 1.35 1.52 1.50 1.77 1.34

1.17 1.63 1.78

1.25 1.28 1.48 1.57 1.51 2.02 2.19 2.45

1.42 1.42 1.52 1.51 1.58 2.16 2.26 2.67

a The R 26 the anion-hydrogen peak positions are estimated from the anion-O peak P gmax radii and Pauling ionic radii Rion are from Marcus; position by subtracting 0.95 Å; the Goldschmidt ionic radii RionG are from Cotton and Wilkinson;14 the experimental solvation free energies ∆GexptM and ∆GexptFK are from Marcus21 and Friedman and Krishnan,11 respectively; ReffC,P ) (Rgmax + RionP)/2; ReffC,G ) (Rgmax + RionG)/2; ReffE,M and RefE,FK are obtained from eq 1 with  ) 78.5 using ∆GexptM and ∆GexpFK, respectively.

8 represents solvation free energies computed using eq 7 with  ) 78.5 and ReffC,P from Table 3 as a function of q2/2ReffC,P. The experimental ∆GexptM values (open circles) as a function of q2/2ReffC,P lie on or close to the solid line, validating the new expression (eq 7) for the solvation free energy. In contrast, hydration free energies computed using the Born formula (eq 1) with R ) RionP (filled diamonds) and R ) Rgmax (filled circles) deviate from the solid line and experimental values, and the deviations magnify with increasing q2/2ReffC,P values, as found for the model systems (Figure 7). In summary, the results shown in Table 3 and Figure 8 establish the generality of eq 7 in predicting solvation free energies. Estimation of Rgmax and Its Relation to Rion. The new expression for the solvation free energy (eq 7) can be used to

estimate unknown Rgmax of ions whose experimental solvation free energies are available. This is shown in Table 4, which lists monatomic cations and polyatomic anions of varying charge (+4e to -2e), whose Rgmax values have not yet been experimentally determined. Using experimental hydration free energies of these ions in eq 1 with  ) 78.5, ReffE radii were computed. Subsequently, RgmaxC radii were determined using eq 4 with Pauling ionic radii (see Table 4). These predicted radii could be verified using molecular dynamics simulations with ionwater potential parameters that reproduced experimental solvation free energies, as employed in this work (see Methods). Since experimental solvation free energies and Rgmax may be difficult to determine for some ions, while Rion are generally known, it is useful to relate Rgmax to Rion. The results of such

Theory of Ionic Hydration

J. Phys. Chem. B, Vol. 103, No. 37, 1999 7965

Figure 8. Plot of solvation free energies for ions in Table 3 vs q2/ 2ReffC,P. The solid line is computed using eq 7. The open circles are experimental ∆GexptM from Table 3. The filled diamonds and filled circles correspond to free energies computed using eq 1 with R ) RionP and R ) Rgmax, respectively.

TABLE 4: Estimation of Rgmax from Experimental Solvation Free Energiesa ion H+

Cu+ Pb2+ V2+ Cr2+ Pd2+ Ag2+ Yb2+ Eu2+ Sm2+ Co3+ Ga3+ V3+ Ti3+ Au3+ Sc3+ Yb3+ Ho3+ Tb3+ Pm3+ Pu3+ Bi3+ U3+ Hf4+ Zr4+ Ce4+ Pu4+ U4+ OHCNSHBF4MnO4CO32S2SiF62PtCl62-

RionP Å

-∆Gexpt, kcal/mol

ReffE, Å

RgmaxC, Å

0.30 0.96 1.18 0.79 0.82 0.86 0.89 1.05 1.17 1.19 0.61 0.62 0.64 0.67 0.70 0.75 0.87 0.90 0.92 0.97 1.01 1.02 1.04 0.71 0.72 0.80 0.93 0.97 1.33 1.91 2.07 2.32 2.40 1.78 1.84 2.59 3.13

251.0 125.5 340.6 436.2 442.2 456.5 445.8 360.9 331.0 328.6 1074.3 1079.1 1008.6 959.6 1056.4 907.0 853.3 829.4 812.6 776.8 773.2 831.7 766.0 1664.7 1622.9 1462.7 1567.9 1520.1 102.8 70.5 70.5 45.4 56.2 314.3 314.3 222.3 163.7

0.65 1.31 1.92 1.50 1.48 1.44 1.47 1.82 1.98 1.99 1.37 1.37 1.46 1.54 1.40 1.63 1.73 1.78 1.82 1.90 1.91 1.77 1.93 1.58 1.62 1.79 1.67 1.73 1.59 2.32 2.32 3.61 2.92 2.09 2.09 2.95 4.00

1.01 2.61 2.67 2.22 2.15 2.01 2.05 2.58 2.79 2.80 2.14 2.11 2.28 2.40 2.09 2.50 2.59 2.66 2.71 2.83 2.81 2.53 2.81 2.44 2.51 2.79 2.41 2.48 1.86 2.74 2.58 4.90 3.44 2.39 2.33 3.31 4.88

Pauling ionic radii and experimental ∆Gexpt are taken from Marcus.21 ReffE are obtained from eq 1 with  ) 78.5 using ∆GexptM. RgmaxC ) 2ReffE - RionP. aThe

RionP

an attempt are shown in Figure 9a, which shows a plot of q2/ RgmaxC vs q2/RionP. The circles, triangles, and squares in Figure 9a are the data points for some mono-, di-, and tripositive ions in Tables 3 and 4, respectively (Ag+, Be2+, Al3+, Cr3+, Fe3+, and Tl3+ showed large deviations and were eliminated from this plot). The solid line in Figure 9a is the best linear fit of all the data points with a regression coefficient of 0.98. The linear fit

Figure 9. (a) Plot of q2/RgmaxC versus q2/RionP for some monocations (circles), dications (triangles), and tripositive ions (squares) in Tables 3 and 4. The solid line is a linear fit to the data points. (b) Plot of solvation free energies ∆Gcalc computed using eq 3 against experimental ∆GexptM for most of the ions in Table 3. The solid line corresponds to ∆GexptM vs itself. The open and filled circles correspond to ∆Gcalc with ∆ in eq 3 equal to the values obtained by Latimer et al.15 and eq 12c, respectively.

allows Rgmax to be expressed in terms of RionP as

q2

q2 ) C1 + C2 Rgmax Rion

(10)

where C1 ) -0.003 and C2 ) 0.4 for cations. A similar fit for negative ions (considering only halide ions in Table 3 since the other anions are polyatomic ions whose ionic radii are not accurately known) gave C1 ) -0.032 and C2 ) 0.87. Note that these C1 and C2 values have been obtained for ions in water and they will be different for ions in nonaqueous solutions or heavy water. Equation 10 can be rearranged to give Rgmax as

Rgmax )

q2Rion C1Rion + C2q2

(11a)

Since C1Rion is generally much smaller than C2q2 for both cations and anions (see Table 3), Rgmax is roughly proportional to Rion:

Rgmax ≈

Rion C2

(11b)

The coefficients C1 and C2 reflect the properties of the particular solvent. They can be related to the empirical constants ∆ obtained by Latimer et al. (see Introduction) for alkali cations and halide anions using eqs 3 and 7; i.e.,

1 ∆ ) (Rgmax - Rion) 2

(12a)

Substituting eqs 11a and 11b for Rgmax yields eqs 12b and 12c, respectively. The radii of the alkali cations range from 0.74 Å

7966 J. Phys. Chem. B, Vol. 103, No. 37, 1999

∆)

Babu and Lim

q2Rion(1 - C2 - C1Rion/q2)

(12b)

2(C1Rion + C2q2)

∆≈

Rion(1- C2) ) cRion 2C2

(12c)

(Li+) to 1.70 Å (Cs+), yielding an average Rion of 1.27 Å. The latter multiplied by c ) 0.75 for cations gives ∆ ) 0.95 Å, close to the value obtained by Latimer et al. (0.85 Å). Similarly, using the average RionP value of the halide anions (1.83 Å) in eq 12c with c ) 0.075 for anions yields ∆ ) 0.14, which is also close to the value obtained by Latimer et al. (0.1 Å). The validity of eq 12c is shown in Figure 9b, which plots solvation free energies ∆Gcalc computed using eq 3 against experimental ∆GexptM for most of the ions in Table 3. The solid line corresponds to ∆GexptM vs itself. The open and filled circles in Figure 9b correspond to ∆Gcalc with ∆ in eq 3 equal to the values obtained by Latimer et al. and eq 12c, respectively. Both sets of ∆Gcalc values fall close to the experimental line except for PO43-. Solvation Entropies and Enthalpies of Ions. Based on the new formula for the solvation free energy (eq 7), the corrected expressions for the entropy and enthalpy are

(

∆SCorr ) -

)

∂∆GCorr ∂T

[

) ∆GCorr

P

( )

∂Rgmax 1 (Rion + Rgmax) ∂T

( )]

-

P

∂ 1 (13) ( - 1) ∂T P ∆HCorr ) ∆GCorr + T∆SCorr ) ∂Rgmax T ∆GCorr 1 + (Rion + Rgmax) ∂T

[

( )

P

-

( )]

∂ T (14a) ( - 1) ∂T P

The contribution of the second term in eq 14a to the solvation enthalpy is generally nonnegligible, as noted by Roux et al.24 in the context of eq 1. Thus, knowledge of the temperature dependence of Rgmax is needed to compute the solvation entropy and enthalpy from eqs 13 and 14a, respectively. In the absence of such information, the hydration enthalpies for all the ions in Table 3 can be estimated from

[

∆HCorr ≈ ∆GCorr 1 -

( )]

∂ T ( - 1) ∂T P

(14b)

with (T/)(∂/∂T)P ) -1.357, the experimental value for water at 25 °C.11 The ∆GCorr are computed using eq 7 with  ) 78.5 and ReffC,P from Table 3. The relative accuracy of eq 14b in predicting solvation enthalpies is illustrated in Figure 10. The solid line in Figure 10 represents hydration enthalpies computed using eq 14b as a function of q2/2ReffC,P. The experimental values39 (open circles) as a function of q2/2ReffC,P lie on or close to the solid line. In contrast, hydration enthalpies based on the Born model (i.e., eq 14b where ∆GCorr is replaced by ∆GBorn) with R ) RionP (filled diamonds) and R ) Rgmax (filled circles) deviate from the respective experimental values and solid line. The trends in Figure 10 are similar to those in Figure 8 for the free energy. Discussion Limitations in Methodology. The calculations were carried out using a rigid and nonpolarizable TIP3P model of water molecules.30 This is a key drawback of the methodology since water molecules in the vicinity of the ions are subjected to strong

Figure 10. Plot of solvation enthalpies for ions in Table 3 vs q2/2ReffC,P. The solid line is computed using eq 14b. The open circles are experimental ∆Hexpt from Rosseinsky.39 The filled diamonds and filled circles correspond to hydration enthalpies based on the Born model (i.e., eq 14b where ∆GCorr is replaced by ∆GBorn) with R ) RionP and R ) Rgmax, respectively.

polarizing fields from the ions.8 However, such polarization effects are, to some extent, implicitly taken into account in the calculations since the potential energy parameters have been parametrized to yield the experimental solvation free energy of the ion in question.27 In other words, polarization and other effects are implicit in the potential parameters. Another potential source of error stems from the truncation of long-range electrostatic forces. Since the solvation of an ion at infinite dilution, without the interference of the electric field from counterions, is of interest, Ewald summation,43 which requires charge neutrality of the simulation system, is not appropriate for handling such long-range forces. Periodic boundary conditions are also not ideal for simulations of ionic solutions since a large number of water molecules is needed to avoid distortions of the electric fields near the boundaries of the central box, and computed solvation free energies have been shown to be sensitive to the radius used in truncating ion-water forces.27,42,44,45 Our choice of a spherical boundary with a potential cutoff at 20 Å, the radius of the simulation sphere, seems reasonably accurate judging from the long-range behavior of the rdfs and solvation energies of mono- and divalent metal ions (see Figures 2 and 4-6). Limitations in Experiments. Experimental determinations of Rgmax have limitations, as discussed by Marcus.26 The uncertainty in Rgmax determined by X-ray or neutron diffraction is estimated to be 0.01 Å or worse. Although the uncertainty in Rgmax determined by EXAFS is reported to be