J. Phys. Chem. 1985,89, 1461-1467
1461
Estimation of Dimer Coulombic Intermolecular Energy and Site Charge Polarization by the Potential-Derived Method Donald E. Williams* and David J. Craycroft Department of Chemistry, University of Louisville, Louisville. Kentucky 40292 (Received: July 16, 1984; I n Final Form: November 12, 1984)
Water dimer was used as a sample system to investigate intermolecularsite charge Coulombic energy and site charge polarization. Optimum site charge models for water monomer and 52 configurationsof water dimer were obtained by the potential-derived method. Electric potentials at points in space around the monomer or dimer were calculated directly from ab initio 6-31G** basis SCF wave functions. A site charge model which included out-of-plane sites gave significantly better agreement to the electric potentials of both monomer and dimer, as compared to a net atomic charge model. The Coulombic intermolecular energies of the dimer configurations were calculated by using site charges optimized for the monomer. The polarization of the site charges in the dimer was found by comparing them with corresponding monomer charges. The purpose of the study was to examine trends in site charge polarization and the intermolecular Coulombic energy component to aid in the design of intermolecular potential energy functions applicable to large molecules.
Introduction To the first approximation for a neutral molecule, the electric potential at points in space outside the boundary of its electron cloud is zero. However, even in such neutral molecules the electron cloud does not distribute itself in such a way as to yield a zero electric potential a t relatively close distances, because of the differing electronegativities of the atoms. It is convenient to regard this electric potential around the molecule as arising from the presence of site charges. A site charge is not an exactly defined physical quantity and approximations must be made in its calculation. A net atomic charge is a special case of a site charge where the nuclear positions are chosen to represent the surrounding electric potential. The energy of intermolecular electrostatic interaction in a molecular dimer may be approximated by application of Coulomb's law for interaction between charge sites of the monomers. The great advantage of using site charges to represent the electrostatic interaction energy is that the energy evaluation is extremely fast and therefore the site charge method can be applied to larger molecules or to clusters of molecules where a more elaborate treatment may not be practicable. For the very large number of intermolecular interactions in a crystal it is highly desirable to have a site charge model which can be rapidly evaluated as a lattice sum to give the Coulombic energy of the crystal. Until recently the population analysis procedure' has been the most widely used method for the determination of net atomic charges in molecules. In population analysis (PA) net atomic charges are obtained from the coefficients of a molecular orbital wave function formed from a linear combination of atomic orbitals (LCAO-MO method). PA net atomic charges are not a definable quantum-mechanical physical property, but they have gained wide use since they show trends which correlate well with electronegativity and chemical reactivity. PA charges are not defined for wave functions which are not of the LCAO-MO type. A number of attempts have been made to obtain net atomic charges by other methods. In contrast to the definition of population analysis charges, the electric potential at points in space around a molecule is a rigorously defined quantum-mechanical property obtainable from any type of wave function. Scrocco and Tomasi2 have reviewed methods of directly calculating the molecular electric potential from the molecular wave function. The most recent development has utilized the idea of calculating the electric potential at a grid of points in space around the molecule directly from its wave function and then finding by a least-squares fitting technique optimum net atomic charges which (1) R. S. Mulliken, J. Chem. Phys., 23, 1833 (1955). (2) E. Scrocco and J. Tomasi, Adu. Quantum Chem., 11, 115 (1978).
0022-365418512089-1461$01.50/0
will best reproduce this surrounding molecular p~tential.~" These site charges are designated potential-derived (PD) to distinguish them from the population analysis (PA) charges. PD charges generally were found to give a more accurate portrayal of the electric potential than PA charges. The O-H- SOhydrogen bond has been studied extensively in the water system7-15and also in organic molecular systems such as alcohols and carboxylic acids.1618 One extreme viewpoint is that the hydrogen bond can be understood solely on the basis of an electrostatic interaction. A site charge analysis for water dimer may be used to give an estimate of the portion of the intermolecular interaction energy which can be viewed as originating from Coulombic interactions between the sites. In order to develop more quantitative and practically useful empirical models for the hydrogen bond, it is desirable to be able to separate out the Coulombic intermolecular energy from other components. In this view, the intermolecular interaction energy is a sum of pairwise additive atom-atom terms such as V(repu1sion)
+ V(dispersion) + V(Cou1ombic) + V(hydrogen bond)
summed over atom pairs from different molecules. The first three terms can be represented by an (exp 6-1) nonbonded potential energy function bk(eXp 6-1) = B exp(-Cr,k) - h , k d
+ 4,qkrjk-l
(3) F. A. Momany, J . Phys. Chem., 82, 592 (1978). (4) P.H. Smit, J. L. Derissen, and F. B. VanDuijneveldt, Mol. Phys., 37, 521 (1979). ( 5 ) S. R. Cox and D. E. Williams, J . Compur. Chem., 2, 304 (1981). (6) D. E. Williams and R. R. Weller, J . Am. Chem. Soc., 105, 4143 (1983). (7) A. Ben-Naim and F. H. Stillinger, 'Structure and Transport Pr0Cese.s in Water and Aqueous Solutions", R. A. Home,Ed., Wiley, New York, 1972. (8) F. M. Stillinger and A. Rahman, J. Chem. Phys., 60, 1545 (1974). (9) E. Huler and A. Zunger, Chem. Phys., 13,433 (1975). (10) H. L. Lemberg and F.H. Stillinger, J. Chem. Phys., 62,1677 (1975). (1 1) G. H. F. Diercksen, W. P.Kraemer, and B. 0. Roos, Theor. Chim. Acra, 36, 249 (1975). (12) R. 0. Watts, Chem. Phys., 26, 367 (1977). (13) W. L. Jorgensen, J . Am. Chem. Soc., 101, 2011 (1979). (14) W. L. Jorgensen and M. Ibrahim, J. Am. Chem. SOC.,102, 3309 (1980). (15) F. T. Marchese, P. K. Mehrotra, and D. L. Beveridge, Phys. Chem. Lert., 85, 1 (1981). (16) J. Kroon, J. A. Kanters, J. G. C. M. van Duijneveldt-van de Rijdt, F. B. van Duijneveldt, and J. A. Vliegenthart, J. Mol. Srrucr., 24, 109 (1975). (17) J. Mitra and C. Ramakrishnan, Inr. J. Protein Res., 9, 27 (1977). (18) C. Ceccarelli, G. A. Jeffrey, and R. Taylor, J . Mol. Srrucr., 70, 255 (1 98 1).
0 1985 American Chemical Society
1462 The Journal of Physical Chemistry, Vol. 89, No. 8, 1985
Williams and Craycroft
where rjk is a nonbonded distance, qj is a net atomic charge, and A-C are adjustable parameters. This function has enjoyed considerable success in modeling the interaction energy of molecules which are not hydrogen bonded in the c r y ~ t a l l and ~ - ~has ~ been used to model molecular clusters.24 The empirical form for V(hydrogen bond) which is compatible with (exp 6-1) functions for non-hydrogen-bonded oxygen and hydrogen atoms is not well established. A clear separation of the Coulombic component would assist in defining the empirical form of V(hydrogen bond). The PD site charge method may be extended to explore site charges in molecular dimers, such as water dimer. In this approach the electric potential surrounding the dimer is found from the dimer wave function in a way that is completely analogous to that used for the monomer. The optimized site charges in the dimer which best reproduce this electric potential can be found by least squares. An approximation to the Coulombic interaction energy of the dimer can then be obtained by evaluating Coulomb’s law directly between the site charges of the different monomers. Further, it is of interest to note whether the charges in the molecular components of the dimer remain essentially unchanged, or if there are significant changes from monomer values as the molecules approach to form the dimer. For our purposes, the change in the site charges going from monomer to dimer is defined as the dimer site charge polarization. In this work we report site charges for water monomer and for water dimer in various configurations. The dimer configurations which were considered included linear, cyclic, bifurcated, and intermediate-type hydrogen-bonded geometries. The reader should consult ref 27 for a detailed definition of these geometries. The net atomic charge site model was modified to utilize symmetrical nonnuclear sites above and below the twofold axis in the plane of the molecule instead of the oxygen site. The use of these charge sites along with foreshortened hydrogen sites gave an improved description of the surrounding electric potential for both the monomer and the dimer. The intent of this study was not necessarily to arrive at the best models for water monomer and dimer. Rather, our intent was to probe this relatively simple system with the idea of obtaining background information on Coulombic site charge interaction and site charge polarization which would be helpful in developing empirical models for nonbonded and hydrogen-bonded energy functions for larger molecules. Methods Wave functions for water monomer and the various configurations of water dimer were calculated by the a b initio self-consistent field molecular orbital method.z5 The geometry of the water molecule was fixed at the experimental values for the bond length, 0.9572 A, and bond angle, 104.52°.26 Fifty-two water dimer geometries were considered. These dimer geometries and their labeling were consistent with those defined by Matsuoka, Clementi, and YoshimineZ7in their configuration interaction calculations of the energy of water dimer. As a compromise between demands for accuracy and the computer time required, the 6-3 1G** Gaussian basis set was s e l e ~ t e d .This ~ ~ ~is ~a ~split valence basis set which includes d-type polarization functions for oxygen and p-type polarization functions for hydrogen. For comparison (19) A. I. Kitaigorodski,“Molecular Crystals and Molecules”, Academic Press, New York, 1973. (20) D. E. Williams and T. L. Starr, Comput. Chem., 1, 173 (1977). (21) L. Y. Hsu and D. E. Williams, Acta Crysfallogr.Sect. A , 36, 277 (1980). . (22) S. R. Cox, L Y . Hsu, and D. E. Williams, Acta Crystallogr., Sect. A , 37, 293 (1981). (23) D. E. Williams, Top. Curr. Phys., 26, 3 (1981). (24) D. E. Williams, Acra Crystallogr., Sect. A , 36, 715 (1980). (25) U.C. Singh and P. Kollman, QCPE, 15, 446 (1982). (26) W. S. Benedict,N. Gailar, and E. K. Plyler, J . Chem. Phys., 24, 1159 (1956). (27) 0. Matsuoka, E. Clementi, and M. Yoshimine, J . Chem. Phys., 64, 1351 (1976). (28) W. J. Hehre, R. F. Stewart, and J. A. Pople, J . Chem. Phys., 51,2657 (1969). (29) P. C. Hariharan and J. A. Pople, Theor. Chim. Acta, 28,213 (1973).
Figure 1. Charge site placement for the MCY and SC models. The
position of the twofold symmetry axis is shown. Distances and angles are given in the text.
purposes, PD net atomic charges for a variety of small molecules are available derived with the use of this basis set.5 The electric potential for a unit positive charge at a point r in the vicinity of the monomer or dimer containing nuclei A at positions R A and with charge ZAis given by
where P,, is an element of the density matrix of the S C F molecular wave function and +, are the atomic wave functions used as a basis set. The first term is summed over the nuclei and the second term is summed over the electrons of the system. The set of points selected for the evaluation of the electric potential was a cubic grid of 1-8, spacing in a 1.2-A thick shell around the molecule. The inner surface of this shell was selected to represent the distance of closest normal intermolecular approach and was taken at the van der Waals radius of the nearest atom plus the van der Waals radius of the smallest approaching atom, hydrogen. The van der Waals radii of oxygen and hydrogen atoms were taken as 1.4 and 1.2 A. This procedure led to the definition of 195 grid points around water monomer and approximately 300 grid points around each dimer configuration. Cox and WilliamsS tested the suitability of this particular grid selection for water monomer and for a number of other small molecules. The electric potential at each grid point was calculated directly from the ab initio wave function. The PD model assumes that the electric potential around the molecule or dimer can be approximated by placing variable point charges at certain sites within the molecule or dimer. A reasonably good portrayal of the surrounding electric potential can be obtained by using nuclear sites only; this is our PD net atomic charge model. Further significant improvement of the fit to the electric potential of water monomer or dimer can be obtained by not using the oxygen site but instead including two symmetrical charge sites out of the plane of the water molecule above and below the twofold symmetry axis: this is our site charge (SC) model. Matsuoka, Clementi, and YoshimineZ7used a model similar to S C except that in their empirical treatment both sites were compressed together onto the twofold axis: this is the MCY model (Figure 1). For both the SC and MCY models the off-nuclear charge sites (which carry a negative charge) are on the hydrogen side of the water molecule. The site charges were determined by fitting them to the molecular electric potentials at the grid points by least-squares. The quantity minimized, for m grid points and n atoms, was given by m
n- 1
R = Cw1[V,O- Z qJr,,-’ + 1
J
ll-1
(C q, - Z)rl,,-llz J
where V,O is the calculated quantum-mechanical electric potential at point i, qJ is the charge at sitej, rlJis the distance between site j and the i-th grid point, and Z is the net charge on the molecule (zero for the water monomer and dimer here considered). The number of independent site charges is decreased by 1 since there is the dependency condition requiring the sum of the site charges be equal to Z . Also, the symmetry of the monomer and of some of the dimer configurations establishes further dependency conditions.
The Journal of Physical Chemistry, Vol. 89, No. 8, 1985 1463
Estimation of Dimer Coulombic Intermolecular Energy
TABLE I: Potential-DerivedNet Atomic Charges and Calculated 6-31C** SCF Energies for 52 Configurations of Water Dimer (PD Model)" net atomic chargeb structure
01
H1
H2
monomer A1 A2 A3 A4 A5 B10 B11 B12 B13 B14 B18 C19 c20 c 21 c22 C23 C27 D28 D29 D30 D3 1 D32 E37 E38 E39 E40 E4 1 F42 F43 F44 F45 F46 G47 G48 G49 G50 G5 1 H52 H53 H54 H55 H56 I57 I58 I59 I60 I6 1 562 563 564 565 566
-786 -820 -808 -799 -792 -792 -937 -907 -883 -868 -846 -789 -829 -818 -806 -801 -795 -784 same as B10 -99 1 -934 -842 -860 -977 -904 -868 -850 -824 -859 -837 -826 -816 -816 -847 -824 -814 -80 1 -799 -864 -834 -817 -806 -792 -812 -81 1 -806 -801 -795 -875 -859 -860 -851 -829
393 424 422 416 41 1 405 49 1 472 456 447 432 395 399 381 369 369 369 396
393 424 422 416 41 1 405 491 472 456 447 432 395 395 406 410 41 1 412 388
02 -840 -854 -840 -823 -804 -97 1 -955 -93 1 -918 -876 -791 -785 -778 -769 -775 -775 -784
475 458 424 427 534 48 1 460 450 426 430 422 413 407 406 434 420 414 404 405 418 398 391 386 386 413 402 395 391 389 450 44 1 43 1 425 415
475 458 424 427 416 412 409 408 403 436 434 430 423 416 416 41 1 407 400 400 416 416 412 410 403 416 424 424 421 413 434 437 438 433 419
-930 -859 -737 -806 -849 -873 -889 -898 -861 -808 -815 -804 -792 -783 -836 -837 -839 -820 -827 -823 -814 -809 -806 -802 -825 -816 -806 -799 -789 -817 -852 -843 -831 -808
H4
ae
sd
EC
377 397 406 410 381 463 459 45 1 446 429 395 405 409 397 399 393 396
435 42 1 40 1 383 405 463 459 45 1 446 429 395 415 400 399 397 396 388
2.8 2.6 2.0 1.9 2.0 2.1 1.5 1.5 1.7 1.8 2.0 2.4 2.4 2.4 2.3 2.4 2.4 2.6
8.3 10.0 8.0 7.8 7.8 8.1 3.2 3.5 3.9 4.4 5.4 8.4 10.0 10.0 9.6 9.7 9.4 9.1
8.24 -16.12 -22.48 -21.16 -13.84 18.64 -9.72 -16.28 -16.28 -1 1.55 -1.05 5.25 -16.59 -19.69 -17.33 -10.24 -0.79
550 49 1 407 423 438 442 444 445 428 421 396 396 397 394 402 395 435 424 423 419 43 1 425 393 391 40 1 390 382 379 377 369 393 399 402 402
421 386 324 389 438 442 444 445 428 380 400 391 38 1 383 43 1 435 397 393 398 434 403 398 423 414 407 41 1 41 1 409 405 439 440 435 422 40 1
1.6 1.9 2.9 2.2 1.5 1.6 1.7 1.8 2.2 2.5 2.5 2.5 2.5 2.6 2.1 2.0 2.1 2.0 2.2 2.2 2.1 2.0 2.0 2.0 2.4 2.2 2.1 2.1 2.1 2.2 2.1 2.0 2.0 2.2
3.3 4.4 6.9 5.7 3.1 3.7 4.1 4.8 6.1 8.7 8.0 7.8 7.7 8.3 6.3 6.3 6.5 6.4 7.1 9.4 8.6 8.2 7.8 7.5 10.4 8.7 8.0 8 .O 7.7 6.8 5.4 5.1 5.4 6.5
-1 1.34 -20.74 -20.19 -1 3.97 26.44 -13.03 -21.99 -21.14 -14.53 4.59 -14.41 -19.04 -17.72 -1 1.82 9.14 -14.41 -18.48 -16.59 -10.71 9.00 -17.14 -21.50 -18.96 -1 1.47 0.58 -17.30 -20.45 -18.14 -11.26 14.86 -13.50 -21.74 -20.82 -14.23
H3
'The numbering of the dimer structures is the same as MCY. bElectronic units X 10'. 'In kJ/mol.
The minimum value of R was found by taking the first partial derivatives of R with respect to the independent charges and solving the resulting set of linear equations by the usual leastsquares method. The statistical weight, w,,for each point was taken as unity, which had been found satisfactory by Cox and Williams.s The goodness of fit for the site charge models was shown by the rms fit parameter, which is defined as
and a relative rms fit, s, expressed as i
Monomer Site Charge Analysis A comparison of population analysis and potential-derived net atomic charges for water monomer has been given previou~ly.~ Using the 6-31G** basis set led to PA charges which were 0.86
percent.
of the PD charges obtained with the same basis set. The PA charges gave a calculated dipole moment of 1.896 D (1 D = 3.336 X C m) while the PD charges gave 2.21 1 D; the observed dipole moment is 1.846.30 The PD vaue is more internally consistent since it is in better agreement with the dipole moment calculated by applying the dipole moment operator directly to the wave function, 2.184 D. An overestimate of the dipole moment by the 6-31G** wave function is typical; the average scale factor for the ratio of calculated to observed dipole moments for various molecules is 0.91 for this basis set.s The relative rms fit of the PA charge model to the electric potential points was 16.6%. Again, the PD model was more internally consistent, since it gave a n improved relative rms fit of 8.3%. Matsuoka, Clementi, and YoshimineZ7fitted their results for the energy of water dimer with an empirical model which included net atomic charges on the hydrogens but no net charges on the (30) G. Birnbaum and S.K.Chatterjee, J . Appl.Phys., 23, 220 (1952).
1464 The Journal of Physical Chemistry, Vol. 89, No. 8,1985
Williams and Craycroft
TABLE II: Potential-DerivedSite Charges for 52 Configurations of Water Dimer Using Out-of-Plane Charge Sites (SC Model)"
structure monom er A1 A2 A3 A4 A5 B10 B11 B12 B13 B14 B18 C19 c20 c21 c22 C23 C27 D28 D29 D30 D3 1 D32 E37 E38 E39 E40 E4 1 F42 F43 F44 F45 F46 G47 G48 G49 G5O G5 1 H52 H53 H54 H55 H56 157 I58 I59 I60 I6 1 562 563 564 565 566
s1 -630 -664 -648 -637 -634 -634 -696 -674 -660 -65 1 -642 -63 1 -673 -662 -654 -649 -642 -63 1
H1 630 673 66 1 650 644 639 686 665 654 646 640 63 1 695 668 654 647 643 632
site chargeb H2 s2 630 -66 1 673 -670 66 1 -66 1 650 -65 1 644 -640 639 -630 686 -636 665 -64 1 654 -64 1 646 -635 640 -632 63 1 -655 63 1 -647 639 -642 642 64 1 -640 637 -639 -63 1 630
H3
H4
Lrc
sd
609 632 638 639 632 640 645 647 646 637 632 650 645 642 640 639 632
695 682 658 643 638 640 645 647 646 637 632 680 666 654 650 643 630
0.3 0.4 0.4 0.4 0.4 0.4 0.8 0.5 0.4 0.4 0.4 0.3 0.5 0.5 0.5 0.4 0.4 0.3
1.o 1.9 1.7 1.6 1.5 1.6 1.7 1.2 1.o 1 .o 1.2 1.2 2.3 2.1 1.9 1.8 1.4 1.1
652 655 652 647 807 721 68 1 66 1 647 664 660 653 647 642 695 671 659 652 642 705 672 655 648 642 675 663 653 647 64 1 677 666 654 648 642
652 655 652 647 653 639 637 633 63 1 625 642 646 646 643 653 645 64 1 639 635 642 640 638 637 635 637 644 647 645 642 654 649 650 648 643
-689 -669 -655 -644 -622 -628 -635 -637 -640 -660 -664 -658 -650 -640 -639 -642 -643 -645 -641 -654 -645 -64 1 -640 -639 -662 -660 -653 -647 -64 1 -640 -675 -668 -655 -642
707 690 668 654 642 645 646 645 644 680 630 635 635 633 656 648 648 649 644 683 646 644 646 642 684 669 658 648 642 58 1 624 632 633 63 1
645 630 628 628 642 645 646 645 644 617 678 664 654 642 638 646 644 646 64 1 650 662 653 643 639 628 638 638 638 637 688 701 684 66 1 646
0.6 0.4 0.4 0.3 0.9 0.6 0.5 0.4 0.4 0.3 0.3 0.3 0.3 0.4 0.5 0.5 0.4 0.4 0.4 0.5 0.5 0.5 0.5 0.4 0.5 0.4 0.4 0.5 0.4 0.6 0.5 0.4 0.4 0.4
1.2 0.9 0.9 0.9 2.0 1.4 1.2 1.1 1.2 1.2 1.1 1 .o 1.1 1.1 1.5 1.4 1.3 1.2 1.2 2.3 2.0 2.1 1.9 1.5 2.0 1.7 1.7 1.7 1.5 1.9 1.2 1.1 1.1 1.1
same as B10 -639 -646 -645 -644 -750 -697 -670 -655 -643 -633 -64 1 -64 1 -64 1 -640 -682 -663 -653 -648 -640 -686 -665 -654 -647 -640 -650 -647 -645 -642 -640 -660 -645 -642 -640 -639
'The 0-H distances were foreshortened by 10%. bElectronicunits oxygens (MCY model). Instead of using the oxygen site, they allowed a third charge site to range along the molecular twofold axis which passes through the oxygen atom. The optimum location for this third site was found to be located on the hydrogen side of the oxygen atom, at a distance of 0.2581 A. The charge at this site was found to be -1.5034e with charges of +0.7517 on the hydrogens. This site charge model, along with exponential atom-atom repulsion functions and an attractive exponential hydrogen-bonding function, was found to give a good fit to the dimer interaction energies. It should be noted that the MCY model did not include normal provision for the attractive dispersion energy. Experimental" and theoretical3*values are available for the dispersion energy which show that it is not negligible in magnitude. It is expected, therefore, that the dispersion energy of the water dimer obtained (31) G. D. Zeiss and W. J. Meath, Mol. Phys., 30, 161 (1975). (32) B. Jeziorski and M. van Hemert, Mol. Phys., 31, 713 (1976).
X lo3.
CInkJ/mol.
percent.
in the CI calculation was incorporated into the other components of this empirical potential function. For the water monomer the MCY charge site model gave a fit of 12.9% to our 6-31G** electric potential grid. If the charge of the MCY model was optimized to the 6-31G** potential points, keeping its position on the twofold axis fixed, the fit improved to 8.3% with a reduced charge of 0.6841e on the hydrogens. However, even this optimized fit was not any better than the fit with the PD charge model. We found that the fit to the electric potential at grid points around the monomer was significantly improved by allowing the charge on the molecular twofold axis to split into two charge sites with mirror symmetry above and below the molecular plane (SC model). Figure 1 shows the charge sites in this arrangement. We also anticipated that the optimum location for the electron density center of the hydrogen atom would be shifted into the bond by 10%. This electron density center shift for a bonded hydrogen atom is well d ~ c u m e n t e d and ~ ~ -leads ~ ~ to a further slight im-
Estimation of Dimer Coulombic Intermolecular Energy provement in the SC model for water monomer. The foreshortened position was more compatible with the repulsion parts of (exp 6-1) nonbonded potential functions. The fit for the SC model improved to 1.0%. The optimum location of the site charges in the SC model was found at a radial distance of 0.2523 A from the oxygen and 48.9O out of the molecular plane. The first lines of Tables I and I1 show values for the PD and SC charges and goodness-of-fit indices for the monomer. In summary, the fits to the surrounding electric potential of water monomer were found as 16.6% for the PA model, 12.9% for the MCY model, 8.3% for the charge-optimized MCY model, 8.3% for the PD model, and 1.0% for the SC model.
Dimer Results Table I shows potential-derived net atomic charges for 52 configurations of water dimer (PD model). The specification of the dimer geometries was consistent with MCY. There are five O...Odistances of 4.5, 5.0, 5.5,6.0, and 7.0 au (1 au = 0.52918 A) for each of MCY’s 10 configurations A-J. We also considered two configurations, B (bifurcated) and C (cyclic), at the maximum distance of 15 au; MCY did not consider the E (linear) model at this longer distance. The rms fit and relative rms fit are given for each case, and also the calculated SCF energies with the 6-31G** basis are shown. The latter were in general agreement with, but were expected to be less accurate than, the CI energies given by MCY. The PD charges gave a rms fit to the potential grid points of 1.5-2.9 kJ/mol or a relative rms fit of 3.2-10.4%. The rms fits did not vary much between the various configurations, although the bifurcated type B and the linear type E were fitted somewhat better. The larger range of the relative rms fit values reflects variations in the magnitudes of the electric potentials for the various configurations; if the magnitudes were small, a larger relative rms fit value was expected. Table I1 shows site charges for the S C model. As was noted above, there was a substantial improvement in the fit to the monomer electric potential. In going from the PD to the S C model, the rmsfit improved from 2.8 to 0.3 kJ/mol and the relative rms fit improved from 8.3% to 1.0%. Similar improvements in fit were obtained by using the SC model for the various dimer configurations. The rms fits ranged from 0.3 to 0.9 kJ/mol and the relative rms fits were 0.9-2.3%. Note that in the S C model we required the out-of-plane sites to have equal charge. Test calculations showed that relaxation of this site symmetry restriction gave little further improvement of fit. The uniformly good fit obtained for the various dimer configurations with the S C model was quite striking, as compared to the PD model. The dimer charges for the PD model showed site charge polarization compared to the monomer charges. Table I shows that configuration B10 has the maximum polarization of the charge on 0 2 , which has increased in magnitude from -0.786e in the monomer to -0.971 in the dimer. For 0 1 , the largest polarizations were for configurations D29 and E37, where its charge changed to -0.991 and -0.977, respectively. These results were reasonable, a -0distances. since all of these configurations had very short H As the monomer separation increased, the polarization decreased as expected. This was most apparent in configurations B18 and C27, where the monomers were farthest apart and there was little change in the PD charges of the dimer as compared to the monomer. For the SC model (Table 11), configurations B10 and D29 showed smaller site charge polarizations than in the PD model. Configuration E37 showed the largest polarization going from the monomer to the dimer, with the charge on the hydrogen involved in hydrogen bonding increasing from 0.630 to 0.807. The very short H...Odistance of 1.53 A in E37 makes this configuration strongly repulsive, as can be seen from its energy in Table
-
(33) R. F. Stewart, E. R. Davidson, and W. T. Simpson, J . Chem. Phys., 42, 3175 (1965). (34) D.E. Williams, J . Chem. Phys., 43, 4424 (1965). (35) T. L. Starr and D. E. Williams, J . Chem. Phys., 66, 2054 (1977).
-
The Journal of Physical Chemistry, Vol. 89, No. 8, 1985 1465
2ol
MCY,PD,SC
1.1 00
I
I
2 .o
2.5
I
I
3 .O
B
I
I
1.5
120
I
1
2 .o
2.6
I
I
3.0
C 100-
D
LINEAR
20
I
1.5
2 .o
I
2.5
3 .O
Figure 2. Magnitude of the Coulombic energy component for the bifurcate, cyclic, and linear geometries of water dimer (kJ/mol), plotted vs. the shortest H.e-0distance (angstroms).
I. At slightly longer, more experimentally expected, distances the S C model generally showed less polarization than the PD model. The B configutation, in particular, showed less polarization in the S C model. We found this interesting because our preliminary calculations using (exp 6-1) nonbonded functions showed that the B configuration of water dimer had an unexpectedly small hydrogen-bond energy. Kitaura and M o r o k ~ m proposed a~~ a method of decomposing the SCF energy into electrostatic, exchange, polarization, and charge-transfer components. DeWit3’ recently reported a decomposition of water dimer SCF energy into components using the same basis set and the same molecular geometry as MCY. The method of energy decomposition used was described by Smit, Derissen, and van D ~ i j n e v e l d t .Table ~ ~ I11 shows the Coulombic (electrostatic) and polarization (induction) energy components for the various configurations of water dimer as reported by DeWit (D model). The intermolecular Coulombic energy was estimated for the PD and S C net charge models by assigning site charges in the dimer equal to the unpolarized monomer PD or SC values and summing directly over the charge sites of the different molecules. When the site charges in the dimer were shifted to their polarized PD or SC values, a larger Coulombic energy was calculated. The additional energy thus obtained gave a rough estimate of the polarization energy. Note that with this definition of polarization (36) K. Kitaura and K. Morokuma, Int. J . Quantum Chem., 10, 325 (1976). (37) H. G. M. DeWit, Ph.D. Dissertation, Universitty of Utrecht, Utrecht, Netherlands, 1983. (38) P. H. Smit, J. L. Derissen, and F. B. van Duijneveldt, Mol. Phys., 37, 501 (1979).
1466 The Journal of Physical Chemistry, Vol. 89, No. 8, 198.5
Williams and Craycroft
TABLE III: Comparison of the Estimated Coulombic and Estimated Polarization Energies of the Site Charge Models with the Estimated Coulombic Energy of MCY and the Coulombic and Polarization Energy of Dewit (D)" D PD sc
structure A1 A2 A3 A4 A5 B10 B11 B12 B13 B14 B18 C19 c20 c21 c22 C23 C27 D28 D29 D30 D3 1 D32 E37 E38 E39 E40 E4 1 F42 F43 F44 F45 F46 G47 G48 G49 G50 G5 1 H52 H53 H54 H55 H56 I57 I58 I59 I60 I61 562 563 564 565 566
Coul
-. -
-/4./
-52.5 -34.1 -22.3 -10.9 -68.4 -38.9 -24.8 -17.2 -10.6 -1.14 -76.0 -42.5 -26.0 -17.2 -9.1 -0.7 same as B10 -55.1 -36.2 -23.8 -12.6 -1 14.9 -63.2 -37.7 -24.6 -12.9 -56.6 -40.5 -27.7 -19.1 -10.4 -71.1 -39.0 -24.1 -16.6 -9.6 -89.9 -49.1 -29.3 -19.1 -10.0 -61.8 -41.0 -26.8 -18.2 -9.7 -62.1 -54.1 -36.0 -23.9 -12.5
PO1 -7.6 -4.3 -2.2 -1.2 -0.4 -7.3 -3.1 -1.6 -0.9 -0.4 0.0 -8.9 -3.5 -1.6 -0.9 -0.3 -0.0 -5.8 -3.4 -1.9 -0.7 -17.7 -6.5 -3.2 -1.8 -0.7 -7.6 -4.5 -2.5 -1.4 -0.5 -8.7 -3.3 -1.5 -0.8 4.3 -10.6 -4.0 -1.8 -1 .o -0.3 -7.8 -3.7 -1.9 -1 .O -0.3 -8.3 -5.7 -3.2 -1.8 -0.6
MCY (Coul) -52.0 -48.8 -36.8 -26.3 -14.0 -29.7 -24.1 -19.5 -15.9 -10.8 -1.29 -66.0 -44.7 -3 1.O -22.2 -12.3 -0.9 -44.4 -36.6 -26.7 -15.3 -82.6 -53.1 -36.9 -26.8 -15.5 -45.9 -42.8 -33.3 -24.6 -14.0 -51.5 -34.9 -24.6 -18.2 -11.0 -72.9 -47.8 -32.6 -23.1 -12.8 -52.8 -42.9 -31.8 -23.2 -12.9 -41.9 -48.7 -37.2 -27.1 -15.4
Coul -3 1.4 -27.0 -20.0 -14.5 -8.2 -31.6 -24.6 -19.3 -15.4 -10.2 -1.15 -36.3 -25.3 -18.1 -13.4 -7.9 -0.7
POI -5.0 -1.9 0.1 0.8 0.4 -27.2 -18.0 -11.6 -8.4 -3.9 -0.1 -1.6 0.6
-32.8 -24.6 -18.1 -10.8 -48.6 -34.6 -25.1 -18.6 -11.1 -27.1 -26.3 -21.2 -16.2 -9.7 -33.3 -24.0 -18.0 -14.0 -9.0 -38.8 -26.7 -19.1 -14.2 -8.5 -3 1.3 -25.6 -19.4 -14.5 -8.5 -25.1 -30.4 -24.0 -17.9 -10.6
-34.1 -16.2 -2.3 -2.6 -18.9 -11.4 -8.6 -6.9 -2.6 -5.1 -0.9 0.2 0.8 0.2 -5.7 -3.8 -2.8 -1.2 -1.1 -3.9 -0.7 0.0 0.3 0.2 -1.5 0.0 0.6 0.6 0.4 -7.9 -5.7 -4.5 -2.5 -0.5
1.O
0.9 0.6 0.0
Coul -38.6 -34.1 -25.2 -18.1 -9.9 -3 1.8 -24.5 -19.2 -15.2 -10.0 -1.12 -45.9 -30.8 -21.4 -15.4 -8.7 -0.6 -33.8 -25.9 -19.1 -11.3 -57.7 -38.0 -26.8 -19.7 -11.6 -32.3 -29.3 -22.9 -17.2 -10.1 -40.1 -27.4 -19.7 -14.8 -9.2 -50.7 -33.2 -22.9 -16.4 -9.3 -37.7 -29.9 -22.1 -16.2 -9.2 -29.2 -33.6 -26.1 -19.3 -11.3
POI
-7.2 -3.0 -0.3 0.4 0.2 0.6 1 .o 0.4 0.4 0.1 0.0 -7.5 -3.0 -1.3 -0.7 -0.3 -0.1 -1.8 -2.1 -0.8 -0.5 -10.8 -1.9 -0.4 0.1 -0.1 -3.7 -2.3 -0.8 -0.4 -0.1 -3.7 -1.5 -0.9 -0.6 -0.2 -8.3 -2.8 -0.9 -0.5 -0.3 -5.9 -2.8 -1.3 -0.6 -0.3 -7.0 -3.6 -1.4 -0.2 -0.1
"All values in kJ/mol.
energy no correction was made for any change in the monomer energy in going from monomer-optimized charges to dimer-optimized charges. Such a correction was not possible because the site charge models are not valid for the calculation of intramolecular Coulombic energies. Nevertheless, these polarization energies could be compared to those obtained by DeWit since the latter also refer to the unpolarized monomer. Table I11 shows the Coulombic and polarization energy components obtained in this way. Examination of the values in the table shows that parallel trends exist between the components of the D model and the corresponding estimates of the MCY, PD, and S C models. Conceptually, the PD and SC models are most closely related, since both were obtained directly from the surrounding electric potential and they differ only in the location of charge sites. The S C model is the more accurate of the two since it gave a better fit to the electric potential. The MCY site charges were obtained from an empirical fit to the CI energy of water dimer, where no explicit provision was made for a dispersion energy component.
Therefore, the MCY Coulombic component may be affected by this neglect of the dispersion component. In the work of DeWit the dispersion energy was indeed evaluated and considered separately. However, the D model was conceptually different from the others, being based on a classification of the S C F integrals into types. This is apparent in Table I11 by noting large numerical differences between values shown for the D model as compared to the others. The D model generally showed the largest magnitudes for the Coulombic and polarization interaction energies. Of the site charge models, MCY gave the best agreement with the D model for the Coulombic energy components. However, no estimate of the polarization energy was available for the MCY model. The PD and SC models generally showed smaller magnitudes of estimated Coulombic and polarization energy. The S C values were closer to the D values than were the PD values. Note that the D model (as well as the MCY model) used a different basis set than the calculations for the S C and PD models.
J . Phys. Chem. 1985,89, 1467-1473
I
1.5
I
2 .o
I 2.5
I
3.0
I
Figure 3. Relative polarization (in percent) of the charge on the nearest hydrogen by the second molecule, vs. the shortest H--0distance (angstroms) in the dimer. The dashed lines show the polarization in the B (bifurcated), C (cyclic), and E (linear) models fitted with PD (net atomic) charges. The solid lines show similar results for the SC (site charge) models.
.
The largest estimated polarization energies for the PD model were a t configuration D29 (-34.7) and at configuration E37 for the SC model (-10.8). Both of these configurations have very short H.. -0distances. In some cases the estimated polarization energy obtained was positive. More of the PD values were positive than were the SC values, which was consistent with the idea that SC was a better model for estimation of polarization energy. In the SC model only the A4, A5, and E40 structures, and the first five B structures, showed apparently positive polarization energies. In general polarization should be a stabilizing energy component. The positive values obtained (which never exceeded 1.0 kl/mol) indicated that this model for the polarization energy was approximate.
1467
Figure 2 compares the magnitudes of the Coulombic component of the dimer energy for the D, MCY, PD, and SC models. The D model showed the largest Coulombic component, except for longer distances in the cyclic and linear dimers. The MCY model showed a Coulombic component closer to the D model in the linear and cyclic structures. For the bifurcated structures the MCY, PD, and SC Coulombic components were very nearly equal and were much less than the D values, especially at short distance. For the linear and cyclic structures the PD model deviated the most from the D model. Figure 3 shows a plot of the relative change in the charge of the nearest hydrogen atom as the second molecule approaches to form the dimer. The dashed curves show the PD model. The three curves shown are for the B (bifurcated), C (cyclic), and E (linear) cases. For a given He. .O distance, the B structure showed the largest site charge polarization, and the C structure the least. The linear E structure was intermediate in its polarization, but since this case was evaluated at the very short H..-O distance of 1.51 A, the polarization reached a maximum of 36%. The solid lines in Figure 3 show the site charge polarization in the SC model. The apparent polarization was smaller, reaching a maximum of 28% at the closest distance for the E structure. The B points no longer showed a larger polarization than the E points; in fact both the B and E points could be plotted on the same curve. In addition, the C points for the SC model are closer to the E and B points than they were for the PD model. The figure showed that the S C model site charge polarization had a more uniform trend vs. the H. SOdistance than the PD model. Acknowledgment. This work was supported by research grant GM16260 from the National Institutes of Health. Registry No. Water, 7732-18-5.
Simulation Study of Rotational and Translational Dynamics of Diatomic Molecules wlth Quadrupolar Interactions Vinayak N. Kabadi and W. A. Steele* Department of Chemistry, The Pennsylvania State University, University Park, Pennsylvania 16802 (Received: August 27, 1984; In Final Form: December 10, 1984)
Molecular dynamics simulationsof two model diatomic fluids have been performed over a range of temperature which includes the solid as well as the liquid at fixed density. Details of the simulation and thermodynamic properties have been reported elsewhere; here, time-correlation functions are reported. Rotational and translational velocity correlations are given as well as a number of reorientational time-correlation functions. The molecular interactions studied were based on site-site models interacting with and without quadrupolar energy. Attempts were made to fit the time correlations for the liquids to simple theoretical results. It is concluded that reorientation in these highly hindered fluids can be described as torsional oscillation with a fluctuating frequency. With the exception of the long-time tails, the translational correlation functions are reasonably well reproduced if one assumes memory functions with Gaussian time dependence.
Introduction In the preceding paper,' a computer simulation study of the thermodynamics of two diatomic systems in their melting regions was reported. In that work, the pairwise potential energies were taken to be those for the well-known sitesite model, with reduced site separation distance L* = L/o held constant at 0.547. For one system, ideal quadrupolar interactions were added, with ) ~1.00. / * Reduced reduced quadrupole moment Q* = Q / ( ~ u ~ = temperatures were varied over a rather wide range (greater than a factor of 2) with reduced density held fmed at p* = p d = 0.625. (1)
V. N. Kabadi and W. A. Steele, J . Phys. Chem., in press.
Here, t and IJ are the well-depth and distance scaling parameters in the Lennard-Jones site-site interaction function. Since the method of simulation was molecular dynamics, time-dependent quantities could also be extracted from the position and velocity data that were generated. In the present paper, some time-correlation functions for single-molecule translation and rotation are reported together with associated quantities such as correlation times and self-diffusion constants. These results have been fitted to theoretical expressions. Before we present the details, it should be pointed out that the values of the mean squared torque and force reported in paper 1 indicate clearly that both rotational and translational motions in these systems are highly hindered, for the temperatures and
0022-3654/85/2089-1467$01.50/0 0 1985 American Chemical Society
