Transferability of ion models - American Chemical Society

7. 0 35. 2/77. 0.00033. Na+. 1.70. 0.040. 1.87. 0.0028. 0 Taken fromref 11. Simulation Procedure ... periodic boundary conditions and a 10-Á nonbond ...
23 downloads 0 Views 663KB Size
J. Phys. Chem. 1993.97, 6524-6529

6524

Transferability of Ion Models Tami J. Marrone and Kennetb M. Merz, Jr.' Department of Chemistry, I52 Davey Laboratory, The Pennsylvania State University, University Park, Pennsylvania I6802 Received: December 14, 1992

Parameters were derived empirically for K+ and Na+ that qualitatively reproduce the absolute free energies of solvation and structural characteristics of these ions in methanol. The procedure used free energy perturbation simulations to obtain the absolute free energies of solvation and molecular dynamics trajectories to examine structural aspects. We find that within the effective two-body model we obtain a reasonable set of compromise parameters. Our results also suggest that we have done as well as we can with an effective two-body potential. Further improvement must come from the use of methods that handle polarization, long-range electrostatic effects, and large simulation sample sizes. The results for these parameters are compared to other parameters (J. Phys. Chem. 1990, 94, 8021) that were derived using a similar strategy (i.e., fitting the free energy of solvation and structural data), but used an alternative simulation protocol. We find that ion parametrization is very dependent on computational strategy, which suggests that ion parameter sets may not be generally transferable.

summation:

Introduction

N Many important biological interactions involve the binding or AG = x[V(A+AA)-V(A)],, (2) transportation of ions. Transport of ions through the cell membrane via carriers or channels affects the polarization of the where Vis the potential energy at the appropriate value of A, N cell membrane which can affect cellular function.' Modeling of is the number of molecular dynamics steps, and AA is 1 / N . these carriers and channels with computer simulations requires During the perturbation simulation, long-range electrostatic reasonable models for both the biomolecules and the ions involved, interactions contribute to this free energy that may not be and much effort has been expended along these lines.24 accounted for in the simulation due to the use of cutoffs.17 One Over the past few decades, several researchers have developed simple way to account for long-range forces is to apply the Born parameter sets to describe the behavior of alkali-metal ions in corre~ti0n.l~This method treats the solvent beyond the cutoff water through molecular dynamics (MD) and Monte Carlo radius ( r ) as a dielectric continuum with a dielectric constant of sim~lations.~-l~ Early efforts centered on using ab initio cale, . In this approximation the free energy contribution beyond culations to garner insight into ion/water interactions and then the cutoff is given as use this data (ion-water distances and interaction energies) to develop ion parameter sets to be used in computer simulation Z2e2 s t ~ d i e s . ~Using J ~ this strategy, attention was paid to obtaining AG = - 8*€# 1(3) a parameter set that reproduced the gas-phase interaction energy and the structural aspects of ion solvation. The weakness in this where 2 is the charge of the ion, e is the electronic charge, and method lies in the lack of attention given to the absolute and €0 is the permittivity of a vacuum. Use of this procedure in ion relative free energies of solvation of these ions and the neglect parameter development was first introduced by Straatsma and of long-rangeelectrostaticinteractionsthat arisein solution. More Berendsen.l7 The Born correction approximates the ion-solvent recently, free energy perturbation (FEP) techniqueshave appeared free energy beyond the cutoff, but it does not consider the free that solve the former weakness by allowing ion parameters to be energy contributions due to solvent-solvent interactions. For developed that reproduce the absolute and relative free energies example, Straatsma and Berendsen have shown that increasing of solvation of the ions in question." the solvent-solventcutoff to a distance larger than the ion-solvent FEP techniques3J3-16have been used by several authors to cutoff makes the calculated free energy of solvation decrease in obtain the absolute free energy of solvation of m ~ l e c u l e s . ~ ~ J ~ -magnitude. ~~ Although large cutoffs for the solvent-solvent The free energy change accumulates as a molecule 'disappears" interactions may be required to include the free energy contrior "appears" in a solvent box to give the free energy of solvation. butions of this interaction, it can be impractical to use this larger The parameters for the initial state (A=O) are perturbed into the cutoff in free energy simulations due to increased computational parameters for the final state (A=l) during the simulation, and time. Thus, the Born correction serves as a crude correction to the free energy is described by the following integral: account for long-range electrostatic forces in the determination of the free energies of solvation of ions.

-(

where AG is the free energy change, V(A)is the potential energy at a particular value of A, and A is the parameter that describes the perturbation. The slow growth method approximates this integral as a summation.3 Assuming the time-dependent perturbation is small, the change in free energy calculated by perturbing one state into another is equal to the following 0022-3654/93/2097-6524S04.~0/0

f)

Herein we present free energy and structural results from ion parameters derived to reproduce the absolute free energies of solvation and structural characteristicsof K+and Na+ in methanol. Furthermore, we compare our ion parameters to those of a recent set derived by Aqvst" using a similar parametrization procedure, but using a different computational approach. Our results indicate some interestingproblems relating to ion parameter development, 0 1993 American Chemical Society

The Journal of Physical Chemistry, Vol. 97, No. 24, 1993 6525

Transferability of Ion Models

TABLE I: Parameter Sets our parameters ion R* (A) e (kcal/mol) K+ 2.17 0.035 Na+ 1.70 0.040 a Taken from ref 11.

TABLE 11: convergence of Absolute Free Energies of Aqvist parameters' R* (A) c(kcal/mol) 2.77 1.87

0.00033 0.0028

Simulation Procedure Thevander Waalsandelectrostaticcomponentsofthe AMBER potential function20were used to model the ion-solvent systems:

where is the total potential energy, the first summation represents the van der Waals interactions, and the second summation represents the electrostatic interactions. Rij is the distance between atoms i and j, cij and Rij* are parameters that describe the shape of the Lennard-Jones potential, qi and qj are the charges of atoms i and j, and e is the dielectric constant. All simulations were done with the AMBER 3.0 revision AZ1 suiteof programs using a constant pressure of 1 atm with a pressure coupling constant of 123 (units: le-O6/bar) and a constant temperature of 298 K.22 The charge used for the ions was +1, and a constant dielectric of one was used throughout. The ions were minimized in a box of 537 methanol molecules23 or 535 TIP3P24water molecules and then equilibrated for 27 ps using periodic boundary conditions and a 10-A nonbond cutoff. All free energy runs were done over 45 ps using periodic boundary conditions. In all cases, the nonbond pair list was updated every 25 time steps and a 1.5-fs time step was used in conjunction with SHAKE.25 To obtain the absolute free energy of solvation, the ion was perturbed into a dummy atom with no charge and no van der Waals interactions. Three runs that were temporally separated by 18 ps of further equilibration were averaged to give the solvation free energy values. The relative free energy of solvation of K+ and Na+ was averaged from a 45-ps perturbation in the forward and a 45-ps perturbation in the reverse direction. After a 63-ps equilibration, coordinates were collected over a 36-ps period for each ion in each solvent. A total of 960coordinate sets were collected in each case. The ion to solvent oxygen radial distribution functions were calculated from the coordinates for each ion in each solvent. The first peak of the radial distribution function was integrated to give the coordination number of the first coordination shell. The parametrization scheme used several production runs of perturbing the ions to dummy atoms in methanol over 45 ps. Parameters were chosen that reasonably reproduced both free energies and structural characteristics of the ions in methanol. In general, we found that parameters that exactly reproduced the free energy of solvation with this method gave poor structure, and parameters that exactly reproduced structure gave poor free energies. Hence, the parameters presented here are a compromise between exactly reproducing the structure and energetics of ion solvation. To test the effect of cut off radius on the calculated free energy of solvation, free energy perturbations were calculated with our K+ model in methanol with a 12.5-A cutoff. The K+ structure equilibrated for 63 ps with a lo-A cutoff in methanol was equilibrated for 18 ps with a 12.5-A cutoff. Two other starting configurations were generated by further equilibrating for 18 and 36 ps, respectively. Free energy perturbations were done over 45 ps to obtain the calculated free energy of solvation for our K+ model in methanol.

Results Table I contains the best set of parameters we found in our parametrization efforts. The other parameter set was derived by

Solvation of E+ in Methanol time of perturbation (pa) 45 90 180

AG (kcal/mol) -67.3 -63.3 -66.6

Aqvistl for K+and Na+ in SPC water with the surface constraint all-atom solvent (SCAAS) model of Warshel and c o - ~ o r k e r s ? ~ ~ ~ ~ In the SCAAS model a spherical boundary is used to model the infinite system instead of the more common periodic boundary condition (PBC) model that usually uses a rectangular boxshape.28 The spherical model is a convenient representation to use, but it suffers from several limitations. First, edge effects can be severe if they are not modeled correctly. The SCAAS model handles this problem effectively with a parametric continuum representation beyond the edge of the simulation sphere. Second, it is presently not possible to do constant-pressure simulations using a spherical boundary model. Finally, PBCs give a very good representation of the solution phase in a nonparametric way (Le., the boundary is handled by explicit solvent-solvent interactions and not by a simplified representation). However, the PBC approach does have its limitations such as an artificially imposed symmetry that can influence dynamics. The PBC approach is very widely used and is the method of choice for the simulation of condensed-phase systems. Therefore, it is very desirable to arrive at ion parameters that employ the usual PBC approach. Theionparametersaregiven inTableI. The Aqvistparameters were converted to Ri* and c values from the Ail and Bit values for the ions. The parameter sets show some substantial differences. The parameters derived using the SCAAS model havevery shallow well depths when compared to the parameters derived here. Furthermore, the spread between the R* values for K+ and Na+ are substantially different from the two models. Hence, clearly the two parameter sets arise from distinctly different minima in parameter space. However, both sets give similar results as seen below. Convergence of free energy simulations is an important issue in simulation^.^^ Therefore, we tested the convergence of the absolute free energy calculations. Table I1 contains these results. Starting from the same initial codiguration, perturbations were done over 45, 90, and 180 ps. The calculation appears to converge at 45 ps since these energies, calculated over various time scales, are within a few kilocalories of each other. The error obtained in our convergence runs is similar in size to the error determined by averaging several runs together at a 45-ps time scale as shown in Table 111. Hence, we feel 45-ps runs not only are computationally convenient but also result in acceptably accurate estimates of the free energy of solvation. The free energy results for Na+ and K+ in methanol and water are given in Table 111. We calculate a AG value using FEP methods and then correct these free energy values for long-range electrostatic effects using the Born equation to arrive at AG-. Comparing AG,, to the average of several derived from experimental data for the free energy of solvation for the ions, AGcxp,both parameter sets overestimate the free energies of solvation. Our parameters overestimatethe AGexpvalue by 11 kcal/mol on average, while Aqvist's do by 16 kcal/mol. Interestingly, the error in energies in water are lower than that seen in methanol. This indicates that the ion parameters should be tested for compatibility in a particular solventprior to use because it may not give similar results in every solvent model. Figures 1-4 contain the ion to oxygen radial distribution functions (rdfs) and the first peak integration for the K+models. Figures 5-8 contain the ion to oxygen rdfs and the first peak integrations for the Na+ models. Comparing the K+ to oxygen rdfs of our parameter set and Aqvist's parameter set (Figures 1 and 2, respectively) shows the latter to have more structure. The first peaks of the rdf are sharper, and the second peaks are more defined. Comparison of the integration of the first peak of the

6526 The Journal of Physical Chemistry, Vol. 97, No. 24, 1993

Marrone and Merz

TABLE III: Solvation Free Energy Resulbs for K+ and Na+ in Methanol and Water ion

solvent

K+

MeOH HzO MeOH

Na+

AG,, (kcal/mol)b -72.7 f 3.3

AG (kcal/mol) -70.0 -65.9 -91.2 -82.9

-75.0 & 3.7 -89.8 f 3.1 -92.4 f 4.0

Hz0

** 0.1 1.5 2.2 * 0.6

AGom (kcal/mol)

AG (kcal/mol) -74.2 i 0.7 - 6 9 0.4

**

-86.1 1.5 -82.3 0.1 -107.3 2.2 -99.3 f 0.6

AG-

(kcal/mol)

-90.3 -85.4 -1 11.2 -108.2

* -95.1 * 2.2 -91.8 * 1.5

* 0.4 1.2 * 2.2 1.5

Reference 11. Methanol values were averaged from values in refs 29 and 31, and water values were averaged from values in refs 29, 30, and 32. 12

1

B

L.

k

0

n

a

I

-Pe C

CI

0

2

-

4

6

a

10

K+ 0 Distance (A) Figure 1. K+ to oxygen rdf for our K+ model in methanol (-) and in water (- -).

-

10

2.0

2.5

3.0

3.5

4.0

4.5

5.0

Distance (A) Figure 3. Integration of the first peak of the K+ to oxygen rdf for our K+ model in methanol (-) and in water (- - -).

a

8-

n

a

Ea

( I

0

6-

I

4-

I )

e

P e

CI

2-

0

2

-

4

6

a

10

2.0

K+ 0 Distance (A)

3.0

2.5

4.0

3.5

4.5

5.0

Distance (A) Integration of the first peak of the K+to oxygen rdf for Aqvist’s

Figure 2. K+to oxygen rdf for Aqvist’s K+ model in methanol (-) and

-4.

in water (- -).

K+ model in methanol

rdfs for K+ (Figures 3 and 4) shows that neither parameter set plateaus in the water case. A plateau indicates a solvation shell. Thus, it is not clear where the first solvation shell ends and the secondone begins. Theintegrationcurverisesslightlymoreswiftly for our parameter set than Aqvist’s parameter set. Although our parameters show poorer structure, they show better energetics. As mentioned previously, we observed similar behavior during the parametrization process. Parameters that gave good free energies gave poor structure, and parameters that gave good structure gave poor free energies. This has led to the compromise parameter sets presented here. Comparison of the rdfs for Na+ (Figures 5 and 6) and the integration of the rdfs (Figures 7 and 8) indicate that both parameter sets show similar structure and energetics in methanol. In water, the structures are similar but the energetics differ by approximately 9 kcal/mol. In this case, our parameters give a better energetic representation. Table IV summarizes the first peak position of the rdfs and the coordination numbers for the ions of both parameter sets in

methanol and water. Experimentally, the first peaks of the rdfs range from 2.6 to 2.95 A for K+ in water and 2.4 to 2.5 A for Na+ in water.34 The coordination numbers range from 6 to 8 for K+ and 4 to 8 for Na+.34 The position of the first peak of the rdf for our parameter set is at the upper limit of the experimental values for water, while the first peak position of Aqvist’s parameter set is in the middle of the range of values. In our case we purposely parametrizedthe ions to give a slightlyhigher value for the position of the first peak, because methanol is bulkier than water and should not get as close to the ion as water. Hence, we obtain first peak positions at the upper limit for that observed for water. However, we find that the first peak position in methanol and water is essentially identical in ail cases. The coordination numbers calculated for both parameter sets agree with the water experimental data. Table V containsthe relative free energyof solvation calculated for each parameter set in methanol and water. Our free energy differences reproduce the experimentally derived relative free

-

(-)

-

and in water (- -).

The Journal of Physical Chemisty, Vol. 97, No. 24, 1993 6527

Transferability of Ion Models 18 I

0

2

6

4

Na+

- 0 Distance (A)

8

2.5

3.0

3.5

Distance

Figure 5. Na+ to oxygen rdf for our Na+ model in methanol (-) and

4.0

4.5

5.0

(A)

Figure 8. Integration of the first peak of the Na+ to oxygen rdf for

' n -

in water (- -). l8

2.0

10

Aqvist's Na+ model in methanol (-)

and in water (- - -).

TABLE n7: Structural Features of K+ and Na+ in Methanol and Water

ion

solvent

K+

MeOH H20 Na+ MeOH H20

grr

n w

@

our parameter set

Aqvist parameter set"

rdf peak (A)

coord

rdf peak (A)

coord

2.95 2.95 2.55 2.55

6 6-7 5.5 6

2.85 2.85 2.45 2.45

6 6-7 5.5 6

experimentb rdf peak

(A)

word

2.6-2.95

6-8

2.4-2.5

4-8

Reference 11. Reference 33.

TABLE V Relative Free Energy Results for K+ Methanol and Water ~~

0

2

6

4

Na+

- 0 Distance (A)

8

~

-

~

-

Na+ io

~~~

our parameter set

Aqvist parameter seta

solvent

AGexp (kcal/mol)b

AAG (kcal/mol)

AAG (kcal/mol)

MeOH H20

-17.2 0.2 -17.4 i 0.4

-17.6 i 0.1 -16.7 k 0.1

-21.1 i 0.4 -21.4 0.1

10

Fgve 6. Na+ to oxygen rdf for Aqvist's Na+ model in methanol (-) and in water (- -).

~

*

Reference 11. Calculated by averaging the differences in the free energy of solvation values in refs 28-31.

~~

98-

165432-

10 2.0

2.5

3.0

3.5

4.0

4.5

5.0

(A) Figure 7. Integration of the first peak of the Na+ to oxygen rdf for our Na+ model in methanol (-) and in water (- - -). Distance

energies of solvation whereas the Aqvist parameter set overestimates the difference in the relative free energies of solvation. Although the absolute free energies are too large, we still obtain excellent agreement with experiment regarding relative free energies.

Discussion Comparing the energetic and structural results of both parameter sets provides insight into the general problems associated with parametrization. Both sets of parameters give slightly different absolute free energiesof solvation and structure using the PBC approach with a lo-A cutoff and the AMBER force field. Comparison of the ability of each parameter set to reproduce the relative free energies of solvation, our parameters doabetterjobinboth thewater andmethanolcases. Theinability of the Aqvist parameters to reproduce the relative free energies within the AMBER force field and PBC approach has previously been recognized by the Kollman group.37 To reproduce the relative free energy of solvation, the potassium R* was adjusted to 2.66 in the AMBER 4.035force field. Free energy results for this adjusted set is found in Table VI. The rdf and integration of the first peak of the rdf for this adjusted set of K+ parameters can be found as Figures 9 and 10, respectively. The relative free energies of solvation calculated using the adjusted K+ set of parameters and Aqvist's Na+ set of parameters agree well the experimental values, and the absolute free energies of solvation of K+in methanol and water are slightly more negative than the values calculated using the original parameter set. The structural features using the adjusted parameter set are similar to the original set, although the adjusted set shows slightly better structure.

Marrone and Merz

6528 The Journal of Physical Chemistry, Vol. 97, No. 24, 1993 TABLE VI: Structural and Energetic Data When R+ Adjusted to 2.66 structural features free energy of solvation K+ .-,N ~ + rdf Peak AG AGm AAG ion solvent (A) word (kcal/mol) (kcal/mol) (kcal/mol) K+ MeOH 2.75 6 -78.4f 1.5 -94.5 h 1.2 -17.7 f 0.1 H20 2.75 6-7 -73.5 f 0.4 - 8 9 . 9 h 0 . 4 -17.5 fO.l ~

12

L

n

i

I

6

II

2

0

I

-

a

6

4

K+ 0 Distance

10

(A)

Figure 9. K+ to oxygen rdf for adjusted K+model in methanol (-) and

in water (- - -). 10

A

i

a

a pc (r

0

2.0

2.5

3.0

3.5

4.0

4.5

5.0

Distance (A) Figure 10. Integration of the first of the K+to oxygen rdf for the adjusted K+ model in methanol (-) and in water (- - -). Interestingly, the parameters determined by Aqvist" using a similar parametrization scheme, but a different simulation protocol and model, gave very good energetic and structural results. However, when transferring these parameters to our computational procedure, the results were not as good. Aqvist argues that periodic boundary conditions, which were used here, will not give reasonable results because of the truncation effects observed by Straatsma and Berendsen in their simulations on ion hydration. However, this may not be the full story, and other factors like polarization effects may be coming into play. Furthermore, the SCAAS model has adjustable parameters that may allow for some of these effects. Therefore, it is difficult to compare these models directly. To address the truncation issue, we examined the effect of varying the cutoff radius during the free energy simulation. The calculated free energy of solvation for our K+model in methanol gave a value of -70.0 1.5 kcal/mol when a 10-A cutoff was

used and a value of -70.9 i 1.2 when a 12.5-A cutoff was used. The Born correction is -16.1 kcal/mol for a lo-A cutoff in methanol and -12.9 kcal/mol for a 12.5-A cutoff. The Born correction predicts a difference of -3.2 kcal/mol for these two systems. The prediction appears to hold true within the errors of our simulation. However, now we have increased the solventsolvent cutoff as well, which is thought to contribute a positive free energy." Since the change in free energy observed by changing the cutoff was only a small fraction of the total energy, we feel other factors such as many-body effects probably play a more crucial role than truncation effects provided the sample size and cutoff are large. Because we observed different results for parameters obtained using different simulation protocols, we think that it is important to define the problem the parameters will be used to described. Wesuggest that the parametrization procedureused to reproduce the free energy of solvation and structure of the ions should be similar to the protocol used in the ensuing simulationsusing these parameters. For example, if the parameter set for the ion is to be used in a simulation that uses PBCs, then the free energy simulation used to parametrize the ion also should use PBCs. Our parameter set was derived to reproduce the absolute free energies of solvation and structural characteristics of K+ and Na+ in methanol within the effective two-body potential. The ions were parametrized using periodic boundary conditions with a 10-A nonbond cutoff. Aqvist's parameter set was derived to reproduce the absolute free energies of solvation and structural characteristics of ions in SPC36 water within the context of the SCAAS Both models used the Born correction to account for long-range electrostatic interactions. The Aqvist parameter set reproduces the energetic and structural properties within the context of its parmetrization. However, this parameter set does not reproduce the absolute free energy of K+ and Na+ in TIP3P water using our scheme with periodic boundary conditions and a 10-A cutoff. This is even more surprising since it has been shown that Aqvist's parameter set reproduces the absolute free energy of Na+ in TIP3P using the SCAASmodel.' Thus, thedifferences cannot be attributed to the water model but the parametrization and simulation scheme. Presently,our ion parameters are the best parameters available for simulations of K+ and Na+ in Jorgensen's methanol model that use periodic boundary conditions within the effective twobody potential. They reasonably reproduce the absolute free energy of solvation and structural characteristics of these ions in methanol. They also reproduce energetia and structural features in TIP3P water. However, their performance would need to be tested before being used in other solvent models, because they may not be transferable. Our parametrization procedure could be improved by using a more rigorous procedure to handle long-range forces. We chose the Born correction with one cutoff because it is very simple to use. Multiple cutoffs for ion-solvent and solvent+lvent interactions could have been used but would have required parametrizationfor the cutoffs. Simulation of these ions with a biological molecule could require parametrization of the biological molecule's nonbond cutoffs. Better schemes could include procedures to handle long-range forces such as the Ewald sum or the reaction field method." Calculationsthat include polarization would raise the level of sophisticationof the calculation but would not include contributions from beyond the cutoff. Free energy perturbations that included polarization and a procedure to include long-range forces could possibly provide a very accurate free energy of solvation of an ion, although this has yet to be proven.

Conclusions We have derived a set of parameters for K+ and Na+ in an effective two-body methanol model. The performance of the ion parameters is reasonable and from our experiences suggests that

Transferability of Ion Models the present parameter sets are the best that can be obtained using two-body potentials and periodic boundary conditions. We also find that the performance of ion parameters is very dependent on the parametrization scheme used. We suggest that ion parameters should be developed that employ a simulation protocol that is planned to be used in future studies. For example, if the ion parameters are going to be used to study ionophores using periodic boundary conditions with a lo-A cutoff, the parametrization protocol should also use these simulation protocols. To obtain further improvement on the present ion parameters, it appears that both polarization and long-range electrostatic interactions need to be better represented. Acknowledgment. This research was supported by the ONR (N00014-90-3-4002) and ACS-PRF (23225-G6,4), and their support is acknowledged. We also thank the Pittsburgh Supercomputer Center for CRAY YMP time and the Center for Academic Computing for IBM 3090 time. References and Notes (1) Dobler, M. Ionophores and Their Structures; John Wiley and Sons: New York, 1981. (2) Pullman, A. Chem. Rev. 1991,91,793. Contribution of Theoretical Chemistry to the Study of Ion Transport Through membranes. (3) Kollman, P. A.; Merz, K.M. J. Acc. Chem. Res. 1990,23,246252. Computer Modeling of the Interactions of Complex Molecules. (4) Eisenman, G.; Alvarez, 0. J . Membr. Biol. 1991, 119, 109-132. Structure and Function of Channels and Channelogs as Studied by Computational Chemistry. (5) Lybrand, T. P.; Ghosh, I.; McCammon, J. A. J . Am. Chem. Soc. 1985, 107, 7793-7794. Hydration of Chloride and Bromide Anions: Determination of Relative Free Energy by Computer Simulation. (6) Chandrasekhar, J.; Jorgensen, W. L. J . Chem. Phys. 1982,77,50805089. The Nature of Dilute Solutions of Sodium Ion in Water, Methanol, and Tetrahydrofuran. (7) Impey, R. W.; Madden, P. A,; McDonald, I. R. J. Phys. Chem. 1983, 87, 5071-5083. Hydration and Mobility of Ions in Solution. (8) Chandrasekhar, J.;Spellmeyer, D. C.; Jorgensen, W. L.J. Am. Chem. SOC.1984, 106, 903-910. Energy Component Analysis for dilute aqueous solutions of Li+, Na+, F-, and Cl-. (9) Kollman, P. A,; Wipff, G.; Singh, U. C. J. Am. Chem. Soc. 1985, 107,2212-2219. Molecular MechanicalStudiesof Inclusionof Alkalications into Anisole Spherands. See also references cited therein. (10) Migliore, M.; Corongui, G.; Clementi, E.; Lie, C. G. J. Chem. Phys. 1988,847766. Monte Carlo study of free energy of hydration for Li+, Na+, K+, F-, and C1- with ab initio potentials. (11) Aqvist, J. J. Phys. Chem. 1990, 94, 8021. Ion-Water Interaction Potentials Derived from Free Energy Perturbation Simulations. (12) Kistenmacher, H.; Popkie, H.; Clementi, E. J. Chem. Phys. 1973,58, 1689. Study of the Structure of Molecular Complexes 11: Energy Surfaces for a Water Molecule in the Field of a Sodium or Potassium ion. (13) Zwanzig, R. W. J. Chem. Phys. 1954,22,1420. High-Temperature Equation of State by a Perturbation Method. I. Nonpolar Gases. (14) van Gunsteren, W. F. Protein Eng. 1988, 2, 5-13. The Role of Computer Simulation Techniques in Protein Engineering. See also references cited therein.

The Journal of Physical Chemistry, Vol. 97,No. 24, 1993 6529 (15) Mezei, M.; Beveridge, D. L. Ann. N.Y.Acad. Sci. 1986,482, 1-23. Free Energy Simulations. See also references cited therein. (16) Jorgensen, W. L. Acc. Chem. Res. 1989,22,184-189. Free Energy Calculations: A Breakthrough for Modeling Organic Chemistry in Solution. See also references cited therein. (17) Straatsma, T. P.; Berendsen, H. J. C. J . Chem. Phys. 1988,89,5876. Free Energy of Ionic Hydration: Analysis of a Thermodynamic Technique to Evaluate Free Energy Differences by Molecular Dynamics Simulations. See also references cited therein. (18) Jorgensen, W. L. J. Chem. Phys. 1988, 89, 3742-3746. Efficient Computation of Absolute Free Energies of Binding by Computer Simulations. Applications to the Methane Dimer in Water. (19) Merz,K.M.,Jr.J.Am.Chem.Soc.1991,113,40~11. COzBinding to Human Carbonic Anhydrase 11. (20) Weiner, S.J.; Kollman, P. A.; Nguyen, D. T.; Case, D. A. J. Comput. Chem. 1986,7,230-252. An All Atom Force Field for Simulations of Proteins and Nucleic Acids. (21) Singh,U.C.; Weiner,P.K.;Caldwell, J. W.;Kollman,P.A.AMBER(UCSF), Version 3.0; Revision A by G. Seibel. AMBER. (22) Berendsen, H. J. C.; Potsma, J. P. M.;van Gunsteren, W. F.; DiNola, A. D.; Haak, J. R.J. Chem. Phys. 1984,81,3684-3690. Molecular Dynamics with Coupling to an External Bath. (23) Jorgensen, W. L. J . Phys. Chem. 1986,90, 1276-1284. Optimized Intermolecular Potential Functions for Liquid Alcohols. (24) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J.; Impey, R. W.; Klein, M. L. J . Chem. Phys. 1983, 79,926. (25) Ryckaert, J. P.; Ciccotti, G.; Berendsen, H. J. C. J. Comput. Phys. 1977,23,327-41. NumericalIntegrationoftheCartesianEquationsofMotion

of a System with Constraints: Molecular Dynamics of N-alkanes. (26) Warshel, A.; King,G. Chem. Phys. Lett. 1985,121,124. Polarization Constraints in Molecular Dynamics Simulations of Aqueous Solutions. The surface constraint all atom solvent (SCAAS) model. (27) King, G.; Warshel, A. J. Chem. Phys. 1989,91, 3647. A Surface Constrained All Atom Solvent Model for Effective Simulations of Polar Solutions. (28) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Clarendon Press: Oxford, 1987. (29) Mitchell, M. J.; McCammon, J. A. J . Comput. Chem. 1991, 12, 27 1-275. Free Energy Difference Calculations by Thermodynamic Integration: Difficulties in Obtaining a Precise Value. (30) Marcus, Y. Ion Soluation; John Wiley: New York, 1985. (31) Gomer, R.; Tryson, G. J . Chem. Phys. 1977, 66, 44134424. An Experimental Determination of Absolute Hall Cell EMFs and Single Ion Free Energies of Solvation. (32) Burgess, J. Metal Ions In Solution; Ellis Horwood Limited: Chichester, 1978; pp 175-219. (33) Noye.s,R.M.J.Am. Chem.Soc. 1%2,84,513-522. Thermodynamica of Ion Hydration as a Measure of Effective Dielectric Properties of Water. (34) Marcus, Y. Chem. Rev. 1988, 88, 1475. Ionic Radii in Aqueous Solution. (35) AMBER 4.0. Pearlman, D. A.; Case, D. A.; Caldwell. J. C.; Seibel, G. L.; Singh, U. C.; Weiner, P.; Kollman, P. A. University of California, San Francisco. (36) Berendsen,H. J.C.;Pastma, J.P.M.;vanGunsteren,W. F.;Hennans, J. In Intermolecular Forces; Pullman, B., Ed.; Reidel: Dordrecht, The Netherlands, 1981; pp 331-342. (37) Kollman, P. A. Personal communication.