A Systematic Appraisal of Density Functional Methodologies for

Hydrogen-Bond Acid/Base Catalysis: A Density Functional Theory Study of Protonated ... Hydrogen Bonding between Amino Acid Backbone and Side Chain ...
0 downloads 0 Views 416KB Size
Download PDF
J. Phys. Chem. 1996, 100, 4781-4789

4781

A Systematic Appraisal of Density Functional Methodologies for Hydrogen Bonding in Binary Ionic Complexes Andrew T. Pudzianowski Computer-Assisted Drug Design, Macromolecular Structure DiVision, Bristol-Myers Squibb Pharmaceutical Research Institute, Box 4000, Princeton, New Jersey 08543-4000 ReceiVed: August 23, 1995; In Final Form: October 31, 1995X

A number of density functional (DF) methodologies were systematically examined for their ability to describe strong hydrogen bonds. A set of 10 ion/molecule systems, 5 cationic and 5 anionic, was chosen on the basis of the availability of experimental and high-level ab initio results against which the DF methods could be compared. All DF models used the Lee-Yang-Parr (LYP) functional for nonlocal electron correlation, and either the 1988 Becke (B) or 1993 Becke “three-parameter” (B3) nonlocal exchange functional. Full geometry optimizations were carried out with the B-LYP combination and Gaussian basis sets, including Pople bases from 6-31G(d,p) to 6-311++G(d,p), and the TZVP and TZVP+ bases, as well as with the numerical DMol DNP and DNPP basis sets, but B3-LYP was used only with 6-311++G(d,p). Vibrational modes and thermodynamic functions were calculated only with the B-LYP/ and B3-LYP/6-311++G(d,p) models. The results support the following major conclusions. Nonlocal DF models require diffuse atomic functions to adequately describe strong H-bonding systems, especially in anionic complexes. With 6-311++G(d,p) the B-LYP and B3-LYP models are slightly inferior, but still quite similar, to MP2. B3-LYP is clearly the better DF model, yielding complexation enthalpies with root mean square deviations of 0.7 and 1.6 kcal/mol, respectively, from the MP2 and experimental values. More satisfactory agreement with experiment can be obtained in individual cases by judicious enlargement of the basis set. Efficient auxiliary basis/model density DF methods are much faster than ex post facto DF schemes based on existing ab initio methodology, which in turn are appreciably faster than MP2 with the same basis set. DF models of good quality will extend to strong H-bonded systems large enough to be of direct bioorganic relevance.

Introduction Hydrogen bonding of an ion with a neutral molecule entails a complexation enthalpy in the range of -15 to -50 kcal/mol, roughly an order of magnitude larger than for an H-bonded complex between neutral molecules. The enhanced energy content is associated with appreciable covalent character, qualities which define the “strong hydrogen bond”.1 Ion solvation has been the historical impetus for a large body of work on the strong hydrogen bond, and certainly remains a major factor. (The reviews by Meot-Ner2 and Deakyne3 survey, respectively, the experimental and theoretical areas through 1987.) More recently the focus has been shifting to biological molecules and endogenous or synthetic ligands. The role of ionic hydrogen bonding in the interactions of amino acid side chains has prompted recent experimental4,5 and theoretical6 work on gas-phase clusters. The importance of strong, “low barrier” hydrogen bonds in enzymatic catalysis is emerging7 and will certainly drive further work in the near future. Implicit in this shift is the need to deal on a fundamental level with molecular systems large enough to be relevant. Unfortunately, hydrogen bonding in general is a difficult challenge for conventional theoretical treatments. Recent ab initio studies8-10 have shown that optimized geometries can vary widely at theoretical levels below MP2/6-31+G(d,p). (The notation specifying the basis set and correlation levels is standard.11) Reliably accurate complexation energies require significantly higher levels of theory. Del Bene has established9 what is probably the most economical computational protocol consistent with high quality. It calls for geometry optimization and evaluation of zero-point energies and thermodynamic X

Abstract published in AdVance ACS Abstracts, March 1, 1996.

0022-3654/96/20100-4781$12.00/0

contributions at MP2/6-31+G(d,p), followed by computation of energies and other static properties at the MP4/cc-pVTZ+ level, which uses the Dunning correlation-consistent valence triple split basis12 augmented with s- and p-type diffuse functions13 on all non-hydrogen atomic centers. Clearly, the demands on computational resources quickly become forbidding, and even with steady, incremental advances in computer capacity it may be decades before accurate ab initio models can be applied to H-bonded complexes of direct biological import. Given this bleak prospect, the advent of density functional (DF) methods14-18 in chemistry has been timely. “Native” DF methods (see below) are faster than conventional RoothaanHall (RH) ab initio procedures at the Hartree-Fock (HF) level, and much faster than correlated RH calculations.19 The local density approximation (LDA) of DF theory includes correlation but overestimates the absolute values of binding energies.18 This effectively limits the LDA to isolated molecules, making it unsuitable for the study of complexes and thermochemistry in general. Fortunately, the LDA can be economically corrected with nonlocal exchange and correlation functionals of the local density and its gradient.20 The Kohn-Sham (KS) formalism15 of DF theory builds its density from one-electron wave functions, so several linearcombination-of-atomic-orbitals (LCAO) methodologies have evolved. Among those using analytical functions the Gaussian basis schemes are most common and further divide into “native” DF auxiliary basis/model density methods,21 and those which adapt existing RH methods ex post facto to do DF calculations.22 The nonanalytical DMol approach23 numerically solves the KS equations for atoms and atomic ions to implement basis sets as three-dimensional grids. © 1996 American Chemical Society

4782 J. Phys. Chem., Vol. 100, No. 12, 1996 As DF methodologies are just reaching an initial plateau of maturity, many questions regarding applicability and performance remain to be answered. The study of hydrogen-bonding problems in particular became possible no more than 5 years ago with conclusive demonstrations that the LDA was unworkable and that nonlocal corrections were effective.21b,24 DF studies on H-bonding published since then are few, although the pace is now increasing.25 Even so, only a small proportion of these results deals directly with strong, ionic hydrogen bonds. In the most comprehensive of the recent studies, Gresh et al.25f included five binary ion/molecule systems. Their results, however, should be considered preliminary because the geometries of complexes were not fully optimized, the basis set employed was double-split and did not include diffuse functions, and neither zero-point energies nor thermodynamic functions were evaluated, so that complexation energies at 0 K were directly compared to experimental complexation enthalpies. A more recent study by Del Bene et al.25g did feature full geometry optimization and thermodynamically corrected energy comparisons, albeit for a series of neutral binary complexes, using one of the DF methods considered here but only with the 6-31+G(d,p) basis. It still remains to be seen whether the quality of high-level ab initio models can be consistently approached by comparably extensive DF methods. Furthermore, even if that proves to be the case it must also be demonstrated that DF methods will deliver such results significantly faster than RH methods can, an obvious prerequisite for calculations on molecular systems of chemically meaningful size. The study reported here was intended to at least partly fill this gap by systematically examining the performance of some of the DF methods mentioned above. The following conditions were deemed important for a useful evaluation: (a) The set of molecular systems studied should be large enough to reveal trends and include both cationic and anionic complexes. (b) Experimental gas-phase complexation enthalpies should be available. (c) RH ab initio calculations at a uniform, high level should be feasible so as to provide optimized geometry, complexation energy, and other direct performance benchmarks for the DF methodologies. (d) At least some of the chosen molecular systems should pose challenges, either as established, difficult computational problems or by presenting choices of complexation schemes for which the relevant experimental results are ambiguous. The following set of 10 ion/molecule systems met the above requirements: H3O+/H2O, NH4+/H2O, CH3NH3+/H2O, NH4+/ NH3, CH3NH3+/NH3, OH-/H2O, CH3O-/H2O, HCOO- (formate)/H2O, CN-/H2O, and HCC- (acetylide)/H2O. The ab initio stage of the study, described in detail elsewhere,10 involved full geometry optimizations and vibrational frequency calculations at the MP2/6-311++G(d,p) level. The DF phase, whose early stages have been described,26 tested a number of Gaussian basis sets up to 6-311++G(d,p), including the TZVP basis designed for DF calculations by Godbout et al.,27 as well as two numerical DMol basis sets.23 Nonlocal corrections included the Lee-Yang-Parr (LYP) correlation functional,28 and the 1988 Becke29a (B) and 1993 Becke29b “three-parameter” (B3) exchange functionals. The more general conclusions include the following. Diffuse atomic functions are absolutely necessary for acceptable results with anionic H-bonded systems, and also appreciably improve results for cationic systems. With 6-311++G(d,p) the B-LYP and B3-LYP models yield results comparable to the MP2 benchmarks in nearly all respects but B3-LYP is noticeably better, agreeing closely with MP2 in cases where B-LYP results are ambivalent and would make a clear-cut interpretation of

Pudzianowski experimental results difficult. Speed comparisons employing the B-LYP/6-311++G(d,p) model show that the auxiliary basis/ model density approach is much faster than an RH method adapted for ex post facto DF calculation, which in turn is faster than MP2. The study supports the overall conclusion that current DF methods can perform well in the demanding area of ionic hydrogen bonding and will extend to systems large enough to be useful in bioorganic applications. Computational Methods Pople basis sets included 6-31+G(d,p), 6-31++G(d,p), and 6-311++G(d,p), derived from the standard, valence doublesplit 6-31G(d,p)30a and triple-split 6-311G(d,p)30b basis sets by addition of a single shell of diffuse s- and p-type functions13 to each non-hydrogen atomic center (+) and an s-type diffuse function to each hydrogen center (++). The “decontracted” TZVP basis, (7/111/411/1*) on C, N, O and (41/1*) on H centers,27 was singly augmented13 for the TZVP+ basis used here. The DMol basis sets were “double-numerical” (DN),23 using KS orbitals for neutral atoms and their +2 atomic ions. Polarization functions were KS orbitals for a one-electron atom related by atomic number Z to a given atomic center: DNP had 3d(Z-1) functions on all C, N, O and 2p(Z)1.3) functions on all H centers, and DNPP had additional 3d(Z+1) functions on C, N, O and 2p(Z)4.0) functions on H centers. Numerical diffuse functions31 became available as most of the DMol work described here was being completed and were not included. The B-LYP functionals were used with all basis sets, but B3LYP was used only with 6-311++G(d,p). Nonlocal corrections were self-consistent so that optimized geometries for H-bonded complexes would not be compromised by the LDA, as in the “perturbative” scheme.24 Local correlation was supplied by the Vosko-Wilk-Nusair (VWN) functional32 for Gaussian basis methods, and by the Hedin-Lundqvist/Janak-Moruzzi-Williams functionals33 for DMol calculations. The Gaussian 92/DFT system34 was used for most of the Gaussian basis set work. Direct SCF was used with the “Tight” convergence criteria. Numerical DF integrations used a modified (75 302) grid, i.e., one with 75 radial shells per atom and 302 points per shell, “pruned” to leave ∼7000 points per atom. A few cases needed a larger (99 302) full grid. Geometry optimizations used the Berny method,35 initially in Cartesian coordinates without symmetry constraints and with default convergence criteria, followed (when relevant) by imposition of symmetry in internal coordinates under “Tight” criteria. Most optimizations started from ab initio geometries,10 and initial Hessian matrices were usually calculated at the HF level, although semiempirical AM136 Hessians also worked well. Harmonic vibrational frequencies were calculated only for the B-LYP/ and B3-LYP/6-311++G(d,p) models. Basis set superposition effects (BSSE) could not be examined here because Gaussian 92/DFT does not evaluate DF quantities with the “ghost” functions of the counterpoise method.37 We note, however, that the counterpoise energies in the companion ab initio study10 averaged about 10% of the complexation ∆E0 with MP2 as opposed to ∼3% at the HF level. A similar relation exists between MP2 and a nonlocal DF method in the results reported by Gresh et al.;25f hence, we anticipate that DF counterpoise energies for the 6-311++G(d,p) basis used here will be similar to those found at the HF level.10 The DGauss 2.3 program system38 was also used for B-LYP calculations with 6-311++G(d,p) and TZVP. The A2 auxiliary basis designed27 for TZVP was used with both atomic orbital (AO) basis sets. High accuracy was specified for the evaluation of integrals and exchange-correlation functionals (INTACC )

Hydrogen Bonding in Binary Ionic Complexes

J. Phys. Chem., Vol. 100, No. 12, 1996 4783

high, XCGRID ) fine). The SCF tolerances were 1.0 × 10-7 hartree for energies and 1.0 × 10-5 electron/bohr3 for densities. The geometry optimization tolerance was 5.0 × 10-4 hartree/ bohr for the largest gradient component. (See Table 1 for definitions of the units.) Calculations with the DMol 2.3.5 program system23 used the XFINE numerical mesh, and multipolar functions to fit the model density were given by LMAX values of 3 and 4, respectively, for H and C, N, O atomic centers. The SCF tolerance was 1.0 × 10-7 hartree, and the geometry optimization tolerance was 1.0 × 10-4 hartree/bohr on the largest gradient component. Atomic calculations for the AO bases included the B-LYP nonlocal corrections. The SYBYL39 and SPARTAN40 systems were the primary graphics interfaces. Normal vibrational modes in particular were analyzed with SPARTAN. Results and Discussion The various results are integrated with discussions in the following way. Part A reports the energies obtained for all optimized species, and complexation energies and geometrical results that pertain directly to hydrogen bonding. These are compared to MP2/6-311++G(d,p) results,10 and the various DF approaches are evaluated and discussed. Part B focuses on B-LYP/ and B3-LYP/6-311++G(d,p), comparing complexation thermochemistry results with MP2 and experimental values. Part C extends this detailed comparison to calculated intermolecular normal vibrational modes in the complexes. Part D examines systems that challenge computational approaches or experimental interpretations. Finally, part E compares a “native” DF and an ex post facto numerical approach. A. Complexation Energies and Optimized Geometries. Table 1 gives the complexation energy ∆E0 and four geometrical parameters for each complex, along with the corresponding MP2/6-311++G(d,p) results.10 A hydrogen bond between an acceptor center A and a donor HD is denoted A- -H-D and described by the structural parameters A- -H, the acceptorhydrogen distance; H-D, the convalent bond distance in the donor; A- -D, the distance between the acceptor and donor atomic centers; and Θ, the angle defined by the H-bond centers AHD. Optimized structures determined with B3-LYP/6311++G(d,p) are pictured in Figure 1. (Full geometries for all DF models considered here are available from the author on request.) Tables S1 and S2 report the total energy E0 for every optimized species, along with covalent donor H-D bond distances in isolated donors. Also given there are the zeropoint energies (Ezp), thermal contributions (Etherm), and entropies obtained with B-LYP/ and B3-LYP/6-311++G(d,p). (See the paragraph at the end of the paper regarding the availability of supporting information.) Table 1 contains all the information needed for a broad appraisal of the DF methodologies examined here. To make this easier, root mean square (rms) deviations of DF results from the ab initio values have been determined. For complexation energies these are ∆(∆E0)rms ) [∆E0(MP2) - ∆E0(DF)]rms. Optimized geometries of the complexes are analyzed in terms of the quantities ∆Xrms ) [X(MP2) - X(DF)]rms, where X is one of the structural parameters defined above. All rms results are assembled in Table 2. The results in Tables 1 and 2 support a number of conclusions. First, diffuse functions are essential for reliable descriptions of strong H-bonding by DF methods. This parallels the results obtained by Del Bene et al. with B3-LYP for neutral H-bonded complexes,25g as well as Del Bene’s ab initio MP2 results for charged complexes.41 Consider the B-LYP/6-31G(d,p) results

Figure 1. B3-LYP/6-311++G(d,p) optimized geometries for cationic and anionic complexes. Final point group symmetries are indicated.

(Table 1) for H3O+/H2O and OH-/H2O, with overly negative ∆E0 values. Augmentation to B-LYP/6-31+G(d,p) improves these energies dramatically, bringing them to within ∼2 kcal/ mol of the ab initio values. Further augmentation to B-LYP/ 6-31++G(d,p) has little effect on ∆E0, a trend evident in the rms values for all complexes in Table 2. This minimal effect of H-centered s-type diffuse functions parallels Del Bene’s results41 for MP2 with augmented basis sets. Optimized geometries also change little except in the two strongest anionic

4784 J. Phys. Chem., Vol. 100, No. 12, 1996

Pudzianowski

TABLE 1: ∆E0 and Optimized H-Bond Geometries complex H3O+/H2O

NH4+/H2O

CH3NH3+/H2O

NH4+/NH3

CH3NH3+/NH3

OH-/H2O

CH3O-/H2O

CN-/H2O

NC-/H2O

∆E0b

H-D, Å

A- -H, Å

A- -D, Å

Θ, deg

MP2c B3 B 6-31++ 6-31+ 6-31 TZVP+ TZVP DNPP DNP MP2 B3 B 6-31++ 6-31+ TZVP+ TZVP DNPP MP2 B3 B 6-31++ 6-31+ TZVP+ TZVP DNPP MP2 B3 B 6-31++ 6-31+ TZVP+ TZVP DNPP MP2 B3 B 6-31++ 6-31+ TZVP+ TZVP DNPP

-35.5 -36.7 -36.5 -37.3 -37.4 -42.4 -36.2 -37.3 -37.5 -39.6 -21.6 -22.0 -21.5 -21.6 -21.6 -21.2 -21.8 -21.8 -19.7 -19.5 -18.9 -18.9 -19.0 -18.6 -19.1 -18.9 -27.6 -28.2 -28.5 -29.6 -29.6 -28.5 -29.2 -29.9 -24.7 -24.3 -24.2 -25.1 -25.1 -24.2 -24.8 -25.3

(a) Cationic Complexes 1.192 1.198 1.212 1.212 1.212 1.215 1.218 1.217 1.202 1.214 1.055 1.064 1.075 1.078 1.078 1.077 1.079 1.064 1.046 1.051 1.062 1.064 1.064 1.063 1.065 1.049 1.120 1.144 1.181 1.202 1.201 1.189 1.193 1.184 1.085 1.093 1.112 1.120 1.120 1.115 1.116 1.105

1.192 1.198 1.212 1.213 1.214 1.215 1.218 1.218 1.204 1.214 1.649 1.637 1.640 1.625 1.625 1.633 1.626 1.599 1.705 1.703 1.709 1.700 1.700 1.705 1.697 1.670 1.574 1.548 1.510 1.470 1.472 1.501 1.493 1.458 1.677 1.683 1.673 1.646 1.646 1.665 1.661 1.621

2.382 2.393 2.420 2.424 2.423 2.426 2.431 2.431 2.402 2.423 2.704 2.701 2.715 2.704 2.704 2.710 2.705 2.663 2.748 2.754 2.771 2.764 2.764 2.768 2.762 2.719 2.694 2.692 2.691 2.673 2.674 2.689 2.686 2.642 2.762 2.776 2.785 2.766 2.766 2.780 2.777 2.726

174.2 175.1 174.2 174.1 174.0 173.4 173.5 173.4 173.3 172.6 179.6 179.8 179.9 179.8 179.8 179.8 180.0 179.9 175.2 178.5 178.9 178.5 178.5 178.5 178.3 178.1 180.0 180.0 180.0 180.0 179.8 180.0 180.0 180.0 180.0 180.0 180.0 179.6 179.7 179.5 179.4 179.1

MP2 B3 B 6-31++ 6-31+ 6-31 TZVP+ TZVP DNPP DNP MP2 B3 B 6-31++ 6-31+ TZVP+ TZVP DNPP MP2 B3 B 6-31++ 6-31+ TZVP+ TZVP DNPP MP2 B3 B 6-31++ 6-31+ TZVP+ TZVP DNPP

-28.3 -29.4 -29.5 -30.3 -30.1 -48.2 -29.6 -36.4 -31.6 -35.0 -26.0 -25.1 -24.3 -25.2 -25.9 -25.7 -27.9 -28.0 -15.7 -16.4 -16.1 -16.2 -16.3 -16.0 -18.0 -22.8 -15.7 -16.3 -16.3 -16.9 -16.9 -16.4 -19.6 -23.7

(b) Anionic Complexes 1.095 1.117 1.160 1.214 1.240 1.252 1.242 1.245 1.196 1.229 1.064 1.062 1.082 1.102 1.109 1.113 1.111 1.095 0.987 0.996 1.010 1.015 1.015 1.015 1.019 1.009 0.992 1.001 1.017 1.024 1.023 1.022 1.030 1.016

1.380 1.363 1.335 1.270 1.243 1.253 1.244 1.245 1.252 1.230 1.426 1.458 1.460 1.420 1.407 1.405 1.411 1.385 1.833 1.802 1.799 1.767 1.769 1.784 1.758 1.679 1.943 1.917 1.897 1.863 1.866 1.881 1.830 1.782

2.474 2.479 2.495 2.483 2.483 2.501 2.486 2.489 2.448 2.459 2.488 2.518 2.539 2.521 2.515 2.517 2.521 2.480 2.807 2.788 2.800 2.776 2.778 2.792 2.772 2.686 2.924 2.908 2.906 2.881 2.883 2.897 2.856 2.796

176.0 176.6 177.2 178.2 178.6 173.1 178.5 178.2 179.4 179.4 175.4 174.8 175.3 176.0 176.5 176.9 176.2 177.5 168.4 169.8 170.7 171.9 171.7 171.5 172.6 176.1 169.3 170.2 171.3 172.7 172.5 172.4 173.2 175.3

modela

Hydrogen Bonding in Binary Ionic Complexes

J. Phys. Chem., Vol. 100, No. 12, 1996 4785

TABLE 1: (Continued) modela

complex HCC-/H2O

HCOO-/H2O (C1)

HCOO-/H2O (C2V)

MP2 B3 B 6-31++ 6-31+ TZVP+ TZVP DNPP MP2 B3 B 6-31++ 6-31+ TZVP+ TZVP DNPP MP2 B3 B 6-31++ 6-31+ TZVP+ TZVP DNPP

∆E0b -18.4 -19.0 -19.1 -19.4 -19.4 -19.2 -21.2 -22.4 -16.9 -16.8 -16.3 -16.8 -16.8 -16.7 -17.8 -19.0 -19.7 -19.2 -18.4 -18.9 -18.9 -18.7 -20.2 -20.9

H-D, Å (b) Anionic Complexes 1.005 1.015 1.033 1.040 1.040 1.041 1.047 1.033 0.996 1.002 1.016 1.024 1.024 1.025 1.026 1.015 0.972 0.976 0.987 0.991 0.991 0.991 0.991 0.977

A- -H, Å

A- -D, Å

Θ, deg

1.874 1.849 1.824 1.794 1.798 1.804 1.780 1.734 1.675 1.680 1.688 1.660 1.659 1.656 1.651 1.588 2.022 2.027 2.045 2.040 2.040 2.040 2.018 1.940

2.869 2.856 2.851 2.831 2.833 2.841 2.823 2.766 2.663 2.673 2.696 2.678 2.678 2.676 2.673 2.600 2.868 2.869 2.899 2.900 2.899 2.895 2.880 2.792

170.0 171.4 172.4 173.7 173.3 173.5 173.7 175.8 170.8 170.2 171.2 172.5 172.5 173.0 173.3 174.8 144.2 143.2 143.6 143.6 143.6 143.3 144.1 144.6

a The first three rows report results for the 6-311++G(d,p) basis with MP2, B3-LYP (B3), and B-LYP (B); subsequent rows give B-LYP results for various basis sets, where Pople basis designations omit G(d,p) for brevity. b Complexation energy reported in kcal/mol. The units used in this study are related to SI units as follows: 1 cal ) 4.184 J; 1 hartree ) 627.51 kcal/mol; 1 Å ) 1.0 × 10-10 m; 1 bohr ) 0.52918 Å. c All MP2 results from ref 10.

TABLE 2: Root Mean Square Deviations of DF Results Relative to MP2/6-311++G(d,p) ∆(AD)rms, Å

∆Θrms, deg

(a) Cationic Systems (n ) 5)a 0.012 0.013 0.033 0.030 0.043 0.051 0.042 0.050 0.038 0.036 0.040 0.040 0.031 0.064

0.009 0.023 0.022 0.022 0.025 0.024 0.037

1.5 1.7 1.5 1.5 1.5 1.5 1.4

0.71 1.03 1.07 0.98 0.76 3.75 4.65

(b) Anionic Systems (n ) 7) 0.011 0.023 0.032 0.037 0.053 0.065 0.062 0.071 0.064 0.066 0.066 0.081 0.044 0.120

0.016 0.029 0.031 0.029 0.023 0.037 0.086

1.0 1.6 2.6 2.5 3.0 2.9 4.8

0.69 0.92 1.16 1.12 0.76 2.96 3.67

(c) Overall (n ) 12) 0.012 0.019 0.032 0.034 0.049 0.060 0.055 0.063 0.054 0.056 0.056 0.067 0.039 0.100

0.014 0.026 0.028 0.026 0.024 0.032 0.070

1.2 1.6 2.2 2.1 2.5 2.4 3.8

DF model

∆(∆E0)rmsb

B3 B 6-31++ 6-31+ TZVP+ TZVP DNPP

0.66 0.74 1.27 1.28 0.76 1.11 1.44

B3 B 6-31++ 6-31+ TZVP+ TZVP DNPP B3 B 6-31++ 6-31+ TZVP+ TZVP DNPP a

∆(HD)rms, Å

∆(AH)rms, Å

n is the number of complexes. b Deviations in kcal/mol.

complexes, OH-/H2O and CH3O-/H2O (Table 1), where the shared hydrogen’s position shifts, as seen in A- -H and H-D, while A- -D is nearly unaltered. The TZVP+ results in Tables 1 and 2 also show improvements relative to TZVP. Diffuse functions have a greater effect on anionic complexes, which is expected given the importance of augmented basis sets in RH methodology for molecular anions.11 In fact, Table 2 shows that for cationic systems all of the DF methods examined, with the exception of B-LYP/6-31G(d,p) and /DNP, are at least marginally acceptable. With anionic systems, however, B-LYP/ TZVP and /DNPP are too inaccurate for dependable work.

Addition of a second set of polarization functions (DNP to DNPP) has a marked effect on DMol results for H3O+/H2O and OH-/H2O: ∆E0 is improved by 2.1 and 3.4 kcal/mol, respectively, and H-bond distances are shortened by ∼0.02 Å. Overall, B-LYP/DNPP gives the shortest A- -H and A- -D distances, and the rms deviations reported in Table 2 reflect this systematic trend. In geometrical terms B-LYP/DNPP resembles the double-split models in Table 2 for cations, but stands apart in its poor performance with anionic complexes. Comparing 6-31++G(d,p) and 6-311++G(d,p) (Table 2), we find that expanding from double- to triple-split yields minor overall improvement in ∆E0 but more significant improvement in H-bond geometry. Again, most of the latter stems from the shared hydrogen’s position in OH-/H2O, which in the B-LYP/ 6-311++G(d,p) model resembles the MP2 results more closely. Note, however, that the overall B-LYP/TZVP+ results, which also involve a triple-split basis, are similar to the double-split B-LYP/6-31+G(d,p) model in geometries but are significantly better in terms of ∆E0 (Table 2). These comparisons serve to emphasize that DF results are quite dependent on basis set details, as is the case with RH methods. Overall, B3-LYP/ and B-LYP/6-311++G(d,p) emerge as the two best DF models relative to MP2/6-311++G(d,p), with a tangible edge to B3-LYP. Hence, we now examine the performance of these two models in more depth. B. Thermochemistry of Complexation for B-LYP and B3LYP Models. All thermodynamic quantities reported here refer to the complexation equilibrium A + HD h A- -HD at 298.15 K and 1 atm of pressure. The complexation energy ∆E0 is adjusted to these standard conditions by addition of the zeropoint energy difference ∆Ezp and the thermal energy term ∆Etherm (including vibrational, rotational, and translational contributions) for complexation, derived from the values reported in Tables S1 and S2 for individual species, and the adjusted energy is symbolized by ∆E°. The complexation enthalpy ∆H° is then evaluated under the ideal gas assumption as ∆H° ) ∆E° - RT, where R is the ideal gas constant and T is the absolute temperature. Table 3 gives ∆H°, ∆S°, and ∆G°

4786 J. Phys. Chem., Vol. 100, No. 12, 1996

Pudzianowski

TABLE 3: Gas-Phase Thermochemistry at 298.15 K for 6-311++G(d,p) Methods complex H3O+/H2O

∆H°(exp)a

model

-35.0c

MP2b

B3-LYP B-LYP NH4+/H2O MP2 B3-LYP B-LYP CH3NH3+/H2O MP2 B3-LYP B-LYP NH4+/NH3 MP2 B3-LYP B-LYP CH3NH3+/NH3 MP2 B3-LYP B-LYP OH-/H2O MP2 B3-LYP B-LYP CH3O-/H2O MP2 B3-LYP B-LYP CN-/H2O MP2 B3-LYP B-LYP NC-/H2O MP2 B3-LYP B-LYP HCC-/H2O MP2 B3-LYP B-LYP HCOO-/H2O (C1) MP2 B3-LYP B-LYP HCOO-/H2O (C2V) MP2 B3-LYP B-LYP

-20.6d -18.4d -24.8d -21.4d -27.0e -23.9e

(-14.6)f (-14.6) -16.2f

(-16.0)e (-16.0)

∆H° -35.8 -37.0 -36.6 -20.4 -20.8 -20.3 -18.3 -18.1 -17.6 -27.1 -27.8 -28.5 -23.6 -23.2 -23.3 -28.6 -29.9 -30.4 -24.7 -23.6 -22.7 -14.7 -15.2 -15.0 -14.9 -15.3 -15.4 -17.4 -18.1 -18.4 -15.2 -15.2 -14.8 -17.8 -17.3 -16.6

∆S° -32.55 -30.54 -32.42 -24.04 -24.01 -24.90 -24.93 -23.72 -25.88 -27.83 -29.63 -29.62 -26.34 -26.81 -26.96 -26.11 -24.80 -25.99 -27.93 -26.54 -27.73 -23.26 -23.24 -24.54 -24.25 -24.01 -25.36 -23.89 -24.63 -25.68 -26.70 -25.12 -26.50 -31.31 -30.66 -31.66

-26.1 -27.9 -26.9 -13.2 -13.6 -21.5 -10.9 -11.0 -9.9 -18.8 -19.0 -19.7 -15.8 -15.2 -15.3 -20.8 -22.5 -22.6 -16.4 -15.7 -14.4 -7.8 -8.3 -7.7 -7.7 -8.1 -16.3 -10.3 -10.8 -10.7 -7.2 -7.7 -6.9 -8.5 -8.2 -7.2

TABLE 4: Root Mean Square Deviations for B3-LYP and B-LYP ∆H° Relative to MP2 and Experimental Values ∆(∆H°)rmsb

∆(∆H°)rms(exp)

(a) Cationic Systems (n ) 5) MP2 B3-LYP B-LYP

0.68 0.80

1.47 1.81 2.03

(b) Anionic Systems (n ) 6) MP2 B3-LYP B-LYP

0.76 1.21

0.95 1.51 1.83

(c) Overall (n ) 11) MP2 B3-LYP B-LYP a

0.73 1.06

model

∆ωrms(cations)a (n ) 31)b

∆ωrms(anions) (n ) 40)

∆ωrms(overall) (n ) 71)

B3-LYP B-LYP

30 27

24 21

27 24

∆G°

a Energies reported in kcal/mol, entropies in cal/(mol‚K). b MP2 results from ref 10. c Hiraoka et al. (ref 42). d Meot-Ner (ref 43). e MeotNer and Sieck (ref 44). f Meot-Ner (ref 45).

modela

TABLE 5: Root Mean Square Deviations of DF Frequencies from MP2 Values

1.21 1.65 1.92

All models used 6-311++G(d,p). b Deviations in kcal/mol.

for the B-LYP/, B3-LYP/, and MP2/6-311++G(d,p) models, along with experimental ∆H°. Table 4 reports rms values for deviations ∆(∆H°) of B-LYP and B3-LYP ∆H° values from MP2 results, as well as ∆(∆H°)rms relative to experimental values. The overall results for ∆(∆H°)rms in Table 4 show that B3LYP is ∼0.4 kcal/mol closer to MP2 and ∼0.3 kcal/mol closer to experimental ∆H° than is B-LYP. These differences can be important in the context of “chemical accuracy” for theoretical models, usually accepted as (1 kcal/mol or about the same as experimental uncertainties in gas-phase ∆H°. The MP2/6311++G(d,p) model is nearly of that quality, with an overall

a

Deviations in cm-1. b n is the number of frequencies.

∆(∆H°)rms of 1.2 kcal/mol relative to experimental values, but small improvements can be costly.10 Its 1.6 kcal/mol figure makes B3-LYP more interesting than B-LYP (1.9 kcal/mol), and also reflects chemical details that B3-LYP is better able to describe. (See Part D below.) Both DF models are somewhat better with the anionic systems because both track the MP2 results rather closely for cations (Table 4), but the MP2 model itself is worse for NH4+/NH3 and CH3NH3+/NH3 than for any other complexes examined.10 With the anionic complexes B3-LYP follows MP2 more closely (Table 4) and MP2 comes within 1.0 kcal/mol rms of experimental ∆H° values. C. Intermolecular Vibrational Modes and Frequencies. Intermolecular modes are those whose nuclear displacements extend over both partners of the complex. (In exceptionally strong, symmetrical complexes like H3O+/H2O where there are no discrete partners, every mode is “intermolecular”, further attesting to the significant covalent character of the bonding.) Table S3 lists the B-LYP/ and B3-LYP/6-311++G(d,p) results for these modes alongside the MP2 results.10 Table 5 gives rms values for the deviations ∆ω of DF frequencies from MP2 values. The two DF models are similar in overall performance and both are slightly better with anionic complexes (Table 5). The overall ∆ωrms values for B-LYP and B3-LYP are, respectively, 24 and 27 cm-1, with about as many frequencies underestimated as overestimated (Table S3). The largest deviations occur in H3O+/H2O and HCC-/H2O, and B3-LYP is somewhat worse in the latter case. An important observation in the MP2 results was a hydrogen bond stretch-compression mode, roughly in the 200-350 cm-1 range, common to all the complexes studied.10 Figure 2 shows the nuclear displacements in the CH3O-/H2O complex as given by B3-LYP/6-311++G(d,p). This mode may play a role in proton transfer, and it is gratifying that the DF models retain this feature (Table S3). Also retained, with the exception of the B3-LYP result for H3O+/H2O, are qualitative relationships noted for the MP2 results between the frequency of this mode and the acid strength of the H-bond donor or the base strength of the acceptor.10 D. Individual Systems. Detailed comparisons and discussions are limited to those systems for which the B-LYP/ and B3-LYP/6-311++G(d,p) results differ appreciably from ab initio models, or which have posed challenges for ab initio treatments or interpretation of experimental results. 1. H3O+/H2O. This system has recently been the subject of high-level ab initio studies8,10,46 and is rather well characterized. The most extensive calculations8,46,47 have essentially converged on ∆H0 ) -35 kcal/mol, the experimental value cited here,42 and on geometrical parameters almost identical to the ab initio values in Table 1. Hence, we are in an excellent position to evaluate the DF results for this system. Neither DF model gives ∆H° to “chemical accuracy”. B-LYP comes 0.4 kcal/mol closer than B3-LYP to the experimental ∆H° (Table 3), but the latter gives H-bond distances substantially closer to the MP2 values (Table 1). There is room for improvement.

Hydrogen Bonding in Binary Ionic Complexes

Figure 2. (a) H-bond vibrational mode in compression direction for CH3O-/H2O complex: B3-LYP/6-311++G(d,p) ω ) 343 cm-1. (b) HCOO-/H2O C2V complex B2 mode: B3-LYP-6-311++G(d,p) ω ) 95 cm-1.

In ab initio work with MP2 and basis sets derived from 6-311++G, expansion of the polarization set from (d,p) to (2d,2p) improved ∆H0 slightly8,10 but the (d,p) result was already acceptable. Still closer agreement with the experimental value required MP4 correlation and extension of polarization to (3df,3pd).8 It may be that larger basis sets could improve the DF models somewhat, perhaps enough for “chemical accuracy”. Some support for this conjecture will be found below. 2. NH4+/H2O and CH3NH3+/H2O. Tables 1 and 3 show that B3-LYP/6-311++G(d,p) matched the MP2 and experimental ∆H° closely for these two systems, but B-LYP was a bit off the mark with CH3NH3+/H2O. Geometry optimizations on both complexes were complicated by the lowest-frequency normal mode (Table S3), which had a negative force constant when “Finegrid” was used. Reoptimization with a (99 302) grid (∼30 000 points per atom) gave positive force constants. Problems attend the low curvature of these energy surfaces near the complexes. For example, B3-LYP with the (99 302) grid gave three conformers of NH4+/H2O (Figure 1 shows an “eclipsed” structure), with gradient norms of 2 × 10-5 hartree/ bohr or less, no negative force constants, and energies less than 1 × 10-7 hartree apart. These conditions demand high numerical accuracy. 3. NH4+/NH3 and CH3NH3+/NH3. The two best DF models produced mixed results relative to MP2. For NH4+/NH3 the B3-LYP ∆H° was 0.7 kcal/mol more negative than MP2 and a full 3.0 kcal/mol lower than the experimental value, but B-LYP

J. Phys. Chem., Vol. 100, No. 12, 1996 4787 was 1.4 kcal/mol lower than MP2 and 3.7 kcal/mol lower than experiment (Table 3). For CH3NH3+/NH3, however, both DF values were nearly equal, ∼0.3 kcal/mol higher than the MP2 result and ∼1.9 kcal/mol lower than experiment. In both cases B3-LYP gave somewhat better H-bond geometries (Table 1). These systems were the most difficult for MP2/6-311++G(d,p).10 ∆H° was ∼2 kcal/mol more negative than experimental values in both cases. Expansion of polarization to (2d,2p) brought ∆H° within 1 kcal/mol of experimental values, and MP4 correlation yielded a further ∼0.5 kcal/mol. Again, it may be that larger basis sets would confer “chemical accuracy” on the B3-LYP model in these cases. 4. OH-/H2O. Here B3-LYP is clearly preferable, giving ∆H° within 1.3 kcal/mol of the MP2 value (but 2.9 kcal/mol lower than experiment) and good agreement with the MP2 geometry (Table 1). Again, however, the molecular system is a challenge for MP2/6-311++G(d,p),10 which yields a ∆H° 1.6 kcal/mol more negative than the experimental value (Table 3). The simple expedients described above do not improve ∆H° to an acceptable extent; in fact, the level of theory must be raised to CCSD(T)/aug′-cc-pVTZ before close agreement with experiment is attained.48 (See ref 10 and references cited therein for an overview of recent ab initio work on this difficult system.) There is some direct evidence in this particular case that more flexible basis sets may yield acceptable improvements in the performance of current DF models. The present study included some experiments with B-LYP in which the 6-311+G(d,p) and TZVP+ basis sets were further augmented with s- and p-type diffuse functions on all hydrogen centers. With geometry optimization these models gave ∆E0 values, respectively, of -26.7 and -27.4 kcal/mol, improvements of 2.8 and 2.2 kcal/ mol relative to B-LYP/6-311++G(d,p) and /TZVP+ (Table 1). 5. Cyanide/H2O. Two complexes are possible: CN-/H2O (N acceptor) and NC-/H2O (C acceptor). These cannot be distinguished by ∆H°, as recent ab initio values10,49 for both complexes are the same within experimental uncertainty and agree well with the experimental ∆H°. The B-LYP and B3LYP models retain these features (Table 3), remaining at least within 0.8 kcal/mol of the experimental ∆H°. Thus, the DF and MP2 results pose the same interpretation, i.e., that the gasphase experimental results refer to a mixture of the two types of complex.10 6. HCC-/H2O. The B-LYP and B3-LYP ∆H° values are both ∼2 kcal/mol lower than the experimental ∆H° (Table 3), but again B3-LYP is closer in geometry to MP2 (Table 1). The MP2 model itself is somewhat marginal for this system,10 giving a ∆H° 1.2 kcal/mol lower than the experimental result (Table 3). All three models would probably benefit significantly from a straightforward expansion of the basis set. 7. HCOO-/H2O. This is another case with two possibilities for complexation: a symmetrical (C2V) complex with double H-bonds involving both carboxylate oxygens as acceptors and both water O-H bonds as donors, and a C1 complex with a single, “linear” H-bond. (See Figure 1.) Recent ab initio studies6,10 indicate that these complexes should be distinguishable on the basis of ∆H°, which is 2.6 kcal/mol lower for the C2V complex according to MP2/6-311++G(d,p) (Table 3). On the other hand, the experimental ∆H° is -16.0 kcal/mol,44 within experimental uncertainty of the MP2 value for the C1 complex (Table III). Along with other compelling experimental results,44 this suggests that the C1 complex is the one actually observed. This situation has been discussed elsewhere;10,44 here, we focus only on a direct comparison between DF and ab initio results.

4788 J. Phys. Chem., Vol. 100, No. 12, 1996

Pudzianowski

TABLE 6: Timing Comparisons for “Native” and ex Post Facto DF Methods (a) CH3NH3+/NH3 at Fixed Geometry: Energy and Gradient modela

CPU time,b s

CPU time (relative)

2195.4 1043.6 448.3 92.0

1.000 0.4754 0.2042 0.0419

MP2 B-LYP HF B-LYP DGauss

(b) Full Geometry Optimization CPU Times for Three Complexes modelc

OH-/H2O

CH3NH3+/H2O

CH3NH3+/NH3

B6 G, s B6 D, s G/D BT G, s BT D, s G/D

1109.0 104.6 10.6

18592.3 1046.8 17.8 8102.4 481.8 16.8

29458.0 1224.3 24.1 20999.0 885.0 23.7

a All with 6-311++G(d,p). b On a single Cray Y-MP/2E processor. Designations are B6, B-LYP/6-311++G(d,p); BT, B-LYP/TZVP; G, Gaussian 92/DFT; D, DGauss 2.3; G/D, ratio of G CPU time to D CPU time. c

Only B3-LYP follows MP2 in putting ∆H° for the C1 complex close to the experimental value; in fact, the B3-LYP and MP2 values are the same (Table 3). On the other hand, B-LYP maintains the relative stability of the C2V complex but puts its ∆H° close to the experimental value. Thus, B-LYP is inconsistent with the relevant experimental observations. The vibrational modes given by both DF models (Table S3) agree well with the MP2 results.10 Even low-frequency details are quite similar; for example, the B-LYP/ and B3-LYP/6311++G(d,p) models both give the same B2 mode found with MP2 for the C2V complex (Table S3), nuclear displacements for which are shown in Figure 2. E. “Native” and ex Post Facto DF Methodologies. A “native” DF method fully exploits the KS formalism by using model densities, derived from auxiliary basis sets, to separate the evaluation of classical Coulomb potential energy, and local and nonlocal exchange-correlation functionals which are densitydependent, from the remaining HF-like calculations involving one-electron orbitals.18 This approach eliminates all four-center AO integrals (and includes correlation) so it is formally expected to be less costly than RH methodology. An ex post facto DF method, on the other hand, does not follow the Kohn-Sham formalism from the outset, but instead starts with existing RH methodology and grafts on DF constructs: RH calculations give a set of one-electron orbitals, and then HF exchange components are removed from the energy and all exchange-correlation functionals are evaluated numerically.22 Thus, an RH method can be adapted to do DF calculations, but the process entails essentially all the work of an HF step along with additional numerical computations, and so is by definition costlier than an uncorrelated RH calculation. It is still, however, less expensive than conventional correlation treatments with a given basis set. The practical aspects of these approaches were explored in comparisons of two implementations, the DGauss system38 representing “native” methods and Gaussian 92/DFT34 representing ex post facto methods. This was done in two stages. First, the energy and its gradient were evaluated with Gaussian 92/DFT for the CH3NH3+/NH3 complex at a fixed geometry, using HF/, MP2/, and B-LYP/6-311++G(d,p). A DGauss B-LYP/6-311++G(d,p) calculation was then done for comparison. The timing results are reported in Table 6a. A second set of calculations was then performed with Gaussian 92/DFT and DGauss, this time involving Cartesian coordinate geometry

optimizations with B-LYP/6-311++G(d,p) on the OH-/H2O, CH3NH3+/H2O, and CH3NH3+/NH3 complexes, and B-LYP/ TZVP on the CH3NH3+/H2O and CH3NH3+/NH3 complexes. Timing comparisons for these calculations are reported in Table 6b. Two things should be recognized before these results are considered. First is the simple fact that geometry optimizations can bring specific software characteristics into play, thus somewhat complicating a comparison of the underlying methodologies. A more subtle consideration involves some slight differences between the Gaussian 92/DFT and DGauss implementations of the TZVP and 6-311++G(d,p) basis sets, but these are quite minor and do not affect our conclusions. The Gaussian 92/DFT results in Table 6a allow us to compare an ex post facto DF approach and MP2 with the same software. The B-LYP/6-311++G(d,p) model needed slightly less than half the CPU time of MP2. Still, the B-LYP calculation took ∼2.3 times longer than HF, so it is clear that such DF approaches will remain tied to some of the computational limitations of RH methodology, as was expected from formal considerations alone. By contrast, the corresponding B-LYP/ 6-311++G(d,p) fixed-geometry calculation with DGauss was ∼11 times faster than with Gaussian 92/DFT, ∼24 times faster than MP2, and even ∼5 times faster than HF. For full geometry optimizations the speed increases reported in Table 6b for DGauss relative to Gaussian 92/DFT ranged from 10.6 times for OH-/H2O to 24.1 times for CH3NH3+/NH3 with B-LYP/6-311++G(d,p) and were similar to those for B-LYP/TZVP. An interesting comparison not made in the table is with the corresponding time for full geometry optimization of OH-/H2O by MP2, which was 2624.3 s; hence, the “native” DGauss B-LYP calculation was ∼25 times faster, while the ex post facto Gaussian 92/DFT B-LYP calculation was ∼2.4 times faster than MP2. Conclusions (a) Like correlated RH methods, LCAO DF methods require diffuse atomic functions to adequately describe ionic H-bonded systems. This requirement is more pronounced for anionic H-bonded systems than for cationic systems. (b) Given (s,p) shells of diffuse functions on all non-hydrogen centers, DF models show little or no improvement on further augmentation with s-type diffuse functions on hydrogens. Again, this parallels the behavior of correlated RH models. (c) Without augmentation, the Gaussian TZVP and numerical DNP and DNPP basis sets developed for DF calculations are inadequate for ionic H-bonded systems. Some effort should be devoted to constructing sets of diffuse functions for these basis sets. In the case of DMol it may also be worthwhile to explore triple-numerical bases, since reliable analytical LCAO models require triple-split basis sets at the least. (d) B-LYP/ and B3-LYP/6-311++G(d,p) models are slightly inferior to MP2/6-311++G(d,p), which approaches “chemical accuracy” for the 10 ionic H-bonded systems studied here. B3LYP was appreciably the better DF model, with complexation enthalpies whose rms deviations from MP2 and experimental values were, respectively, 0.7 and 1.6 kcal/mol. The B3-LYP model was similar to the MP2 model in nearly all respects, whereas B-LYP/6-311++G(d,p) was less consistent. (e) There is evidence that judicious additions to the 6-311++G(d,p) and TZVP+ basis sets will improve the DF models, perhaps to “chemical accuracy”. In particular, the use of diffuse functions of higher angular momentum than occupied valence shells, e.g., p-functions on hydrogens, may be a fruitful approach.

Hydrogen Bonding in Binary Ionic Complexes (f) Efficient auxiliary basis/model density DF methods are much faster than correlated RH models at the same basis set level. Hence, current DF models can already extend theoretical treatments of good quality to ionic H-bonded systems several times larger than those to which RH MP2 models can be applied. Acknowledgment. Gaussian basis set calculations were done on Cray Y-MP and C94 computers in the BMSPRI HighPerformance Computing Group, and DMol calculations were done on a Silicon Graphics Onyx server in the Macromolecular Structure Division. I thank Messrs. Richard Shaginaw, John Stringer, and Stewart Samuels and Dr. Robert Bruccoleri for their advice and help with these resources, and Dr. Joseph Villafranca for maintaining a stimulating research environment. Thanks are also due to Drs. Jan Andzelm and George Fitzgerald for much helpful advice. Supporting Information Available: Table S1 which reports electronic and zero-point energies, thermal energies, entropies, and donor H-D bond lengths for isolated ions and molecules, Table S2 with the same information for the optimized complexes, and Table S3 with intermolecular vibrational modes and frequencies (10 pages). This material is contained in many libraries on microfiche, immediately follows this article in the microfilm version of the journal, can be ordered from the American Chemical Society, and can be downloaded from the Internet; see any current masthead page for ordering information and Internet access instructions. References and Notes (1) Kollman, P. A.; Allen, L. C. J. Am. Chem. Soc. 1970, 92, 61016107. (2) Meot-Ner (Mautner), M. In Molecular Structure and Energetics; Liebman, J. F., Greenberg, A., Eds.; VCH: New York, 1987; Vol. 4, Chapter 3. (3) Deakyne, C. A. In Molecular Structure and Energetics; Liebman, J. F., Greenberg, A., Eds.; VCH: New York, 1987; Vol. 4, Chapter 4. (4) Meot-Ner (Mautner), M. J. Am. Chem. Soc. 1988, 110, 3071-3075. (5) Meot-Ner (Mautner), M. J. Am. Chem. Soc. 1988, 110, 3075-3080. (6) Zheng, Y.-J.; Merz, K. M. J. Comput. Chem. 1992, 13, 11511169. (7) (a) Cleland, W. W.; Kreevoy, M. M. Science 1994, 264, 18871890. (b) Frey, P. A.; Whitt, S. A.; Tobin, J. B. Science 1994, 264, 19271930. (8) Frisch, M. J.; Del Bene, J. E.; Binkley, J. S.; Schaefer, H. F., III. J. Chem. Phys. 1986, 84, 2279-2289. (9) Del Bene, J. E. Int. J. Quantum Chem.: Quantum Chem. Symp. 1992, 26, 527-541. (10) Pudzianowski, A. T. J. Chem. Phys. 1995, 102, 8029-8039. (11) Hehre, W. G.; Radom, L. A.; Pople, J. A.; Schleyer, P. v. R. Ab Initio Molecular Orbital Theory; Wiley-Interscience: New York, 1986. (12) Dunning, T. H., Jr. J. Chem. Phys. 1989, 90, 1007-1023. (13) Clark, T.; Chandrasekhar, J.; Spitznagel, G. W.; Schleyer, P. v. R. J. Comput. Chem. 1983, 4, 294. (14) Hohenberg, P.; Kohn, W. Phys. ReV. B 1964, 136, 864-870. (15) Kohn, W.; Sham, L. J. Phys. ReV. A 1965, 140, 1133-1138. (16) Parr, R. G.; Yang, W. Density-Functional Theory of Atoms and Molecules; Oxford: New York, 1989. (17) Density Functional Methods in Chemistry; Labanowski, J. K., Andzelm, J. W., Eds.; Springer: New York, 1991. (18) Ziegler, T. Chem. ReV. 1991, 91, 651-667. (19) Dixon, D. A.; Andzelm, J.; Fitzgerald, G.; Wimmer, E.; Jasien, P. In ref 17, pp 33-48. (20) (a) Langreth, D. C.; Mehl, M. J. Phys. ReV. B 1983, 28, 1809. (b) Perdew, J. P.; Harbola, M. K.; Shani, V. In Condensed Matter Theory; Arponen, J. S., Bishop, R. F., Manninen, M., Eds.; Plenum: New York, 1988; Vol. 3, p 235.

J. Phys. Chem., Vol. 100, No. 12, 1996 4789 (21) (a) Dunlap, B. I.; Connolly, J. W. D.; Sabin, J. R. J. Chem. Phys. 1979, 71, 3396. (b) Salahub, D. R.; Fournier, R.; Mlynarski, P.; Papai, I.; St-Amant, A.; Ushio, J. In ref 17, pp 77-100. (c) Andzelm, J. In ref 17, pp 155-174. (22) (a) Gill, P. M. W.; Johnson, B. G.; Pople, J. A.; Frisch, M. J. Int. J. Quantum Chem.: Quantum Chem. Symp. 1992, 26, 319-331. (b) Johnson, B. G.; Gill, P. M. W.; Pople, J. A. J. Chem. Phys. 1992, 97, 78467848. (23) Delley, B. J. Chem. Phys. 1990, 92, 508-517. The DMol program is available commercially from BIOSYM, San Diego, CA. (24) Hill, R. A.; Labanowski, J. K.; Heisterberg, D. J.; Miller, D. D. In ref 17, pp 357-372. (25) (a) Sim, F.; St-Amant, A.; Papai, I.; Salahub, D. R. J. Am. Chem. Soc. 1992, 114, 4391-4400. (b) Kim, K.; Jordan, K. D. J. Phys. Chem. 1994, 98, 10089-10094. (c) Kieninger, M.; Suhai, S. Int. J. Quantum Chem. 1994, 52, 465-478. (d) Wei, D.; Salahub, D. R. J. Chem. Phys. 1994, 101, 7633-7642. (e) Barone, V.; Orlandini, L.; Adamo, C. Chem. Phys. Lett. 1994, 231, 295-300. (f) Gresh, N.; Leboeuf, M.; Salahub, D. In Modeling The Hydrogen Bond; Smith, D. A., Ed.; ACS Symposium Series 569; American Chemical Society: Washington, DC, 1994; Chapter 6. (g) Del Bene, J. E.; Person, W. B.; Szczepaniak, K. J. Phys. Chem. 1995, 99, 10705-10707. (26) Pudzianowski, A. T. Presented at the 8th International Congress of Quantum Chemistry Satellite Symposium, Thirty Years of Density Functional Theory: Concepts and Applications, Cracow, Poland, June 1994. (27) Godbout, N.; Salahub, D. R.; Andzelm, J.; Wimmer, E. Can. J. Chem. 1992, 70, 560-571. (28) Lee, C.; Yang, W.; Parr, R. G. Phys. ReV. B 1988, 37, 785-788. (29) (a) Becke, A. D. Phys. ReV. A 1988, 38, 3098. (b) Becke, A. D. J. Chem. Phys. 1993, 98, 5648-5652. (30) (a) Hariharan, P. C.; Pople, J. A. Chem. Phys. Lett. 1972, 66, 217. (b) Krishnan, R.; Binkley, J. S.; Seeger, R.; Pople, J. A. J. Chem. Phys. 1980, 72, 650. (31) Guan, J.; Duffy, P.; Carter, J. T.; Chong, D. P.; Casida, K. C.; Casida, M. E.; Wrinn, M. J. Chem. Phys. 1993, 98, 4753-4765. (32) Vosko, S. J.; Wilk, L.; Nusair, M. Can. J. Phys. 1980, 58, 12001211. (33) (a) Hedin, L.; Lundqvist, B. I. J. Phys. C. 1971, 4, 2064-2083. (b) Moruzzi, V. L.; Janak, J. F.; Williams, A. R. Calculated Electronic Properties of Metals; Pergamon: New York, 1978; p 26. (34) Gaussian 92/DFT, Revision G.1; Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Gill, P. M. W.; Johnson, B. G.; Wong, M. W.; Foresman, J. B.; Robb, M. A.; Head-Gordon, M.; Replogle, E. S.; Gomperts, R.; Andres, J. L.; Raghavachari, K.; Binkley, J. S.; Gonzalez, C.; Martin, R. L.; Fox, D. J.; Defrees, D. J.; Baker, J.; Stewart, J. J. P.; Pople, J. A. Gaussian, Inc.: Pittsburgh, PA, 1993. (35) Schlegel, H. B. J. Comput. Chem. 1989, 3, 214. (36) Dewar, M. J. S.; Zoebisch, E. G.; Healy, E. F.; Stewart, J. J. P. J. Am. Chem. Soc. 1985, 107, 3902-3909. (37) (a) Boys, S. F.; Bernardi, F. Mol. Phys. 1970, 19, 553-566. (b) van Duijneveldt, F. B.; van Duijneveldt-van de Rijdt, J. G. C. M.; van Lenthe, J. H. Chem. ReV. 1994, 94, 1873-1885. (38) Andzelm, J.; Wimmer, E. J. Chem. Phys. 1992, 96, 1280. DGauss is part of the UniChem system of computational chemistry programs available from Cray Research, Inc. (39) SYBYL Version 6.04; TRIPOS Associates, Inc.: 1699 South Hanley Road, Suite 303, St. Louis, MO 63144, 1993. (40) SPARTAN Version 3.1; Wavefunction, Inc.: 18401 Von Karman Ave., No. 370, Irvine, CA 92715, 1994. (41) (a) Del Bene, J. E. J. Comput. Chem. 1987, 8, 810-815. (b) Del Bene, J. E. J. Phys. Chem. 1988, 92, 2874-2880. (42) Hiraoka, K.; Takimoto, H.; Yamabe, S. J. Phys. Chem. 1986, 90, 5910-5914. (43) Meot-Ner (Mautner), M. J. Am. Chem. Soc. 1984, 106, 1257-1264. (44) Meot-Ner (Mautner), M.; Sieck, L. W. J. Am. Chem. Soc. 1986, 108, 7525-7529. (45) Meot-Ner, M. J. Am. Chem. Soc. 1988, 110, 3858-3862. (46) Xie, Y.; Remington, R. B.; Schaefer, H. F., III. J. Chem. Phys. 1994, 101, 4878-4884. (47) Pudzianowski, A. T. J. Chem. Phys. 1995, 102, 7761. (48) Del Bene, J. E.; Shavitt, I. J. Mol. Struct. (Theochem) 1994, 307, 27-34. (49) Del Bene, J. E. Struct. Chem. 1989, 1, 19-27.

JP9524688