J. Phys. Chem. B 1998, 102, 5673-5682
5673
Continuum Corrections to the Polarization and Thermodynamic Properties of Ewald Sum Simulations for Ions and Ion Pairs at Infinite Dilution Shinichi Sakane,† Henry S. Ashbaugh,‡,§ and Robert H. Wood*,† Department of Chemistry and Biochemistry, Department of Chemical Engineering, and Center for Molecular and Engineering Thermodynamics, UniVersity of Delaware, Newark, Delaware 19716 ReceiVed: January 20, 1998
Ewald sum simulations of the solvation thermodynamics of charged solutes employ an unphysical uniform neutralizing background charge to ensure convergence of the periodic lattice sum. Continuum corrections to the solvation free energy, entropy, and volume of single ions and ion pairs for the effect of the neutralizing background and periodic boundary conditions at finite cell size are derived in order to allow efficient calculations of the ionic properties at infinite dilution. The derivation presented in this paper shows the physical origin of the effects and can be easily extended to multiple charge sites. Corrections are small for high dielectric constant solvents but become increasingly important as the ion size is increased, the dielectric constant is decreased, or the unit cell size is decreased. An alternative way of calculating the thermodynamic properties from Ewald sum simulations is proposed for which the corrections are small for low dielectric constant solvents. Tests for low and high dielectric constant water show that, after the appropriate continuum corrections are applied, the free energy of charging an ion using Ewald sum simulations agrees with the results for potential truncation simulations and spherical boundary simulations (when corrected for truncation effects). The corrected free energy of hydration is not sensitive to the system size even for low dielectric constants. The continuum model correctly predicts that Ewald sum simulations yield a solvent polarization at large distance from an ion that is smaller than the polarization of a truly infinite system.
I. Introduction One of the long-term goals of our laboratory is to develop methods for the prediction of the thermodynamic properties of dilute aqueous solutions under extreme conditions where experimental measurements are either difficult or impossible. Presently, we are exploring the use of computer simulations of classical Hamiltonian models to make these predictions. This method requires (1) accurate calculation of the properties of the classical Hamiltonian model, (2) refinement of the Hamiltonian model to fit both experimental data and quantum mechanical calculations, and (3) estimation of the uncertainties of the calculation due to model errors. For nonpolar solutes it has been shown that the proposed method works very well.1 However for ions at infinite dilution, there are questions about the accuracy of the calculated properties from a classical Hamiltonian model due to the long-range nature of the electrostatic forces and the approximations to these forces that are made. Commonly used methods for calculation of electrostatic interactions for the simulation of ions at infinite dilution are Ewald summation (ES),2-4 reaction field (RF),4 and potential truncation (PT)4,5 with periodic boundary conditions (PBC) and spherical boundary conditions, e.g., the surface-constrained allatom solvent model (SCAAS).6,7 In principle all the methods are exact in the limit of infinite system size with the appropriate * Corresponding author. E-mail:
[email protected]. Fax: 302-831-6335. † Department of Chemistry and Biochemistry, and Center for Molecular and Engineering Thermodynamics. ‡Department of Chemical Engineering, and Center for Molecular and Engineering Thermodynamics. § Current Address: Department of Physical Chemistry 1, Chemical Center, S-221 00 Lund, Sweden.
treatment of the surface effects.3,4,8 Because of the long range of the forces involved, however, convergence to the infinite dilution limit can be very slow unless corrections for finite box size are made. In a previous paper,9 Wood derived continuum corrections for free energies of ionic hydration to PT simulations. These corrections are for the effects of truncated electrostatic forces with ion-solvent interactions cutoff at Ris, solventsolvent interactions cutoff at Rss, and PBC. The simulations of aqueous single ions of Straatsma and Berendsen5 with varying Ris and Rss were corrected and give the same free energies for different cutoff radii. Hummer et al. used the ES method with the Hamiltonian including the “self-interactions” (i.e., the interactions of the ion with its images and with the neutralizing background) and showed the free energy of charging was independent of system size if more than 16 water molecules are used.10 However, significant disagreement existed between Hummer et al.’s results and the corrected PT results9 of Straatsma and Berendsen.5 Ashbaugh and Wood derived a simple correction for the surface effects in the PT method produced by ion-water interactions with truncation based on the ion-oxygen distance.8 By application of these corrections, the results of Straatsma and Berendsen agree well with those of Hummer et al. To explore the utility of Hummer et al.’s method10 as a function of solvent and solute, we derive continuum corrections to the solvation free energies for the ES simulations of ions and ion pairs. These corrections are similar to the energy corrections for point defects in ionic crystals of Leslie and Gillan.11 Figueirido et al.12 and Lynden-Bell and Rasaiah13 proposed first-order corrections to the ES simulations for single ions, which become significant in low dielectric constant solvents. To the first order, our continuum corrections agree
S1089-5647(98)00822-0 CCC: $15.00 © 1998 American Chemical Society Published on Web 06/26/1998
5674 J. Phys. Chem. B, Vol. 102, No. 29, 1998
Sakane et al.
Figure 1. Schematic illustration of the “infinite lattice correction”, ∆AIL-Correction. The gray shading denotes the dielectric medium surrounding the ion. The square lattice denotes the periodic boundary conditions. The diagonal lines indicate the uniform neutralizing background. The negative charges outside the ion are the bound charges produced by the electric field of the ion. The bound charges produced by the uniform neutralizing background are not shown.
with those proposed by Figueirido et al. and Lynden-Bell and Rasaiah. After the present paper was written, Hummer et al. derived the same corrections for single ions including higher order terms by solving Poisson’s equation.14 Their correction for ionic solvation free energies (eq 12 in ref 14) is the same as eq 4 presented here (discussed later). In this paper, we (1) present an alternative derivation to that of Hummer et al.14 for single ions, which shows the physical origin of the corrections, (2) use this method to derive corrections for ion pairs, (3) derive corrections for ∆U, ∆S, and ∆P, (4) explore the magnitude of the corrections for different conditions, (5) derive a new correction that is more accurate for low dielectric constant solvents, (6) demonstrate that, even for rather small system sizes, corrected ES simulations give the same results as PT simulations and SCAAS simulations when appropriate corrections are applied to each method, and (7) show that the polarization of the solvent around an ion is lower in ES simulations than it is at infinite dilution. II. Continuum Correction Methodology A. “Infinite Lattice Correction” Method. We begin with the correction to the simulated free energy of charging a single ion or ion pair with infinite lattice ES (with conducting boundary conditions)3,4 involving the interactions with a neutralizing background charge and periodic images due to PBC. This is the correction to ES simulations with the Hamiltonian including the self-interactions, hereafter referred to as infinite lattice correction (IL), ∆AIL-Correction. The correction to the Helmholtz free energy from simulations at constant volume is derived here, and the correction for constant pressure will be discussed in a later section. Figure 1 shows the schematic illustration of the continuum correction. The correction is just the free energy of charging ion(s) at infinite dilution minus the free energy of charging ion(s) with the ES Hamiltonian with both free energies calculated with the continuum approximation,
∆AIL-Correction ) ∆A∞cont - ∆AES cont
(1)
where ∞ indicates infinite dilution. The first term of eq 1 for a single ion is simply the Born equation,15
∆A∞cont
( )( )
1 z2e2 1 )14π0 2Rion r
(2)
For a pair of ions (ion A and B) separated by rAB (rAB . RA + RB), the first term of eq 1 is approximated as
∆A∞cont )
[
(
)( )]
2 zA2e2 zB2e2 1 zAzBe 1 + 14π0 rABr 2RA 2RB r
(3)
where zRe and RR are the charge and the radius of ion R,
respectively, and r is the relative dielectric constant of the medium. Note that eq 3 and the following continuum correction equations for ion pairs become less accurate when rAB gets close to RA + RB because of the low dielectric constant inside RA and RB.16 Calculation of ∆AES cont for both single ions and ion pairs is derived in the Appendix. The correction to the free energy of charging a single ion for the ES Hamiltonian including the selfinteractions is then
∆AIL-Correction SingleIon
(
)( )
2 2 2 z2e2ξself 1 πz e Rion 1 )+ 1(4) 3 2r 4π0 3L r
where ξself is the self term3,10,17-20
ξself ) lim [ψES[r] - 1/|r|] |r|f0 and ψES is an effectiVe pair potential of the ES,3,4,10,17,19
ψES[r] )
[
1 4π0
∑n
erfc(η|r + n|) |r + n|
∑ k*0
4π
L3|k|2
+
(
|k|2
exp -
4η2
) ]
+ ik‚r -
π
L3η2
(5)
The sums extend over all periodic cells with the real space vector n and Fourier space vector k ) 2 πn/L2 with a convergence parameter η. The last term in eq 5 should be added for the system with a net charge.19 The self term is ξself ) (1/(4π0)) × (-2.837297/L) for a cubic lattice.10,18-20 The corrections of Figueirido et al.12 and Lynden-Bell and Rasaiah13 are just the first term of eq 4. The correction of Hummer et al.14 is the same as eq 4. For an ion pair, the correction is
[ ( )
2 1 zAzBe - zAzBe2ψES[rAB] 4π0 rAB
) ∆AIL-Correction IonPair
]( )
(zA2 + zB2)e2ξself 1 + 2 r
[
]( )
2 2 1 πe (zA + zB)(zARA + zBRB ) 1 1(6) 3 4π0 3L r 2
Simulations with fixed rAB are simpler for the practical application of the equations for ion pairs because ψES is a function of the vector rAB. The correction to the entropy and the pressure change of solvation of ion(s) is calculated from ∆AIL-Correction,
(
∆SIL-Correction )SingleIon
[ ( ( [ (
-
)
)
]( ) ( )
V
2 2 2 z2e2ξself ∂r 1 πRion z e + 4π0 2 ∂T 3L3
)∆PIL-Correction SingleIon -
)
∂∆AIL-Correction SingleIon ∂T
)
∂∆AIL-Correction SingleIon ∂V
)
T
V
1 (7) r2
)
]( ) ( )
2 2 2 z2e2ξself ∂r 1 πRion z e + 4π0 2 ∂V 3L3
T
1 (8) r2
Since the temperature and volume derivatives of the dielectric
Ewald Sum Simulations
J. Phys. Chem. B, Vol. 102, No. 29, 1998 5675
constant of model waters are difficult to obtain, empirical equations of state for water can be used to estimate the magnitude of these corrections. Following directly from the equations above, the corrections for ∆U is given by ∆UIL-Correction ) ∆AIL-Correction + T∆SIL-Correction. B. “Central Cell Correction” Method. Here, we derive an alternative method of correcting ES simulations, the central cell correction (CC). In this method, we use simulations with an ES Hamiltonian to generate an ensemble of configurations representing the water structure inside a sphere of radius R0 around the ion (R0 e L/2). The free energy of charging the ion is calculated using full 1/r potentials inside R0 (not Ewald potentials) and then corrected for the fact that the polarization inside R0 in the simulation is different from what it would be for an ion at infinite dilution, Cin. Then the usual Born correction, Cout, is added for the region outside R0. The correction for the inaccurate polarization around an ion inside R0 is approximated as the continuum calculation of the free energy of interaction of the solvent molecules inside R0 at infinite dilution minus the continuum calculation of the free energy of interaction inside R0 under the ES Hamiltonian with PBC,
Cin ) ∆A∞cont[R0] - ∆AES cont[R0]
(9)
where ∆A∞cont[R0] is the Born free energy of charging an ion that only interacts with solvent between Rion and R0,
∆A∞cont[R0] ) -
( )(
)( )
1 z2e2 1 1 1 14π0 2 Rion R0 r
(10)
∆AES cont[R0] is calculated from the polarization at r[r, θ, φ] from the ion in the radial direction in an ES lattice in a continuum dielectric, PES[r]. Assuming the polarization is proportional to the electric field, P ) 0(r - 1)E, the polarization at r in the radial direction in a dielectric of r is
( )
r 1 PES[r] ) Evac[r]‚ 0 1 |r| r
(11)
where Evac[r] is the electric field in an ES lattice in vacuum. By use of equations defined by Heyes,21 Evac[r] is obtained for a point charge of ze at the center of each cell of an infinite ES lattice without any solvent molecules. ∆AES cont[R0] is calculated from numerically integrating PES[r] over Vin, the volume of the sphere between Rion and R0:
∆AES cont[R0] ) -
ES
( )∫ P r [r] dV
1 ze 4π0 2
Vin
2
(12)
We have found by several numerical tests that an approximate angle-averaged polarization, PES avg[r], suggested in section V.C can be used to simplify the integration and give very accurate results. The angle-averaged polarization PES avg[r] is approximated by assuming that the neutralizing background inside r reduces the electric field and the polarization, and the neutralizing background outside r and in neighboring cells has no effect on the average:
PES avg[r] ≈
( )(
)( )
1 ze0 4πr3 1 112 4π0 r r 3L3
Following eq 12, we have
(13)
∆AES cont[R0] ) -
ze 20
∫RR PESavg[r] dr ) 0
ion
( )[
]( )
1 2π 2 1 1 z2e2 1 - (R - Rion2) 1 (14) 4π0 2 Rion R0 3L3 0 r
We next add the interaction that the ion would have at infinite dilution with the solvent molecules outside R0 in a continuum dielectric, Cout. This is just the usual Born correction for truncated interactions,
Cout ) ∆A∞cont[∞] - ∆A∞cont[R0] ) -
( )( )
1 1 z2e2 1(15) 4π0 2R0 r
The CC method is similar to the PT correction of Wood.9 In the previous nomenclature,9 Cin ) (Css + Ccell) and Cout ) Cis ) CBorn. The total correction is then a sum of Cin and Cout:
) ∆ACC-Correction SingleIon -
( )[
]( )
2π 2 1 1 z2e2 1 + (R - Rion2) 1 (16) 4π0 2 R0 3L3 0 r
C. Correction to the Gibbs Free Energy. We have derived the corrections to the Helmholtz free energy of charging single ions or ion pairs obtained from ES simulations at constant volume. When simulations are performed at constant pressure, these correction equations are no longer accurate because the system box length, L, fluctuates. The exact Gibbs free energy, for instance, for the IL correction of a single ion is calculated as
) ∆GIL-Correction SingleIon
〈
∫01 -
(
)( )〉
2 2 2 λz2e2ξself[L(λ)] 1 2πλz e Rion 1 + 13 r 4π0 3L (λ) r
dλ
λ
(17) where λ is the coupling parameter varying from 0 to 1 along the charging path, and 〈...〉λ denotes the ensemble average at λ. We can approximate this by using the average L at each λ:
≈ ∆GIL-Correction SingleIon
[
∫01 -
(
)( )]
2 2 2 λz2e2ξself[〈L〉λ] 1 2πλz e Rion 1 + 13 r 4π0 3〈L〉λ r
dλ
(18)
As expected, we found from simulations that this approximation introduces a very small error even under conditions where the volume fluctuation is large (discussed in section V.A), but this approximation will not be accurate near the critical point. The corrections to the Gibbs free energy of ion pairs and the CC method can be calculated in the same manner. III. Survey of Continuum Corrections A. Magnitude of Continuum Correction. To explore the magnitude of the continuum correction, we first calculate the corrections to the free energy of charging single ions with an ES Hamiltonian for a number of different conditions with the IL and CC method. The corrections are functions of the dielectric constant of the solvent, the ion charge, the ion radius, and the size of the unit cell (plus the arbitrary cutoff radius R0 for CC). From the free energy corrections for several example
5676 J. Phys. Chem. B, Vol. 102, No. 29, 1998
Sakane et al. TABLE 2: Infinite Lattice Correction for Ion Pairs, IL-Correction (kJ/mol) ∆AIonPair (a) qA ) 1 e, qB ) -1 e b ∆AIL-Correction IonPair
r ) 1.1
r ) 15
r ) 80
L (Å)
RA ) RBa (Å)
rAB ) 5Å
rAB ) 10 Å
rAB ) 5Å
rAB ) 10 Å
rAB ) 5Å
rAB ) 10 Å
10 20 60 100
1.0 1.0 1.0 1.0
106.3 15.4 2.4 1.3
53.1 3.4 1.5
7.8 1.1 0.2 0.1
3.9 0.2 0.1
1.5 0.2 0.0 0.0
0.7 0.0 0.0
(b) qA ) qB ) 1 e b ∆AIL-Correction IonPair
r ) 1.1
Figure 2. Continuum corrections to the simulated free energy of charging a single ion (q ) 1 e, Rion ) 1 Å) with ES for (a) “infinite lattice correction” and (b) “central cell correction” method with a sphere radius R0 ) L/2.
TABLE 1 (a) Infinite Lattice Correction for Single Ions (q ) 1 e), ∆AIL-Correction (kJ/mol) SingleIon L (Å)
Rion (Å)
10 10 20 20 60 60 100 100
1.0 3.0 1.0 3.0 1.0 3.0 1.0 3.0
r ) 1.1
∆AIL-Correction SingleIon r ) 15
r ) 80
179.3 180.4 89.6 89.7 29.9 29.9 17.9 17.9
14.5 25.4 6.7 8.1 2.2 2.2 1.3 1.3
3.9 15.4 1.4 2.8 0.4 0.5 0.2 0.3
L (Å)
RA) RB (Å)
rAB ) 6Å
10 10 20 20 60 60 100 100
1.0 2.0 1.0 2.0 1.0 2.0 1.0 2.0
559.0 560.6 344.9 345.1 119.0 119.0 71.6 71.6
r ) 15
rAB ) 10 Å
rAB ) 6Å
311.6 311.8 118.2 118.2 71.4 71.4
46.4 62.7 26.0 28.0 8.8 8.8 5.3 5.3
r ) 80
rAB ) 10 Å
rAB ) 6Å
rAB ) 10 Å
23.5 25.6 8.7 8.8 5.2 5.3
13.4 30.7 5.5 7.6 1.7 1.7 1.0 1.0
5.0 7.2 1.7 1.7 1.0 1.0
a R and R are ion radii of A and B. Calculations with R ) R ) A B A B 2.0 Å for qA ) 1 e, qB ) -1 e showed that the corrections are not b sensitive to the ion size. rAB is the distance between ions A and B.
(b) Central Cell Correction for Single Ions (q ) 1 e,R0 ) L/2), ∆ACC-Correction (kJ/mol) SingleIon L (Å)
Rion (Å)
r ) 1.1
∆ACC-Correction SingleIon r ) 15
r ) 80
10 10 20 20 60 60 100 100
1.0 3.0 1.0 3.0 1.0 3.0 1.0 3.0
-15.8 -14.7 -8.0 -7.8 -2.7 -2.7 -1.6 -1.6
-162.3 -151.4 -81.6 -80.3 -27.3 -27.2 -16.4 -16.3
-171.7 -160.2 -86.4 -84.9 -28.8 -28.8 -17.3 -17.1
cases (Table 1a), we find that the IL corrections to the free energy are negligible in high dielectric constant solvents except for small L and large Rion; however, the correction is much larger in a low dielectric constant solvent. In contrast to the IL corrections, the CC method yields corrections that are very large for solvents with high dielectric constant but very small for low dielectric constants solvents (Table 1b). The dielectric constant and system size dependencies of IL and CC corrections for a single ion (q ) 1 e, Rion ) 1 Å) are shown in Figure 2. The results show the opposite trends in the IL and CC corrections versus dielectric constant. The correction values of IL and CC cross over around r ) 2; therefore, the IL method is preferred for r > 2. In effect, this means that the CC method is only preferred in very low dielectric constant solvents. Table 2 show the IL corrections to the free energy of charging ionic pairs in high and low dielectric constant solvents. The corrections are negligible as expected for the neutral ionic pair
Figure 3. Continuum corrections (IL method) to (a) the free energy and (b) internal energy of hydration of a sodium ion (Rion ) 1 Å) from ES simulations with 256 water molecules.
(qA ) 1 e, qB ) -1 e) in high dielectric constant solvents and are small even for very low dielectric constant solvents unless the system size gets too small (Table 2a). The corrections are much larger for the charged ion pair of qA ) qB ) 1 e (Table 2b). The system size and dielectric constant dependencies of the corrections are very similar to that of single ions but larger. The results are not sensitive to the distance between ions. Figure 3 shows the IL corrections to the free energy and the internal energy of hydration of a sodium ion using the dielectric constant of water calculated from an equation of state22 as a function of density and temperature. The curves are cut at the saturation density for the temperatures below the critical temperature. The corrections for the free energy and the internal energy have the same trend and become significant above the
Ewald Sum Simulations critical temperature and at low density. In contrast to the internal energy, the correction to the enthalpy (not shown in the figure) becomes very large near the critical point. In principle, the corrections to ∆H, ∆V, and ∆S at constant pressure go to infinity at the critical point of water since (∂r/∂T)P f ∞ and (∂r/∂P)T f ∞ at the critical point. B. Accuracy of Continuum Model. It is well-known that the continuum approximation is inaccurate in the first solvation sphere and that the Born equation is only accurate with empirical ion radii. As a result, the present corrections are only accurate to the extent that the corrections are not sensitive to the ionic radius used. As a test, we calculated the corrections with different ion sizes (Tables 1 and 2). The results are not sensitive to the ion size except for the IL corrections for single ions and non-neutral ion pairs at high dielectric constants with very small system size (in contrast, a 10% change in ion radius produces an approximately 10% change in the Born hydration energy). We conclude that the continuum correction models are accurate enough in most situations because the solvent contribution within the first solvation shell is small. Continuum models are based on the assumption that the charging free energy has a quadratic dependence on the ionic charge as is predicted in the Born model. In compressible solvents, however, electrostriction near the ion changes the quadratic dependence on the ionic charge because the solvent density and dielectric constant near the ion increase as the ion is charged. This is one of the effects that makes the Born equation inaccurate in the first solvation shell. Therefore, the continuum correction is expected to be less accurate where the relationship between the charge and the electrostatic potential at the ion is nonlinear. However, for the test case in dilute steam where electrostriction is very large, we find that the error in the corrected free energy of hydration due to the nonquadratic behavior is very small (see section V.A). IV. Tests of Corrections Using Simulations To assess the accuracy and convergence properties of the various methods of calculating corrections to free energies of solvation of ions, we have performed simulations of single ions in water at two extreme conditions. We have calculated the properties for a sodium ion in water at ambient conditions where the dielectric constant is very high and for a chloride ion in water at low density where the dielectric constant is very low for several different system sizes. The molecular dynamic simulations have been performed with the DL_POLY program.23 The system examined contained a single ion and N water molecules. We studied sodium ion in SPC water24 in a canonical (NVT) ensemble at 298 K and 0.997 g/cm3 by controlling the temperature with a Nose´-Hoover thermostat.25-28 TIP3P water29 was used for the simulation of chloride ion at 623.15 K and 27.36 bar for comparison with previous results for this Hamiltonian using PT.30 The simulations at the low-density condition were performed in the NPT ensemble with Berendsen’s thermostat31 to avoid incorrect bulk densities due to a large cluster of water molecules around the ion under these conditions. The ion-water and water-water interactions are represented by Coulomb interactions between all the charge sites including the partial charges of oxygen and hydrogen on water, and by Lennard-Jones interactions between the ion and water oxygens and between the water oxygens. The Lennard-Jones parameters of ion-water interactions are those of Straatsma and Berendsen5 (Table 3a). The parameters for the Cl--SPC water interaction are used directly for the Cl--TIP3P water interactions. A ° qvist32
J. Phys. Chem. B, Vol. 102, No. 29, 1998 5677 TABLE 3: Ion-SPC Water Interaction Potential Parameters (a) Lennard-Jones Parameters of Berendsen5 x
o,x (kJ/mol)
σo,x (Å)
Na+ Cl-
0.200 546 0.537 866
2.85 3.75
(b) Lennard-Jones Parametersa of A° qvist32 x
o,x (kJ/mol)
σo,x (Å)
Na+
0.086 849
3.247
a
This is calculated for ion-SPC water interactions using the mixing rule as in ref 32.
has shown that his sodium ion-SPC water model can be directly transferred to a sodium ion-TIP3P water model to give very similar results. Electrostatic forces were calculated using the Ewald lattice summation with optimized parameters.23 PBC with a simple cubic lattice and cell length L was used. Real space interactions were truncated at around 0.45L using a Verlet neighbor list. Simulations were performed for an ion for N ) 32, 64, and 128 to check the system size dependence of the results. At ambient conditions, the results were averaged over 100 ps for each λ with a time step of 0.5 fs, and for low densities the results were averaged over 200 ps at λ ) 0 to 2 ns at λ ) 1 with a time step of 2 fs (the autocorrelation time is shorter for smaller λ, where λ is the coupling parameter). Besides the PT results of Wood and Thompson,30 we performed PT simulations of Cl- with 64 TIP3P water molecules at 623.15 K and 27.36 bar in an NPT ensemble with Berendsen’s thermostat. The procedure was the same as that of ES simulations at the same conditions. The free energies were calculated by thermodynamic integration with the polynomial approximations.33 The first and second derivatives of the free energy with respect to λ from the ES simulation at a few values of λ were fitted simultaneously (with appropriate weights) by a polynomial. The Helmholtz or Gibbs free energy of charging the ion, ∆Fsim(λ ) 0 f 1) ) F(1) F(0), is approximated with the nth-order polynomial F(λ) ≈ n anλn. The coefficients, a1, a2, ..., an, are obtained by ∑i)1 fitting the first and second derivatives of F(λ), which are related to the mean and fluctuation of the potential at the ion, (∂F(λ)/ ∂λ)λ ) 〈φ〉λze and (∂2F(λ)/∂λ2) ) -β〈(φ - 〈φ〉)2〉λz2e2, respectively. We performed the simulations for λ ) 0, 1/2, and 1 at ambient conditions and for λ ) 0, 1/4, 1/2, 3/4, and 1 at 623.15 K and 27.36 bar. Third- and fifth-order polynomials were used to fit the simulation results at ambient temperature and 623.15 K, respectively. For the IL method, the mean and fluctuation of φ with the self-interaction term10 are
〈φ〉λ ) 〈φES〉λ + λze〈ξself〉λ
(19)
and
1 〈(φ - 〈φ〉)2〉λ ) 〈φ2ES〉λ - 〈φES〉λ2 - 〈ξself〉λ + β 2λze(〈φESξself〉λ - 〈φES〉λ〈ξself〉λ) + λ2z2e2(〈ξ2self〉λ - 〈ξself〉2λ) (20) where φES is the Ewald potential at the ion and is the sum of the effective pair potential, ψES, defined in eq 5 over all water oxygen and hydrogen sites, 3N
φES )
zieψES[ri] ∑ i)1
(21)
5678 J. Phys. Chem. B, Vol. 102, No. 29, 1998 where zie is the charge of the atom site on water molecules at ri from the ion and ξself is the self term (section II.A). Note that ξself is constant and the last two fluctuation terms of eq 20 are zero for the NVT ensemble. For the CC method, the Coulomb potential at the ion, φ, is obtained by 1/r full electrostatic calculations from the equilibrium configurations of the ES simulations inside R0. In this calculation, an appropriate boundary treatment should be used to avoid surface effects due to molecule-based summation (e.g., the radial integration of charge densities proposed by Hummer et al.34 or molecule-based summations with a correction for the surface effect.8) Here, we integrated over the O-H bond dipole interactions based on the center of the bond dipole to eliminate the surface potential, assuming the dipoles are randomly oriented at a distance R0 away from the ion.35 The dipole-based sorting method was also used to calculate the polarization around a solute to avoid any errors caused by the oxygen-based sorting. In order to correct for the long-range electrostatic effects, it is important to know the dielectric constant of the model water. The dielectric constant of SPC water at ambient condition is 67.8.36 For the low-density high-temperature simulation, the dielectric constant of TIP3P water is estimated by using the Kirkwood equation37 with a Kirkwood g factor of 1.0268 at a density of 0.01 g/cm3 from Pitzer’s estimate for water.38 It is expected that the g factor for TIP3P water in the gas phase will be similar to that of real water in the gas phase. In fact, the calculated dielectric constant is quite insensitive to the value of g. Thus, the dielectric constant of water with a density of 0.01 g/cm3 is obtained as 1.096 48 when g ) 1.0268 and as 1.0939 when g ) 1.00. V. Results of Simulations A. Free Energy of Hydration. Here, we discuss the test results of the continuum corrections to the free energy of hydration of single ions at ambient conditions and at 623.15 K and 27.36 bar with different system sizes (Table 4). The results of IL and CC calculations show excellent agreement (with small errors for CC with N ) 32) even when the correction is as large as 150 kJ/mol. All the corrected results of ES (except for CC with N ) 32) yield the same value, -401 ( 3 kJ/mol at ambient condition and -83 ( 2 kJ/mol at 623.15 K and 27.36 bar. At ambient conditions, the IL corrections for the hydration free energy of sodium ion are about 1% of the free energy of charging from the ES simulations even for N ) 32 (Table 4a). Hummer et al.10 noted that the results with the self-interaction term were essentially independent of the number of water molecules, so they correctly concluded that the corrections are very small. The values for the CC corrections are much larger, but these corrections still give accurate results with N g 64. The Gibbs free energies of charging Cl- in low dielectric constant water obtained by the constant pressure simulations are shown in Table 4b. The continuum corrections for this system are calculated by using the ensemble average of the box size, L, at λ ) 0, 1/4, 1/2, 3/4, and 1 (see section II.C). The estimated errors in the corrections due to the use of average L are very small (less than 1 kJ/mol). The results show that the IL corrections are as large as 35% of the simulated free energies, while the CC corrections are very small. Therefore, CC is the preferred method for these conditions. Unlike the situation at ambient conditions, ∆Gsim for IL with the self-interaction term is system-size-dependent at conditions with low dielectric constant. After the correction, however, the results are accurate even for N ) 32.
Sakane et al. TABLE 4 (a) Hydration Free Energy of Na+ in SPC Watera at Ambient Conditions (kJ/mol) N
∆Asimd
128 64 32 128 64 32 216 512 512
-402(2) -401(2) -403(2) -282(2) -253(2) -197(2) -421(2)e -424(4)e -461(4)e
methodb ES ES ES ES ES ES PT PT PT
IL IL IL CC (R0 ) 6.8 Å) CC (R0 ) 5.2 Å) CC (R0 ) 3.95 Å) Ris ) Rss ) 9 Åc Ris ) Rss ) 9 Åc Ris ) 12 Å, Rss ) 9 Åc
∆Acorrection ∆Acorrected -400(2) -398(2) -399(2) -399(2) -404(2) -391(2) -400(2) -400(4) -396(4)
2.2f 3.1f 4.4f -117.4f -150.8f -194.4f 20.5g 24.2g 64.6g
(b) Hydration Free Energy of Cl- in TIP3P Watera at 623.15 K and 27.36 bar (kJ/mol) methodb ES ES ES ES ES ES PT PT
N IL IL IL CC (R0 ) 30 Å) CC (R0 ) 23 Å) CC (R0 ) 17 Å) Ris ) 12 Å, Rss ) 18 Åb Ris ) Rss ) 23 Åh
∆Gsimd
128 -108(2) 64 -116(2) 32 -131(2) 128 -80(2) 64 -81(2) 32 -85(2) 214 -74.1(0.4)i 512
-79(2)
∆Gcorrection ∆Gcorrected 25.4(0.2)j 32.7(0.2)j 42.6(0.2)j -2.4(0.2)j -2.7(0.2)j -3.6(0.2)j -5.9k
-83(2) -83(2) -88(2) -82(2) -84(2) -89(2) -80.0(0.4)
-3.7k
-83(2)
a
The Lennard-Jones parameters of Straatsma and Berendsen (Table 3a) are used for both Na+-SPC and Cl--TIP3P interactions (see text). b The abbreviations are ES, Ewald sum; PT, potential truncation; IL, infinite lattice correction; CC, central cell correction. c The same values of Ris and Rss are used for the simulations, the free energy calculations, and the continuum corrections. d The free energy of charging calculated from the simulation results. The statistical errors estimated from block averages are given in parentheses. e These are the Gibbs free energy of charging, ∆Gsim, from PT simulations of Straatsma and Berendsen5 at constant pressure. These Gibbs free energies are directly compared with the Helmholtz energies from the ES simulations above, but the difference should be negligible. f The continuum corrections to ∆Asim with Rion ) 1 Å and r ) 67.8. g ∆Gcorrection for PT are the results of Wood9 subtracting the surface energy of Ashbaugh and Wood8 (-79.5 kJ/mol). h The values of Ris and Rss stated were used for the simulation, and Ris ) 20 Å, Rss ) 23 Å were used for the free energy calculations and the continuum correction in order to avoid extraneous effects produced near Ris by the truncation. i The results of Wood and Thompson.30 j The continuum corrections to ∆Gsim with Rion ) 1 Å and r ) 1.096 48. The statistical errors due to the volume fluctuations estimated from block averages are given in parentheses. k ∆Gcorrection for PT were calculated by using continuum correction method of Wood9 and subtracting the surface energy8 (0.7 kJ/mol) at this density.
Figure 4. Average ES potential at Cl- as a function of z ()q/e) including the self-interaction term of Hummer et al.10 (623.15 K and 27.36 bar).
It should be noted that we found a nonquadratic behavior in the charging process at 623.15 K and 27.36 bar (Figure 4). As we stated in section III.B, the continuum corrections are less accurate in compressible solvents. The compressibility effect of solvents, however, becomes negligible as L f ∞. We find that the corrected results at this condition do not depend on L for N g 64, and the three different methods (IL, CC, and PT)
Ewald Sum Simulations TABLE 5: (kJ/mol)
J. Phys. Chem. B, Vol. 102, No. 29, 1998 5679
Hydration Free Energy of Na+ in SPC Water at Ambient Conditions with the Ion-Water Parameters of A° qvista
boundary
method
PBC PBC PBC PBC PBC SCAASb SCAASc SCAASc
ES ES PT PT PT PT PT
IL CC (R0 ) 6.8 Å) Ris ) 10 Å, Rss ) 10 Å Ris ) 10 Å, Rss ) 10 Å/∞d Ris ) 11.5 Å, Rss ) 10 Å/∞d Ris ) 11.5 Å, Rss ) 10 Å/∞d Ris ) 10 Å, Rss ) 10 Å Ris ) Rss ) ∞
∆Asim
∆Acorrection
-358(2) -235(2) -458(1)e -438(2)e -437(3)e -440(2)e -452(1)e -417(2)e
2.2 -117.4 20.2(2.2)f 3.2(0.4)f 2.5(0.1)f 2.8(0.1)f 19.7(2.2)f 0f
∆Asurfaceg
∆Acorrected
79.5 79.5 79.5 79.5 79.5 59.3
-356(2) -352(2) -358(2) -355(2) -355(3) -358(2) -353(2) -358(2)
Table 3b (ref 32). b RSCAAS ) 20.13 Å (ref 39). c RSCAAS ) 14.89 Å (ref 39). d In addition to a solvent-solvent cutoff of 10 Å, all solvents within Ris interact with each other (ref 39). e The simulations results of Åqvist (ref 39). f The continuum correction for solvent-solvent interactions cutoff at Rss (ref 9). g The corrections due to the surface (extrinsic) potential8 (see text). a
give very accurate results ((2 kJ/mol). Therefore, we can conclude that the error in the continuum correction due to the nonquadratic behavior is small under these conditions. B. Consistency with Other Methods. With the agreement of the IL and CC method for ES simulations established, let us consider the other methods: reaction field (RF), potential truncation (PT) with PBC, and surface-constrained all-atom solvent (SCAAS). Hummer et al.’s generalized RF gives the same results as ES,10 so it is correct for high dielectric constant solvents. It seems likely it will have errors similar to those of ES for low dielectric constants and large ions because it is based on a similar neutralizing background. The corrected hydration free energies from ES simulations are compared with the PT results in Table 4. At ambient condition, the PT simulation results of Straatsma and Berendsen5 with the continuum corrections9 and the surface potential correction of 79.5 kJ/mol8 agree with the ES results within 4 kJ/mol for different system sizes. The corrected ES results at 623.15 K and 27.36 bar also agree with the corrected PT simulations for several system sizes within 3 kJ/mol. To compare with the PBC (with PT) and SCAAS results of A° qvist,39 the hydration free energies from ES simulations at constant volume with the sodium ion-SPC water parameters of A° qvist32 (Table 3b) were performed at ambient conditions. Table 5 shows the ES simulations results with IL and CC corrections, and the PBC and SCAAS results39 with the continuum corrections of Wood.9 Since this Hamiltonian model is different from what was used for the previous discussions (different ion-water parameters), the results in Table 5 should not be compared with the free energies in Table 4. After the surface potential corrections estimated below are applied to the PBC and SCAAS simulation results, all the results in Table 5 agree with each other. The surface potential for the PT with PBC and SCAAS calculations in which the ion-solvent interactions are truncated at Ris with oxygen-based sorting is 〈φ〉0 ) -79.5 kJ mol-1 e-1.8 We have assumed that the potential calculations for PBC and SCAAS of A° qvist39 used a cutoff based on the oxygen distance and calculated 〈φ〉0 accordingly. For the SCAAS with no truncation (Ris ) Rss ) ∞), the unknown surface potential on the spherical boundary can be estimated from the difference in 〈φ〉0 for the SCAAS results and for the ES calculations that have no surface potential.8 A° qvist and Hansson40 found that with no truncation of ion-solvent interactions (Ris ) ∞) and local reaction field (LRF) for solvent-solvent interactions outside 10 Å, the potential at the uncharged sodium ion is 〈φ〉0 ) -16.3 kJ mol-1 e-1. They state that LRF gives essentially identical results to no truncations, so it seems likely their 〈φ〉0 for LRF is also the SCAAS result without truncations. We obtained 〈φ〉0 ) 43.0 ( 1.2 kJ mol-1 e-1 from ES simulations
Figure 5. (a) MD simulation and continuum angle-averaged polarization around Na+ at 298 K and 0.997 g/cm3: (- - -) Born polarization; (-) continuum polarization prediction of the ES (r ) 67.8, Rion ) 1 Å); (O) polarization from ES simulation with 512 SPC waters in NVT ensemble. The error bars show the statistical errors estimated from block averages. The error bars for r >10 Å are the size of the symbol. (b) MD simulation and continuum polarization around Cl- at 623.15 and 27.36 bar: (- - -) Born polarization; (-) continuum polarization of ES environment (r ) 1.096 48, Rion ) 1 Å) ; (O) polarization from ES simulation with 256 TIP3P waters in NVT ensemble.
with the Na+-SPC parameters of A° qvist32 (Table 3b). The resulting surface potential for the SCAAS with no truncations is ∆〈φ〉0 ) -16.3 - 43.0 ) -59.3 kJ mol-1 e-1. After these surface potential corrections (∆Asurface ≈ -ze*∆〈φ〉0, assuming a linear charging process) to the SCAAS simulation results are applied, the ES, PT, and SCAAS results agree within 3 kJ/mol for this Hamiltonian. Assuming this estimate of the surface potential is correct, we have agreement among all the methods (ES, RF, PT, and SCAAS) of calculating ions at infinite dilution. C. Polarization. The angle-averaged polarization around an ion from molecular dynamics simulations with ES is shown in Figure 5 and compared with the continuum calculation of the ES polarization and the Born polarization. The ES polarization is obtained by averaging PES[r] of eq 11, and the Born polarization at r is simply obtained from the electric field in a dielectric of r,
PBorn[r] )
( )( )
1 ze0 1 14π0 r2 r
(22)
The polarization around Na+ in the simulation at ambient conditions oscillates around the continuum polarization above
5680 J. Phys. Chem. B, Vol. 102, No. 29, 1998 9 Å (Figure 5a). The continuum model correctly predicts that the simulation polarization is lower than the Born polarization for large r. Figure 5b shows the polarization above 20 Å from Cl- at 623.15 K and 27.36 bar. The simulation polarization approaches the continuum polarization for r > 20 Å. In both simulations, the ES calculations give increasingly lower polarization as the distance from the ion increases and becomes about 1/ of the Born model at r ) L/ , so they should be used with 2 2 caution to explore solvent structure near L/2. Of course PT simulations also produce incorrect polarization near Ris. As mentioned earlier, a very accurate estimate of the average polarization at distance r in the ES calculation can be obtained by reducing the Born polarization by the neutralizing back3 3 Born[r] (eq 13). ground inside r, PES avg[r] ≈ (1 - 4πr /(3L ))P This equation assumes that the neutralizing background outside r and in neighboring cells has no average effect. Numerical tests show that this approximate angle-averaged polarization, PES avg[r], differs by
