J . Phys. Chem. 1994,98, 8253-8255
8253
COMMENTS Comment on “Transferability of Ion Models” Johan Aqvist Department of Molecular Biology, Uppsala University, Biomedical Centre, Box 590, S - 75124 Uppsala, Sweden Received: October 4, 1993 In a recent paper,’ Marrone and Merz derive force field parameters for the interaction of Na+ and K+ ions with water and methanol. They employ a parametrization procedure, involving molecular dynamics (MD) free energy perturbation (FEP) simulations of these ions in solution, which is similar to that earlier introduced by US.^ In this type of approach the ion parameters are adjusted to reproduce the absolute free energy of solvation as well as relevant radial distribution functions (rdfs), particularily the first ion-oxygen peak. The solvation free energies are calculated as a sum of two contributions, one from the explicit system that is simulated microscopically and an additional continuum (or Born) contribution due to the use of a finite ionsolvent interaction range. The main difference between the procedure used by Marrone and Merz and ours is that their computations are carried out with periodic boundary conditions (PBC) in the (N,T,P)ensemble using a 10-8, nonbonded cutoff, while we used a restrained spherical boundary, the so-called SCAAS model.3 In this model there is usually no need for an interaction cutoff radius since the finite size of the system limits the number of intermolecular interactions. Here, we will discuss those simulations of ref 1 that were carried out with T I P 3 P a water as the solvent. Although most of our calculations reported in ref 2 used the SPC water model,$we also reported some results with TIP3P. We found that these two water models gave very similar results for AGh,d,(Na+) as well as for the corresponding rdf.2 Marrone and Merz obtain, however, significantly different results in their simulations using our potentials, for both the absolute and the relative AG~dr)s(between Na+ and K+) as well as the rdfs, and they conclude that these discrepancies must be due to different simulation schemes and that ion parameter sets thus may not be generally transferable between such schemes. While absolute AG‘s are likely to be affected by boundary conditions, a conclusion that will be reached below, we find the reported differences in the relative Na+ K+ free energy and the first rdf peaks (Marrone and Merz’s Tables IV and V) surprising since these quantities should be mainly determined by the nearest solvation shells. In order to shed some light on these issues, we have carried out a number of MD/FEP simulations using the SCAAS boundary model as well as PBC in the ( N ,V,T) ensemble. (We will refer to all free energies below as AG‘s although the Helmholtz’ quantity is in some cases the appropriate one.) These calculations which are summarized in Table 1 all employ the same protocol: 10 ps equilibration + 70 ps data collection distributed over 35 different values of the FEPcoupling parameter, X. (The X points wereunevenly distributedfor optimal sampling efficiency, and at each point the first 0.4 ps was discarded.) For each perturbation process two independent trajectories of the above length were calculated, in the forward and reverse directions, in order to assess the convergence. Thus, each entry in Table 1 reflects a total simulation time of 160 ps. For the free energy difference between Na+ and K+ using TIP3P water, PBC, and the same cutoffs as Marrone and Merz, we obtain AAG = 18.1 f 0.1 kcal/mol, which is close to our earlier result using the SCAAS/SPC model with no cutoffs,* AAG = 17.6 f 0.5 kcal/mol. (Marrone and Merz say that they obtain a value of 21.4 kcal/mol using our parameters’.) The rdf
-
peaks that we get for Na+ and K+ are 2.40 and 2.74 A, respectively, again in agreement with ref 2. (Marrone and Merz’s Table IV give the corresponding peaks a t 2.45 and 2.85 A.) Our only explanation for the differing results reported in Tables IV and V of Marrone and Merz is that they did not use our ion-water potentials, but instead used those that would result from applying an additive combination rule, Rij* = ‘/2(Ri* + Rj*), with the Lennard-Jones term expressed in terms of tij and Rij* parameters. While the self-interaction parameters for our ion models are correctly given in Marrone and Merz’s Table 1, the application of an additive (arithmetic) R* combination rule for the ionwater interactions would yield incorrect potentials. We say incorrect here, since our multiplicative (geometric) combination rule is explicitly given in eq 4 of ref 2. It turns out that if the additive rule is applied for the ion-TIP3P interaction, one ends up with “effective” ( A I ,BI)parameters (eq 4 of ref 2) of (144.34, 3.90) and (707.28,5.07) forNa+andK+, respectively. Comparing these values of Table 1 of ref 2, it can be seen that the Na+ parameters are not much altered while those for K+ change significantly and in such a way as to make the ion look bigger. This would thus explain why the both the AAG and the K+ rdf peak are overestimated by Marrone and Merz. As a check, we did a calculation of the relative hydration energy using this incorrect potential, which gave AAG = 21.8 f 0.3 kcal/mol (Table 1, second entry), in agreement with the results of Marrone and Merz. It is also said in ref 1 that “the inability 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”. If this is so, we suggest that maybe also the Kollman group did not use our potentials as specified by the combination rule in eq 4 of ref 2. In fact, the apparent revision of our K+ parameters in the AMBER 4.0 force field6 quoted by Marrone and Merz yields a potential very close to our original one if the additive combination rule is used. Marrone and Merz also show that with this “revision” our original results2 for the free energy difference and K+-water rdf are recovered. This is fine, but to pick out parameters from a set of potentials and subject them to a different mixing rule than they were designed for (Tables IV and V of Marrone and Merz), thereby changing the interaction functions, is not a meaningful test of the original potentials. In this context, one may wonder whether Marrone and Merz are really using the original OPLS potential for methanol, since it also is based on the geometric mixing rule.4b Here, use of the arithmetic rule instead would, e.g., yield a 7% difference in the 0.-CH3 repulsive LennardJones term. Marrone and Merz also claim that the PBC model “is the method of choice for simulations of condensed-phase systems”. This is clearly true for systems possessing real symmetry, such as crystals, and perhaps for uncharged liquids, but for ionic solutions with long-ranged electric fields one may question the adequacy of PBC in their usual implementation (i.e., with relatively small box size and cutoffs). In Table 1 we report the results of a number of different simulations that illustrate the behavior of both spherical boundary conditions and PBC models with respect to some key parameters. First, we can see that for the process of creating/annihilating an uncharged atom with the Na+ Lennard-Jones (LJ) parameters the SCAAS and constantvolume PBC models give very similar results. This is as expected since the corresponding “cavitation free energy” is mainly determined by “packing” interactions (Le., relatively short-ranged solute-solvent van der Waals and solvent-solvent dipole interactions) and does not involve any long-range solvent polarization
0022-3654/94/2098-8253$04.50/0 0 1994 American Chemical Society
8254
Comments
The Journal of Physical Chemistry, Vol. 98, No. 33, 1994
TABLE 1: Summary of Free Energy Calculations. Process model system size (A) box = 21.74 PBCJTIP3P Na+ Kt box = 21.74 PBCJTIP3P Na+ Kt R s c m = 11.34 nil NaO SCAASJSPC box = 21.74 nil NaO PBC/SPC R s c =~11.34 SCAASJSPC NaO Na+ R s c m = 14.89 SCAAS JSPC Nao Na+ box = 34.15 NaO Na+ PBCJSPC R s c m = 14.89 NaO Na+ SCAASJSPC R s c m = 14.89 NaO Na+ SCAASJSPC box = 34.15 PBC/SPC NaO Na+ box = 34.15 PBCJSPC Nao Na+ SCAASJSPC R s c =~20.13 Nao Na+ SCAASJSPC R s c m = 14.89 Nao Na+ R s c m = 1 1.34 NaO Na+ SCAASJSPC R s c m = 11.34 SCAAS/SPC nil Nao
---------
cutoffs (A)
AG (kcal/mol)
comment
+18.1 f 0.1
+21.8 f 0.3 +2.8 0.2 +3.1 f 0.2 -99.7 k 0.3 -99.7 f 0.3 -109.4 f 0.2 -108.0 f 0.1 -108.0 f 0.2 -104.6 f 0.3 -104.5 k 0.7 -105.1 f 0.4 -100.5 f 0.3 -91.3 f 0.9 +2.9 f 0.3
add. mix. rule no cutoffs no cutoffs no cutoffs no pol cons dual R L dual R L dual G,no pc no pol cons Marrone and Merz params Marrone and Merz params
a Each entry corresponds to a total simulation time of 160 ps. All calculations were carried out with a modified versionZof the ENZYMIX programg at 300 K using an MD time step of 2 fs. The ion interaction functions of ref 2 are used in all cases except the second entry, which employs the erroneous potential(s) resulting from the additive mixing rule and the last two entries that correspond to the potentials of ref 1. For all entries denoting ion charging the continuum contributions2are included. The small difference for AGh,d,(Na+) compared to ref 2 probably reflects the longer trajectories used here and possibly also the minor increase in the SCAAS radial force constant in the revised model) compared to the original one.lo The notation R', = 1OJm means that, in addition to a regular water-water cutoff of 10 A, all waters within the ion cutoff sphere (RR) interact with each other. The notation Rfw= m for the SCAAS model is equivalent to Rfw = R ~ C Msince the system is finite. The last two entries, which represent calculations done with the Marrone and Merz parameters (no cutoffs), show that these give an underestimation of the solvation energy together with the SCAAS model (the most often cited experimental value for AGhydr(Nat)appears to be -98.2 kcal/mol).ll
effects. Also, since the radial restraints on surface layer of the SCAAS sphere acts so as to restrict volume changes one would expect similar results to (N,T,V) PBC simulations for this type of process. (Here, one would probably expect (N,T,P)calculations to give a slightly smaller value for the cavitation energy, but this is likely to be a second-order effect compared to the electrostatic energies considered below.) Second, as shown earlier,2,3 the SCAAS model is virtually size independent. This result is again illustrated by comparing the Na+ charging energy obtained for two different sizes of the SCAAS sphere, namely R ~ C A=A1 1.34 ~ A and RSCAAS = 14.89 A, in both cases with no cutoffs (cf. Table 1). Third, it can be seen from the table that a PBC calculation using R:, = R:, = 10 A gives a charging energy for Na+ of AGQ = -109.4 kcal/mol. Thisoverestimation is thus in thesame range as that reported in Table I11 of Marrone and Merz for AGhydrof Na+ using our parameters and TIP3P water. If we now reduce the water-water cutoff, R$w,to 10 A also in the SCAAS model, a similar value of AGQ = -108.0 kcal/mol is obtained regardless of whether the ion-water cutoff is set to R:, = 10 A or R:w = RSCAAS.(In the former case the surface polarization restraints on the water sphere are turned off since the surface molecules are beyond the ion cutoff.) Thus, we conclude that the SCAAS and PBC models appear to converge if a solvent-solvent cutoff is introduced also in the spherical case. An important reason for why AG then becomes more negative is that waters relatively close to the ion do not interact with each other, resulting in an overpolarization of the ion cutoff sphere (e.g., with R L = 10 A, two waters situated diametrically opposite to each other with respect to the ion and slightly more than 5 8, from it will not see each other!). This type of effect has been observed earlier2*'and is, in fact, also reported in ref 1 upon changing R:, from 10 to 12.5 A but regarded by Marrone and Merz as being within their convergence errors. Fourth, it is interesting to see whether the spherical boundary and PBC approaches also converge in some limit closer to the actual SCAAS model with no cutoffs. One can, e.g., easily examine our conjecture concerning overpolarization within the ion cutoff sphere by introducing dual water-water cutoffs so that in addition to a regular R:, of, say, 10 A we let all waters with the ion cutoff sphere interact with each other. We carried out two such simulations using PBC and ion-water cutoffs of 10 and 1 1.5 A, respectively. The corresponding results are given in Table
1, and one now finds that AGQ has risen by about 5 kcal/mol, as expected. As a further check an SCAAS calculation was carried out with similar setup, Le., R L = IO/=, R:w = 11.5 A, and RSCM = 20.13 A (no polarization constraints since the ion does not see the boundary in this case). This simulation yielded AGQ= -1 05.1 kcal/mol, in close agreement with the corresponding PBC results (cf. Table 1). Compared to the results with the standard SCAAS model with no cutoffs, in which all atoms within the SCAAS sphere interact, we however still observe an overestimation of the charging free energy by some 5 kcal/mol. In order to exclude the possibility that it is the surface polarization restraints of the SCAAS model that are responsible for this difference, a calculation using RSCM = 14.89 8, without these restraints was also carried out and yielded AGQ = 100.5 kcal/ mol. Hence, the effect of the surface polarization restraints on the energetics of ion charging seems negligible provided the system size is large enough. We must thereforeconclude that in addition to the overpolarization within the ion cutoff sphere there is another effect that originates from the interaction of solvent molecules inside this sphere with those outside of it (that do not see the ion). Let us now recall the rationale behind the Born correction that is used to account for solvent polarization outside the ion cutoff sphere. This contribution to the free energy is obtained in terms of the integral
representing the electrostatic energy stored in the medium outside the cavity of radius a. Thus, all polarization of the medium outside the microscopic cavity-is taken into account by the dielectric constant relating the E and D fields to each other. It is therefore fundamentally incorrect to both apply the Born correction and at the same time let explicit solvent molecules outside of the ion cutoff sphere interact with waters inside it, since their polarization contributes to the free energy. That is, even though they do not interact directly with the central ion they do see waters within the sphere that are polarized toward the ion. The result of using a solvent-solvent cutoff that lets waters inside the ion cutoff sphere interact with those outside of it, in combination with a Born correction based on R$ thus leads to some extent of double bookkeeping. This was already pointed out in ref 2, and the present results indicate that this effect is
Comments quite important. Although it is true that interactions with the surroundings are also to some extent included, albeit parametrically, in the SCAAS model, we have already noted above that the effect of the polarization restraints on the calculated free energies is rather small. Finally, it might again (see ref 2) be of interest to compare the minimum energies for gas-phase Na+ and K+ monohydrates with the experimental values reported by Kebarle.8 The TIP3P water model together with our ion parameters yields minimum energies of -23.2 and -18.2 kcal/mol for Na+ and K+, respectively. The corresponding values obtained with the Marrone and Merz potentials are -20.8 and -1 5.9 kcal/mol while the experimental result is -24.0 kcal/mol for Na+ and -17.9 kcal/mol for K+.s Here it can also be noted that the energy minimum for K+increases to -17.2 kcal/mol if our parameters are used with the wrong mixing rule. In summary, we find that some of the results reported by Marrone and Merz using our ion potentials are inconsistent with the results presented here and in ref 2, and we suggest that this is due to their using another LJ combination rule than that specified in our original paper. The effect of this is that the K+-water potential changes significantly, and consequently, the Marrone and Merz calculations overestimate AAGh,d,(Na++K+). There is, however, also a real transferability problem between the different boundary treatments which is illustrated, e.g., by the absolute values of AGhydr for Na+. In this case, the use of a different mixing rule has little effect on the resulting potential, and we have shown here that the discrepancies are likely due to two factors related to long-range electrostaticinteractions,namely, the PBC implementation with a small water-water cutoff and an inconsistent continuum correction, both of which contribute with negative terms to the hydration energy. We think that this provides an explanation for the reported difficulty in ref 1 to simultaneously reproduce energetics and structure, in particular since these problems do not appear in the spherical boundary
The Journal of Physical Chemistry, Vol. 98, No. 33, 1994 8255 case where the parameters also are in good agreement with gasphase data. We should perhaps point out that the arguments presented above do not provide a basis for dismissing PBC in general. They do, however, elucidate some rather serious problems with standard PBC implementations for evaluating absolute ion solvation free energies. More sophisticated PBC treatments involving lattice summation or reaction field methods might well solve these problems. But, it is our feeling that the “standard” PBC model (with small box and cutoffs) is often accepted without questioning its physical soundness.
References and Notes (1) Marrone, T. J.; Merz, K. M., Jr. J. Phys. Chem. 1993, 97, 6524. (2) Aqvist, J. J. Phys. Chem. 1990, 94, 8021. (3) King, G.; Warshel, A. J. Chem. Phys. 1989, 92, 3647. (4) (a) Jorgensen, W. L.; Chandrasekhar,J.; Madura, J. D.; Impey, R. W.; Klein, M. L. J. Chem. Phys. 1983,79,926. (b) Jorgensen, W. L. J . Phys. Chem. 1986, 90, 1276. (5) Berendsen,H.J. C.;Postma,J. P. M.;vanGunsteren,W. F.;Hermans, J. In Intermolecular Forces; Pullman, B., Ed.; Reidel: Dordrecht, The Netherlands, 1981; pp 331-342. (6) 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. (7) Straatsma, T. P.; Berendsen, H. J. C. J. Chem. Phys. 1988,89,5876. ( 8 ) Kebarle, P. Annu. Rev. Phys. Chem. 1971, 28,445. (9) Warshel, A.; Creighton,S.In Computer SimulationofBiomolecular Systems; van Gunsteren, W. F., Weiner, P. K., Eds.; ESCOM: Leiden, 1989; p 120. (10) Warshel, A.; King, G. Chem. Phys. Lett. 1985, 222, 124. (1 1) (a) NBS 500, Selected Values of Chemical Thermodynamic Properties;National Bureau of StandardsCircular No. 500;US.GovernmentPrinting Office: Washington, DC, 1952. (b) Rosseinsky, D. R. Chem. Rev. 1965,65, 467. (c) Noyes, R. M. J. Am. Chem. Soc. 1962,84,513. (d) Friedman, H. L.; Krishnan, C. V. In Water, a Comprehensive Treatise; Franks, F., Ed.; Plenum: New York, 1973; p 1. (e) Burgess, M. A. Metal Ions in Solution; Ellis Horwood Ltd.: Chichester, 1978. (12) Valleau, J. P.; Whittington,S. G. In Modern Theoretical Chemistry; Berne, B. J., Ed.; Plenum: New York, 1977; Vol. 5, Part A, p 15. (13) Brooks, C. L. 111; Karplus, M. J . Chem. Phys. 1983, 79, 6312.
