460
J . Phys. Chem. 1990, 94, 460-466
I ncorporat ing Electric Polarizabilities in Water- Water Interact ion Potentials Satoru Kuwajima* and Arieh Warshel* Department of Chemistry, University of Southern California, Los Angeles, California 90089- 1062 (Received: March I , 1989; I n Final Form: July 18. 1989)
Water-water interaction potentials that include many-body corrections through classical electric polarizabilities and the corresponding induced dipoles are investigated. Significant attention is given to the transferability of gas-phase dimer potentials to condensed phases. The potential functions used are parametrized by using information about the gas-phase dimer potential surface from ab initio CI calculations (the raw data used for the parametrization of the MCY model, modified for its incorrect behavior at long distances). The potential thus obtained is examined and tuned by extensive calculations of the crystal structures and lattice energies of proton-ordered polymorphs of ice (ices 11, IX, and VIII). These crystal calculations involve a novel energy-minimizationprocedure for the treatment of the induced dipoles and their energy contributions. Some of the popular water-water potentials (MCY, TIPS2,and SPC) are also examined for comparison. Our potentials reproduce fairly well experimental crystal structures and lattice energies, demonstrating that the transferability problems suffered by pairwise-additive potentials such as the MCY model can be removed by explicitly including polarizabilities and induced dipoles.
Introduction
by using experimental results. Nevertheless, while pursuing a combined refinement using experimental and theoretical conModel potentials that mimic microscopic molecular interactions straints, it is possible to examine the transferability problem using are an essential element in computer simulation of chemical and the available surfaces. Thus we demonstrate in this work that biological processes. While the majority of current simulation the transferability of a potential, based on a reasonable gas-phase studies employ pairwise-additive potentials (which are hereafter a b initio surface (constrained to have its correct asymptotic form), called pair-interaction models) the interest in using more versatile can be significantly improved once the effects of induced dipoles models is steadily growing, owing in part to the rapid increase are included. in computer power. A class of such improved potentials can be Among the available pair-interaction potentials, we consider obtained by taking into account electric polarizabilities of molthe raw a b initio C I data used for the construction of the Matecules and the corresponding induced dipoles.’V2 The corresponding suoka-Clementi-Yoshimine (MCY) model8 to be quite appromodels (referred to in this paper as polarizable-atom models) priate for our purpose despite its limitations (see below). Alrepresent more than a formal exercise; first, they include some ternative high-quality calculation^^^ do not provide a convenient of the three-body effects missing in pair-interaction models, which dimer potential surface but only few points on this surface. might be important in consistent treatment of ground-state Although the MCY model has been successful in reproducing properties (although such effects may be compensated for by various quantities of liquid water,9J0it predicts a density too low empirical parametrization procedures). Second, the inclusion of (reportedly” by about 24%) at atmospheric pressure or, correinduced dipoles effects is a crucial element of consistent and spondingly, a pressure toO high (3000-7000 atom9vi2)at the normal reliable calculation of excited-state processes’ and may also be density. The average interaction energy is also reported to be too essential for consistent studies of charge-separation processes. high by about 1 kcal/mol ( 14%).9These deficiencies are a typical Thus it is important to obtain more refined polarizable-atom manifestation of many-body effects, which cause poor transfermodels, and a systematic step in this direction is reported in the ability of model potentials from one phase to another (studies of present work. the many-body effects in water interactions that included trimer Polarizable-atom models of various forms have already been and used analytical three-body potentials were re* ~ ~ ~ * calculations ~~~ proposed‘” and used in simulations of ions in s o l ~ t i o n s ~and ported recently). 3 ~ 1 4 in simulations of liquid However, it appears that the As to the gas-phase surface of the water dimer, although the inclusion of induced dipoles and the resulting decrease in perMCY surface is not perfect, it seems to do quite a good job; the manent dipoles makes it harder to maintain the correct short-range potential might underestimate the 0-0distance by 0.1 A but this structural features of liquid water using simple potentials with is far from being certain since the observed distance might reflect only three interaction centers.’ While developing improved pozero-point and thermal extension associated with the shallow dimer larizable-atom potentials we take in this work an intermediate potential. At any rate the 0-0 distance for the MCY gas-phase step, examining whether polarizable-atom models can retain an surface appears to be better than that given by current potential interaction potential that consistently describes the gas-phase dimer functions (see Table 11). Regardless of the possible improvement (e.g., correct monomer dipole moment) while providing reasonable of the MCY surface what we try to show here is that the results for the condense phase. transferability of this (and any other reasonable gas-phase dimer I n examining the above transferability problem one may like surface) can be improved very significantly by using polarizato take the pair-interaction potential directly from gas-phase ab ble-atom models. initio calculations. However, although such calculations are now approaching quantitative level for systems of few atoms, they are not sufficiently reliable to be taken as the ultimate source of (8) Matsuoka, 0.;Clementi, E.; Yoshimine, M. J . Chem. Phys. 1976, 64, information. Instead i t seems preferable to consider calculated 1351. (9)Lie, G. C.; Clementi, E.; Yoshimine, M. J . Chem. Phys. 1976,64, 2314. surfaces as only part of the data base that should be augmented -
( 1 ) (a) Warshel, A.; Russell, S. T. Q.Reu. Biophys. 1984, 17. 283. (b) Warshel, A,; Levitt, M. J . Mol. Biol. 1976, 103, 227. (2) Stillinger, F. H.; David, C. W. J . Chem. Phys. 1978, 69, 1473. (3) (a) Warshel, A. J . Phys. Chem. 1979, 83, 1640. (b) Warshel, A. J . Phys. Chem. 1982, 86, 2218. (4) Lybrand, T. P.; Kollman, P. A. J . Chem. Phys. 1985, 83, 2923. ( 5 ) Barnes, P.; Finney, J . L.; Nicholas, J. D.: Quinn, J. E. Nature 1979, 282, 4.59. (6) Rullmann, J. A . C.: van Duijnen, P. Th. Mol. Phys. 1988, 63, 451. (7) King, G.and Warshel, A . J . Chem. Phys., in press.
(10) Impey, R. W.; Klein, M. L.; McDonald, I. R. J . Chem. Phys. 1981, 74, 647. ( I I ) Owicki, J . C.; Scheraga, H. A. J . Am. Chem. SOC.1977, 99, 7403. (12) Mezei, M.: Beveridge, D. L. J . Chem. Phys. 1982, 76, 593. ( 1 3) (a) Clementi, E.; Kolos, W.; Lie, G. C.; Ranghino, G. Int. J. Quantum Chem. 1980, 17, 377. (b) Clementi, E.; Corongiu, G . I n t . J . Quanrum Chem. B i d . Symp. 1983, I O , 31. (14) Yoon, B. J.; Morokuma, K.: Davidson, E. R. J . Chem. Phys. 1985, 83, 1223. ( 1 5 ) Berendsen, H. J. C.; Postma, J. P. M.; van Gunstren, W. F.; Hermans, J. I n intermolecular Forces; Pullman, B., Ed.; Reidel: Dordrecht, Holland, 1981; pp 331-342.
0022-3654/90/2094-0460$02.50/00 1990 American Chemical Society
The Journal of Physical Chemistry, Vol. 94,No. 1, 1990 461
Water-Water Interaction Potentials
between the interaction centers, and the summations on H’ runs over the two H atoms of each molecule. The induced dipole moment of each water molecule is taken to be a point dipole located on the oxygen atom (no dipole on the H atoms). The induced dipole & of the kth molecule is given by b k = d E Q , k + Ep,k) (3)
0
H
H
Figure 1. Water monomer in the polarizable-atom models; the 0-H distance is 0.9572 A and the H-0-H angle is 104.52O.
The present work does not engage in liquid-water simulations. Instead, proton-ordered polymorphs of ice are used as the representative test for the performance of water-water interaction potentials in condensed phases. As reported by Morse and Rice16 (MR), the transferability problems of the MCY model in ice are quite similar to those in liquid water. Moreover, molecular crystals can provide detailed information about intermolecular forces and have been used extensively in the optimization of a wide variety of force fields ( e g , ref 17 and 18). Further extension of such approaches to the refinements of force fields for charged and highly polar molecules should include induced dipoles, and the present study provides a systematic step in that direction. The paper starts by outlining our polarizable-atom potential and describing the parametrization procedure (which involves a b initio and experimental information about the gas-phase dimer surface). This is followed by a description of the nevel procedure used for evaluations of the first derivatives of our potential, which are needed for effective energy minimization. The resulting model is then tested by extensive calculations of crystal structures and lattice energies of proton-ordered ices. This is used for the fine tuning of a single shielding parameter, keeping the gas-phase dipole moment and the atomic polarizabilities (which are based on experimental information) at their gas-phase values. Finally, the performance and validity of the model are discussed and the importance of including induced dipoles in consistent simulations of various forms of water and aqueous solutions is emphasized. Model Potentials Our model potentials are composed of pairwise-additiveenergies (electrostatic and van der Waals interactions) and the energy associated with the induced dipoles
U = C cpair(k.k3 k>k’
+ Ccind(k) k
(1)
where k and k’ refer to water molecules. Following the MCY model, we adopt the four-center model for the pairwise-additive interactions, placing a negative charge at a point (M) on the bisector of the HOH angle (Figure 1). Unlike the MCY model we use standard exponential-6 type potentials for the van der Waals interactions. For simplicity, the interaction center of the van der Waals attractions is placed only on the oxygen atom. The water molecule is treated as a rigid body, keeping the atoms at their experimental geometry19 (see Figure 1). The pairwise-additive part of the potential then takes the form
where am, boo, etc., are constants, roo, etc. denote the distance (16) Morse, M. D.; Rice, S. A . J . Chem. Phys. 1982, 76, 650. (17) (a) Warshel, A.; Lifson, S. J . Chem. Phys. 1970, 53, 582. (b) Huler, E.; Sharon, R.; Warshel, A. QCPE 325 Quantum Chemistry Program Exchange, Indiana University, 1976. (c) Warshel, A. In Modern Theoretical Chemistry (Semiempirical Methods); Segal, G . A,, Ed.; Plenum Press: New York, 1977; p 133. (18) Hagler, A. T.; H u h , E.; Lifson. S. J. Am. Chem. Sac. 1974, 96, 5319. (19) Benedict, W. S.;Gailar, N.; Plyler, E. K . J . Chem. Phys. 1956, 24, 1 139.
where 01 denotes the (isotropic) polarizability of the water monomer and &,k EP,kgives the electric field at the oxygen position, with E , , and EP,krepresenting the contributions from the permanent (partial) charges and induced dipoles, respectively. The anisotropy of the polarizability is neglected here and the fields produced in each molecule by its own partial charges are not included in the field E , , of eq 3. An important feature of our model is the modification of the field &J, from the usual Coulomb field. As discussed later, it is necessary to introduce a distance-dependent shielding factorfs(r) in order to reduce the field strength. The EQ,k fields in our model are expressed as
+
where rOHt and TOM, denote the distance between the 0 atom of the kth molecule and the H and M centers of the k‘th molecule, respectively, with rOH,and rOM,representing the corresponding vectors. We employ the shielding factor of the form where rs is a constant. While we consider that the use of the shielding factorsf,(r) has a physical basis (see ref 2 in which a similar factor was employed), they might as well be regarded as a practical invention for the time being. In order to avoid further complexity, we do not use any sheidling factor for the EP,kfields, which are given by Ep,k = [-Wkr/roo> + 3 ( b k ~ ’ r o o ) ~ 0 o ’ / ~5 1o o (6) k’#k
The induced-dipole energy
tind
in eq 1 is given by’
Note that the EP,kfield is not included in the right side of eq 7 . As stated in the introduction to this paper, we would like to start the refinement of the above model using information about the gas-phase dimer. In doing this it is important to keep those parameters which can be deduced uniquely from experimental studies as nonadjustable parameters. Thus we fixed the charges (qH and qo) and the parameter rOMusing the observed gas-phase dipole moment2’ (1.855 D) and used the experimental polarizability a = 1.42 A3 for the induced dipoles of eq 3. Furthermore, the constant cm of the van der Waals attractions in eq 2 is fixed to the value obtained by the well-known London approximation, cm = (3/4)Zpa2, with the experimental polarizability and ionization potential Zp = 12.59 eV. The remaining parameters, excluding the shielding distance r, (that will be considered below), were subjected to a least-squares optimization with respect to the a b initio potential surface of the gas-phase dimer designated as tinter in the MCY paper. This surface, however, was modified to correct for some of its deficiencies. That is, the long-range part of the MCY surface corresponds to a dipole-dipole potential with a monomer moment of 2.2 D as compared to the experimental value of 1.855 D. Since weforce our surface to have the correct asymptotic form (using the experimental dipole moment) there are unavoidable deviations at large dimer separation between our model and the original ab initio surface (this surface contains only negative energy configurations at large separation). Furthermore, the a b initio surface is likely to overestimate the binding energy by some tenths of kcal/mol (Table 11). For these reasons, we modify the MCY surface used in the fitting procedure, shifting (20) Zeiss, G.D.; Meath, W. J . Mol. Phys. 1977, 33, 1155. (21) Clough, S. A,; Beers, Y.; Klein, G . P.; Rothman, L. S. J . Chem. Phys. 1973, 59, 2254.
moments with different r, values. The 2.2 D potentials were refined by using 49 dimer configurations (out of the 65 consisting of the ab initio surface) with an 0-0distance of less than or equal to 7 bohr. A good fit to the rest of the surface was automatically achieved owing to the fixed monomer moment. As the fit was found to be poor for the 10 configurations of the shortest 0-0 distance, a smaller weight of 0.1 was given to those configurations (see the considerations above). As seen from Table I1 the 2.2 D model with different r, values yields very similar dimer properties and root mean square error to the a b initio surface. As expected, those results are also close to the corresponding results of the MCY model.
TABLE I: Parameters for Polarizable-Atom Models" P/ 1.855 P12.2 qH 0.63760* 0.65576* a00 2516.8 358.40 boo 2.6975 2.2859 aOH 16.834 2.5806 bo" 2.5688 1.8422 QHH 15.771 2.1625 ~ H H 2.2612 1.7828 coo 32.270* 32.270: 0.53485* 0.44721* 1.429* 1.429* $3;: 1.45 1.45 r5. A I n atomic units unless otherwise noted. The following unit-conversion factors were used: 1 bohr = 0.529177A le1 = 4.80325 X esu, and 1 au = 627.503 kcal/mol. b T h e 0 - M distance in Figure I . The parameters designated by * are not treated as adjustable parameters (see text). (I
TABLE 11: Optimized Dimer Structure, Dimer Binding Energy, and Rms Deviations from the ab Initio MCY Potential Surface Dotential P/ 1.855 MCY
TIPS2
sPC
Kuwajima and Warshel
The Journal of Physical Chemistry, Vol. 94, No. I , 1990
462
rM.
2.836 2.872 2.788 2.752
A
O."der
48.7 37.6 46.5 22.2
-Au 5.60 5.87 6.20 6.61
Potentials with F~~~ = 2.2 Dd rs = 0.9 2.860 41.7 5.87 rs = 1.1 2.862 41.2 5.88 rr = 1.3 2.866 40.6 5.86 39.7 5.85 P12.2 ( r 5 = 1.45) 2.867 rs = 1.6 2.866 39.1 5.84 expt 2.98 f 0.04c 60 f IO' 5.44f 0.7f
rmsC
0.82(0.42) 0.19 (0.13)
0.21 (0.09) 0.21 (0.09) 0.21 (0.09) 0.21 (0.10) 0.21 (0.10)
" T h e angle between the 0-0 line and the HOH bisector of the hydrogen-accepting molecule. *Binding energy in kcal/mol. ' I n kcall mol. The rms error to the partial surface in which the dimer configurations of the shortest 0-0 distance (4.5bohr) are excluded is shown in parentheses. Polarizable-atom models with the monomer dipole moment of 2.2 D; the rs parameters are in A. 'Dyke, T. R.; Muenter, J. S . J . Chem. Phys. 1974,60, 2929. fCurtiss, L. A.; Frurip. D. J . ; Blander, M. J . Chem. Phys. 1979. 7 1 . 2703.
it upward in a uniform way by 0.33 kcal/mol. Finally, since we used the convenient exponential-6 potential, which limits the flexibility of our fitting procedure, we found it essential to place only a small weight on the points of the MCY surface with the shortest 0-0distance. This procedure may also be viewed as an attempt to stay close to the relevant experimental information (Table 11) with a simple atom-atom van der Waals potential. The Powell methodZZwas used in the optimization procedure which includes the 39 dimer configurations with an 0-0distance of less than or equal to 6 bohr, with the configurations of the shortest 0-0 distance (2.38 A) given the weight of 0.01. The shielding distance, rs, which turns out to be insensitive to the resulting mean-square error, was examined with several values, ranging from 0.9 to 1.6 A, and left to be tuned by the calculations of the ice phases. The resulting potential that is referred to as the P/ 1.855 is given in Table I. The performance of this model with regard to the dimer properties is considered in Table 11. As seen from this table the root mean square error between the P/ 1.855 potential and the modified ab initio surface is fairly small: 0.22 kcal/mol when the configurations of the shortest 0-0 distance are excluded (the corresponding error to the overall surface is 0.66 kcal/mol). In addition to the P/ 1.855 potential we examine and refined a potential with 2.2 D monomer dipole moment, which corresponds to the MCY moment. The potential designated as P/2.2 in Table I is roughly the best among the potentials of 2.2 D monomer (22) Press, W. H.; Flannery, B. P.; Teukolsky, S. A,; Vetterling, W . T. Numericol Recipes; Cambridge University Press: New York, 1986. (23) Frisch, M . J.; Del Bene, J. E.; Binkley, J. S.; Schaefer, H. F., 111 J . Chem. Phys. 1986. 8 4 . 2279.
Incorporation of Induced Dipoles in Lattice-Energy Minimization Although energy minimizations in crystals are mostly routine work these days," full energy minimization that includes induced dipoles has not been reported to the best of our knowledge. In this section we briefly explain basic aspects of the computation, focusing on the treatment of induced dipoles. The main features of our computational procedure are analytical evaluation of the first derivatives of the lattice energy, noniterative calculations of induced dipoles (see ref 1 for discussion of the iterative model), and the thorough use of the Ewald method for electrostatic interactions including the induced dipoles. The novel aspect of our approach is the calculation of the induced-dipole energy and its first derivatives. The induced-dipole energy per unit cell is given by n
where the summation is over the atoms in a central unit cell, pn is the induced dipole on the nth atom, and EQ,, is the Q-field (Le., field produced by permanent charges) on the nth atom. We consider here a general case in which crystals are composed of polarizable atoms with given permanent charges. The induced dipoles are calculated by solving the linear equation' (9) where pni and EQ,,i denote the component i (i = x, y , z ) of the p, and E , , vectors. The matrix (Sni,mj) is defined by
where a, is the (isotropic) polarizability of the nth atom, 6 , , is the Kronecker 6, and the dipole-dipole interaction matrix (TdPj) is given by
(rmj - rnj + l j ) / I r m - rn + 1151 (11) where r,, is the position vector of the nth atom with r, = (rm,rny,rZ) and the summation is over the lattice vectors I, I = I,a Ibb LCcwith the basic translation vectors a, b, and c. When n and m are equal in eq 11, the zero vector I = (O,O,O) is excluded from are evaluated by the standard the summation. The quantities TniJnj Ewald (or Ewald-Kornfeld) method. The first derivatives of the induced-dipoleenergies are evaluated as follows. We start by rewriting eq 9 as
+
,IO
+
X.Y.2
The derivative of Uind with respect to a variable X is obtained by using eq 12 along with the matrix relationship 6S1/6X = S-'(dS/dX)S-'; the resulting expression is
auind/ax= n
r
m
J
I n the case that .Y = rnrrmanipulation of eq 13 leads to
Water-Water Interaction Potentials
The Journal of Physical Chemistry, Vol. 94, No. 1, 1990 463 xy,z
auind/arni
= -qnEP,ni -
pnjFij('n)
J
(14)
where EP,niis the component i of the P field (Le., field produced by induced dipoles) on the nth atom and Fij(rn) denotes the field gradient at the position r, that includes both the permanent-charge and induced-dipole contributions. In the case that X i s a component of the basic translation vectors a, b, and c, eq 13 is rewritten as
TABLE 111: Dependence of the Results of Ices I1 and IX on the Shieldine Parameter ( r e ) 1.1
1.3
1.45
1.6
expt
2.584 1.36 18.37
2.700 0.78 16.63
Ice I1 2.765 0.51 15.30
2.799 0.65 14.51
2.821 0.86 13.89
2.799
2.560 ~ 0 0 0 ~7.87 -AUc 18.26
2.676 6.28 16.59
Ice IX 2.741 4.84 15.32
2.773 3.80 14.58
2.794 2.96 14.00
2.770
0.9
rooa umb
-AUc rooa
where the derivative of Ep,niis taken with the induced dipoles fixed. The derivatives of Q or P fields appearing in the right side of eq 14 and 15 are evaluated by the Ewald method in a similar manner to the evaluations of similar quantities in the usual minimization procedure. Additional complexity arises from the shielding factor introduced by eq 4. After evaluating the quantities with the normal Coulomb interactions by using the Ewald method, correction terms have to be incorporated. The correction of the Q field is AEQ,ni = -
[ I -fs(r)]qn~(~ni - rmi)/r3
(16)
m#n
where r = Irn - rml and the summation is over the atoms of the whole lattice. Spherical truncation is employed in the computation of eq 16. The corresponding correction to the field gradient is
where r = Ir, - rml andf,'(r) denotes the derivative off,(r). A similar correction term has to be applied to the derivative dEQ*/aX in eq 15. Furthermore, the term q,,EP,niin eq 14, which represents the reaction to the forces which act on induced dipoles from the permanent charges, has to be modified. The correction term for EP,niin eq 14 is given by
where r = (r, - rml. Although EP,niin eq 14 is modified by eq 18, the P field contribution to the field gradients Fijin eq 14 and the term aEp,,,/dX in eq 15 remain unchanged because they represent induced-dipole-induced-dipole interactions that are assumed here not to be subjected to the shielding modification. The implementation of our computation procedure includes additional details which are given below. (i) The van der Waals interactions are calculated by the usual spherical truncation combined with a simple continuum correction AUvw= ( 2 r / u ) E E n a n a l s m U a a t ( r )dr 2r a a'
Rc
(19)
where n, denotes the number of the ath kind of atoms in the unit cell, uad(r)is the interaction potential between the ath and a'th atoms, R, is the cutoff radius, u is the unit-cell volume, and the summations are over the atom kinds. The truncation correction AU,, depends only on the unit-cell parameters through 1/ u . (ii) The energy minimization is performed by the conjugate gradient method.22 (iii) Spatial symmetries other than the translational one are not incorporated in the computation; as a result, the minimization procedure does not preserve the symmetry and calculated structures retain traces of broken symmetry. Results and Discussion Ices 11, IX, and VI11 are the known proton-ordered phases of ice. Crystal structures determined by neutron diffraction are
14.1
14.0
"Average nearest-neighbor 0-0 distance, in A. Root mean square error of calculated nearest-neighbor 0-0-0 angles relative to the experimental ones, in degrees. Lattice energy, in kcal/mol.
available for ices I1 and IX.24925 For ice VI11 detailed neutron or X-ray diffraction studies have not been published but key information from which the full structure can be deduced are cited in the literature.26 A small amount of disorder (about 4%) in H atoms was reported in ice IX but it is neglected here. Let us first examine the effects of the shielding factor f , ( r ) introduced in eq 4 and 5. Table 111 shows the comparison of the calculated and experimental structures and lattice energies for different r, parameters in ices I1 and IX. In striking constrast to dimer properties (Table 11), the calculated results of ice strongly depend on the rs parameter. With increasing rs value the nearest-neighbor 0-0 distance becomes longer and the lattice energy decreases in magnitude. Furthermore, the ice I1 crystal at r, = 0.8 A appears to lose meaningful equilibrium structures, with some of the water-water distances apparently collapsing to zero. Therefore, the use of the shielding factor is unavoidable within the framework of the present model potentials. We choose r, = 1.45 A as a standard value by using the average 0-0 distance as a criterion. This value of r, yields thefs factors of 0.786 and 0.973 for a typical hydrogen-bonded 0-H distance of 1.8 A and 0-0 distance of 2.76 A, respectively. Thus, the primary effect of the shielding is to reduce the fields produced by the partial charges of the nearby hydrogen-bonded H atoms. The shielding factors are fairly close to unity for the 0-0 distances, and this justifies our omission of the shielding factor in the evaluation of the induced-dipole-induced-dipole interactions. Tables IV, V, and VI summarize results for ices 11, IX, and VIII, respectively (the quantities examined in these tables are the same as those of the M R paper16). Results of the MCY and two other pair-interaction models, which are computed here, are also shown for comparison. The MCY model reported here is the version designated tinterin the MCY paper. (In the MR paper another version of the MCY model was examined, along with the ST2 and other models.) Experimental lattice energies are estimated from the lattice energy of ice Ih and the heats of transformation among the ice polymorph^.^^.*^ As for computational details, the truncation radii used in the energy minimizations were 10-12 A for the direct-lattice interactions; 50-100 terms of the reciprocal-latticevectors were retained for the calculations of the electrostatic interactions (including induced dipoles). An estimate of the errors involved in the results of Tables IV-VI is obtained by the mutual discrepancy in quantities which should give a single value by symmetry. This estimate indicates that most of the results are correct to the last digit but some of them might be in error by about *2 at the last digit. We examine the performance of the MCY model in Tables IV-VI. The MCY model consistently gives too long nearest(24) Kamb, B.; Hamilton, W. C.; LaPlaca, S. J.; Prakash, A. J . Chem. Phys. 1971, 55, 1934. (25) LaPlaca, S. J.; Hamilton, W. C.; Kamb, B.; Prakash, A. J . Chem. Phys. 1973, 58, 567. (26) Whalley, E. In The Hydrogen Bond; Schuster, P . , Zundel, G., Sandorfy, C., Eds.; North-Holland: Amsterdam, 1976; Vol. 111, Chapter 29. (27) Bertie, J. E.; Calvert, L. D.; Whalley, E. Can. J . Chem. 1964, 42, 1373.
The Journal of Physical Chemistry, Vol. 94, No. I , 1990
464
TABLE IV: Calculated and Observed Structures and Energies for Ice II quantity’ P / l . 8 5 5 P/2.2 MCY TIPS2 S P C expt rWqb8,
--
1-1
I II I1
II II I
-
1-11
2.790 2.785 2.774 2.793 3.175
2.777 2.805 2.775 2.838 3.253
2.908 2.928 2.863 2.959 3.372
2.772 2.789 2.770 2.8 16 3.214
2.739 2.759 2.735 2.780 3.168
2.803 2.768 2.780 2.844 3.240
10.2 9.8 4.1 1.8
9.4 10.2 4.6 2.6
12.1 9.9 6.9 3.1
8.5 7.9 3.6 4.1
8.3 9.1 2.5 4.2
9.1 8.9 1.4 7.8
114.4 82.9 115.8 128.2 87.7 130.1 1 1 8.9 121.6 81.4 100.8 107.8 126.8 0.99 7.72 113.0 1.203 13.96
115.5 85.3 115.6 125.9 87.2 130.1 118.6 122.2 81.3 98.5 109.2 128.3 0.65 7.76 1 1 3.0 1.177 14.51
112.8 85.9 116.4 128.1 86.7 129.1
116.6 85.7 114.5 124.9 89.9 127.6 119.6 126.6 79.0 98.5
OH,‘ deg
--
1-1
I II
I1
II
-
11 1
0 ~ 0deg , ~
1-1-11
I-1-11 1-1-1
I-1-11
-- --- --- --
I-1-11
II I 11 II II II II II I II II I II II I II I1 1 I-11-1 ,Jowle a=b=d a=p=yj Pg
-AUh
118.8 124.0 82.2 98.4 108.4 125.9 1.67
8.10 113.1 1.044 11.35
120.8 80.5 111.7 123.4 93.4 128.7 119.9 124.6 76.1 103.2 105.0 102.3 128.0 128.7 3.85 2.28 7.74 7.59 113.1 112.9 1.199 1.255 13.15 13.73
116.0 83.8 115.1 126.2 87.9 130.4 118.7 121.9
81.0 99.6 108.9 127.7 7.78 113.1 1.180 14.1
“The arrows designate the direction of hydrogen donation. The Roman numerals designate molecules according to the notation of ref 24. See also ref 16. *Nearest-neighbor0-0 distance except for 1-11 that depicts a non-nearest-neighbor pair. ‘Angle between the hydroxyl 0 - H line and the corresponding hydrogen-bond 0-0 line. dNearestneighbor 0-0-0 angle. Rms error of the calculated nearest-neighbor O-O-O angles to the experimental ones. [Unit-cell parameter, length in 8, and angle in degrees. gDensity, in g/cm3. *Lattice energy, in kcal /mol. TABLE V: Calculated and Observed Structures and Energies for Ice I X auantitv‘ P/1.855 Pi2.2 MCY TIPS2 SPC expt
--
roo. A I 2’ I 1 ”’
2-1
nonbonded
2.770 2.761 2.808 3.31 3.64 3.56
2.772 2.752 2.795 3.81 3.84 3.91
2.873 2.845 2.900 4.13 3.95 4.26
2.769 2.757 2.762 3.47 3.65 3.76
2.742 2.117 2.751 3.41 3.64 3.59
2.163 2.750 2.797 3.45 3.65 3.73
8.6 9.3 5.2
9.2 7.2 8.1
9.6 9.7 5.7
7.3 6.4 3.7
7.6 6.3 7.3
9.7 8.3 3.5
92.4 98.8 144.1 94.1 110.5 102.6 129.5 91.5 108.2
93.8 103.4 140.5 95.1 113.0 100.1 134.6 94.2 103.9 92.4 3.80 6.86 6.88 90.0 1.108 14.58
101.6 100.9 138.2 89.2 117.0 98.1 131.0 102.2 102.2 92.6 3.35 6.90 7.61 90.0 0.983 11.75
97.7 99.4 143.8 91.7 112.7 98.4 127.8 101.7 105.5 91.5 0.54 6.71 6.85 90.0 1.163 13.30
92.2 99.4 142.3 94.5
96.8 99.2 143.7 91.8 112.7 99.2* I281
88.0
2.61 6.82 6.54 90.0 1.178 13.85
111.0
103.6 131.7 96.5 107.4 87.3 3.17 6.74 6.51 ~~
~
100.8
106.1 90.9 6.73 6.83
900
900
1212 1361
1160 140
’The numbers designate molecules according to the notation of ref 25. See Table I V for the meanings of the other symbols. *This angle is quoted differently in both ref 16 and 2 5 but the value given here is derived from the primary data of ref 25.
neighbor 0-0 distances by about 0.1 A and too low density by 10-15%. The lattice energy is too high by 2-4 kcal/mol. Nearest-neighbor O-O-O angles, on the other hand, are fairly
Kuwajima and Warshel well reproduced. The polarizable-atom models, on the other hand, are free from the above-mentioned deficiencies of the MCY model while mostly sustaining its successful results. Thus, the polarizable-atom models, which were built from almost the same a b initio surface as the MCY model, represent substantial improvement over the MCY model. When we compare the polarizable-atom models with liquidwater-based pair-interaction models (Le., the TIPS2 and S P C models), no overall improvement is observed in ices I1 and IX. While the polarizable-atom models, particularly the P/ 1.855 model, appear to give better results than the S P C model,15 the TIPS2 is comparable or even slightly better for ices 11 and IX. These observations mean that the transferability between liquid water and ice is generally good in pair-interaction models. The situation is somewhat different in ice VIII, where the P/ 1.855 model is clearly better than the TIPS2 model. Ice VI11 is abnormal in that it involves non-hydrogen-bond water-water contacts.’6,26 (Note that the shortest 0-0contact in Table VI is not a hydrogen bond.) Therefore, ice VI11 provides a test of model potentials for the unusual non-hydrogen-bond contacts. The relatively successful results of the P/1.855 model in ice VI11 may indicate that this model carries beter transferability than the liquid-water-based models in which interactions infrequent in liquid water are not built in. However, as will be argued in the Concluding Remarks section, the agreement in the ice-phase case (which is not much different than water) is not the main issue of the present paper. Table VI1 lists induced dipole moments and energies for the ices and the dimer at their calculated equilibrium structures. The induced dipoles in the ices are almost parallel to the permanent dipole of the water molecule. In ices I1 and IX the sum of the permanent and induced dipoles amounts to about 2.8 and 3.3 D according to the P/ 1.855 and P/2.2 models, respectively. These numbers are much larger than the monomer dipole moments of 2.19 (MCY), 2.24 (TIPSZ), and 2.27 (SPC) D in pair-interaction models. The results here are in accord with the notion that the monomer dipole moment in ice and liquid water is substantially larger than the gas-phase value due to the polarization of water molecules.29 Since dielectric constants are sensitive to the monomer dipole moment, it is reasonable to expect that the polarizable-atom model will be of some importance in simulations of dielectric properties. The MCY and TIP4P (a variant of TIPSZ) models in fact have been reported to give a static dielectric constant too low for liquid water,30though good agreement with experiment has been reported with a modified S P C model in which water molecules are flexible3’ (in nonrigid water models the molecular flexibility simulates a part of the effect of the electronic polarizability but of course with incorrect frequency dependence). An interesting point emerges from a closer examination of the energy results, which can be analyzed in terms of the hydrogenbond energy (the hydrogen-bond energy of ice is identified here with the corresponding lattice energy divided by 2 ) . That is, from the energy results in Tables I1 and IV-VI, the difference between the hydrogen-bond energy of ice and that of the gas-phase dimer (i.e., the binding energy) is experimentally about 1.6 kcal/mol, which is the average of ices I1 and IX, while model potentials yield -0.10 (MCY), 0.41 (TIPSZ), 0.23 (SPC), 1.35 (P/1.855), and 1.42 ( P / 2 . 2 ) kcal/mol. Though the same quantity is of less physical significance for ice VI11 (because of the existence of the non-hydrogen-bond contacts), it is formally given as 1.3 kcal/mol experimentally and -1 .OO (MCY), -0.54 (TIPS2), -0.78 (SPC), 0.61 (P/1.855), and -0.10 (P/2.2) kcal/mol by model potentials. Thus, the hydrogen-bond energy is generally enhanced in ice as compared to the dimer but the pair-interaction models badly underestimate this enhancement (the reasonable energy provided by such models is obtained at the expense of having too large dimer (28) Jorgensen, W. L.; Chandrasekhar, J.; Madura. J. D.; Impey, R. W.; Klein M. L. J . Chem. Phys. 1983, 7 9 , 926. (29) Eisenberg, D.; Kauzmann, W. The Srruclure and Properlies of Wafer;Oxford University Press: New York. 1969. (30) Neumann, M. J . Chem. Phys. 1985, 82, 5663; 1986, 85. 1567. (31) Anderson. J.; Ullo, J. J.; Yip, S. J . Chem. Phys. 1987, 8 7 . 1726.
The Journal of Physical Chemistry, Vol. 94, No. 1, 1990 465
Water-Water Interaction Potentials TABLE VI: Calculated and Observed Structures and Energies for Ice VI11
auantitv"
P I 1.855
Pl2.2
MCY
TIPS2
SPC
expt
2.938 2.852 3.028 3.340 3.346 1.3 110.6 107.2
3.033 2.968 3.100 3.625 3.383 0.2 112.2 104.1
3.130 2.920 3.375 3.129 3.676 3.4 108.6 111.3
2.975 2.966 2.983 3.624 3.329
2.969 2.799 3.159 3.187 3.408
1 12.0 104.6
2.943 2.952 2.934 3.717 3.235 3.7 113.3 102.0
70.6 70.4 177.6 75.1 68.2 177.8 126.4 53.6
68.6 74.3 178.4 77.5 67.0 178.4 128.0 52.0
74.7 62.2 173.4 74.4 68.7 174.3 124.4 55.6
68.2 75.2 179.8 75.6 67.9 179.8 127.7 52.3
66.6 78.2 179.8 77.8 66.8 179.8 129.0 51.0
72.4 67.0 174.9 76.6 67.5 175.5 126.1 53.9
111.9 102.7 110.5 56.0 128.7 2.55 4.73 6.98 90.0 0.042 1.535 12.41
107.3
124.4 99.9 107.5 62.2 130.1 2.76 5.17 7.07 90.0 0.115 1.268 9.75
105.0 104.2 112.0 52.5 127.9 5.18 4.71 7.28 90.0 0.004 1.484 11.33
101.6 102.5 113.3 50.8 128.8 6.50 4.57 7.41 90.0 0.004 1.544 11.67
118.1 98.9 109.5 59.0 130.6
roo, 8,
H bonded
type I b type I l b type I l l b type IV OH, deg
deg
Om:
e)-$
deg type 1 type I type I type I1 type I I type 11 type 111 type .. 111 1-1 11-11 1-1 I 1-111 11-111
~ooo
a=b
c Ly=p=y K/
P
-AU
101.0 112.2 53.7 129.5 4.26 4.78 7.46 90.0 0.029 1.402 11.51
0.0
110.3 107.9
4.80 6.988 90.0 0.0888 1.486 13.5
"The designation of bond and angle types follows that of ref 16. See Table IV for the meanings of the other symbols. bNon-hydrogen-bond0-0 contact. '0-0-0 angle for hydrogen-bonded neighbors. dO-O-O angle involving one non-hydrogen bond. cO-O-Oangle involving no hydrogen bond. /Structural parameter which describes the relative position of the two ice I, lattices composing the ice VI11 lattice (see ref 16). ZFrom ref 26. A slightly different value (0.087) is adopted in ref 16. TABLE VII: Induced-Dipole Moments and Induced-Dipole Energy" Ice I I b P/l.855 pj2.2
PI
PI1
0.975 (3.8) 1.106 (3.9)
0.980 (1.7) 1.114 (1.3)
-Uind 3.29 4.25
PI
P2
-uind
0.930 (5.7) 1.088 (3.7)
0.940 (0.1) 1.071 (0.0)
3.1 1 4.16
Ice I X b ~~
Pl1.855 Pl2.2
Ice VI11 0.679 (0.0) 0.695 (0.0)
PJl.855 Pl2.2
2.28 2.35
DimerC Pll.855 Pl2.2
PI
112
-uind
0.262 (34.7) 0.286 (26.8)
0.221 (48.0) 0.249 (46.3)
0.52 0.64
" Dipole moment in debyes and energy in kilocalories per mole. The angle (in degrees) between the directions of the induced dipole and the bisector of the HOH angle is given in parentheses. bThe Roman- or Arabic-numeral subscripts of the dipole moments correspond to the designation of the water molecules in Table IV or V. c p , and cc2 designate the induced dipoles on the hydrogen-accepting and hydrogendonating water molecules, respectively.
binding energy). Furthermore, the major part of the energy enhancement comes from the induced-dipole energy in the polarizable-atom models. From Table VI1 the induced-dipole contribution to the hydrogen-bond-energy enhancement for ices I1 and IX on average is 1.08 and 1.46 kcal/mol by the P/1.855 and P/2.2 models, respectively. Thus, it is suggested that the induced dipoles may be the major factor in the variation of the hydrogen-bond energy in different phases. The induced-dipoledriven variation is estimated to be about 1 kcal/mol between the
dimer and ices I1 and IX. This is about 15% of the total hydrogen-bond energy and is not small enough to be ignored. Concluding Remarks
The present work examines polarizable-atom models, giving special attention to the transferability of gas-phase dimer potential to condensed phases. The P/1.885 model developed here retains the gas-phase monomer dipole moment and a good fit to the properties of the gas-phase dimer, while giving encouraging results for the properties of ices. This demonstrates that the poor transferability suffered by the pair-interaction models is mostly removed by using the polarizable-atom model. The incorporation of induced dipoles appears to provide a consistent and effective way of extrapolating dimer potentials to condensed phases, but it should be kept in mind that an adjustable parameter, the shielding constant of rs in eq 5, has played an essential part in the refinement of our polarizable-atoms models. The use of condensed-phase results in tuning parameters for studies of condensed phases is clearly an appropriate and, at present, highly recommended procedure. However, it should be challenging and instructive to find a way to evaluate the proper value of rs and the effects is represents without using condensed-phase information. Although we have emphasized the advantage of the polarizable-atom model, it should be noted that the liquid-water-based pair-interaction potentials have proved to work fairly well for the polymorphs of ice. In fact, it is possible that the inclusion of three-body charge-transfer interactions and consideration of the quantum mechanical effects associated with the nuclear motions (in addition to the induced dipoles effect) is needed for making a polarizable-atom model for pure liquid water or ice, superior to a well-parametrized pair-interaction potential. However, the inclusion of induced dipoles is an essential element in consistent calculations of charge-transfer reaction^'.^^.^^ and excited-state (32) Hwang, J.-K.; Warshel, A. J . Am. Chem. SOC.1987, 109, 715.
J . Phys. Chem. 1990, 94, 466-471
466
processes’s2 and should be included in treatment of systems containing ions. These effects might also be important in consistent treatment of the bulk dielectric constant and are essential for treating its frequency dependence. One may argue that in view of the success of pair-interaction models in describing bulk properties, there is no need for polarizable-atom models. However, it is important to mention that the present model provides better transferability from gas phase to condensed phases than pair-interaction models. That is, the P/ 1.855 does not use the monomer dipole and atomic polarizabilities as adjustable parameters and therefore reproduces the correct long-range behavior for the gas-phase dimer. This is not the case for pair-interaction models that must adjust the monomer dipole moment in order to get a good fit for condensed phase properties. The implication of this fact is apparent if one considers the difference between the hydrogen-bond energy in ice and the gas-phase dimer (see discussion in the previous section). Pairinteraction models can, of course, be adjusted to give good energies in condensed phases but this is usually done at the expense of overestimating the gas-phase energy. The reason is simple, about 15% of the total hydrogen-bond energy comes from induced dipoles and this real effect cannot be reproduced consistently without (33) Hwang, J.-K.; King, G.; Creighton, S.; Warshel, A . J . Am. Chem.
including induced dipoles in the given model. As to structural aspects of the bulk, it is important to point out that the use of a smaller permanent dipole moment in polarizable-atom models leads in general to more “diffuse” forces, making it harder to fit for correct short-range orientational features.’ Thus it is no surprise that previous liquid-water simulations employing polarizable-atom models, which are not very carefully paramet r i ~ e d ,have ~ . ~ been less successful than those of the best pairinteraction models. Obviously the present P/1.855 model is not considered by us as our “best” potential for simulations of solution processes. We intend to continue refining polarizable-atom potentials by adding to the current fitting procedure gas-phase experiments, improved a b initio surfaces, and the experimental information about liquid water. Such studies as well as investigation of polarizable-atom potential in studies of hydrated ionic crystals are planned for the near future. In addition we intend to continue our studies of charge-transfer reactions, where the inclusion of induced dipoles is quite important, replacing our previous polarizable-atom water models3~7~32~33 by more refined potentials.
Acknowledgment. We thank Dr. F. Sussman for valuable discussions. This work was supported by Grant No. GM-24492, National Institutes of Health. Registry No. H,O. 7732-18-5.
Soc. 1988, 110. 5297.
Polymer Swelling. 8. Sorption of 1- Iodoalkanes by Poly(styrene-co-divinylbenzene) L. A. Errede 3M Corporate Research Laboratories, 3M Center, Bldg. 201 -2N-22, St. Paul, Minnesota 551 33 (Received: December 13, 1988; In Final Form: June 19, 1989)
The results obtained in this study of poly(styrene-co-divinylbenzene) swelling in liquids that comprise the homologous series I(CH,),H confirm that the number, a,, of adsorbed molecules per accessible phenyl group of polymer at liquid saturation is a function of three parameters: cyo, the intrinsic dynamic packing efficiency when steric hindrance owing to the polymethylene chain is not a factor; A, the incremental increase in steric hindrance to adsorption per added methylene group; and B, the entropic or enhanced steric hindrance factor owing to self-association of the liquid when n > 6. Thus, the logarithm of a, decreases linearly with increasing n, where n < 8, and also with the corresponding difference (6p, - 6,)*, where 6,, is the solubility parameter of the polymer and 6, is the calculated solubility parameter of the test liquid. The extension of the straight line obtained thereby was used to establish the effective solubility parameters (S’,) for those liquids with n > 6. The difference ( 6 , - a’,) varies with the magnitude of self-association, which increases with (n - 6) for liquids with n > 6.
Introduction I reported’s2 that the volume, S, of liquid sorbed per gram of poly(styrene-co-divinylbenzene) [hereafter referred to as poly(Sty-co-DVB) or (Sty),,(DVB),] at thermodynamic equilibrium with excess liquid is a function of the average number, X [ A = 1 / x at small x ] , of backbone carbon atoms between cross-link junctions as given by
s = C(X1/3 - X01/3) (1) where C (in milliliters of liquid per gram of polymer) is the relative swelling power of the test liquid and 1/X, is the critical cross-link density of the polymer above which S (in milliliters of liquid per gram of polymer) is zero, Le., immeasurably small. Since as,the number of adsorbed molecules per accessible phenyl group in the polymer, varies directly with C and inversely with the molar volume of the sorbed liquid, eq 1 can be recast to give aG
= asA
where A = is the macrostructural factor proportional to the uloosenessnof the cross-linked polymer network such that the product of A and a, gives the molecular sorption capacity, aG,which is the total number of sorbed molecules per phenyl group in the liquid-saturated gel supported by polymer with cross-link density 1/X. The presence of two types of sorbed molecules in such gels was established kinetically by monitoring evaporation from the enmeshed poly(Sty-co-DVB) particles saturated with volatile liq u i d ~ .The ~ ~nonadsorbed ~ molecules in the gel escape first, and the process follows zero-order kinetics, the rate constant of which indicates association with molecules of its own kind. The adsorbed molecules escape thereafter, and the process follows first-order kinetics, the rate constant of which indicates association with the adsorbent polymer. A typical pattern of such desorptions is shown in Figure 1, which records residual a, as a function of time for evaporation of
(2)
( I ) Errede, L. A . Polymer Swelling. 2 . J . Appl. Polym. Sci. 1986, 19, 1749. (2) Errede, L. A . Polymer Swelling. 7. J . Phys. Chem. 1989, 93. 2668.
0022-3654/90/2094-0466$02.50/0
(3) Errede, L. A.; Kueker, M. J.; Tiers, G . V. D.; Van Bogart, J . W . C. Polymer Drying. 1. J . Polym. Sci., Part A-1 1988, 26, 3375. (4) Errede, L. A.; Van Bogart, J. W. C. Polymer Drying. 2. J . Polym. Sci., Parr A-1 1989, 27, 2015.
0 1990 American Chemical Society
