J. Phys. Chem. 1995,99, 8013-8016
8013
Ab Initio Calculations of ProtiumDeuterium Fractionation Factors in OzHs+ Clusters Arthur S. EdisonJJ John L. MarkleyJ3s and Frank Weinhold*>f?o Graduate Biophysics Program, Department of Biochemistry and National Magnetic Resonance Facility ai Madison, and Theoretical Chemistry Institute and Department of Chemistry, University of Wisconsin at Madison, Madison, Wisconsin 53706 Received: January 13, 1995; In Final Form: February 28, 1995@
Fractionation factors have been calculated for the H/D isotopic equilibrium reactions of O2H5+ clusters. The calculations used ab inirio geometries and force constants to evaluate the reduced isotopic partition functions for each isotopic species. We find good agreement with experimental gas-phase measurements [Graul, S . T.; Brickhouse, M. D.; Squires, R. R. J. Am. Chem. SOC.1990,112,631-6391 at higher levels of Hartree-Fock theory. By contrast, Mgller-Plesset second order perturbation theory (MP2) yielded somewhat poorer agreement with experiment, presumably because the calculated hydrogen bond distances are underestimated at this level.
Introduction Equilibrium protium/deuterium isotope effects are sensitive probes of molecular structure and hydrogen bonds in a variety of molecules.'-'0 For the exchange reaction A-H
+ B-D
t A-D
+ B-H
(1)
the equilibrium constant (known as the fractionation factor) will, in general, not be equal to unity. Deviations from unity occur as a result of the tendency of deuterium and protium to fractionate in a nonstatistical distribution, in order to minimize the Gibbs free energy of the reaction. Fractionation factors, measured for many compounds (for example, refs 2, 8, 11, and 12), are often less than unity for strongly hydrogen-bonded A-H complexes and can be as low as 0.3 for very strongly hydrogenbonded symmetric clusters of the form A- -H- -A,* Recent experimental data have shown that amide hydrogens in proteins (not often thought to be involved in symmetric, strong hydrogen bonds) have ranges of fractionation factors from 0.3 to 1.7.'~~ In the Bom-Oppenheimer approximation, the shapes of the electronic potential energy surfaces are invariant with respect to isotopic substitution. Using this approximation, Bigeleisen and Mayer showedI3 that the equilibrium constant for eq 1 can be estimated by the ratios of reduced isotopic partition functions
where
' Graduate Biophysics Program.
*
Department of Biochemistry and National Magnetic Resonance Facility at Madison. 5 Theoretical Chemistry Institute and Department of Chemistry. Current address: Department of Zoology, University of Wisconsin, Madison, WI 53706. Abstract published in Advance ACS Abstracts, April 15, 1995. @
Here, ui = hvi/kT, vi are the harmonic vibrational frequencies, s are the symmetry numbers, and the prime refers to the light isotope. Although the symmetry numbers often play no role in the reaction and are neglected,I3-l5 here they must be included. Previous studies have used eq 2 with harmonic vibrational frequencies obtained from ab initio molecular orbital calculation^.^^'^-'^ Hout and co-workersI4 were able to reproduce experimental trends (often with a high degree of accuracy) for isotope exchange reactions involving non-hydrogen-bonded molecules with calculations that used 3-21G and 6-3 lG* basis sets at the Hartree-Fock level and 6-31G* at the MP2 level of theory." By using HF/6-3 lG* without complete geometry optimization, Saunders and co-workers were able to obtain fractionation factors proportional to experimental values for bicyclobutonium and cyclopropylcarbinyl ions.I5 Weil and Dixon have used frequencies from HF/4-3 1G calculations with some empirical adjustments to successfully calculate the fractionation factor for hydrogen-bonded methoxide anion^.^ Williams has carried out an empirical calibration against experimental data of calculated fractionation factors for carbonhydrogen bonds; he concluded that calculations at HF/4-3 1G provide reasonable correlations with experiment (but are consistently too large). l6 Graul and co-workers have reported gas-phase deuterium equilibrium isotope data from protonated water clusters obtained from collision-induced dissociation experiments.I2 Owing to the small cluster sizes, these gas phase data provide an ideal test for the theory of equilibrium isotope effects. The O2H5+ cluster has been studied extensively theoretically. Am0 ng the more noteworthy studies have been Del Bene's analyses of the effects of basis set and electron correlation on hydrogenbonding energies21.22and Frisch and co-worker's calculations at very high levels of theory (MP4/6-31l++G(3df,3pd) //MP2/ 6-31 l++G(2d,2p)).*O We report here the results of calculations of fractionation factors for O2H5+ and provide a direct comparison with the experimental results of Graul and co-workers.I2 These calculations have employed a wide range of basis sets at both the Hartree-Fock and the second-order Mgller-Plesset perturbation levels of theory, so that the effects of basis set and electron correlation (at MP2 level) could be evaluated. The highest level used in this study was MP2/6-3 1l++G(2d,2p)//MP2/631 l++G(2d,2p). After briefly describing our methods, we
0022-3654/95/2099-8013$09.00/0 0 1995 American Chemical Society
8014 . I Phys. . Chem., Vol. 99, No. 20, 1995
Edison et al.
TABLE 1: Total Energies and Distances between the Oxygen Atoms (R)and the Central Hydrogen to Oxygens (r) in the 02H5+ Cluster at Different Levels of Theow Rb(&
reaction
#(A)
basis
E(hartrees)
3-21G 6-31G* 6-3 1+G* 6-311++G* 6-311++G** 6-311++G(2d,2p)
RHF -151.56094 -152.352 49' - 152.357 50 - 152.406 60 -152.435 49 -152.439 38
2.380 2.457f 2.482 2.473 2.417 2.431
1.190 1.4068 1.444 1.452 1.362 1.386
1.190 1.051 1.038 1.022 1.056 1.046
6-31G* 6-31+G* 6-3 1 1++G*' 6-3 1 1++G*' 6-311++G(2d,2p)
MP2 -152.732 93" -152.743 16 -152.822 99 -152.82241 -152.918 43k
2.415 2.418 2.436 2.387 2.389'
1.210 1.211 1.361 1.194 1.196
1.210 1.211 1.076 1.194 1.196
F(A)
TABLE 2: Equilibrium Constants for the Isotopic Equilibrium Reactions of OzH5+
"Each level of theory used geometries optimized (no imaginary frequencies unless specifically mentioned) at that level. 01- 0 2 distance. 01-H distance, where the H is the atom in the central hydrogen bond. 02-H distance, where the H is the atom in the central hydrogen bond. e Reference 23 reports - 152.352 43 hartrees, but ref 19 reports -152.352 49, both with 6-31G(d). f Reference 23 and ref 19 both report 2.455 A with, 6-31G* and 6-31G(d), respectively. g Reference 23 reports 1.404 A. Reference 23 reports -152.737 99 hartrees, but our geometry is in agreement. IC, structure with no imaginary frequencies. 1 C2 structure with one imaginary frequency (-615.8 cm-I). Refereye 20 reports -152.961 025 hartrees. Reference 20 reports 2.387 A, but the 0-H distances are in agreement.
'
present calculated equilibrium constants for single, double, triple, and quadruple deuterium substitutions in 02Hs+. These results were obtained at the restricted Hartree-Fock (RHF) level by using the 3-21G, 6-31G*, 6-31+G*, 6-311++G*, 6-311++G**, and 6-311++G(2d,2p) basis sets and at the Moller-Plesset second-order perturbation (MP2) level by using the 6-31G*, 6-31+G*, and 6-31 l++G*, and 6-31 l++G(2d,2p) basis sets. All of the theoretical levels reproduce the qualitative experimental trends, with the higher level RHF calculations performing best. In accordance with a previous study which indicated that MP2 calculations tend to overestimate the strengths and underestimate the lengths of hydrogen bonds,I7we find the MP2 level gives poorer results for isotope fractionation factors in these strongly hydrogen-bonded complexes.
Methods Equation 2 was calculated with the QUIVER computer program,15 which uses the Cartesian coordinates and force constant matrix from Gaussian 92.25 All of the calculations at both the RHF and MP2 levels were fully geometry optimized with no symmetry constraints, and additional MP2 calculations were obtained using C2 ~ y m m e t r y . 'Analytical ~ ~ ~ ~ frequencies for each set of optimized coordinates were calculated, and no imaginary frequencies were observed, with the exception of a single large frequency (-615.8 cm-I) which resulted from an MP2/6-3 11++G* calculation with C2 symmetry; this indicated that the equilibrium geometry is not C2 at this level. This result was unexpected, because all previously reported MP2 structures studied have been found to have CZ~ y m m e t r y . ~ ~ ,However, *~,*~ the MP2/6-3 11++G* structure optimized with no constraints was found to be of C, symmetry, had lower total energy than the C2 structure (Table l), and yielded only real positive frequencies. Interestingly, the calculated fractionation factors for the MP2/6-31 l++G* structure are in better agreement with experiment than other basis sets at the MP2 level (Table 2). These results indicate that the equilibrium geometry is severely distorted with even slight changes in theoretical level and that anharmonic corrections may be unusually important in this
basis
A
B
C
D
E
F
3-21G 6-31G* 6-3 1+G* 6-311++G* 6-311++G** 6-311++G(2d,2p)
0.48 0.56 0.60 0.64 0.48 0.56
RHF 0.97 0.97 0.97 0.97 0.97 0.97
0.46 0.54 0.61 0.63 0.49 0.54
0.48 0.56 0.64 0.66 0.50 0.56
0.46 0.54 0.61 0.64 0.49 0.54
0.46 0.55 0.62 0.64 0.50 0.55
6-31G* 6-3 1+G* 6-311++G*b 6-311++G(2d,2p)
0.32 0.32 0.48 0.32
MP2 0.97 0.97 0.97 0.97
0.31 0.30 0.48 0.31
0.32 0.30 0.50 0.32
0.32 0.30 0.49 0.32
0.32 0.30 0.49 0.32
Experiment' 0.85 0.63
0.16
0.63
0.65
0.74
A: DHO- -H- -OH* 2 H20- -D- -OH2 (sa/sb = 1/4d) B: DzO- -H--OHz Z DHO- -H- -0HD (sa/sb = 2/ld) C: D2O- -H- -OH2 2 DHO- -D- -OH2 (sa/sb = 2/ld) D: D2O- -H- -OHD 2 DzO- -D- -OH2 (sa/sb = 1/2d) E: D20--H- -0HD 2 DHO- -D- -0HD (sa/sb = l/ld) F: Dz0- -H- -OD2 2 D20- -D- -0HD ( s J s ~= 4/1') a The reactions, labeled A-F, are shown at the bottom of the table. C, structure. Experimental equilibrium constants are from Table 4 of ref 12 and were corrected by dividing by the ratios of symmetry numbers given for each reaction. The symmetry numbers, sa and s b , are for the reaction a 2 b.
SCHEME 1 H
H,r
\
species. We found some discrepancies with previous studies in energies and geometries; these are all indicated as footnotes in the appropriate tables. The rotational symmetry numbers are determined by the point group of the clusters. However, elucidation of the point group is not completely straightforward. Del Bene and co-workers found that differentAlevels of theory and different basis sets yielded different equilibrium structures for the 02H5+ clusters; this demonstrates the flat, anharmonic nature of the potential energy s ~ r f a c e . ' Scheme ~ 1 (adapted from ref 19) shows the symmetry of the structures found at different levels of theory. With the exceptions of 3-21G and 6-311++G* for RHF and MP2, respectively, RHF predicts a C, and MP2 predicts a C2 structure. Regardless of the symmetry of the calculated complexes, the effective experimental symmetry number will be 4 for the fully protonated or deuterated species as the result of rapid interconversion between different static configurations of the four non-H-bonded hydrogens. This is consistent with Graul's assignments of statistical expectations for the experimental distributions of isotopes in the experimental 02H5+ clusters.'* Analogous statistical symmetry numbers used for other species are given at the bottom of Table 2. Thus, for comparison with our calculated fractionation factors, the experimental data shown in Table 2 were divided by the ratios of the symmetry numbers. 9320
Results and Discussion Table 1 gives the calculated total energies, oxygen-oxygen distances, and the distances between the central hydrogen and
Clusters
J. Phys. Chem., Vol. 99, No. 20, 1995 8015
TABLE 3: Average Values of the Symmetry ,Uncorrected Equilibrium Constant Calculated for the General Reaction LzO- -H- -LDO 2 LzO- -D- -LHO, Where L Is Either H or D basis (s/s')f (avg) basis (s/s')f (avg)
Although the MP2 geometries appear to be the cause of the poor calculations of fractionation factors in 02H5+, we cannot rule out the possibility that the apparent success of the better RHF calculations is due to a cancellation of errors. Vibrational anharmonicities might be substantial and might oppose the effects of electron correlation. Without oversimplifying the problem to a single vibrational mode,2 it is difficult to incorporate the effects of anharmonicities in the vibrational partition function, and we do not attempt to estimate this correction in the present work.
ProtiumDeuterium Fractionation Factors in
~~~
02H5+
~
3-21G 6-31G* 6-31+G*
0.46 0.55 0.62
RHF 6-311++G* 6-311++G** 6-31 l++G(2d,2p)
0.64 0.50 0.55
6-3 lG* 6-31+G*
0.32 0.30
MP2 6-311++G* (Cs, 6-31 l++G(2d,2p)
0.49 0.32
Conclusions
each oxygen atom for various theoretical levels. Table 2 gives the calculated and experimental equilibrium constants, the ratios of symmetry numbers, and experimental equilibrium constants for the isotopic reactions of the O2H5+ cluster. The frequencies for all our calculations are similar to previously published values and are available from the authors. From Table 2, it is clear that, for any given theoretical level, the calculated equilibrium constant was remarkably consistent for the general reaction LZO- -H- -LDO(+)
2 LZO-
-D- -LHO(+)
(3)
where L is either H or D. The average value of the equilibrium constant for eq 3 for each basis set and level of theory is shown in Table 3. The average values shown in Table 3 can be compared with the experimental values for
'/,H,O+
+ '/,D20 * '/,D,Of + '/,H,O
(4)
Reported experimental equilibrium constants for eq 4 range from 0.69 & 0.02 at 25 "C in the liquid phase26-28to 0.79 f 0.01 in the gas phase29 or for lyonium ion in a ~ e t o n i t r i l e . ~ The ~ calculations presented in this work represent gas phase values; in a companion paper on peptide fractionation factors, we have converted the results to liquid phase to provide a direct comparison with experiment (Edison, A. S.; Weinhold, F.; Markley, J. L. J. Am. Chem. SOC.,submitted). Graul and coworkers found experimental equilibrium constants that ranged from 0.66 to 0.80 as the number of total deuterons increased.I2 Our calculations showed no similar increase in fractionation factor when the number of deuterons was increased. The insensitivity to the basis set found for the calculations for reaction B (Table 2) results from a near-zero Gibbs free energy difference for the reaction. It appears from Tables 2 and 3 that the best theoretical results are obtained by using RHF/6-31+G* and RHF/6-31 l++G* levels of theory. All other levels of theory lead to frationation factors that are somewhat too low with respect to experiment. This trend can be compared to the calculated oxygen-oxygen bond lengths (Table l), which ranged from 2.38 to 2.48 A. The levels of theory that had the best agreement with experimental fractionation factors had 0-0 distances of about 2.47-2.48 A. On the other hand, the levels with the worst agreement with experiment had 0-0 distances between 2.38 and 2.42 A. Other studies on the neutral water dimer showed that RHF/6-31G* reproduced the experimental 0-0 bond length (2.98 A) whereas MP2/6-31G* calculations yielded a shorter (2.92 A) distance.I7 In an extensive set of calculations of the hydrogen bond energies of O*H5+ and other clusters, Del Bene found that the experimental binding enthalpy could be reproduced quite accurately by using the highest level of theory (MP4/6-311+G(2d,2p)) for the electronic association energy at a Hartree-Fock C, structure with the RHF/6-31G(d) zero-point and vibrational energy corrections.21 Del Bene's results seem to be consistent with ours.
We have calculated the equilibrium constants for a series of isotopic species of the 02H5+ cluster and have compared our results to gas phase experimental results.I2 The fractionation factor was found to depend sensitively on H bond geometry. We found quite good agreement with experiment for 6-3 1+G* and 6-311++G* basis sets at the Hartree-Fock level. RHF/ 3-21G and 6-3 lG* reproduced the qualitative experimental trends (although with lower quantitative accuracy) and may well be useful in studies of larger complexes where highly extended basis sets are impractical. The MP2 level of theory appeared less reliable for these calculations, most likely because it underestimates hydrogen bond lengths.
Acknowledgment. We thank Profs. W. Wallace Cleland, Maurice Krevoy, and Robert Squires for helpful discussions and Dr. Stewart Loh for providing the stimulus for pursuing this study. Prof. Martin Saunders kindly provided the source code for QUIVER which was used in calculations of the reduced isotopic partition functions. Ab initio calculations were performed with Gaussian 92 installed on a KubotdStardent TITAN in the Chemistry Department and a Cray Y-MP at Cray Research, Inc., and with GAMESS 31 installed on a Silicon Graphics Indigo computers in the Biochemistry Department. A.S.E. was a trainee of an NIH Molecular Biophysics Training Grant (GM08293). J.L.M. has support from NIH Grant RR02301. References and Notes (1) Kreevoy, M. M. In Isotopes in Organic Chemistry; Buncel, E., Lee, C. C., Eds.; Elsevier Scientific Publishing: Amsterdam, 1976; Vol. 2, pp 1-31. (2) Kreevoy, M. M.; Liang, T. M. J. Am. Chem. SOC.1980,102,33153322. (3) Weil, D. A.; Dixon, D. A. J. Am. Chem. SOC. 1985, 107, 68596865. (4) Kresge, A. J.; More O'Ferrall, R. A,; Powell, M. F. In Isotopes in Organic Chemistry; Buncel, E., Lee, C. C., Eds.; Elsevier Scientific Publishing: Amsterdam, 1987; Vol. 7, pp 177-273. ( 5 ) Hibbert, F.; Emsley, J. In Advances in Physical Organic Chemistry; Bethell, D., Ed.; 1990; Vol. 26, pp 255-379. (6) Cleland, W. W. Biochemistry 1992, 31, 317-319. (7) Loh, S. N.; Markley, J. L. Biochemistry 1994, 33, 1029-1036. (8) Loh, S. N.; Markley, J. L. In Techniques in Protein Chemistry; Angeletti, R., Ed.; Academic Press: New York, 1993; Vol. 4, pp 517524. (9) Cleland, W. W.; Kreevoy, M. M. Science 1994,264, 1887-1890. (10) Frey, P. A.; Whitt, S. A.; Tobin, J. B. Science 1994, 264, 19271930. (11) Jarret, R. M.; Saunders, M. J . Am. Chem. SOC. 1985, 107, 26482654. (12) Graul, S. T.; Brickhouse, M. D.; Squires, R. R. J. Am. Chem. SOC. 1990, 112, 631-639. (13) Bigeleisen, J.; Mayer, M. G. J. Chem. Phys. 1947, 15, 261-267. (14) Hout, R. J., Jr.; Wolfsberg, M.; Hehre, W. J. J. Am. Chem. SOC. 1980, 102, 3296-3298. (15) Saunders, M.; Laidig, K. E.; Wolfsberg, M. J. Am. Chem. SOC. 1989, 1I 1 8989-8994. (16) Williams, I. H. J . Phys. Org. Chem. 1990, 3, 181-190. (17) Hehre, W. J.; Radom, L.; Schleyer, P. v. R.; Pople, J. A. Ab Initio Molecular Orbital Theory; John Wiley & Sons: New York, 1986. ~
8016 J. Phys. Chem., Vol. 99, No. 20, 1995 (18) Kollman, P. A.; Allen, L. C. J. Am. Chem. SOC. 1970, 92, 61016107. (19) Del Bene, J. E.; Frisch, M. J.; Pople, J. A. J. Phys. Chem. 1985, 89, 3669-3614. (20) Frisch, M. J.; Del Bene, J. E.; Binkley, J. S.; Schaefer 111, H. F. J. Chem. Phvs. 1986. 84. 2279-2289. (21) Del Bene, J. E. J. Comput. Chem. 1987, 8, 810-815. (22) Del Bene, J. E. J. Phys. Chem. 1988, 92, 2874-2880. (23) Lee, E. P. F.; Dyke, J. M. Mol. Phys. 1991, 73, 375-405. (24) Peeters, D.; Leroy, G. J. Chim. Phys. 1991, 88, 411-420. (25) Frisch, M. J.; Trucks, G. W.; Head-Gordon, M.; Gill, P. M. W.; Wong, M. W.; Foresman, J. B.; Johnson, B. G.; Schlegel, H. B.; Robb, M. A.; 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 92, 1992 Revision A ; Gaussian Inc.: Pittsburgh, PA, 1992.
Edison et al. (26) Kresge, A. J.; Allred, A. L. J. Am. Chem. SOC. 1963, 85, 15411541. (27) Heinzinger, K.; Weston, R. E., Jr. J. Phys. Chem. 1964, 68, 744751. (28) Salomaa, P.; Aalto, V. Acta Chem. Scand. 1966, 20, 2035-2042. (29) Larson, J. W.; McMahon, T. B. J. Am. Chem. SOC. 1986, 108, 1719-1720. (30) Kurz, J. L.; Myers, M. T.; Ratcliff, K. M. J. Am. Chem. SOC. 1984, 106, 5631-5634. (31) Schmidt, M. W.; Baldridge, K. K.; Boatz, J. A,; Jensen, J. H.; Koseki, S.; Gordon, M ,S.; Nguyen, K. A.; Windus, T. L.; Elbert, S. T. QCPE Bull. 1990, 10, 52-54. JP9501660
