J. Phys. Chem. 1994,98, 8256-8257
8256
Reply to Comment on “Transferability of Ion Models” Tami J. Marrone and Kenneth M. Merz, Jr.’ Department of Chemistry, I52 Davey Laboratory, The Pennsylvania State University, University Park, Pennsylvania 16802 Received: December 7, 1993
We recently reported1 a set of parameters for K+ and Na+ which were designed to be used in conjunction with Jorgensen’s methanol (MeOH) model as solvent.2 At the same time we compared our MeOH parameters to a set derived by Aqvist for these ions in SPC3 and TIP3P4 water.5 We compared the performance of our parameter set to that of Aqvist in TIP3P water and in MeOH using periodic boundary conditions6 (PBCs) instead of the SCAAS7 model employed by Aqvist. From this study we found that there was a transferability problem for ion models between simulation protocols and possibly between solvent models. In Aqvist’s comment he states that we areusing ”incorrect potentials” and, furthermore, he suggests that PBCs do not give a realistic representation of ions in solution. We address these comments below. We first address the assertion that we are using “incorrect potentials”. By directly converting Aqvist’s A and E to R* and t, we obtain the self-interaction terms given in Table 1. These represent a direct conversion of Aqvist’s original parameters and are correct as given, but problems arise when the Lennard-Jones (L-J) interaction between the ions and TIP3P (or SPC) water is evaluated. Indeed, the mixing rule employed by Aqvist was a multiplicative one, while AMBER8 employs an additive mixing rule for the evaluation of Ri,*. ( t i , is evaluated using a multiplicative mixing rule within AMBER.8) By directly applying both of these mixing rules to ion/TIP3P Lennard-Jones interactions, we obtain the parameters given in Table 2. One can see that the Na+/TIP3P water parameter is nearly identical regardless of whether an additive or multiplicative mixing rule is used. The problem arises with the K+/TIP3P water parameters which give significantly different Rij* values when an additive or multiplicative mixing rule is used. Hence, using this set, we obtained free energies and solvent structure, which were in poor agreement with Aqvist’s values.’ However, it was realized by us and by Kollman that since the mixing rules differed, we needed to modify the original parameters to take this into account. This lead to the use of an R * / 2 of 2.66 instead of 2.77. Using this new R * / 2 and an additive mixing rule, we obtain the new parameter values given in Table 2, which now are very close to the values obtain using Aqvist’s parameters and a multiplicative mixing rule. Hence, excellent agreement with Aqvist’s original solvation free energies and solvent structure was obtained.’ Are we using the “incorrect potentials” as averred by Aqvist? It is true that there are small differences in the parameters given in Table 2 for Ri,* and ti,, but we feel that these differences are minor and only had a small influence on the outcome of our simulations. Therefore, we have not used “incorrect potentials” as suggested by Aqvist. However, it appears confusion has occurred because of our choice of headings in two of the tables of the original manuscript.’ Thus, Tables IV and V should of labeled the“Aqvist parameter set” as being the“unmodified Aqvist parameter set” or as being the “directly converted Aqvist parameter set”. This would of made it explicitly clear that the data presented in these tables did not employ the Aqvist ion potential functions, but ones in which a different set of mixing rules was used. The OPLS MeOH model uses a multiplicative mixing rule for R* (or u)2+9to determine the solvent/solvent interactions. Since we were interested in studying ionophores using an all-atom 0022-3654/94/2098-8256$04.50/0
TABLE 1: Ion Parameters Sets Na+ K+ K+’ 3.74 5.54 5.32 R’ (4 0.0028 0.00033 0.00033 c (kcal/mol) a Kollman’s modification of Aqvist’s original K+ parameter set. TABLE 2: Ion/TIP3P Water Parameters Sets Evaluated Using Multiplicative and Additive Mixing Rules additive multiplicative Na+ K+ K+* Na+ K+ R$ (A) 3.637 4.538 4.428 3.636 4.426 trj (kcal/mol) 0.02053 0.00706 0.00706 0.02053 0.00706 a Kollman’s modification of Aqvist’s original K+ parameter set. AMBER representation,lOJl we had the difficulty dealing with the inconsistent mixing rules employed by OPLS and AMBER for the solute/solvent interactions. The only difference that arise when the additive mixing rule is used to determine the R,* (or uii> term is for the C-0 Lennard-Jones interaction. (The 0-0 and C - C give identical numbers regardless of which mixing rule is used.) Again AMBER and OPLS use the same mixing rule for ci1.2,9J0 In terms of Rij* (or aij) we obtain values of 3.84 8, (3.42 A) and 3.82 A (3.40 A) for the additive and multiplicative mixing rules, respectively. Hence, the location of the minimum in the Lennard-Jones profile is only shifted by 0.02 A using the additive mixing rule. Furthermore, the value of the repulsive and attractive Lennard-Jones terms vary by 6% and 3%, respectively. From test M D simulation we found that the properties of the MeOH model employing the additive mixing rule were very similar to those obtained using the multiplicative one.2 Hence, instead of developing a new MeOH solvent model, we used the OPLS parameter set employing an additive mixing rule in our studies.’J1 Other schemes have been suggested by Kollman and co-workers to solve this difficulty.l2 The second major issue raised by Aqvist is the use of PBCs versus the SCAAS model developed by Warshel and employed by Aqvist in his studies. Problems encountered through the use of PBCs with ions in solution have been recognized for yearsI3-’5 and will only be summarized here. The first is the over polarization of the ion cutoff sphere, which was originally observed by Straatsma and Berendsen.14 We agree that this is an important issue, but the test simulations reported in our original paper gave relatively large error bars (86.1 f 1.5 kcal/mol with a lo-A cutoff and 83.8 f 1.2 kcal/mol with a 12.5-A cutoff‘). Therefore, we were unable to make a firm decision regarding this issue. We do not view this as an attempt to ignore this issue but as an attempt to be careful in the analysis of our results. Indeed, in the original description of the SCAAS method’ with cutoffs of 8,9, and 10 A the charging free energy for Na+ 4 was -102.4, -98.6, and -99.8 kcal/mol, respectively. The magnitudes of these free energy differences are similar to the ones we observed, and it was concluded that the SCAAS method is size consistent. Later simulations using the SCAAS model showed reduced free energy differences.5 The next issue regards the use of the Born correction in molecular simulations. This model was originally developed to give the solvation free energies of ions using a continuum representation, but in recent years it has been used to correct free energy simulations when long-range interactions are neglected for charged solutes.14 This approximation recovers the free energy resulting from long-range electrostatic interactions, but in our implementation (as well as Aqvist’s) the forces on the ion and water molecules resulting from the long-range interactions are neglected. It is true that in the PBC model there is an interaction between the water molecules within the cutoff with those outside
-
0 1994 American Chemical Society
Comments
The Journal of Physical Chemistry, Vol. 98, No. 33, 1994 8257
of it, which will have an effect on the computed solvent structure and free energies. However, in our case we only used the Born correction to estimate the free energy arising from the neglected long-range interactions and never stated otherwise. Moreover, the SCAAS model does include interaction between water molecules within the sphere with the continuum outside of the sphere.’ This is required to obtain the correct density using the SCAAS model as well as to eliminate surface effects introduced by the SCAAS model itself.’ Hence, it is likely that the SCAASmodel also has an effect on the computed solvent structure and free energies when it is employed due to averaged interaction between water molecules within the sphere and without it.
the quality of results reflected in Aqvist’s original paper. However, it should be stressed that Aqvist never indicated that this was the case.5 In summary, the statement that we are using “incorrect potentials” is not accurate, since the shape of our potentials are for all intents and purposes are identical to those developed by Aqvist, but with an additive mixing rule and not a multiplicative one. Finally, PBCs nor the SCAAS model is perfect, but the former cannot be ruled out until long-range interactions are accounted for using an approach that is better than simply using the Born correction.
Another important consideration regarding the SCAAS model is it ability to handle the dynamics of ion/water systems. King and Warshel’ nicely show that a neat water SCAAS simulation can recover the dynamics of the system, but since the SCAAS model restrains the ion to be in the center of the simulation sphere, quantities like self-diffusion of the ions, velocity autocorrelation functions, etc., cannot be readily evaluated. Another aspect of the SCAAS model is the requirement to carry out a constantvolume simulation. This is not a major obstacle in many cases, but when a constant-pressure simulation is desired PBCs offer the best alternative.
(1) Marrone, T. J.; Merz, Jr., K. M. Transferability of Ion Models. J. Phys. Chem. 1993, 97,6524-6529. (2) Jorgensen, W. L. Optimized Intermolecular Potential Functions for Liquid Alcohols. J. Phys. Chem. 1986, 90, 1276-1284. (3) Berendsen, H. J. C.; Grigera, J. R.; Straatsma, T. P. The Missing Term in Effective Pair Potentials. J. Phys. Chem. 1987, 91, 62896271. (4) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J.; Impey, R. W.; Klein, M. L. Comparison of Simple Potential Functions for the Simulation of Liquid Water. J . Chem. Phys. 1983, 79, 926. (5) Aqvist, J. Ion-Water Interaction Potentials Derived from Free Energy Perturbation Simulations. J . Phys. Chem. 1990, 94, 8021. (6) Allen, M. P.;Tildesley, D. J. Computer Simulation of Liquids; Clarendon Press: Oxford, 1987. (7) King,G.; Warshel, A. ASurfaceConstrainedAll AtomSolvent Model for Effective Simulations of Polar Solutions. J . Chem. Phys. 1989,91,3647. (8) Weiner, S. J.; Kollman, P.A.; Case, D. A.; Singh, U. C.; Ghio, C.; Alagona, G.; Profeta, S.;Weiner, P. A New Force Field for Molecular Mechanical Simulation of Nucleic Acids and Proteins. J . Am. Chem. Sac.
Can we retain the use of PBCs, while correcting for both the long-range energy and force interactions? Through the use of Ewald sums or reaction field methods we can recover both of these components of the long-range interactions without employing spherical models which require the useof a parametrized boundary interaction.6 These methods are still subject to long-range ordering of the simulation *b0xn,l4 but through the use of large sample sizes in the simulation box this can be mitigated to some extent. Studies along these lines appear warranted prior to dismissing the use of PBCs. Finally, it is not surprising that our gas-phase ion-water interactions are low, since our parameters were designed to be used with methanol and not water. Clearly, as we state in the paper, Aqvist has developed a set of ion water parameters that are very good, but care should be taken when they are used outside the SCAAS model. For example, it is clear from the previous discussion that the ion parameters developed by Aqvist will not give the SCAAS results when they are used with PBCs and a cutoff. Hence, the use of these potentials in a simulation of a protein with its associated counterions using PBCs will not give
References and Notes
1984,106, 765-784. (9) Jorgensen, W. L.; Briggs, J. M.; Contreras, M. L. Relative Partition
Coefficients for Organic Solutes from Liquid Simulations. J . Phys. Chem. 1990, 94, 1683-1686. (10) Weiner, S. J.; Kollman, P. A.; Nguyen, D. T.; Case, D. A. An All Atom Force Field for Simulations of Proteins and Nucleic Acids. J. Comput. Chem. 1986, 7, 230-252. (1 1) Marrone, T. J.; Merz, K. M., Jr. Molecular Recognitionof Potassium Ion by the Naturally Occurring Ionophore Nonactin. J. Am. Chem. SOC. 1992, 114, 7542-7549. (12) Thomas IV, B. E.; Kollman, P. A. Free Energy Perturbation
Calculations of the Relative Binding Affinities of an 8-Subunit Cavitand for Alkali Ions in Methanol. J. Am. Chem. Sac. 1994, 116, 3449-3452 and references cited therein. (13) Bader, J. S.; Chandler, D. Computer Simulation of the Mean Force between Ferrous and Ferric Ions in Water. J. Phys. Chem. 1992,96,64236427. (14) Straatsma, T. P.; Berendsen, H. J. C. Free Energy of Ionic
Hydration: Analysisof a Thermodynamic Technique to Evaluate Free Energy Differences by Molecular Dynamics Simulations. J . Chem. Phys. 1988,89, 5876. (15) Valleau,J. P.; Whittington,S. G. AGuide toMontecarloforStatistica1 Mechanics: 1. Highways. In Statistical Mechanics. Part A: Equilibrium Techniques; Berne, B. J., Ed.; Plenum Press: New York, 1977; Vol. 5; pp 137-166.
