J. Phys. Chem. 1995, 99, 5187-5195
5187
Apparent Local Dielectric Response around Ions in Water: A Method for Its Determination and Its Applications Jin-Kee Hyun Department of Physics, Washington State University, Pullman, Washington 99164-2814
C . Satheesan Babu and Toshiko Ichiye* Department of Biochemistry/Biophysics and of Chemistry, Washington State University, Pullman, Washington 99164-4660 Received: September I , 1994; In Final Form: December 20, 1994@ The solvation properties of ions in aqueous solution are determined by the unique structure of water and its dielectric properties. The orientational structure of the water molecules is governed in part by the dipolar nature of water but displays a much richer behavior due to hydrogen bonding. Here, properties of the orientation angle 8 of a water dipole as a function of distance r from the ion, namely, the mean value (cos e), and the distribution &cos 8;r),are calculated from molecular dynamics simulations of C1- in water. The orientational behavior is compared to that of a point charge and a point dipole in a dielectric continuum, which can be solved analytically. More specifically, the microscopic field due to a point charge experienced by the dipole is different than the macroscopic field due to a point charge so that the problem must be viewed as a dipole in a cavity immersed in a dielectric continuum, using arguments similar to those of Kirkwood in relating the dielectric constant E of a fluid to the dipole moment p of molecules in the fluid. By comparing (cos e), calculated from simulation with the analytic expressions, an apparent “local” dielectric response E( r ) dependent on r is calculated; E(r) is termed apparent because it is not a position-dependent dielectric constant to be used in the Poisson or Poisson-Boltzmann equations but rather reflects the cumulative effects of such a positiondependent dielectric constant. The variation of e ( r ) from the bulk dielectric constant allows one to quantitate the deviations from continuum behavior. In the first two shells, E(r) varies significantly with r, and a reduction in e(r) in the first shell is clearly seen. Furthermore, e(r) appears to approach continuum bulk behavior by about 8 8, despite the considerable orientation even at r = 10 8,. In addition, the simple expression for P(cos 8 ; r ) in the continuum approximation using e ( r ) is used to calculate solvation properties of ions. It is also demonstrated that, in the limit where E(r) is assumed constant (Le., that of bulk) and the density is assumed constant (Le., a continuum), the Born solvation energy expression is a limiting case of our e,xpression for the solvation energy.
I. Introduction The solvation structure and energetics of ions in aqueous solution are important to many subjects such as protein stability and electron transfer. However, water structure around ions is complex. For instance, the Born solvation energy equation is accurate only when a special radius, often called the Born radius, is used to describe an ion rather than the crystal radius. While some of the change in radius can be attributed to an actual change in “size” of an ion upon solvation, the standard Born radii are in fact parameters that also correct for the two major assumptions of the Born equation, namely, that the water is a continuum and the dielectric permittivity has the same value even at short distances from the ion.’ Thus, the Born model a priori will not differentiate between cations and anions of the same size and charge magnitude, even though Born radii generally show a larger increase over crystal radii for cations than anions. Rashin and Honig have explained this by noting that the radius of the cavity formed by the ion in water is the relevant radius to the Born equation, not the crystal radii2 Further, t h e ~ r y , ~ and experiment6 all show that the first shell of water around an anion prefers a hydrogenbonding orientation whereas that around a cation has essentially dipolar interactions. Of the two assumptions of the Born equation, water can be easily shown not to be a continuum at separation of a few
* To whom correspondence @
should be addressed. Abstract published in Advance ACS Abstracts, March 15, 1995.
0022-3654/95/2099-5 187$09.00/0
molecular diameters by examining the numerous radial distribution functions of water around ions from t h e ~ r ysimulati0n,42~ ,~ and experiment.6 However, although a reduction in the dielectric permittivity is thought to occur around ions,’ the apparent local dielectric response at distance of a few molecular diameters has only been studied indirectly.’ This reduction is often referred to as dielectric saturation in electrochemical literature, but true dielectric saturation has not been demonstrated since it would imply that the reduction is due to a complete orientation of water molecules near ions by the high electric field of the ions so that the waters can no longer respond to further changes in the field. The concept of a local dielectric response is not new.’ It can be thought of as a combination of several effects including possible dielectric saturation and the specific water structure around an ion of specific size, charge sign, and charge magnitude. It is obviously dependent on the orientation of water around an ion, but although work has been done studying the orientation of water around ions4q5and between pairs of water m o l e ~ u l e s the , ~ ~connection ~ with the local dielectric response has not been made. A simple estimate of the degree of polarization or orientation of water around an ion can be made by considering a point charge and a single point dipole both immersed in a dielectric continuum. Account must also be taken of the fact that the microscopic field at the water dipole due to the ion is not equal to the macroscopic field. By Boltzmann averaging, an equation can be derived that relates the distribution of water molecule
0 1995 American Chemical Society
Hyun et al.
5188 J. Phys. Chem., Vol. 99, No. 14,1995
orientations as a function of distance r from a charged particle to the dielectric permittivity of the continuum. However, the dielectric continuum model will break down at small separations between the charged ion and the water molecule of interest. To study the deviations from continuum behavior, we have evaluated an apparent local dielectric response function ~ ( r by ) solving the analytic expression for the dipole in the continuum using simulation data for the water dipole orientations. We use the term apparent because this ~ ( ris)not the position-dependent dielectric constant needed in the Poisson equation or the Poisson-Boltzmann equation, which is the basis of such programs as DelPhi.lo This is somewhat reminiscent of the work of Dang and Pettit6 in which they fit a potential of mean force between two ions from simulation to obtain an estimate of the bulk dielectric constant, although they did not examine the dielectric permittivity at close separations. This provides a simple means of quantitating the structure-makingor structurebreaking properties of an ion by the deviation of r(r) from the continuum value. Also, this could provide a simple means of replacing solvent via a quasi-continuum approximation, for instance, in calculating the solvation energy due to the polarization of the solvent for charges in aqueous solution. Finally, it can be used to measure the degree of orientation at boundaries in simulations of ions in solution using cutoffs, periodic boxes, etc. This paper is organized as follows. In section 11, we outline our theoretical methodology and the details of the simulations performed. We also demonstrate the connection between our work and the Born solvation energy. In section 111, we present our results for the apparent local dielectric response and the solvation energy. The connection to other continuum theories is also discussed. In the last section, we give our conclusions.
II. Methods A. Theory. Consider a system consisting of a single charged ion with charge 41 immersed in water, in which the water is modeled as discrete particles (see Figure la). This will be referred to as the full system and may be modeled by a simulation of an ion in water with potential parameters such as TIPS," SPC,12 ST2,13 etc. The orientation of any given water dipole will be influenced by the presence of the ion and may be described by 8,defined as the angle between the dipole vector of the water and the internuclear vector between the water and the ion. Here, 8 = 0" corresponds to a bifurcated hydrogen bond and 8 = 180" corresponds to the oxygen pointing to the ion. The probability P(cos 8;r) of a given value of 8 as a function of r, the separation between the ion and the water molecule, may be calculated from a simulation of an ion in solution. Note that this probability is defined so that it is independent of the number of molecules at r and that the complete molecular radial distribution function g(r) may be written as
where 4 is the angle between the plane of the water molecule and the plane formed by the dipole vector and the internuclear vector and P '(#;I-) is the corresponding probability of 4 at r. Both P(cos 8;r) and P'(4;r) are normalized to a value of 1 at a given value of r, i.e., S i l d cos 8 P(cos 8;r) = 1 and I r d 4 P '(4;r)= 1. The orientational effects of the ion on a particular water molecule can be considered by the following model: replace all of the water with a dielectric continuum except for the water molecule whose dipole orientation is being examined and model the ion as a point charge q and the water as a point dipole of moment ,u (see Figure lb). Thus, the energy is independent of
(a>
Figure 1. Models of ionic solvation: (a) full and (b) reduced point dipole.
4 and P'(4;r) = 1/(2n). For simplicity, we will consider a nonpolarizable model of water here. The energy for a point dipole in an electric field E is given by in which cos 8 is the angle between the dipole and the electric field. The field due to an ion in a dielectric continuum with a dielectric constant E is dependent on the distance from the ion and is given by E(r) = &r2
(3)
However, this does not take into account the detailed forces between the solvent molecules. For pure liquids, the Kirkwood formula for the relation between the molecular dipole moment and the bulk dielectric constant takes into account such forces. The original derivation by Kirkwood involves the calculation of the total dipole moment of a spherical specimen of solvent embedded in a dielectric continuum in the presence of an external field.14 However, another interpretation is to consider an Onsager type model15 in which a particular solvent molecule is considered to be a point dipole in a spherical cavity of radius Q, corresponding to the radius of a water molecule. However, the field inside the cavity E differs from the external field EO(such as eq 3) and is given by
E=-----3r 2r
+
lE0
(4)
(considering a nonpolarizable model of water although Onsager and Kirkwood-Frohlich relations are extended to polarizable models and thus these results are easily extended). This relation
Local Dielectric Response around Ions in Water is true when EOis a uniform field (as used in the Onsager and Kirkwood derivations) and also when Eo is due to a point charge outside of the cavity of the dipole, but the latter is true only for the internal field at the center of the cavity, which means the dipole must be located at the center of the cavity. In the discussion so far, we consider the external EO to be due to a point charge without a cavity, since it can be shown that the field outside a point charge located in the center of a cavity (the ion cavity) is simply that of a point charge,16and further it can be shown that, in the limit of infinite separation, the cavitycavity effects are negligible. Of course, when the dipole cavity and the ion cavity are close, cavity-cavity effects become important as17 will be discussed in section IIID. To arrive at the Kirkwood formula, the solvent molecule is considered to have an effective dipole moment p*, where
in which gk = (kf)/Np2, known as the Kirkwood g factor.14 (Note: this is the infinite system gk factor, so the finite system values calculated from simulation18 must be modified appropriately.) Thus, the deviation from point dipole behavior is described by a spherical cavity containing a dipole of moment p*, in which the dipole moment has been modified to take into account the effects of polarization of the environment and the detailed interactions with the neighboring molecules. However, at separations of a few molecular diameters between the ion and the water molecule of interest, one may expect that the continuum dielectric approximation breaks down so that the dielectric response becomes r-dependent and will be denoted as Ewl(r). The subscript I denotes the ion type since the r dependence will depend on the ion, and w refers to water as EwI(r) will approach the bulk dielectric constant for water, E,, as r 00. The utility of the Onsager type interpretation of the Kirkwood formula described in eqs 4 and 5 is that one can then describe the interaction between a solvent molecule and an ion as that between a solvent molecule (point dipole) of moment p* in a spherical cavity and the field E(r) in which
-
J. Phys. Chem., Vol. 99, No. 14, 1995 5189
with B = l / k ~ Tk~ , being the Boltzmann constant and T being the temperature. Furthermore, the average value of cos 8 at a given value of r is then (cos e), = coth@E(r)] - l/wpE(r)l
(10)
The apparent local dielectric response including cavity effects cfI(r) can now be evaluated by calculating (cos e), from a simulation of an ion in water and solving eqs 8 and 10 for cfI(r). Alternatively, P(cos 8;r) may be calculated at each r from simulation and dI(r)obtained by fitting eqs 8 and 9 to the data. In either case, c,I(r) can be obtained from €’ using ](eq I) 7 if gk is known. B. The Solvation Energy and Connection to the Born Solvation Energy. P(cos 8;r)combined with g(r) may be used to calculate other average properties of the system via eq 1, assuming that P ’(q5;r)is a constant. For instance, the average solvation energy U may be calculated by
(U(Rl,R2))= -4n~h:dr
f,d
COS
0 p q COS 8 g(r) x P(cos 0;r) (1 la)
= -4n~p4.h: dr g(r)(cos
e),
(1 1b)
where e is the number density and where the limits of the integration are set at R1 and Rz, both distances from the ion, so that the energy of any shell around the ion can be calculated. In practice, g(r) and P(cos 8;r)may be obtained from simulation for small values of R2, but for large R2 where they are not known may be estimated by
in which R is the distance of closest approach of water to the ion. and (Of course, by the Kirkwood formulation one may invoke an r-dependent gk factor, but it would be very difficult to calculate.) Thus, the system becomes equivalent to an ion and a point dipole of moment p* (not in a cavity) in a dielectric continuum, with an r-dependent dielectric response
+
efI(r)= (2cw1(r) 1)/3 Another interpretation is that it is equivalent to an ion and a point dipole of moment p (not in a cavity) in a dielectric continuum with an r-dependent response
in which e, is the bulk dielectric constant for the model of water being studied. It is interesting to note that the formula for the Born solvation energy can be derived from eq 11. If g(r) is approximated by the continuum value, i.e., eq 12, and c’I(r) is approximated by a constant value, Le., eq 13, then eq l l b becomes
where we used eq 10 for (cos e),. Using a Taylor series expansion of coth(x), this is approximated by
For simplicity, we use the latter interpretation, and we will refer to cfI(r)as the apparent local dielectric response including cavity effects. Thus, the field is Now, using the Kirkwood formula, The energy for this system is then given by eq 2. Finally, by Boltzmann averaging, one can determine the probability of cos 8 at a given value of r, P(cos 8;r).
5190 .I Phys. . Chem., Vol. 99, No. 14, 1995
Hyun et al.
TABLE 1: Cutoffs and Boundary Conditions for Simulations simulation rOno roff rliStLI box lengthb P 39.9 25.0 10.0 10.5 TIPS(1) 8.0 13.0 13.5 39.9 25.0 TIPS(2) 11.0 16.0 17.5 39.9 50.0 TIPS(3) 14.0 16.0 17.5 39.9 50.0 SPC(1) 14.0 r,,, and r,ff are values where nonbonded interactions switch are turned on and off, respectively, and r1istis cutoff for the neighbor lists. Values are given in A. For truncated octahedral box, the distance between parallel square faces in A. Length of simulation (after equilibration) in ps. and eq 13 in eq 15, one obtains
(17) Since the electrostatic potential at the water 4 = (U)/q in this approximation is linear in charge, AG = V2(U),19so that this formula gives the Born solvation energy equation if R is replaced by rion. This may be viewed as part of the origin of the discrepancies between the Born and crystal radii, since R is determined by both the size of the ion and the approach of water. Alternatively, one may interpret as meaning that, in the Bom approximation, the water has zero radius since it has no structure and thus R rion. Thus, the Born equation may be derived from eq 11 in the limit that the solvent is a continuum (eq 12), E(r) is a constant (eq 13). and pqp/e‘ 0. The connection of our work with the Born model shows that it has the proper limiting values. However, we stress that in the full form, i e . , using eq 1 for g(r). eqs 8 and 9 for P(cos 8;r), and an r-dependent €‘, as( defined I),in eq 7, we describe the solvation of an ion without treating the ionic radius as an adjustable parameter as such as is done to obtain Born radii. Of course, both g(r) and e(r) must be calculated, but this may be obtained from a simulation where the ionic radius is well-defined. C. Simulation. Simulations of chloride ion in water were performed using the macromolecular mechanics/dynamics program, CHARMM (version 22g2).20 All simulations were in the microcanonical ensemble at a temperature of approximately 298 K. The Verlet algorithm was used to propagate the dynamics. The TIP3P water1’ as implemented in the CHARMM parameters (version 19) and SPC12water models were used and the interaction parameters for the chloride ion with water were taken from Brooks.21 In each simulation, a single chloride ion was simulated with 1062 water molecules in a truncated octahedral box of length 39.9054 8, using periodic boundary conditions.22 All van der Waals and electrostatic interactions were switched smoothly to zero, and neighbor lists were used. Although several works23have shown cutoffs to be less accurate than Ewald sums for charged systems in aqueous solution, for the work presented here it is necessary to simulate a system consisting of a chloride ion at very low dilution, and so our simulation box is net charged and Ewald sums cannot be used (unless a jellium model is used to neutralize the charge, which would change the physics of the system). To examine the effects of cutoffs on our results, several different cutoffs were employed (Table 1). These points are addressed further in the Results and Discussions section. For the TIPS(1) and TIPS(2) simulations (see Table 1 for definitions), initial coordinates for the water were taken from a previously equilibrated box of pure water at 1 gm/cm3. Equilibration consisted of reassigning velocities from a MaxwellBoltzmann distribution at 298 K every 0.2 ps for 5 ps, rescaling velocities to 298 K every 0.2 ps for 15 ps. and then allowing
-
-
the system to evolve unperturbed for another 5 ps. For the TIPS(3) simulation, a simulation was started from the equilibrated configuration of water and ion obtained from the TIPS(2) simulation with different cutoff distances so that the system came to equilibration immediately. Finally, for the SPC( 1) simulation, the equilibrated coordinates of water and ion of the TIPS(2) simulation were taken as a initial coordinates with SPC potential parameters, and the same cutoff as that of TIPS(2) was used. After minimization for 0.5 ps, equilibration was obtained by reassigning velocities from a Maxwell-Boltzmann distribution at 298 K every 1 ps for 50 ps and then rescaling velocities to 298 K every 1 ps for 45 ps. After equilibration with the cutoff of TIPS(2), the cutoff was extended to that of TIPS(3), and equilibrium was obtained by rescaling velocity to 298 K every 0.2 ps for 25 ps and then allowing the system to evolve unperturbed for 10 ps. In all cases, a time step of 0.002 ps was used during equilibration but was decreased to 0.001 ps during the analysis portion of the simulation. The lengths of each simulation after equilibration are given in Table 1.
III. Results and Discussions Simulation results for two models of water are shown here: the TIP3P model” and the SPC model.12 Both are used in commercial software packages, and the latter has many studies of its dielectric properties. Results for the TIPS(3) and the SPC( 1) simulation are shown in Figures 2 and 3, respectively. The radial distribution functions, g(r), for the SPC(1) simulation matched other simulations of C1- ion in SPC water21 and the Cl--oxygen g(r) is shown for reference in Figure 3a. The g ( r ) of the TIPS(3) simulation is similar, but shows slightly more pronounced peaks, and is shown in Figure 2a. In what follows, the distance of a water molecule from the ion is defined as the distance of the oxygen from the ion. Although other choices are possible, this choice gives better agreement with the energetics of the full charge model than using, for instance, the center of charge. A. Orientation of Water and the Apparent Local Dielectric Response. The average orientations (cos e), for the TIPS(3) and SPC(1) simulations are shown in Figures 2b and 3b, respectively. The value of cos 8 = -1 corresponds to the oxygen pointing to the central ion, cos 8 = 1 corresponds to a bifurcated hydrogen bond orientation, and cos 8 x 0.65 corresponds to a linear hydrogen bond. Also shown are (cos e), calculated from eqs 8 and 10 evaluated assuming two continuum values for E’I(r); namely, e,.,. the bulk dielectric constant, and e’,,, = ( 2 ~1)/3gk ~ (see eq 7), the value including microscopic forces, Le., the cavity effects. For SPC water, e’,,, x 17-18 (from ref 9, E,,, x 68 and gk x 2.6 and from ref 23, e,,, = 57.3 and gk = 2.27). The value of E for TIP3P water is assumed to be close to that of SPC because of the similarity of the models and because a fit to the long tail region of the potential of mean force between C1- and C1- from a simulation in TIPS water gives a value of 72.4 From Figures 2b and 3b, it is obviou! that there is orientational dependence even at large r (near 10 A) in the simulations, and the orientation at large r is comparable to that calculated from eqs 8 and 10 with dI(r) = e’,,,. This is important when determining values for cutoffs and for periodic box sizes, regardless of whether cutoffs are used. Moreover, the degree of orientation is much higher than when calculated from eqs 8 and 10 using e’l(r) = e,.,, indicating that the (dipole) cavity effects, including Rk, are important. An apparent local dielectric response including cavity effects may also be calculated from simulation data by solving for the value of dl(r) that gives the observed value of (cos e),.in the simulation at each value of r using eq 10. This cavity apparent local dielectric response is plotted as a function of r for the
+
Local Dielectric Response around Ions in Water
J . Phys. Ckem., Vol. 99, No. 14, 1995 5191
(a)
4,
4
3
n L
Y
2
e l 1
.
0
5
10
15
r ( A )
0.2
0.0
0
5
10
15
15
100 80 80
60
.
-
- 260 h
-
1.n 40
-
*
Y
m
h
v L
w
60
260
W
40
- 130
U 1
’
m
I
h U 1
130
20
20
’
0 I
Figure 2. TIPS(3) simulation. (a) CI--oxygen radial distribution function. (b) (cos e), as a function of r: the solid line is from
simulations, and dotted line and dash-dotted line are from eq 10 evaluated using e’ = 18 and cw = 68, respectively. (c) Apparent local dielectric response as a function of r: the solid line is calculated from (cos e), by fitting to eq 10, and the dotted line with solid triangles is calculated from P(cos 0 ; r )by fitting to eq 9. TIPS(3) and SPC(1) simulations in Figures 2c and 3c, respectively, to quantitate the deviation from continuum behavior. In Figures 2c and 3c, the left-hand axis corresponds to c’I(r), the cavity apparent local dielectric response for a point dipole of moment p not in a cavity, and the right-hand axis corresponds to ewI(r),the apparent local dielectric response for a point dipole of moment p in a cavity, with the two related by eq 7. Although E , I ( ~ ) is more physically relevant than E’I(r) since it approaches the bulk dielectric constant cw at larger r, both are reported since
0
5
10
0 15
r ( A ) Figure 3. SPC(1) simulations. (a) Cl--oxygen
radial distribution function. (b) (cos e), as a function of r: the solid line is from simulations, and dotted line and dash-dotted line are from eq 10 evaluated using e’ = 18 and ew = 68, respectively. (c) Apparent local dielectric response as a function of r: the solid line is calculated from (cos e), by fitting to eq 10, and the dotted line with solid triangles is calculated from P(cos 8;r) by fitting to eq 9. c’I(r) is obtained directly from eq 10 whereas cwI(r)is dependent on the value of R k used. Figures 2c and 3c show that cwI(r) approaches 60-70, the value expected for SPC and TIP3P water (and that c’I(r) approaches 15-20, the value expected for the microscopic continuum value via eq 7), after the second shell. Since the details of the distribution are lost in looking at (cos e),, the orientational probabilities P(cos e;r) were also calculated for shells of water around the chloride ion, which were defined approximately by the minima in the g(r). To obtain sufficient statistics to calculate P(cos 0;r)for a range of
5192 J. Phys. Chem., Vol. 99, No. 14, I995
Hyun et al.
’ i Oms
I
0.7
-
fi
a
a
8 0.5
u)
8
u
W
W
n
n 0.3
t
0.1-1 .o
-0.5
0.0
0.5
cos e
cos e
0.7 h
a
8 0.5 U
0.3
- 1 .o
0.0 0.5 1 .o cos e Figure 4. P(cos 8;r)as a function of cos 8 for the first shell (‘4 A) around the chloride ion (a) from the TIPS(3) simulations and (b) from the SPC(1) simulation. -0.5
8, it was necessary to use fairly wide shells. These results for the first two shells are shown in Figures 4 and 5. For the first shell between 0 and 4 8, (Figure 4), P(cos 0;r) from the simulation deviates completely from the behavior of eqs 8 and 9. This is because of the presence of the linear hydrogen bond orientation, which cannot be mimicked by a point dipole. In the second shell between 4 and 6 8, (Figure 5 ) , P(cos 0;r)from the simulation is fit well by a least-squares fit to eq 9, where the value of e ’ ~ ( r is ) varied, as do the results for subsequent shells, which are not shown here. The apparent local dielectric response calculated from P(cos 0;r) by fitting to eqs 8 and 9 are also shown in Figures 2c and 3c. In what follows, we use e(r) from (cos e), in preference to that from fitting P(cos 0;r) because it involves simply solving an equation rather than a least-squares fit and because it is easier to collect accurate data for a single value of (cos e), at every value of r rather than a histogram of cos 6 at each r. The first two shells show large deviations of cwl(r)from the bulk cw. While the first peak in g ( r ) corresponds to a maximum in (cos e), and thus a minimum in cwl(r),the second peak in g(r) corresponds to a minimum in (cos e), (see Figures 2 and 3) and thus a maximum in cwl(r). This behavior will be discussed further in section 1II.D. B. Effects of Simulation Details on E @ ) . The determination of cW1(r) by comparison with simulation requires accurate treatment of the long-range electrostatic forces in the simulation, which are most commonly treated via cutoffs, but also by Ewald sums or reaction fieldsz2 For pure water systems, the orientational ~ t r u c t u r eand ~ ~computed ~ ~ ~ , ~ dielectric ~ constant25are greatly affected when Coulombic interactions were cut off at 8
.o
0.0 0.5 1 .o cos e Figure 5. P(cos 8;r) as a function of cos 0 for the second shell (4-6 A) around the chloride ion (a) from the TIPS(3) simulations and (b) from the SPC(1) simulations. The solid line is from simulations and the dashed line is the best fit to eq 9. -1
-0.5
8, although the energetics and radial structure are not significantly affected.26 On the other hand, simulations that used either Ewald sums8 or reaction fields25,27,28 to effectively include the neglected interactions significantly improved the orientational structure. Similarly, the computed dielectric constant has been found to be substantially smaller when cutoffs are used than when a reaction field9 or an Ewald sum* is used. In simulations of peptides in water, several workers have shown that large cutoffs of at least 10-15 8, are needed to adequately simulate these system^,^^^^^ although other simulations of a-helical peptides in water have shown that even with cutoffs of 14 8, there can be large differences from when an Ewald sum was used23 and a simulation of a phospholipid monolayer-water system showed that cutoffs of 30 8, were necessary.31 For studies of ions in water, the problem is more complicated if the system is net charged as Ewald sums cannot be used (a jellium model may be used to neutralize the net charge, but the physics of such a system is different), although a reaction field method for such systems has just been reported.32 Neglect of long-range interactions produced a 15% difference in the solvation free energy when going from 7.5 to 8.5 8, cutoff was used.33 Thus, at the present time, cutoffs are the only practical means of calculating simulations at very low or infinite dilution. The average orientations (cos e), were calculated from TIPS simulations with various cutoffs (Figure 6) to investigate the effects of cutoffs on our results. The behaviors of all three simulations with different cutoffs for the first and the second shell are similar (for the first shell, they are indistinguishable) with maxima near 3 and 6 8, and minima near 5 and 7 8,.
J. Phys. Chem., Vol. 99, No. 14, 1995 5193
Local Dielectric Response around Ions in Water r n
I
I
I
0.8
0.2
0.0
0
10
15
C. Solvation Energies. The average solvation energies of different shells for the TIPS(3) simulation calculated using eq 11 are shown in Table 2. For reference, “exact” values were calculated from the simulation by averaging the solvation energy over the entire trajectory
in which r, denotes the distance of atom a from the C1-, and the sum extends over all water molecules with the oxygen distance between R1 and R2. In addition, the dipole approximation is tested by calculating the average energy using the simulation data when all the waters are replaced by point dipoles located at the oxygen of each water, i.e.,
e), as a function of r in TIPS simulations (see Table 1): TIPS(1) (dotted line), TIPS(2) (dash-dotted line), and TIPS(3) (solid line).
Figure 6. Cutoff dependence of (cos
implying that cutoffs do not affect the orientation at such close in which r, now denotes the distance of the oxygen of water ion-water distances. The most accurate simulation, TIPS(3), molecule “a” from the C1-. As noted in the beginning of section shows relatively good agreement with that calculated from eqs 111, the energetics are best when the point dipole was placed at 8 and 10 with c’I(r) = clWbetween 8 and 12 8, (see Figure 2b), the center of the oxygen rather than at the points such as the center of charge. Table 2 shows that the dipole approximation indicating continuum behavior is reached after the second shell. However, there is a minima at 13 8, and a subsequent maxima breaks down at small distances as expected, decreasing the magnitude by about 20 kcaymol. However, beyond the first at 14 A, the latter distance corresponding to the beginning of the potential energy switching region. The SPC( 1) results are shell, the dipole approximation is in very good agreement with the “exact” charge distribution. Thus, beyond about 4 A, water similar. Moreover, both the TIPS( 1) and TIPS(2) simulations also show oscillations before the beginning of their respective behaves essentially as point dipoles in terms of the solvation switching regions, with the TIPS(2) showing oscillations at 8 energy. 8, even though the switching region begins at 11 8,. In sum, The solvation energies using eq 11 are calculated with g(r) the results for (cos e), up to 7 8, appear consistent between the from simulation, denoted [g(r)]-exact~~, or using the continuum three TIPS simulations, and the behavior after 8 8, appears to estimate in g(r) 1 r R = 2.9 A, denoted [g(r)lCont, and with be like a continuum based on the TIPS(3) simulation, although (cos e), from simulation, denoted (cos f3)r,-exact~~, or using the this is not definitive. continuum estimate (eq 13) with cw= 68 and cfw= 18, denoted The oscillations at the beginning of the switching function The (U(R1,RZ)) calculated using [g(r)]“exact*. and (COSe),,,,,t. region are due to waters that experience less than the full (cos f3),,-exact~~, referred to as I, may be viewed as a test of the Coulomb interaction with the ion due to the switching function integration, since they are directly from the simulation and thus but also interact with water molecules that are fully oriented should equal to that from eq 19 (dipole). The values are in by the ion. A similar oscillation in g(r) in the cutoff region good agreement with those of dipole approximation. has been analyzed34 and can be seen in the g(r) shown here The (U(R1,Rz))calculated using [g(r)]~sexact., and (COSe)r,cont, (Figures 2a and 3a). The peak in (cos e), of TIPS(2) at 8 8, is referred to as 11, shows the contribution of the deviations from most likely due to the propagation of the anomalous orientation continuum behavior of the orientation of the dipoles. For the at the start of the switching region, and both TIPS(1) and first shell, the average solvation energy is 20 kcdmol less in TIPS(2) show slightly greater peaks at 7 8,. Since this magnitude than that of I, and for the second shell, it is 7 kcal/ anomalous orientation is less in the TIPS(3) simulation, due to mol greater in magnitude than I, which indicates that there are the smaller Coulomb force at greater distance, there is apparently large deviations from continuum behavior in the orientation of less propagation of this effect to closer distances. However, it the dipoles in this region (Figure 2a). However, for the outer is also possible that the propagated effect cancels a real shells, the results of I1 are within about 0.8 kcdmol of those of I, indicating that the continuum estimate of c’I(r) = dw is oscillation in orientation. reasonable for outer shells. Another practical difficulty which affects the calculation of dielectric properties from simulation is the length of the The (U(R1,Rz))calculated using [g(r)lCont and (cos f3),;LeXaCt*, simulation. In recent years, many evaluations of the bulk referred to as III, may be compared with I to see the contribution dielectric constant of various models of water s i m u l a t i ~ n ~ ~ *of~ variations ~ with g(r) while [P(cos e;r)],, is used in both cases. have noted that simulation times of nanoseconds are necessary For the first shell, the energy is about 25 kcaumol smaller in because one needs to calculate the fluctuations in the dipole magnitude since the real g(r) has a large peak that affects the moment of the simulation box, (Ai@), which is a very slowly calculated (U(R1,Rz)) considerably in this region. However, converging property. The Debye relaxation time, t ~ .is values become much closer as shell becomes large, and then measured experimentally to be 9.3 ps at 298 K.35Although eventually, the values are in ood agreement within 0.05 k c d we are not attempting to quantitate the bulk dielectric constant mol for the outer shells (’6 ) since g(r) becomes constant at (and, in fact, we use literature values for gk), the dielectric large distances. This may be considered as the effect of the properties measured here are expected to converge faster as they number density of water molecules around the ion on the are dependent on the relaxation time of the individual water average solvation energy. molecules, which has been experimentally measured as 3.6 ps A completely continuum estimate, referred to as IV, in which at 293 K.36 However, this is potentially an easier means of the (U(R1,Rz)) is calculated using [g(r)lContand (cos 6’)r,cont, calculating the bulk dielectric constant. shows very large deviations from I in the first shell and large
-
K
5194 J. Phys. Chem., Vol. 99, No. 14, I995
Hyun et al.
TABLE 2: Average Solvation Energy (kcaymol) sheli (A) exactO dipoleb I’ IId 111‘ IVf Born -67.26 f 2.39 -67.60 0-4 -86.72 f 3.25 -46.74 -41.25 -29.42 -31.04 4-6 -19.98 f 1.53 -19.58 f 1.46 -19.56 -26.93 -20.05 -28.06 -27.28 -13.62 f 1.68 -13.42 f 1.70 -13.37 6-8 -14.00 -13.92 -14.36 -13.64 8-9 -5.58 f 1.32 -5.57 f 1.38 -5.57 -4.78 -5.59 -4.81 -4.54 9-10 -3.57 f 1.00 -3.56 f 1.05 -3.62 -3.89 -3.58 -3.85 -3.64 -3.44 f 0.83 10-11 -3.37 f 0.88 -3.32 -3.17 -3.32 -3.15 -2.97 11-12 -3.19 f 0.68 -3.21 f 0.67 -3.26 -2.61 -3.27 -2.63 -2.48 ‘Calculated from simulation using eq 18. Calculated from simulation using eq 19. Calculated using e 11 with [g(r)]-eracr~~ and (cos 0 Calculated using eq 11 with [ (r)]-exact~. and (cos 0 ),.,. e Calculated using eq 11 with [g(r)],,,,, R = 2.9 , and (cos O)r,-exz,-, f Calculated using eq 11 with [gfr)],,,,, R = 2.9 and (cos 0),,,,,,.
1,
deviations from I as well in the second shell. Overall, from Table 2, the dipole approximation, the oscillations in the real g ( r ) and the oscillations in ~ ( r individually ) decrease the magnitude of the calculated (U(Rl,R*)) values by 20-25 kcaY mol in the first shell, but only the oscillations in € ( I ) significantly affect the second shell. This can be interpreted as that the noncontinuum behavior of the structure of the solvent, namely, electrostriction, is important only in the first shell whereas the variations of the dielectric response is important in the first and second shell. Furthermore, the comparison of IV with I and the “exact” values indicates that overall IV is in reasonable agreement for the outer shells, and so it may be used to estimate the solvation energy beyond the cutoff region. Finally, we compare IV with the Bom solvation energy, the difference being that the limit Pqpk‘ .- 0 has been taken in the former (eq 15). The results are similar and for large r both approach exact values, indicating either may be used to approximate solvation energies beyond cutoffs in simulations. However, although the Bom equation is simpler for a simple ion, the method demonstrated by IV is readily extendible to arbitrary charge distributions of the solute(s). The radii needed to get the correct total solvation energies of the simulation ( i e . , Bom-type radii) are 2.0 8, to match the “exact” results and 2.3 8, to match the “dipole” results. We stress that these radii are different from cavity radii introduced by Rashin and Honig2 since we have not considered water as a dielectric continuum. D. Connections to Continuum Theory. The work described here makes a connection between microscopic and continuum pictures of aqueous solvation of ions, much in the spirit of works of Debye, Onsager, and Kirkwood, in relating in the bulk dielectric constant E of a fluid to the dipole moment p of molecules of that fluid. However, some care must be taken in the interpretation of the apparent local dielectric response given here. One major point is that this ~ ( ris) not a positiondependent value for dielectric screening to be substituted into Poisson’s equation or a Poisson-Boltzmann description as found in programs such as DelPhi,lo but rather the cumulative screening specifically between an ion and a water molecule. We also note that, by assuming a point dipole, some of the variation in ~ ( ris) not due purely to electrostatics but rather to the fact that there is some structure in the charge distribution (Le., three-point charges in SPC or TIPS). Although it is also possible to solve for P(cos 8;r) and (cos e), for a three-point (or more) charge model of the water, it must be done numerically and thus becomes less useful for more complicated solutes. In addition, as mentioned earlier, the ion cavity has been ignored, which is responsible, at least in part, for the reduced value of c’(r) near the ion. The field E the dipole experiences is affected somewhat at close separations, and although one can derive analytical expressions for this field, the relationship between E and the ion field &(r) (eq 3) is much more complicated than the simple expression given by eq 4 and is dependent on 8. This again makes it less useful for more complicated solutes.
x
)r:exaccw.
A simple explanation for the behavior of the first shell is that, at very close separations,one may expect the local dielectric permittivity to be very small, because there is less dielectric material between the particles, Le., an ion cavity-dipole cavity effect. This reduction in the apparent local dielectric response is often termed dielectric saturation in the electrochemical literature; however, it is not truly dielectric saturation unless it can be shown that the water molecules have become so strongly oriented by the field of the ion (perhaps in part due to the lower apparent dielectric response due to the cavity-cavity effects) that they cannot respond. Thus, ~ ’ ( rmust ) be reduced not only because there is less dielectric material but also because the position-dependent dielectric constant must be less near the ion. Comparisons to the expected values for (cos e), when cavitycavity effects are included, with the dielectric media having a single value, are necessary to prove or disprove actual dielectric ~aturati0n.l~ Moreover, for water, the situation is more complicated since for the preferred orientation of the first shell around C1- is the linear H bond (Figure 4) whereas it would be cos 8 = 1 for a dipole. In fact, Figures 2b and 3b show that the first shell waters have an average (cos e), that is higher than that expected for a point dipole (in a cavity) in a continuum with dielectric. This means that while the most probable value of cos 8 for water is less than that for a dipole in a dielectric continuum as seen in Figure 4, the average value is higher because it is more ordered and thus has a smaller range of values for cos 8. In other words, for the first shell waters around negative ions, it is a competition between the lower most probable value of cos 8 due to hydrogen bonding favoring a higher apparent local dielectric response and the greater ordering due to cavity-cavity effects favoring a lower apparent dielectric response, with the latter making a larger contribution. The behavior of the second shell also shows a large deviation from continuum behavior in that cW1(r)is much greater than cw in Figures 2c and 3c. This is a demonstration of the structurebreaking properties of a negative ion, in that the apparent local dielectric response is higher than the bulk dielectric constant near the ion. This may be due in part to the hydrogen-bonding orientation of the first shell waters, which means that the nonhydrogen bonding hydrogen points toward the second shell, thus reducing the orienting effects of the C1-. Since the water is less ordered, (cos e), is reduced so that ewl(r)is increased. Another very interesting consideration arises from the works of Rashin and H o n i g , * ~Pratt ~ ~ rather than the free energy AA, the calculations are much less computationally intensive, allowing for a much extensive testing of the effects of different charges and ionic radii. More importantly, by examining the contribution in shells around the ion, we can both show that the large r behavior approaches the correct continuum and examine where the deviations from continuum occur. Furthermore, the separation of the full orientation dependent g(r) into g(r) and P(cos 0 ; r )(eq 1) is also a separation into a density term and dielectric term.
-
-
A2338
IV. Conclusions This paper presents a way of obtaining a distance-dependent apparent local dielectric response from simulation. By relating the problem to that of a single ion and a single point dipole in a cavity immersed in a dielectric continuum, a distancedependent apparent local dielectric response d ( r ) may be calculated by comparing the analytic expressions for the orientation of water dipoles around an ion to simulation data. In particular, we solve the analytic expression for (cos e), using the simulation values for (cos e), to obtain €’(I). The results shown here are for C1- in both TIP3P and SPC water, which show dielectric saturation (in terms of a reduction of an apparent local dielectric response) and hydrogen bonding play important roles in determining the r-dependent apparent dielectric response, and which imply that C1- is a structure-breaking ion. Also, while there is considerable orientational order at large r, the orientational dependence appears to have continuum behavior beyond 8 A. Finally, expressions for the solvation energy are given, which are demonstrated to lead to the Born solvation energy equation in the limit where the solvent is a continuum, the dielectric response is a constant, and Pqpld 0, and are shown to be good estimates of solvation energies.
-
Acknowledgment. T.I. acknowledges Mickey Schurr, B. Montgomery Pettitt, Lawrence R. Pratt, Bany Honig, and Peter J. Rossky for helpful discussions. This research was supported in part by grants from the (U.S.) National Science Foundation
(MCB-9118085) and the American Chemical Society-Petroleum Research Foundation (23983-G6) and in part by funds provided by Washington State University. We also thank the VADMS Laboratory at WSU for computational resources.
References and Notes (1) Conway, B. E. Ionic Hydration in Chemistry and Biophysics; Elsevier Scientific Publishing: New York, 1981. (2) Rashin, A. A.; Honig, B. J. Phys. Chem. 1985, 89, 5588. (3) Pettitt, B. M.; Rossky, P. J. J. Chem. Phys. 1986, 84, 5836. (4) Dang, L. X.; Pettitt, B. M. J. Phys. Chem. 1990, 94, 4303. (5) Impey, R. W.; Madden, P. A.; McDonald, I. R. J. Phys. Chem. 1983, 87, 5071. Mezei, M.; Beveridge, D. L. J. Chem. Phys. 1981, 74, 6902. Bounds, D. G. Mol. Phys. 1985.54, 1355. Heinzinger, K.; Vogel, P. C. Z. Nahqforsch. 1974,29A, 1164. Vogel, P. C.; Heinzinger, K. Ibid. 1975, 30A, 789. Heinzinger, K.; Vogel, P. C. Ibid. 1975, 31A, 463. Vogel, P. C.; Heinzinger, K. Ibid. 1976, 3IA, 476. Bopp, P.; Heinzinger, K.; Jancso, G. Ibid. 1977,32A, 620. Palinkas, G.; Reide, W. 0.;Heinzinger, K. Ibid. 1977, 32A, 1137. (6) Narten, A. H; Vaslow, F.; Levy, H. A. J. Chem. Phys. 1973, 58, 5018. Hertz, H. G.; Radle, C. Ber. Bunsen-Ges Phys. Chem. 1973, 77,521. (7) Booth, F. J. Chem. Phys. 1951, 19, 391, 1327, 1451; 1955, 23, 453. See also ref 1. (8) Belhadj, M.; Alper, H. E.; Levy, R. M. Chem. Phys. Lett. 1991, 179, 13. (9) Alper, H. E.; Levy, R. M. J. Chem. Phys. 1989, 91, 1242. (10) Gilson, M. K.; Sharp, K. A.; Honig, B. J. Comput. Chem. 1992, 97, 7716. Gilson, M. K.; Honig, B. Proteins 1988, 4, 7. (11) Jorgensen, W. L. J. Am. Chem. SOC. 1981, 103, 335, 341, 345. (12) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; Hermans, J. In Intermolecular Forces; Pullman, B., Ed.; D. Reidel: Dordrecht, 1981. (13) Stillinger, F. H.; Rahman, A. J. Chem. Phys. 1974, 60, 1545. (14) Kirkwood, J. G. J. Chem. Phys. 1936, 4, 592. (15) Onsager, L. J. Am. Chem. SOC.1936, 58, 1486. (16) Bottcher, C. J. F. Theory of Electric Polarization; Elsevier Scientific: New York, 1973; Vol. I. Stratton, J. A. Electromagnetic Theory, 1st ed.; McGraw-Hill: New York, 1941. (17) Ichiye, T.; Hyun, J. K. Unpublished results. (18) Neumann, M. J. Chem. Phys. 1986, 85, 1567; see also refs 8 and 23 and references therein. (19) Rick, S . W.; Beme, B. J. J. Am. Chem. SOC. 1994, 116, 3949. (20) Brooks, B. R.; Bruccoleri, R. E.; Olafson, B. D.; States, D. J.; Swaminathan, S . ; Karplus, M. J. Comput. Chem. 1983, 4 , 187. (21) Brooks 111, C. L. J. Phys. Chem. 1986, 90, 6680. (22) Allen, M. P.; Tildesiey, D. J. Computer Simulation of Liquids; Clarendon Press: Oxford, 1987. (23) Schreiber, H.; Steinhauser, 0.Biochemistry 1992,31,5856. Smith, P. E.; Pettitt, B. M. J. Chem. Phys. 1991, 95, 8430. (24) Smith, P. E.; van Gunsteren, W. F. J. Chem. Phys. 1994, 100,3169. (25) Steinhauser, 0. Mol. Phys. 1982, 45, 335. (26) Pangali, C. R.; Rao, M.; Beme, B. J. Mol. Phys. 1980, 40, 661. (27) Barker, J. A.; Watts, R. 0. Mol. Phys. 1974, 26, 789. (28) Watts, R. 0. Mol. Phys. 1974, 28, 1069. (29) Loncharich , R. J.; Brooks, B. R. Proteins: Srruct., Funct., Genet. 1989, 6, 32. (30) Kitson, D. H.; Hagler, A. T. Biochemistry 1988, 27, 5246. (31) Alper, H. E.; Bassolino, D.; Stouch, T. R. J. Chem. Phys. 1993, 98, 9798. (32) Alper, H.; Levy, R. M. J. Chem. Phys. 1993, 99, 9847. (33) Chandrasekhar, J.; Spellmeyer, D. C.; Jorgensen, W. L. J. Am. Chem. SOC.1984, 106, 903 (34) Brooks 111, C. L.; Pettitt, B. M.; Karplus, M. J. Chem. Phys. 1985, 83, 5897. (35) Hasted, J. B. In Water: A Comprehensive Treatise; Franks, F., Ed.; Plenum Press: New York, 1973; Vol. 11. (36) Jonas, J.; DeFries, T.; Wilbur, D. J. J. Chem. Phys. 1976,65,582. (37) Jean-Charles, A.; Nicholls, A.; Sharp, K.; Honig, B.; Tempczyk, A.; Hendrickson, T. F.; Still, W. C. J. Am. Chem. SOC.1991, 113, 1454. (38) Jayaram, B.; Fine, R.; Sharp, K.; Honig, B. J. Phys. Chem. 1989, 93, 4320. (39) Pratt, L. R.; Hummer, G.; Garcia, A. E. Biophys. Chem. 1994, 51. 47. JP942350B
