The Journal of
Physical Chemistry
~~~
Registered in US.Patent Office
0 Copyright, 1981, by the American Chemical Society
VOLUME 85, NUMBER 1
JANUARY 8,1981
LETTERS A Revised Potential Function for the Water Dimer in the EPEN/2 Form Francis T. Marchese, P. K. Mehrotra, and David L. Beverldge” Chemistry Department, Hunter College of the City University of New York, New York, New York 10021 (Received: August 27, 1960; In Final Form: November 6, 1960)
A revised potential function for water-water interactions is presented in the EPEN form based on the CI potential of Matsuoka, Clementi, and Yoshimine. Monte Carlo calculations performed with this function produce a radial distribution function that is in good agreementwith experiment and comparable to distribution functions derived from other potentials.
Introduction Theoretical studies of molecular assemblies require the rapid evaluation of large numbers of intermolecular energies from a potential function representation of the intermolecular interaction. There have been diverse means proposed for constructing these potentials. Considerable work on model potentials,l on empirical potentials functions based on experimental data,2-12and on theoretical (1)J. 0.Hirshfelder, C. F. Curtiss, and R. B. Bud, “Molecular Theory of Gases and Liquids”, Wiley, New York, 1954. (2)F.H. Stillinger and A. Rahman, J. Chem. Phys., 60,1545 (1974). (3)M. Weissman and L. Blum, Trans. Faraday SOC.,64,2605(1965). (4)P. H. Smit, J. L. Derissen, and F. B. Van Duijneveldt, Mol. Phys., 37,521 (1979). (5)F.A. Momany in “Environment Effects on Molecular Structure and Properties”, Reidel, Dordrecht, Holland, 1976,pp 437-58. (6)D. E.Williams, J. Chem. Phys., 45,3770 (1966). (7)D.E.Williams, J. Chem. Phys., 47,4680 (1967). (8)L. Hsu and D. E. Williams, Znorg. Chem., 18,79 (1979). (9)S.Lifson and A. Warshel, J. Chem. Phys., 49,5116 (1968). (10)A.T.Hagler, E. Huler, and S. Lifson, J. Am. Chem. SOC.,96,5319 (1974). ( I 1) A. T.Hagler and S. Lifson, J. Am. Chem. SOC.,96,5327(1974). 0022-3654/81/2085-0001$01.00/0
potentials representative of quantum mechanical calculat i o n ~ ~ “has ’ ~ been reported. Recently, Snir, Nemenoff, and Scheraga (SNS)2“-24have proposed empirical potential (12)S.Lifson, A. T. Hagler, and P. Dauber, J. Am. Chem. SOC.,101, 5111 (1979). (13)E.Clementi, “Determination of Liquid Water Structure, Coordination Numbers for Ions, and Solvation for Biological Molecules”, Springer-Verlag, New York, 1976. (14)S.Swaminathan, R. J. Whitehead, E. Guth, and D. L. Beveridge, J. Am. Chem. SOC.,99,7817(1977). (15)S. W. Harrison, S. Swaminathan, and D. L. Beveridge, Znt. J. Quant. Chem., XIV,319 (1978). (16)W. L. Jorgensen and M. E. Cournoyer, J. Am. Chem. SOC.,100, 4942 (1978). (17)W. L. Jorgensen, J. Am. Chem. Soc., 101,2011 (1979). (18)W. L.Jorgensen, J . Chem. Phys., 71,5034 (1979). (19)W. L. Jorgensen and M. Ibrahim, J . Am. Chem. SOC.,102,3309 (1980). (20)J. Snir, R. A. Nemenoff, and H. A. Scheraga, J.Phys. Chem., 82, 2497 (1978). (21)R. A.Nemenoff, J. Snir, and H. A. Scheraga, J. Phys. Chem., 82, 2504 (1978). (22)R. A.Nemenoff. J. Snir. and H. A. Scheraea. J . Phvs. Chem... 82.. 2513 (1978). (23)R. A. Nemenoff, J. Snir, and H. A. Scheraga, J. Phys. Chem., 82, 2521 (1978). I
0 1981 American Chemical Society
.
2
The Journal of Physical Chemistry, Vol. 85,No. 1, 1981
function methodology whereby parameters for molecular fragments could be readily transferred from one system to another such that intra- and intermolecular potentials for large molecules can be constructed from parameters determined for prototype systems. The resulting EPENI2 (empirical potential based on the interactions of electrons and nuclei) achieved quite notable success for predicting dipole moments, dihedral angles, rotational barriers, and crystal structures for diverse systems. Still problematic, however, was the description of the water-water interaction central to the study of aqueous solutions of importance in chemistry and biochemistry. The original EPEN parameterization for water was based entirely on experimentally determined quantities. The resulting function incorrectly characterized some repulsive regions of the dimer surface. Subsequent reparameterization was founded on the small quantum mechanical data base of Matsuoka, Clementi, and Yoshimine (MCY).25 Monte Carlo simulations on liquid water based on this function led to an overestimation in the average number of nearest neighbors and the placement of the second maximum of the radial distribution function. These results indicated that the EPEN/2 potential for the water dimer did not possess enough tetrahedral character. Either the EPEN/2 model for water is not flexible enough to provide a correct account of the entire dimer surface, or a larger data base is required to develop a potential function which adequately represents these features. Our recent experiences1*J5with potential functions derived from quantum mechanical calculations indicate that about 200 dimer configurations, well distributed over the dimer surface, are sufficient to define a potential function that usually is satisfactoy for liquid-state computer simulation. We have thus proceeded to develop a potentid function for the water dimer in the EPEN/2 form, representative of nonempirical quantum mechanical calculations. This function, along with the preexisting EPENI2 information, makes a larger number of aqueous solvation problems accessible within only a slightly enlarged EPEN/2 philosophy. We report herein the details and preliminary characterization studies of this function.
Methodology Within the EPEN/2 modelm all electrons and nuclei are treated as point charges. Therefore, water consists of four sets of negative charges, each with values of -2 au, corresponding to four electron pairs, and three positive charges of +6, +1, and +1 au corresponding to oxygen and hydrogen valence charges, respectively. The interaction of these charges may be divided into a coulombic term, Emd, represented as sums over pairwise interactions between all charges QiQj
Ecoul= 332.0719Xij
rij
kcal/mol
where qi and q, are charges on molecules i and j , and r is the distance between qi and qj in angstroms; and an nonbonded interaction, ENB,given as sums over pairwise interactions between electrons only, where ENB= Cq.1%.(A..e-Bgf4J 51 - Cijrij") kcal/mol (2) ij
Sets of A, B, and C parameters are assigned for each type of electron encountered (e.g., water possesses two sets of ~~
~
(24)R.A. Nemenoff, J. Snir, and H. A. Scheraga, J.Phys. Chem., 82, 2527 (1978). (25)0.Matsuoka, E.Clementi, and M. Yoshimine, J. Chem. Phys., 64,1351 (1976).
Letters
+1
H
Figure 1. EPEN/P model for the water molecule.
TABLE I: Revised EPEN/2 Parameters for Water Potential Functiona type
.
.~
lone pair, e-
bonding pair, e-
A B
9379.63 3195.14 5.01708 4.46344 C 0.001742 10.5947 a Oxygen-lone pair distance = 0.2724 A ; oxygen-bonding pair distance = 0.5505 A ; lone pair angle = 107.94".
parameters corresponding to bonding and lone pair electrons). For the interaction between electrons these parameters obey the following combination rules as given by SNS: A, = (A..A..)1/2 [I 11
Bjj = (Bii + Bjj)/2 c.. = (C..C..)1/2 11 81 11 The angle between lone pair electrons and the distance of electronic centers along bonding or pseudobonding axes are also treated as parameters. This results in the EPENI2 model for water depicted in Figure 1.
Calculations The revised EPENI2 water potential was generated from a data base consisting of 295 water dimer configurations. Configurational energies were calculated from the MCY-CI potential for 229 geometries used in the construction of an earlier quantum mechanical potential function by Popkie, Kistenmacher, and Clementi.26 Added to this were the 66 dimer configurations and energies originally used in the construction of the MCY-CI potential.25 The latter configurations are predominantly attractive in nature and the attractive region of the hypersurface should be improved by this addition. Fitting the sum of eq 1and 2 was performed with the GAUSHAUS nonlinear least-squares program which combines the Gauss (Taylor series expansion) method with the method of steepest decent. The resultant potential function was used in a Monte Carlo (MC) simulation in the canonical (T,V,N) ensemble consisting of 100 particles with simple cubic boundary conditions at 25 "C. An equilibrated configuration obtained from a previously reported extended MC calculation with the MCY-CI potential served as an initial guess for the simulation. Complete details of the method may be found in ref 27. Configurational averages were formed (26) H. Popkie, H. Kistenmacher, and E. Clementi, J. Chem. Phys., 1325 (1973).
The Journal of Physical Chemistry, Vol. 85, No. 1, 1981 3
Letters
TABLE 11: Radial Distribution Functions
exptl revised EPENI2 MCY-CI ST2 STO-3G
/ - s.m Figure 2. Isoenergy contour surface for the water dimer computed from the revised EPENIP potential. The minimum on this surface corresponds to an interaction energy of -6.5 kcal/mol. Successive contours represent energy increments of 0 3 kcal/mol. The dimer configuration which corresponds to the minimum is also depicted. 3.2
r
n
2.81
04 00
-
I I
,
I
1
1
,
I
I
I
I
(27) S. Swaminathan and D. L Beveridge,J. Am. Chem. SOC.,99,8392 (1977). (28) M. Mezei, S. Swaminathan,and D. L. Beveridge, J.Chem. Phys., 71, 3366 (1979). (29) W. Jorgensen, J. Am. Chem. SOC.,101, 2016 (1979).
R g(R) R g(R) R g(R) R g(R) R g(R)
R(max 1 )
R(min 1)
R(max 2 )
R(min 2 )
2.85 2.33 2.75 3.05
3.50 0.85
4.65 1.13
5.70 0.86
3.45 0.80 3.45 0.96 3.40 0.62 3.30 0.63
4.60 1.12
5.70 0.90
4.25 1.10
5.60 0.90
4.60 1.24
5.70 0.90
4.47 1.19
5.53 0.88
2.85 2.62 2.80 3.36 2.70 3.40
TABLE I11 : Coordination Number Distribution 1
2
3
0.0
0.02
0.15
0.46 0.28 0.07 0.01
MCY-CI 0.0 0.02
0.16
0.45 0.29 0.08 0.01
revised EPEN/2
4
5
6
7
and -5.84 kcal/mol. (HZO), Computer Simulation. The oxygen-oxygen radial pair distribution function, g(R), from a Monte Carlo simulation, based on the revised EPEN/2 potential, is given in Figure 3. The g(R) computed from a previously reported simulation based on the MCY-CI potential is included for comparison. The overall agreement with experiment is good and the quality of results are comparable with distribution functions produced via other potentials (Table 11). A noticeable difference is found in the positions and heights of the first maxima (Table 11, Figure 2), even though the EPEN/2 function has been parameterized against the MCY-CI potential. However, the distribution of coordination numbers in Table I11 reveals only slight changes in the number of nearest neighbors signifying a retention of the essential tetrahedral characteristics of the original potential. Thus, the shift in position and the increase in the height of the first peak does not signal a change in the computed water structure. The second peak in the g(R) is displaced toward larger water-water separations relative to the MCY-CI results, but as seen from Table 11, the position of the maximum compares favorably with experiment. Conclusions The above results suggest that the EPEN/2 functional form is flexible enough to reproduce structural and thermodynamic properties of liquid water, and is parameterized on a large number O(200-300) of dimer configurations can lead to a satisfactory description of liquid water. Since EPEN/2 parameters are transferable to a good extent, it appears attractive to use the EPEN/2 functional form in conjunction with quantum mechanical calculations to produce functions for other prototype molecules that may be combined to produce intermolecular potentials for larger systems. We are actively pursing this possibility at the present time. Acknowledgment. We thank Dr. Mihaly Mezei for furnishing the water dimer coordinates and Professor William Jorgensen for providing the program GAUSHAUS. This research was sponsored by the General Medical Sciences Division of the National Institutes of Health, Grant GM24149.
