J. Phys. Chem. 1994, 98, 4695-4701
4695
Sensitivity Analysis of a Polarizable Water Model Sheng-Bai Zhu IBM Research Division, Almaden Research Center. 650 Harry Road, San Jose, California 951 20-6099 Chung F. Wong' Department of Physiology and Biophysics, Mount Sinai School of Medicine of the City University of New York, New York, New York 10029-6574 Received: October 8, 1993; In Final Form: February 21, 1994'
Recently, we have extended the combined sensitivity analysis/molecular dynamics to identify the most important parameters of the flexible S P C and TIP3P water models in determining the thermodynamical and structural properties of liquid ~ a t e r . ~Since , ~ the S P C and the TIP3P models do not include molecular polarizability explicitly, it is difficult to examine the role of molecular polarizability in determining the properties of liquid water. In this paper, we carry out a sensitivity analysis on a newly developed polarizable water model. We have systematically investigated the possible responses of the structural and thermodynamical properties of this polarizable water model when the potential parameters of the model are perturbed. Discovered is the particularly important role of molecular polarization in describing the properties of liquid water. On the other hand, the use of a point-charge, rather than a dispersed-charge, model appears to be adequate for simulating liquid water. In conjunction with the results obtained from our previous ~ t u d i e s ,these ~ , ~ findings provide useful hints for improving water models.
1. Introduction
It has long been recognized that nonadditive interactions due to many-body and electric polarization effects can play an important role in characterizing liquid water. Quantum chemical calculations5-7 have indicated that at least 10%of the total intermolecular potential energy of a water trimer is resulted from three-body contributions. Clementi and coworkersclo were the first to explicitly incorporate three-body terms into an analytical potential function for molecular dynamics simulations of liquid water. On the other hand, Zhu and Robinson" have introduced a potential function, determined by the relative configuration of the first-nearest neighbors, to better represent intermolecular hydrogen bonding in liquid water. Much efforts have also been spent to incorporate the effects of electronic polarization into water models. It is well-known that water is a polar and highly polarizable molecular liquid. Under normal thermodynamic conditions,the average molecular dipole moment is raised by about 30% from the vapor to the liquid phase and over 40% to the ice-I state.12 Often, the polarization effects are treated in an average way by using effective atomic partial charges in a water model. This treatment inevitably assumes that the instantaneous induced dipole moment of all the molecules in the condensed phase are identical. A recent study by Zhu and Wong4 has shown that this may be a poor assumption. Moreover, the approximation cannot reflect the variation of the average dipole moment with temperature and pressure, not to mention the intensity of an external electric field. Such a defect of the effective charge water models may become problematic in studying the properties of water in an anisotropic environment. This may be the case for water molecules near an interface,13in the interior of a protein,140rnear the surface of a protein or DNA molecule. In principle, the induced dipole moments of water molecules in a system containing N water molecules can be obtained by solving a set of 3N simultaneous linear equations.I5 More often, this set of equations is solved by a time-consuming iteration procedure, starting from an initial guess of the induced dipole *Abstract published in Advance ACS Abstracts, April 1, 1994.
0022-3654f 94 f 2098-4695%04.50f 0
moments pi% and terminated when self-consistency is achieved. A number of schemes have been proposed to reduce the computational effort for determining the induced dipole moment in molecular dynamics simulations. For example, van Belle et al.I4have carried out only one iteration in each time step, starting from the values given in the previous step. Ahlstriim and coworkers16 have devised a predictorxorrector algorithm in which an iterative solution of the induced dipole moments is carried out every five time steps. This was based on the fact that the electric field produced from neighboringwater molecules varies relatively slowly in time. An alternative approach, proposed by Sprik and Klein,I7views polarization as an explicit degree of freedom and treats the fluctuating polarizability as a Drude oscillator.I8 On the other hand, Zhu and c o - ~ o r k e r s ~ ~take , * 0 polarization effects into account by allowing the values of the point charges associated with a water molecule to vary according to the instantaneous local electric field produced by its environment. The variable charge scheme proposed by Zhu and Robinsonl9.20 is simple in computation yet provides stable solutions. It can be used in parametric sensitivity analysis of polarizablewater models. In this paper, we first develop a three-site polarizable water model and then perform a constant (N,V,T) molecular dynamics simulation to investigate the sensitivity of the structural and thermodynamical properties of liquid water to the choice of the potential parameters of the polarizable model. As discussed before,14 there are several reasons for systematically studying the sensitivity of the properties of a molecular system to perturbations of the parameters of a potential model that determines the properties: (1) It can help to identify the parameters of a potential function that play negligible roles in determining the properties of a molecular system and therefore can be removed from the potential model to simplify the model. (2) It guides the refinement of the parameters of a potential function by suggestinghow the parameters of a potential function should be adjusted to get better agreement with a set of experimental data. This can be done in a least-squares fashion in which all the parameters of a potential are changed simultaneously and consistently to best fit the experimental data. This was usually not done in refining the parameters of a water model. Very often, only one or two properties of water were calculated 0 1994 American Chemical Society
4696
The Journal of Physical Chemistry, Vol. 98, No. 17, 1994
as a function of one or two parameters of a potential function. The values of the parameters that gave the best agreement with the one or two experimental properties would then be chosen as the parameters of the potential model for future simulations. This approach had the drawback that only a few experimentallydetermined properties were used and only a few model parameters were adjusted to try to get these few properties right. The use of these model parameters may not predict other properties of water correctly. It is therefore not uncommon that one water model may give reasonable description of one property of liquid water but not another. A better water model can be constructed by refining all the potential parameters simultaneously and consistently to best fit a large amount of experimental data. A sensitivity analysis provides useful information, in the form of sensitivity coefficients, for the consistent refinement of the potential parameters. (3) One can elucidate the properties of a molecular system in atomic and energetic terms by studying how differentcomponentsof interaction forces determine the properties of the system. Such an analysis can also yield useful insights into designing new molecules that present certain desired structural and thermodynamical properties. The focus of this paper, however, is not on refining the parameters of a water model. Instead, we construct a new polarizable water model by using parameters transferred from several existing water models. By examining how small perturbations of these parameters affect different properties of liquid water, one can get a feeling on the relative importance of different components of the potential model in determining the thermodynamical and structural properties of water. From this analysis, one can also examine whether certain features of the polarizable water model can be removed to simplify the model. In addition, we are particularly interested in examining the role of molecular polarizabilityin affecting the simulated propertiesof liquid water. This is harder to study by using the more commonly used effective point-charge models in which the effects of molecular polarizability are implicitly accounted for in an average sense by increasing the magnitudes of the atomic partial charges on the water oxygens and hydrogens. This study naturally extends our previous studies on the SPC and TIP3P water models.3~~
2. Polarizable Water Model In this section, we describe the characteristics of the new polarizable water model. The new model is a hybrid of several existing water models. It is relatively simple, containing only three sites but including molecular polarizability and flexibility explicitly. The equilibrium geometrical parameters are taken from those of an isolated water molecule. In other words, the equilibrium 0-Hbond length equals 0.9572 A and the equilibrium H-0-H bond angle is 104.52'. The partial charges carried by the oxygen and hydrogens are composed of two parts: a constant component, q(O), determined by the permanent dipole moment of the water monomer, and a variable component, q("),which is a function of the instantaneous electric field acting on the water monomer from the surrounding water molecules. These charges are distributed according to the squared Lorentzian function,
taking the oxygens and hydrogens as charge centers. Here, the constant a measures the broadening of the distributed charges. At long distances,these charges interact with each other according to the familiar Coulomb's law. As the separation of the charge centers becomes small, Le., alr - rd
