2815
J. Phys. Chem. 1983, 8 7 , 2815-2820
which appears to be reasonable. In applying these formulas to few-electron systems, complications of the type discussed in section IV arise. VI. Conclusion In practice, then, how should one determine the chemical potential of the ground state of an electronic system of interest? The simplest way, and surely this is a great tribute to Robert Mulliken, is simply to apply eq 39-41, that is straight line interpolation between two or three
integral-N E(N) values, calculated or taken directly from experiment. To obtain the chemical potential from a single calculation, one may employ a transition-state method." Or, the explicit formulas from a particular density functional model may be used, as for example the formulas for Hartree-Fock theory given by Donnelly,lo or the formulas of classical Thomas-Fermi theory. We commend the separation of density into number and shape factors, as in eq 8 or eq 54, for further use in the density functional theory of few-electron systems. Acknowledgment. We have written this paper partly in response to a paper by Nguyen-Dang, Bader, and Ess6n,20with which we disagree. We have benefited from discussions with many people, most particularly Richard Bader, Me1 Levy, Norman March, and John Perdew. R.G.P. records with gratitude that Henry Eyring has been an inspiration to him through his entire career. What a theoretical chemist! This research was aided by grants to the University of North Carolina from the National Institutes of Health and the National Science Foundation.
(19) N. H. March, J . Phys. Chem., 86,2262 (1982). Equations 4.12 and 4.13 are trivially in error.
(20) T. T. Nguyen-Dang, R. F. W. Bader, and H. Essgn, Int. J . Quantum Chem., 22, 1049 (1982).
follows lim ( 6 T 8 / 6 p ] = -I r-a lim [6BZ,/6p] = I + p = Y.(I - A)
?-+a
(75) (76)
where I and A are the exact ionization potential and electron affinity. The functional derivatives here are not at constant N. Taken together with some results of March,lg eq 75 and 76 would imply (77) lim ( [ 6 T / 6 p ] - [ 6 T 8 / 6 p ] ]= 0 ?.-+a
A New Two-Body Water-Water Potential E. Clementi" and P. Habltz IBM Corporation, IS/TG,Poughkeepsie, New York 12602 (Received: January 10, 1983)
The potential energies for the water dimer in various geometrical configurations have been calculated with the direct configuration interaction method, including all single and double excitationsout of the Hartree-Fock reference. The basis set includes two d sets on the oxygen atom and describes the experimental dipole moment within 0.5%. We have selected 169 geometries to scan the minimum and repulsive area of the potential surface. The results match the experimental minimum geometry and energy for the dimer. Comparison with previous work, Matsuoka et al., show appreciable deviations in the repulsive area. We have fitted the quantum mechanical result and calculated the second virial coefficient,which is, unfortunately, a factor of 2 off the experimental values.
Introduction Analytical potentials to describe the interaction between two or more molecules can be obtained by using either laboratory or ab initio data. As known, there have been several attempts to obtain the interaction between water molecules with the final aim to describe water in the condensed phase; among the empirical potentials, we recall those of ref 1-4 as some of the earliest works. We are interested in ab initio potentials. The standard expansion series to describe the interaction energy for a system of N molecules can be written as follows: (1) 3. D. Bernal and F. D. Fowler, J . Chem. Phys., 1, 515 (1933). (2) J. S. Rowlinson, Tram. faraday Soc., 47, 120 (1951). (3) J. A. Barker and R. 0. Watts, Chem. Phys. Lett., 3, 144 (1969). (4) A. Rahman and F. H. Stillinger, J. Chem. Phys., 55,3336 (1971); F. H. Stillinger and A. Rahman, ibid., 57, 1281 (1972).
E(1, ..., N) = C[EO(ij) + €'(i,j)] + C[E0(ij,k)+ ~'(ij,k)+ ] ... + [to(l, ..., N)
+ ~'(1,..., N)]
(1)
where we have collected between square brackets the two-, three-, and N-body terms, each one decomposed into a Hartree-Fock contribution, 6, and an electronic correlation energy contribution, E'. In (1) we have assumed that each one of the N molecules is kept rigid at the equilibrium geometry (no vibrations). In this work we continue our analyses on the first few terms in the series of eq 1. We recall that the first water-water ab initio two-body potential of near Hartree-Fock (HF) quality was the one by Popkie et al. in 19735 and that other papers6" did quantitatively point out ~~
(5) H. Popkie, H. Kistenmacher, and E. Clementi, J . Chem. Phys., 59, 1325 (1973).
0022-3654/83/2087-2815$01.50/00 1983 American Chemical Society
2816
The Journal of Physical Chemistry, Voi. 87, No. 15, 1983
the relative importance of d ( i j ) relative to co(ij):indeed about 34% of the total two-body energy is due to the correlation energy correction (at the dimer equilibrium geometry). The three-body terms in eq 1 were also analyzed a number of times in the past 10 years. Since, in general, only a few geometrical configurations for the water trimers were analyzed, no three-body analytical potential could be derived. For a survey on this problem we refer to a monograph8 and a recent numerical studyg where the relations between to(i,j,k) and c’(i,j,k) are quantitatively analyzed. The main conclusions of these studies clearly point out (1)that the three-body correction is not negligible and not much smaller than ~ ’ ( i j(2) ) , that the nonadditivity contribution can be obtained by computations within the Hartree-Fock approximation, or (3) that the nonadditivity correction can be reasonably approximated by induction energies relations at large distances.8 We recall that the three-body terms can be defined through the two-body terms, namely, the three-body correction for the trimer (1,2,3) is given by c0(1,2,3)+ c’(1,2,3) = E(1,2,3) - C[to(i,j)+ €’(id1
(2)
Therefore it is essential to have a reliable two-body potential i n order to discuss the effects of three-body corrections. In 1976 Matsuoka et al.IOpresented an ab initio twobody potential which included the electronic correlation energy, obtained with configuration interaction (CI) methods using the ALCHEMY computer program;” this potentiallo is often referred to as the MCY potential. In the present work we present a new two-body potential obtained (a) with a more flexible basis set, (b) for a large number of geometrical conformations, and (c) with more extended CI interactions. The increase in the basis set flexibility is needed in order to obtain a more accurate dipole moment; we recall that the previously computed onelowas about 20% larger than the experimental value. As is known, the dipole-dipole interaction is a main contribution to the water-water interaction energy. The increased sampling on the potential surfaces is needed especially in the repulsive region, where relatively too few configurations were analyzed.1° Finally a more complete inclusion of the interacting configuration was advised, since the original CI computation had a rather large uncertainty’O which leads to two quite different analytical fits whose superiority could not be decided. The computation reported below has been performed with a direct CI computer program kindly made available to us by P. Siegbahn.12 A quantitative assessment of the MCY pair-potential shortcoming (relative to an “exact” two-body potential) includes an error in the fitting which can be as large as 0.5 kca/mol (the standard deviation is smaller,l0however) and (6) G. F. H. Dierckson, Theor. Chim. Acta, 21, 335 (1971). (7) D. Neumann and J. W. Moskowitz, J. Chem. Phys., 49, 2056 (1968). See also D. Hankins, J. W. Moskowitz, and F. Stillinger, ibid., 53, 4544 (1970). (8) E. Clementi, ‘Computational Aspects for Large Chemical Systems“, Springer-Verlag, Heidelberg, 1980. (9) P. Habitz, P. Bagus, P. Siegbahn, and E. Clementi, Int. J.Quantum Chem., 23, 1801 (1983). For MC simulations with two- and threebody potentials, see E. Clementi and G. Corongiu, Int. J . Quantum Chem., in press. (10) 0. Matsuoka, E. Clementi, and M. Yoshimine, J . Chem. Phys., 64, 1351 (1976). (11) See, for example, the review paper by P. S. Bagus and A. R. Williams, IBM J . Res. Develop., 25, 793 (1981). (12) B. Roos and P. E. M. Siegbahn, “Modern Theoretical Chemistry”. H. F. Schefer 111, Ed., Plenum Press, New York, 1977.
Clementi and Habitz
an inherent error (limitations in the method are extensively described in ref 10) of about 0.4 kcal/mol. This brings about an uncertainty which can be as large as 1 kcal/mol in the use of eq 2; alternatively stated, the uncertainty is sufficiently large as to limit the reliability of conclusive discussions of three-body corrections and therefore, of conclusive statements concerning computer experiments on water in the condensed state when based on ab initio two-body potentials. The MCY fit was not pushed to the limit; indeed a very accurate fit of a not too accurate set of values did not appear as a terribly rational goal. Before concluding this section we recall that we have used the Hartree-F~ck,~ the Hartree-Fock perturbati~n,’~ and the Hartree-Fock and CI’O potentials in Monte Carlo (MC) simulations of liquid water;I4-l6for each improvement in the potential, several features of the MC simulation did correspondingly improve, relative to experimental data. This is particularly true for the pair correlation functions g(OO),g(OH),and g ( H H ) and for the scattering inten~ities.’~As it was pointed out in the 1976 Monte Carlo simulation,16other features remained unsatisfactory, for example, the internal energy is too small and the pressure is too large;I6equivalent criticism were later advanced by considering the agreement with the experimental second virial coefficient.lsa None of these drawbacks was surprising since we deal with studies of liquid water performed with interaction potentials described with a series expansion truncated to the first term, the twobody contribution. Several papers published in the past few years pointed out these and similar We might note that there has been a somewhat larger effort in literature at scrutinizing the MCY potential than at recomputing it from more accurate ab initio computations.
Numerical Method As outlined in the Introduction, one target for this paper is a more complete scan of the potential surface by calculating more points than previously. However, this goal reduces the possibilities of using very large basis sets. But nevertheless, additional flexibility for the basis set compared with the previous worklo is required because the interaction of the dipole moments of the two water molecules gives the main structure of the water dimer potential surface; therefore, it is important that the dipole moment of the water molecule is well described by the calculation. According to the experience reported in ref 22-24, a diffuse d set on the oxygen atom of the water molecule and a diffuse p set on the hydrogen atom are most important to describe the correct properties of the molecule. Therefore, we have used the basis set described previously in ref 9 (13) H. Kistenmacher. H. Pookie. E. Clementi. and R. 0. Watts. J . Chem: Phys., 60, 4455 (1974); B. Jeziorski and M.van Hemert, Mol. Phys., 31, 713 (1976). (14) H. Popkie, H. Kistenmacher, and E. Clementi, J . Chem. Phys., 59, 1325 (1973). (15) G. C. Lie and E. Clementi, J. Chem. Phys., 62, 2195 (1975). (16) G. C. Lie, E. Clementi. and M. Yoshimine, J . Chem. Phys., 64, 2314 (1976). (17) A. H. Narten and H. A. Levy, J . Chem. Phys., 55, 2263 (1971); J. C. Dore, Faraday Discuss., Chem. SOC., 66, 82 (1978). (18) (a) G. C. Lie and E. Clementi, J . Chem. Phys., 64, 5308 (1976); (b) G. S. Kell, G. E. McLaurin, and E. Whalley, ibid., 48, 3805 (1968). (19) M. D. Morse and S. A. Rice, J . Chem. Phys., 74, 6514 (1981). (20) R. W. Impey, M. L. Klein, and I. R. McDonald, J . Chem. Phys., 74, 647 (1981). (21) R. W. Impey, P. A. Madden, and I . R. McDonald, Chem. Phys. Lett., 88, 589 (1982); D. G. Bounds, ibid., in press. (22) H. J. Werner and W. Meyer, Mol. Phys., 31, 855 (1976). (23) G. D. Zeiss, W. R. Scott, N. Suzuki, D. P. Chong, and S. R. Langhoff, Mol. Phys., 37, 1543 (1979). (24) (a) E. Clementi and G. Corongui, Chem. Phys. Lett., 90, 359 (1982); (b) E. Clementi and H. Popkie, J . Chem. Phys., 57, 1077 (1972).
The Journal of Physical Chemistry, Vol. 87, No. 15, 1983 2817
A New Two-Body Water-Water Potential
TABLE I: Basis Sets Used for the Calculations type basis set exponent coefficient Oxygen S S S
S S S
S S S S
S
P P P P P P P d d S S S
S S S
P
3 1 195.6 4 669.38 1062.62 301.426 98.515 3 35.460 9 13.617 9 5.386 1 8 1.538 73 0.605 50 0.220 64 114.863 26.876 7 8.320 77 2.972 37 1.128 48 0.423 60 0.150 74 1.0 0.15
Hydrogen 115.691 17.373 4 3.953 4 1.116 74 0.361 27 0.125 85 0.40
0.000 210 0.001 628 0.008 450 0.034 1 9 1 0.110 3 1 1 0.269 493 0.423 545 0.283 038 0.238 7 0 1 0.591 286 1 0.002 266 0.017 1 9 2 0.075 3 4 1 0.212 775 0.370 478 0.398 523 1
S
II
-. i
H u-H
Flgure 1. The different geometries, labeled with capital letters, which are considered In this paper. The reference molecule, I in the view from the top of the molecular plane ( z direction), I1 in the molecular plane (y dkectbn), is identified by the chemical symbols for oxygen and hydrogen. The arrows mark the direction in which the potential surface is scanned.
1 1 0.000 9 1 0 0.005 442 0.031 700 0.111 320 0.155 788 1
TABLE 11: Contributionsto the Minimum Potential Curve (millihartree)
1
E = (1 - C2)Ecorr
Tlh
K* c
with an additional d set on oxygen with exponent I = 0.15, and instead of the p set with exponent ( = 1 on hydrogen, we used a p set with exponent I = 0.4, which represents a compromise between an energy-optimized and a property-optimized exponent. The full basis set is given in Table I. The HF energy in this basis set for a single water molecule in experimental geometry is -76.055 318 au (to be compared with a near Hartree-Fock limit value24of -76.06587 au). The CI calculation with all single and double excitations generated from the HF reference configuration, the 1s core of oxygen frozen and the ls* deleted from the virtual space, yields an energy of -76.267 725 au. (The calculations include 108 000 configurations for the case of the dimer.) The dipole moment in the H F approximation is 0.782 au and in CI is 0.735 au. This is sufficiently near the experimental value of 0.730 au.25 A CI calculation including all single and double excitation is not size consistent. In the asymptotic limit for infinite distance (100 au) between the two water molecules, the CI energy (-152.515 172 au) does not approach the sum of the CI energies of two water molecules, which is lower by 0.020278 au (12.7 kcal/mol) (see Table 111). This well-known effect can be understood by considering the contribution of the quadruple excitations, which are included in the product of the wave functions of the two molecules, but not in the asymptotic limit of our calculation. Therefore, we choose for reference the calculation at 100 au distance. Davidson%has suggested a simple relation which, on the basis of perturbation arguments, allows an estimate of the energy contribution E of the quadruple excitations, namely (3)
(25)T.R.Dyke and J. S. Muenter, J. Chem. Phys., 59, 3125 (1973). (26) S.R.Langhoff and E. R. Davidson, Int. J. Quant. Chem., 8,61 (1974).
p'
I
~
~
R, au
SCF
CI
4.50 4.75 5.00 5.25 5.50 5.75 6.00 6.50 7.00 7.50 8.00
12.440 3.070 - 2.059 -4.682 - 5.847 -6.184 -6.062 - 5.237 -4.248 - 3.400 -2.674
6.545 -1.953 -6.271 - 8.205 -8.814 -8.709 -8.210 -6.798 - 5.348 -4.200 -3.180
E,,, - 5.895 -5.023 -4.212 -3.523 -2.967 -2.525 -2.148 -1.561 -1.100 -0.800 -0.506
Davidson
SCFSUP
final
-0.480 -0.453 -0.428 -0.394 -0.361 -0.330 -0.290 -0.231 -0.173 -0.135 -0.089
0.66 0.57 0.50 0.45 0.41 0.38 0.36 0.31 0.26 0.20 0.15
6.725 -1.836 -6.199 -8.149 -8.765 -8.659 -8.140 -6.718 -5.261 -4.135 -3.119
where C is the coefficient of the HF determinant in the total wave function and E,,,, the correlation energy. By inclusion of the quadruple excitations with (3), the CI theory is lowered for the monomer to -76.279 068 au and for the dimer to -152.551 587 au. The size consistency effect in this approach is a factor of 3 smaller, but still 4.1 kcal/mol (see Table 111). This result, on one hand, shows the approximate character of the Davidson correction, eq 3, but, on the other hand, it points out that it represents a notable step in the right direction. Therefore, we included in our calculations the Davidson correction and took the corrected dimer at 100 au distance as a reference. We chose a total of 165 geometries; 34 of these are done only in HF approximation with an interpolated correlation contribution added to the interaction energy. This interpolation is reliable and easily applicable since, for a given geometry, the correlation contribution shows, to very good accuracy, an exponential dependence from the oxygen-oxygen distance. The geometries are chosen so (a) that the second water molecule approaches the rigidly held first molecule from many directions, (b) that all types of interaction are represented, and (c) that as much symmetry as possible is maintained to reduce computer time. The different geometries are shown in Figure 1. In this way we have scanned not only the minima of the interaction potential, as stressed by Matsuoka et a1.,l0but also many repulsive interactions with the closest approach at two oxygen atoms, two hydrogen atoms, and oxygen and hy-
2818
The Journal of Physical Chemistry, Vol. 87,No. 15, 1983
Clementi and Habitz
TABLE 111: Total Energies in t h e S t a n d a r d (Old) and Increased Basis S e t (New) ( a u ) monomer dimer a t R = 1 0 0 dimer a t m i n i m u m SCF CI
CI"
" With Davidson
old new old new old new
-76.055 319 --76.055 745 - - 7 6 . 2 6 77 2 5 - 76.268 8 6 3 -76.279 068 - - - 7 6 , 2 8 03 3 0
-152.110 637 -152.111 492 -152.515 172 --152.517 271 -152.551 687 -152.553963
-152.116 484 -152.117 260 -152.523988 -152.525 889 -152.560 1 6 2 - 152.562 925
interaction energy -0.005 -0.005 -0.008 -0.008 -0.009 -0.008
847 768 814 618 175 962
correction.
drogen atoms. However, in order to save computer time, we have excluded from the selection of the geometries cases without some symmetry; this selection might be questionable. Table I1 details our results for the minimum potential curve (see Figure 1, case T) with the oxygen-hydrogenoxygen atoms on a straight line, the two planes defined by the three atoms of each molecule intersecting a t a right angle and the line through the hydrogen bond forms a 60" angle with the intersection of the two planes. In the first column in Table I1 we show the oxygen-oxygen distance, then the SCF and CI interaction energy, the correlation contribution for the single and double excitations, and the contribution of the Davidson correction to the interaction energy. The sixth column shows the SCF superposition error in the given basis set. These numbers show the change in the SCF interaction energy if the reference is not calculated at infinite distance, but taken as the sum of the two systems each calculated in the basis set of the dimer at given geometry. In this way, we take care of the additional flexibility for the description of the monomer caused by the basis set of the other molecule. By including the SCF superposition error and the Davidson correction, we get the final interaction energy which was used for the fitting of the dimer potential. Obtaining an estimate of the CI superposition error is more difficult t h a n for the SCF case. Of course it is possible to perform CI calculations for the monomer in the dimer basis set at different geometries. This we have done for geometries B, C, and E, and we show the result in Figure 2. We have plotted the CI superposition error vs. the correlation contribution of the interaction energy. We observe, with some surprise, a linear correlation. The computed values for the CI superposition error are of questionable validity, however. We recall that by computing the CI superposition error as above outlined, the sum of two monomers, each in the dimer basis set, is selected as the reference system. This includes an increase of the size consistency error because the doubled basis set will increase the number of quadruple excitations by a factor of about 100. This size consistency effect, the neglect of the unlinked quadruple excitation, in the language of the cluster expansion, raises the reference and reduces the CI superposition error appreciably. T o obtain a feeling for this calculation, we increased the basis set by an additional diffuse s and p basis set on the oxygen atom of the water molecules with exponents { = 0.07 and { = 0.05, respectively. In this basis set, the SCF-superposition error is reduced by a factor of 3, so it can be expected that the CI superposition error will also be very small. The results are shown in Table 111. The new basis set yields higher (less attractive) energy at the minimum geometry as expected from the sign of the superposition error, but the effect is much smaller than expected from Figure 2. Of course, the new basis set does much more than only reduce the superposition error-the change in the SCF interaction energy is smaller than ex-
Figure 2. The C I superposition error plotted vs. the correlation contribution of the interaction energy. Details are outlined in the text.
pected from the superposition error in Table 11-but it shows that the calculated CI superposition error in Figure 2 is too large. Therefore, we have decided to neglect the CI superposition error; from Table I11 we have some reason to expect that a cancellation of errors will favor this decision.
Results and Discussion For 16 of the 20 geometries of Figure 1, the computational results are shown in Figure 3 as points. Full listing of all coordinates and energies are available from the authors by request. The minimum of this calculated potential energy surface is a t 5.5, kcal/mol a t Roo = 2.97 A (case T) and 0 = 50 (5)' (case N). This compares very well with experimental data obtained by molecular beam electric resonance s p e c t r o s ~ o p y :Roo ~ ~ = 2.97 A, 0 = 51 (6' for the proton accepting water axis and 0 = 58 (6)' for the proton-donating water axis with respect to Roo. Because of the expectation of only a small energy differences, we did not optimize the position of the hydrogen bonding proton. The binding energy is given by Curtiss as 5.5 f 0.5 kcal.2s For the "bifurcated" geometry, type C, which has also a minimum geometry as shown in Figure 3, our result for the minimum is Roo = 2.96 A, E = 3.61 kcal. This compares well with the results a t other calculations as they are collected by Van Heusbergen et al.29in their Table I. We recall that this structure is almost 2 kcal higher than the minimum geometry. We have fitted our results with an Ansatz similar to the one used in ref 10, and we have weighted and energies near the minima, heavily. As explained in ref 10, we define a fitting point N in the plane of the water molecule on the (27) T. R. Dyke, K. M. Mack, and J. S. Muenter, J.Chem. Phys., 66, 498 (1977); T. R. Dyke and J. S. Muenter, ibid., 60, 2929 (1974). (28) L. A. Curtiss, D. J. Frurip, and M. Blander, Chem. Phys. Lett., 54, 575 (1978). (29) B. Van Hensbergen, R. Block, end L. Jransen, J. Chem. Phys., 76, 3161 (1982).
The Journal of Physical Chemistty, Vol. 87, No. 15, 1983 2819
A New Two-Body Water-Water Potential
TABLE IV: Parameters for Eq 6 and 7
i
WW
Ai
1
0-0
2 3 4 5
0-H 0-H H-H H-H
7 2 317.219 416 026 4 -7.969 812 258 0 9 363.718 1 0 6 377 3 1 6 0 3 . 4 1 8 394 1 0 1 9 7.952 022 940 0 0.471 636 359 2
P
ai
1.892 938 049 4 0.585 972 034 0 2.347 935 855 1 1.896 609 625 2 0.373 522 098 4
symmetry axis, which contains the negative point charge (-2q). Equal positive point charges g are on the hydrogen atoms. The charges g and the position of N on the symmetry axis are coupled, to ensure the correct dipole moment, by the following equation: g = 0.331927/(1 - P) au (4) and the distance R(0-N) of N from its oxygen atom 0 is
(5) R(0-N) = 1.107168P au Equations 4 and 5 define the fitting parameter p . For a correct description of the repulsive part of the potential it was necessary to introduce a second exponential for the hydrogen-hydrogen interaction, and for numerical stability it was necessary to multiply the attractive oxygen-hydrogen exponential with R(0-H). The final form is E = El(R(0-0’)) + E,(R(O-Hi’)) + E,(R(O-H,’)) + E,(R(O’+H,)) + E,(R(O’+H,)) + E3(R(H1-Hl1) + Es(R(H1-H;)) + E3(R(HZ-Hl’)) + E3(R(Hz-H,’) + g2[4/R(N-N’) l/R(Hl-Hl’) + l/R(H,-H;) + l/R(H,-Hl’) + l/R(H,-H2’) 2/R(N-H1’) - 2/R(N-H2’) - 2/R(N’-HJ 2/R(N’-H2) (6)
+
L
where
El(R(O-O)) = Al exp(-a,R(O-O)) EZ(R(0-H)) = A,R(O-H) exp(-a,R(O-H)) + A3 exp(-a3R(O-H)) (7) E3(R(H-H)) = A4 exp(-a,R(H-H)) + A , exp(-a&(H-H))
-3
c
-3
c
Flgure 3. The 16 insets are labeled according to the geometries explained in Figure 1. The crosses corresponds to the computational results, the solid line to the fit given in Table IV, and the dotted line to the fit I1 of Table V I 1 of ref 10.
The parameters Al, ..., A,, ayl,..., a,, and P are given in Table IV in atomic units. The resulting energies are in milliiartree. The standard deviations for all points smaller than 1, -1, and -2.8 kcal are 0.026, 0.22, and 0.17 kcal, respectively. The different standard deviations result from the stronger weight for lower energies in the fit. The results of this fit are shown in Figure 3 as solid lines together with the second fit of Matsuoka et al.1° (dotted line). These pictures show appreciable differences. In the repulsive area the old calculation is more repulsive, with a shift to larger distances and too deep minima. While the latter difference is due to the too large dipole-dipole interaction in the old calculation, the former is due to the interpolation procedure for the fit. Indeed the old calculation stressed the mininum area in the selection of the dimer’s geometries; for each direction of approach the angular position of the water molecule was selected to yield only curves with minima. The description of repulsive curves was obtained by extrapolation which cannot be sufficiently accurate. The too repulsive potential also partly explains why the liquid structure calculations, done with the old fit, show too high pressure. Figure 4 shows a two-dimensional scan of the potentials. In these pictures the reference molecule is held fixed in positions I and I1 of Figure 1, and the second water molecule is systematically moved through the x-y plane. For each position of the oxygen atom the Eulerian angles,
2820
The Journal of Physical Chemistry, Vol. 87, No. 15, 1983
B
C
Figure 4. Isoenergy maps (B) and its threedimensional representations (A) for two cross sections of the interaction energy hypersurface. The interactions are computed for given ozygen positions, with angles optimized for minimum energy, shown in C. Units in millihartree.
TABLE V : Computed and Experimental Second Virial Coefficient (cm3/mol)
a
T,"C
exptla
this work
MCY potlb
100 200 300 400
-580 -212 -117 -73
-1124 -433 -233 - 147
-827 -289 -141 -78
Reference 18b.
Reference 18a.
corr scaledC -497 -216 -120 -73
See text.
defining the angular position of the second water molecule in space, are optimized to minimize the interaction energy. The contour to contour energies are given in the figure in millihartree. Parts A and B show this minimum energy dependence of the position of the oxygen atom in threedimensional and contour line representation. Part C shows resulting geometries. A striking feature of this figure is the nearly circular shape of the minimum with the oxygen atom at center; alternatively stated, the potential is not as structured as in ref 10. Part I shows the strongly preferred direction along the hydrogen bond and part I1 the weaker preferred direction of the lone electron pair. In Table V we report the results for the second virial coefficient at four different temperatures. The computation was performed with the technique and program of ref 18a. Unfortunately, the virial calculated with our new fit is a factor of two off the experimental value,lsb constant over the whole temperature range. We refer to ref 18a for a comparison of several two-body potentials and their corresponding second virial coefficients. There it is shown that the virial is also too high in absolute value and notably different for the different fits, but not parallel (on the
Clementi and Habitz
logarithemic scale) to the experimental values; however, the overall agreement of the MCY potential is markedly better than the one found with this new potential. We have made another fit through the result of ref 10 giving more weight to the already overstressed minimum points; with this new fit which had the same standard deviation as previouslylO the second virial was notably worse (about a factor of 2.5 off). This result shows that the virial is very sensitive to the fitting technique (namely, both Ansatz selection and interpolation). Therefore, we have tried a few more flexible analytical representations for the fit of our two-body interactions in order to achieve a better description of the second virial coefficient. Unfortunately our trials were not successful. In other trials, the correlation contribution to the water-water interaction was decreased by a factor 0.5; such a factor is not unreasonable if one accepts as valid the estimates (the linear dependence) of the CI superposition correction. In this way, we obtained a potential with a correct second virial (Table V, column 5 ) , but the water dimer minimum energy is now only -4.7 kcal, a value notably small. After a number of trials we decided to neglect the corrections for the CI superposition error and we posponed the match to the experimental second virial coefficients for further investigations. We realize that, by so doing, we have restricted our original goals simply to a more complete scan of the potential energy surface (relative to the one of ref 10) and to a better basis set. Conclusions The new ab initio two body-water interaction potential, above presented, is rather close to the exact one at least near equilibrium; we feel that the maximum error might be about 0.6-0.9 kcal/mol. However, a substantially larger basis set24bwith inclusion off functions (oxygen atom) and d functions (hydrogen atoms) might lead to unexpected and large changes. The CI computations of the dimer are more accurate than those by Matsuoka, Yoshimine, and Clementi,lo which up to now was the most accurate ab initio energy surface for water-water. We recall that this two-body potential can be further improved by using available experimental data on the interaction of two water molecules. An improvement from quantum mechanical computations seems, however, less straightforward if carried out with current computational methods and/or facilities, since a substantially large basis set24would be required. Therefore we are now working at a new two-body potential using a density matrix functional approach30 which has none of the above drawbacks. Molecular dynamics simulations on liquid water using the potentials given in this paper are reported e l ~ e w h e r e . ~ ~ Acknowledgment. It is a pleasure to acknowledge the permission to use a number of computer programs, such as a direct CI program,12a second viral program,lsb and a fitting program, the latter made available from M. Gratarola (Institute G. Donegani, Novara, Italy). Registry No. Water, 7732-18-5. (30) V. Caravetta and E. Clementi, to be published. (31) P. Habitz and E. Clementi, to be published.
