J. Phys. Chem. 1991, 95, 10457-10464 the extent of the preferential water solvation of the proton is yet to be determined. Finally, we point out that the generality of our conclusions can be tested even in cases where the water cluster models fail to predict the behavior of the R O H type acids. In particular, it was found that the dissociation rates from various naphthol derivatives decrease in the series water-MeOH, water-EtOH, and water-PrOH in a way which parallels the decrease in their c, values.24 A similar trend in the same cosolvent series was found in the dissociation rate of the H20molecule which dissociates according to H 2 0 Q OH- H+.79 In addition, the proton dissociation rates of several naphthol derivatives were found to be unchanged or even to increase in water-urea solutions,lW which is a strong gas-phase base, and its aqueous mixtures have a very large static dielectric constant. All these observations are consistent with a general structure-reactivity law which dictates the proton dissociation dynamics and point out the importance of the water cosolvent in determining the proton-transfer rate.
+
VI. Summary In this study we have demonstrated that photoacid dissociation rates obey a general structurereactivity law which correlates the proton dissociation rate with the pK, values of the precursor acid. By doing so we have shown that photoacids conform to the same empirical general reactivity laws which dictate the behavior of the overwhelming majority of ground-state acids. At the same time we do not find any support to the idea that these rates are determined in water-organic solutions by the kinetic availability
10457
of large water molecule clusters. An alternative view has been expressed in recent correspondence with Professor G. Robinson. H e proposes that cationic acids prepolarize the water structure around them and as such they are suitable neither for a direct comparison with the ROH type acids nor for the critical test of the validity of the various water-cluster models. Clearly more work is required to establish the generality of the two types of approaches. Our observations also indicate that the proton-transfer rate is greatly influenced by the chemical and physical properties of the water cosolvent. In particular, the basicity of a water molecule which is hydrogen bonded to an organic molecule increases when the gas-phase proton affinity of the organic molecule is larger than that of water. Such induced basicity in water was very recently reported for the water-DMSO system from direct spectrophotometric measurements of the coordination constant of 1:l water-DMSO dimer to a Ni2+ metal complex in nitrobenzene.lM As for the proton solvation shell, we consider the question of its precise nature in binary solution to be largely open. Both direct cosolvent participation and inductive modification of the water basicity should be considered in formulating a complete model. Acknowledgment. This work was supported by an NSF grant. Registry NO. 1 - A M P , 1606-67-3. (106) Iwamoto, E.; Nishimoto, J.; Yakoyama, T.; Yamamoto, K.; Kumamaru, T. J. Chem. Soc., Faraday Trans. 1991,87, 1537.
Tunneling Paths in Intramolecular Proton Transfer Norihiro Shida,+Jan Almlof,* and Paul F. Barbara* Department of Chemistry, University of Minnesota, 207 Pleasant Street S.E., Minneapolis, Minnesota 55455 (Received: April 16, 1991)
We recently introduced a procedure for determining effective tunneling paths in nonthennal symmetric proton-transfer reactions. The procedure involves a “maximum probability path” (MPP) based on the ground-state vibrational wave function, which can be determined within the reaction surface Hamiltonian framework by ab initio electronic structure methods. In this paper we further define the MPP concept and explore this path using the malonaldehyde molecule and the formic acid dimer as examples, as well as calculations on simpler model systems. The results verify our qualitative expectations for tunneling in heavy-light-heavy systems. Furthermore, a comparison of the MPP to other more conventional paths, such as the IRC, helps elucidate certain physical aspects of multidimensional tunneling. Finally, simple extensions of the MPP procedure suggest possible directions for new methods for calculating tunneling splitting in complex molecules.
splittings of molecules such as, e.g., malonaldehyde or ammonia, 1. Introduction are analogous to the tunneling splitting of carboxylic acid dimers: Nuclear tunneling between isomers in polyatomic molecules which involve a double proton transfer, as shown in Figure 2. plays an important role in the static and dynamic spectroscopy One characteristic feature of all of these systems is that several of inter- and intramolecular hydrogen bonded molecules in their vibrational degrees of freedom with very different effective masses ground and excited electronic states. Tunneling, for example, can are involved in the tunneling motion. In addition to the obvious be the dominant mechanism of irreversible excited-state intralight hydrogenic coordinate, more massive degrees of freedom are molecular proton transfer, and tunneling is thus essential to the strongly coupled to the tunneling motion. The stretching disultrafast kinetics of this reaction A less common but placement of the oxygen atoms bridged by the tunneling proton simpler role for nuclear tunneling of hydrogen-bonded systems is observed in smaller symmetrical systems like malonaldehyde (Figure l ) , for which nuclear tunneling is responsible for the (1) Many of the papers in this special issue deal with the role of nuclear splitting of the lowest energy vibrational wave function of these tunneling in kinetics and spectroscopy. Also, see the special issue of Chemical systems into two levels that can be detected spectros~opically,~~~Physics (1989, 136, No. 2) devoted to “Spectroscopy and Dynamics of Elein analogy with the splitting observed in an one-dimensional mentary Proton Transfer in Polyatomic Systems”. (2) For a recent review of excited-state proton-transfer research, see: double-minimum potential. Barbara, P. F.; Walsh, P. K.; Brus, L. E. J. Phys. Chem. 1989, 93, 29. For polyatomic systems the tunneling process is complex due (3) Rowe, W. F.; Durest, R. W.; Wilson, E. B. J. Am. Chem. Soc. 1981, to the large number of vibrational degrees of freedom. Tunneling 103, 6292. Present address: Division of Theoretical Studies, Institute for Molecular Science, Myodaiji, Okazaki, 44 Japan. To whom correspondence should be addressed.
0022-3654/91/2095-10457$02.50/0
(4) Firth, D. W.; Beyer, K.; Dvorak, M. A.; Reeve, S. V.; Grushow, A.; Leopole, K. R. J . Chem. Phys. 1991, 94, 1812. (5) Rambaud, C.; Oppenllnder, A.; Pierre, M.; Trommsdorff, H. P.; Vial, J. C. Chem. Phys. 1989, 136, 335.
0 1991 American Chemical Society
The Journal of Physical Chemistry, Vol. 95, No. 25, 1991
10458
-
r H--..f2 2.. *..,*.::.-o
\
r3
II
r1__--H+
0 ::: ..........0
H
I\
tunneling
H/‘\C/‘\H
HO‘\
A
r3
F H
I
/‘\H
Figure 1. Illustration of intramolecular proton tunneling in the malonaldehyde molecule.
p
,,,
.., , ., ., .H-
,O-H
0 \-H
4
H-C
‘0-H
............0
H
L
Tunneling
...........0
\\
H-C
‘0,., , .. . .. ..H
-0/c-H
Figure 2. Illustration of intermolecular double-proton tunneling in the formic acid dimer.
is especially important, since the barrier for hydrogenic motion decreases with decreasing 0-0 distance. These attributes relate the intramolecular tunneling splitting problem to the quantum mechanical calculation of chemical rates for hydrogen exchange reactions between heavy atoms, such as C1+ HCl ClH C1.6 Thus all of these proton- and hydrogen-transfer reactions fall into the ‘heavy-light-heavy” classification of chemical reactions. Theoretical progress on multidimensional nuclear tunneling has advanced greatly in the last decade. Calculations on even medium-size systems can still be prohibitively expensive, and less costly approaches continue to be developed. Both tunneling splitting c a l c ~ l a t i o n s and ~ - ~ the ~ related rate constant calculations6-” (see above) are made difficult by the now well-known feature of heavy-light-heavy systems, namely, that the tunneling occurs far from the so-called intrinsic reaction coordination, IRC, defined as a steepest-descent path in a mass-weighted, Cartesian coordinate system, starting from the saddle point (transition state) and descending to the reactant or product.I2 The IRC is a specific type of ‘reference” reaction path. The reaction path concept has played a central role in the calculation of classical and quantum reaction rates in small polyatomic systemsI3-I5 and recently also
-
+
(6) For recent reviews on semiclassical methods for calculation of reaction rates, see: Truhlar, D. G.; Garret, B. C. Annu. Reu. Phys. Chem. 1984, 35, 159. Miller, W. H. J . Phys. Chem. 1983, 87, 381 1. (7) Carrington, T.; Miller, W. H. J . Chem. Phys. 1986, 84, 4364. (8) Truhlar, D. G.;Steckler, R.; Gordon, M. S. Chem. Reu. 1987,87,217, Truhlar, D. G.; Brown, F. B.; Steckler, R.; Isaacson, A. D. In The Theory of Chemical Reaction Dynamics; Clary, D. C., Ed.; D. Reidel: Dordrecht, Holland, 1986. (9) Shida, N.; Barbara, P. F.; Almlof, J. J . Chem. Phys. 1989, 91, 4061. (IO) Shida, N.; Barbara, P. F.; Almlof, J. J. Chem. Phys. 1991,94, 3633. ( 1 1) Substantial progress has recently been made on quantum rate theory in complex systems. Some of this work has been described in: Hanggi, P.; Talkner, P.; Borkovec, M. Reo. Mod. Phys. 1990,62,251. Some other recent references include: Warshel, A,; Chu, Z. T. J . Chem. Phys. 1990, 93, 4003. Voth, G. A,; Chandler, D.; Miller, W. H. J . Chem. Phys. 1989, 91, 7749. Bongis, D.; Hynes, J. T. J. Chem. Phys. 1991, 94, 3619. (12) Fukui, K. J. Chem. Phys. 1970, 74,4161. Fukui, K. In The World of Quantum Chemistry; Daudel, R., Pullman, B., Ed.; Dordrecht: Netherlands, 1974. Fukui, K. Acc. Chem. Res. 1981, 14, 363. (13) Hofacker, G. L.; Seiler, R. J . Chem. Phys. 1969,51, 3951. Marcus, R. A. J. Chem. Phys. 1969, 45, 4493, 4500. Yamashita, K.; Yamabe, T.; Fukui, K. Chem. Phys. Lett. 1981, 84, 123. Yamashita, K.; Yamabe, T.; Fukui, K. Theor. Chim. Acta 1982, 60, 523. Yamabe, T.: Koizumi, M.; Yamashita, K.; Tachibana, T. J. Am. Chem. SOC.1984, 106, 2255. Truhlar, D. G.; Isaacson, A. D.; Skodje, R. T.; Garrett, B. C. J. Chem. Phys. 1982, 86, 2252. Garrett, B. C.; Truhlar, D. G.; Wagner, A. F.; Dunning, Jr., T. H.; J . Chem. Phys. 1983, 78, 4400. Skodje, R. T.; Schwenke, D. W.: Truhlar, D. G.; Garrett, B. C. J . Am. Chem. Soc. 1984,88,628. Brown, F. B.; Tucker, S. C.; Truhlar, D. G. J. Chem. Phys. 1985,83, 4451. Garrett, B. C.; Abusalbi, N.; Kouri, D. J.; Truhlar, D. G. J. Chem. Phys. 1985,83, 2242. Garrett, B. C.; Truhlar, D. G.; Bowman, J. M.; Wagner, A. F. J. Phys. Chem. 1986,90, 4305. Hancock, G. C.; Rejto, P. R.; Steckler, R.; Brown, F. B.; Schwenke, D. W.; Truhlar, D. G. J. Chem. Phys. 1986, 85,4997. Joseph, T.; Steckler, R.; Truhlar, D. G . J. Chem. Phys. 1987, 87, 7036. Gray, S. K.; Miller, W. H.; Yamaguchi, Y.; Schaefer 111, H. F. J. Chem. Phys. 1980, 73, 2733. YOsamura, Y.; Schaefer 111, H. F.; Gray, S. K.; Miller, W. H. J . Am. Chem. SOC.1981, 103, 1904. Waite, B. A.; Gray, S. K.; Miller, W. H. J . Chem. Phys. 1983, 78. 259. Hougen, J. H.; Bunker, P. R.; Johns, W. C. J . Mol Spectrosc. 1970, 34, 136. Truhlar, D. G.; Garrett, B. C. Acc. Chem. Res. 1980, 13, 440. Truhlar, D. G.; Garrett, B. C. Annu. Reu. Chem. 1984, 35, 159.
Shida et al. for condensed-phase reactions.l’ Qualitatively, it is known that tunneling in heavy-light-heavy systems occurs in a swath along an effective path that can deviate significantly from the IRC, involving a lighter effective mass but a higher effective barrier than tunneling along the IRC. The tunneling path especially avoids the saddle-point region which is accessed by motion of the heavy particles in a heavy-light-heavy system. Very few papers have directly dealt with the calculation and representation of “tunneling paths” in proton-transfer reactions. One early example is the calculation by Kuppermann et al. of “stream-lines” of the nuclear wave function for the prototypical H H2 reaction, leading to a clear pictorial repreH2 H sentation of the deviation from the IRC valley of the thermally assisted tunneling process in this reaction.I6 Another early example is due to Marcus and Coltrin.” Indeed, some of the methods for calculating quantum-mechanical rate constants implicitly rely on an effective path concept, although not explicitly in the sence of a representation of that path.6J1J8 In two recent papers on tunneling splitting calculations employing ab initio electronic structure methods for malonaldehydeg and formic acid dimer,I0 we found it informative to calculate effective tunneling paths for these systems using a procedure that determines a “maximum probability path” (MPP) based on the ground-state vibrational wave function. In this paper we further define the MPP concept and explore this path using the malonaldehyde and formic acid dimer examples, as well as calculations on simpler model systems. The results verify our qualitative expectations for tunneling in heavy-light-heavy systems. Furthermore, a comparison of the MPP to other more conventional paths, such as the IRC, helps elucidate certain physical aspects of multidimensional tunneling. Finally, simple extensions of the MPP procedure suggest possible directions for new methods for calculating tunneling splitting in complex molecules.
+
-
+
2. Maximal Probability Path In analogy to the IRC, which is a curvilinear path on the potential energy surface, a maximal probability path MPP can be defined as a curvilinear path on the negative of the nuclear probability density surface (i.e., -i$(q)12) of the ground vibrational stationary state. We are specifically concerned with a symmetric proton exchange reaction with a large barrier, such as in malonaldehyde (Figure 1) or in the formic acid dimer (Figure 2). The nuclear probability density of the ground vibrational state is clearly highest near the two equilibrium geometries. Consequently, the surface defined by -I$(q)12 possesses global minima near the equilibrium geometries. These global minima in the negative of the probability are connected by a MPP passing through an effective saddle point in -1$(q)12. The location of the effective saddle point (ESP) is most conveniently defined for symmetric reactions in terms of a dividing surface, as follows. By invoking symmetry a 3N - 7 dimensional dividing hypersurface can be defined as a midpoint along all nontotally symmetric coordinates. The point of maximal nuclear probability on this surface represents a saddle point separating the two wells on the -I$(q)12 surface. This saddle point on the negative probability is topologically analogous to the saddle point on the potential energy surface. Specific examples are given below. The MPP is defined as the steepest descent path on the negative of the (14) Miller, W. H.; Handy, N. C.; Adams, J. E. J. Chem. Phys. 1980, 72, 99. (15) Shida, N.; Almlof, J.; Barbara, P. F. Theor. Chim. Acta 1989, 76, 7. Hoffman, D. K.; Nord, R. S.; Rudenberg, K. Theor. Chim. Acta 1986, 69, 265. Panicir, J. Collect. Czech. Chem. Commun. 1975, 40, 1 1 12. Basilevsky, M. V.; Shamov, A. G. Chem. Phys. 1981, 60, 347. Jsrgensen, P.; Jensen, H. J. A,; Helgaker, T. Theor. Chim.Acta 1988, 73, 55. (16) Kuppermann, A.; Adams, J. I.; Truhlar, D. G. In Electronic and Afomic Collisions: Abstracts of Papers; VI11 ICPEAC, Beograd 1973; Kobic, B. C., Kanepa, M. U., Eds.: Institute of Physics: Beograd, 1973. (17) Marcus, R. A.; Coltrin, M. E. J . Chem. Phys. 1977, 67, 2609. (18) A recent example is the least action path: see, e.g.: Garrett, B. C.; Truhlar, D. G. J . Chem. Phys. 1983, 79, 4931.
Tunneling Paths in Intramolecular Proton Transfer probability density surface beginning at the effective saddle point, descending to the global minima (near the equilibrium geometries). A limitation of the steepest descent topological definition of the MPP and IRC is that these paths are not defined past the minima on the surface on the side opposite to the saddle point. It is useful, however, to have a path defined in this region for various types of reduced dimensionally vibrational calculations; see below. In this work we have used the concept of the gradient extremall5to determine the appropriate paths in the region that is not accessible by the steepest descent method. An important practical issue is that quite often, very little probability is found around the ESP. Accordingly, the MPP calculation, which requires at least first derivatives of the nuclear probability distribution function, becomes numerically unstable. The use of the absolute value of the vibrational wave function I\k(q)l is then more practical than using a nuclear probability distribution function. It can be easily shown that the steepest descent path or the gradient extremal path of the nuclear probability distribution function is identical to that of the corresponding wave function itself. Both steepest descent paths and gradient extremal paths are strongly dependent on the choice of coordinate system. We have adopted mass-weighted orthogonal coordinates (mass-weighted Cartesian coordinates in a rigid sense) for MPP calculations analogous to IRC which is a steepest descent path of a potential energy surface in mass-weighted Cartesian coordinates. In this coordination system, we can more directly discuss the relationship between MPP and IRC (see below). The ground vibrational wave function for the systems we examine is highly localized near the reactant and product equilibrium geometries. Since the potential energy surface is approximately harmonic in this region, a rough approximation for I\k(q)I is available in the superposition of normal vibrational wave functions of the reactant and product. Employing this approximation, a rough estimate for the MPP can be obtained as follows: (1) Generate normal vibrational wave functions at the reactant and product equilibrium geometries. (2) Use the sum of the two normal vibrational wave functions as an approximate ground vibrational wave function for the total system. (3) Calculate the MPP using the approximate wave function. We have named this approximated MPP a normal vibrational path (NVP). We explore the deviation of this path from the accurate MPP for various molecular tunneling systems in the following sections. We also employ the NVP as an approximate scheme for making computationally inexpensive estimates for tunneling splittings in polyatomic proton-transfer examples.
3. Intramolecular Proton Transfer in Malonaldehyde In section 2, we defined the MPP in the entire space, using mass-weighted Cartesian coordinates. Here, we present a somewhat different approach, using a “reaction surface” subspace to describe the reaction. The calculational scheme used can be summarized as follows: (1) Three internal coordinates, shown as r l , r2, and r3 in Figure 1 were used as surface variables. (2) The remaining degrees of freedom were treated as dependent variables, adjusted to minimize the potential energy (minimum energy surface). (3) A three-dimensional reaction subspace was determined by (1) and (2). (4) The dynamic couplings to the bath were treated adiabatically (zero-point energy correction). ( 5 ) The dynamic calculations were carried out within the reaction surface Hamiltonian f r a m e ~ o r k . ~ In addition, for the purpose of the MPP calculation, we introduced a new set of coordinates, q, where the kinetic energy expression is orthonormalized:
q is essentially the set of mass-weighted Cartesian coordinates.
(For an extensive discussion of the precise relation between q and r, see Appendix 3 in ref 9.) Figure 3 shows the contour plot of the potential energy surface with the zero-point energy correction, together with the four different reaction paths. q1 and q2 correspond to the mass-scaled parameters r1-r2,and r3 in the diagram
The Journal of Physical Chemistry, Vol. 95, No. 25, 1991 2.0
h
Each contour line :
-
10459
0.002 (a.u.)
1.4 -
Ii d
0.8
-
$ W
s
0.2 -
-0.4 -
-1.0
-1.5
-0.75 q,
0.0 ( a.u. x
0.75
1.5
.IXiiT)
Figure 3. Contour plot of the potential energy surface and the various reaction paths of malonaldehyde in the ortho-normalized mass-weighted internal coordinate system. The coordinates q1 and q2 are discussed in the text (more details are given in ref 9). The four reaction paths correspond to the following: MEP, minimum energy path; MPP, maximal probability path; NVP, normal vibrational path; SLP, straight-line path between the reactant and the product.
above, while q3is adjusted to minimize the potential energy. Due to the difference in mass between the tunneling proton and the oxygens, the angle between the entrance and the exit channels of the potential energy surface in Figure 3 becomes fairly small in this mass-weighted coordinate system. The geometrical saddle point is far from the two equilibrium geometries, and the curvature of the MEP is therefore quite significant. (Note again that the MEP cannot be defined for the direction of the exit channel. For this portion, we have used the gradient extremal path.15) The initial motion of the MEP (from the reactants) is mainly described by q2 (r3),which is interpreted as a heavy motion of the 0-0 stretch. Then, the path sharply curves and the character of the motion changes drastically to q1 (rl-r2),which is the coordinate of the proton movement. The straight line path (SLP) in Figure 3 is a least action path between the reactant and the product, which never contains the heavy motion described by q2. On the SLP, the pure hydrogenic motion drives the proton-transfer reaction. The SLP is the shortest path between the reactant and the product. However, the reaction barrier of 15 kcal/mol is about twice as high as that of the MEP. The eigenfunction of the reaction surface Hamiltonian and the approximated wave function are plotted in Figures 4 and 5 , respectively, together with the MPP or the NVP. In these figures, the variable q3 is adjusted to maximize the probability for each given value of q1 and q2. In Figure 4, the vibrational wave function is highly localized both at the reactant and the product and is fairly harmonic around reactants and products. ESP lies on the point (ql, q2) = (0.0,0.56),far from the geometrical saddle point. In Figure 5 , the shape of the wave function around reactants and products is quite similar to that of the exact wave function in Figure 4. However, the position of ESP shifts some from that of the exact wave function in Figure 4 and it is closer to the geometrical saddle point. The distinction of the four reaction paths is clearly illustrated in Figure 3. The MPP has a definition similar to the EVP (expectation value path) in Figure 6 of our previous work.9 These two paths are, however, quite different. The EVP calculation in the previous work was determined in internal coordinates, whereas the present MPP calculation refers to massweighted internal coordinates. With the obvious exception of the SLP, the initial directions of the three paths are virtually identical and mainly describe the heavy motion of the oxygen atoms. However, the MPP and the NVP take shortcuts from the MEP
-
10460 The Journal of Physical Chemistry, Vol. 95, No. 25, 1991
Shida et al.
+ r2 - r3 - r4 p2 = rl - r2 - r3 + r4 p I = rl
~3
= RI + R2
(4.1) (4.2) (4.3)
Using the above set of coordinates p, the three orthonormalized mass-weighted coordinates q are introduced for the reaction path calculations: 9 = T(P)P
(4.4)
where T(p) is a transformation matrix from p to q. Figure 6 shows the potential energy surface in q1 and q2 with the zero-point energy correction and the four different reaction paths. The variable q3 is optimized to minimize the potential energy for any given value of q1 and q2. q1 and q2 correspond to p1 and p3 in 4.1 and 4.2, with appropriate mass scaling factors. Physically, q1 represents the synchronized double proton transfer (light motion) and q2 is the relative motion of the two monomers -1.0 , (heavy motion). Because of the double proton movement, the reduced mass ratio -1.5 -0.75 0.0 0.75 1.5 between the light motion and the heavy motion is about twice that 9, ( a.u. x i . of the ) previous example. Therefore, the skewed angle of the potential energy surface becomes very small and the difference Figure 4. Contour plot of the ground vibrational wave function and the MPP of malonaldehyde. The axis system is identical to that in Figure between the four paths becomes even more striking than in , I . The contours were drawn at equal spacing in arbitrary units of wave malonaldehyde; The MPP deviates significantly from the MEP function amplitude. or the SLP, and the NVP is very close to the MPP in the double proton-transfer reaction of formic acid dimer.
1
TABLE I: Parameters Used To Construct the Model Hamiltonian in Eo 5.5 N m, 1 2 3 4 5 6
2 4 6 0 2 4
0 0 0 2 1 1
-0.035 0.080 -0.010 0.010 -0.060 0.020
TABLE II: Tunneling Splitting, Zero-Point Energy, and the Reaction Barrier Height Obtained for the Two-Dimensional Model System Using Various Calculational Schemes" zero-point barrier ht, tunneling
MEP MPP NVP
splitting, cm-I
energy, cm-'
kcal/mol
0.007 0.044 0.037 0.17 2.1
290.2 308.1 349.5 1667.0 1926.5
9.5 12.7 15.9 22.9 9.5
SLP exact "The degrees of freedom transversal to the reaction path are ignored completely. and quickly change their characters to the hydrogenic (light atom) motion, consistent with the model previously proposed for the heavy-light-heavy mass combination ~ y s t e m s . ~ J ~ 4. Double Proton Transfer in Formic Acid Dimer Our second example is the double proton transfer reaction in the formic acid dimer (FAD; see Figure 2). Details of this calculation were reported previously.1° The procedure employed in the study of FAD is essentially equivalent to that used for malonaldehyde. For the potential energy surface calculation, high level ab initio SCF and MCPF calculations were carried out. The surface variables used are symmetry-adapted internal coordinates (p) as follows:
5. Vibrational Paths in a Two-DimensionalModel System We have calculated the vibrational wave functions in a twodimensional model system within the reaction path Hamiltonian f r a m e w ~ r k ~to~assess * ' ~ the differences between the different paths. The model Hamiltonian we have chosen is 2
a(qi9q2)= -XEd2/dq: i= I
N
+ i=ECjq1"q2" I
(5.1)
where the parameters C,, li, and mi are listed in Table I. These parameters were chosen to approximately reproduce the general features of the potential energy surfaces of malonaldehyde and FAD in our previous investigations. A contour plot of this model potential is shown in Figure 7,together with the MPP, the NVP, the MEP and the SLP. The contour plots of the wave functions used to evaluate the MPP and the NVP are shown in Figures 8 and 9 respectively. On the MEP and the SLP, very little amplitude is found around the midpoint of the reaction (Le., q1 = 0). The shape of the approximate wave function in Figure 9 is similar to the exact wave function in Figure 8, especially in the vicinity of the reactant or product. We have calculated the lowest pair of the eigenstates in this model system. The results are summarized in Table 11. The tunneling splitting, the zero-point vibrational energy, and the classical barrier height of this system are 2.1 cm-', 1926.5 cm-', and 9.5 kcal/mol, respectively. In Table 11, the one-dimensional calculations on the MPP, the NVP, the MEP, and the SLP are also shown. In these calculations, the degrees of freedom transversal to the reaction paths are completely ignored. In these four different one-dimensional calculations, the SLP gives the closest prediction of the exact tunneling splitting. However, even in that case, the predicted tunneling splitting is 1 order of magnitude smaller than the "exact" value. The effect of the bath modes is presented in Table 111. In these calculations, the contribution from the transversal degrees of freedom was included adiabatically to the potential energy on the reaction paths. In a rigid sense, the adiabatic separation cannot be carried out except for a MEP, as discussed in section 4. But, in these calculations, we have neglected the linear terms of the potential energy part with respect to the transversal coordinates. The results in Table 111 show that the tunneling splittings increase dramatically when the effects from the bath are incorporated (see the differences of the tunneling splittings in Tables I1 and 111). The MPP gives the best prediction of the tunneling splitting in these four calculations.
The Journal of Physical Chemistry, Vol. 95, No. 25, 1991 10461
Tunneling Paths in Intramolecular Proton Transfer
Each contour line : 0.002 (a.u.)
2.0
1:
h
2.5
1.4
F
d
ed
0.8
-7
lS9 1.3
i W
e7
0.2 W
i
-0.4
N
0.1
1 /
-1.0
-1.5
0.7
U
-0.75 q,
0.0 ( a.u. x
0.75
1.5
-0.5
dXii'iX )
-1.5
Figure 5. Contour plot of the approximate ground vibrational wave function used to evaluate the NVP of malonaldehyde. The axis system is identical to that in Figure 1. The contours were drawn at equal spacing in arbitrary units of wave function amplitude. Each contour line :
4.5
h
q,
0.0 ( a.u. x
0.75
1.5
.I a.m.u.
Figure 7. Contour plot of the potential energy surface and the four reaction paths in the two-dimensional model system. 2.5-
0.065 (a.u.)
-i
3.5
1.9-
E d
Ii d
7
-0.75
7-
x 1.3-
2.5
i
j
v v
cI
1.5
c1
U
U
0.5
0.7'
0.1-
MEP -0.5 -2.5
- 1.25 q,
-0.5 -
0.0 ( a.u. x
1.25
2.5
GXii.)
-1.5
-0.75
0.0
q1 ( a.u. x
4
0.75
1.5
a.m.u.1
Figure 6. Contour plot of the potential energy surface and the four reaction paths of formic acid dimer in the ortho-normalized massweighted internal coordinate system. The ql and q2 coordinates are defined in the text.
Figure 8. Contour plot of the ground-state vibrational wave function and the MPP fo the two-dimensional model system. The axis system is identical to that in Figure 3. The contours were drawn at equal spacing in arbitrary units of wave function amplitude.
TABLE 111: Tunneling Splitting, Zero-Point Energy, and the Reaction Barrier Height of the Various Calculational Schemes in the Two-Dimensional Model Svstem"
TABLE I V Tunneling Splitting and tbe Zero-Point Energy of the Various Calculational Schemes in the Two-Dimensional Model System"
MEP MPP NVP SLP exact
tunneling splitting, cm-'
zero-point energy, cm-'
0.029 0.78 0.65 0.17 2.1
275.4 284.9 283.6 1667.0 1926.5
barrier ht, kcal/mol 5.7 9.1 12.1 22.9 9.5
MEP MPP NVP SLP exact
tunneling splitting, cm-' 0.093 1.o 0.40 2.2 2.1
zero-point energy, cm-' 1974.4 1919.4 2008.2 1813.7 1926.5
"The degrees of freedom transversal to the reaction path are treated adiabatically.
"The degrees of freedom transversal to the reaction path are treated diabatically.
The diabatic calculations in the two-dimensional system are summarized in Table IV. For these calculations, the reaction path Hamiltonian formalism was emplooyed. When off-diagonal couplings are taken into account, the tunneling splittings change
remarkably from the adiabatic calculations. (One would expect the diabatic calculation on the SLP to be exact for this system, since the bath mode of the SLP contains only the q2 coordinate and the potential surface is completely quadratic for q2. The minor
10462 The Journal of Physical Chemistry, Vol. 95, No.
Shida et al.
i, 1991
2.5
:1
-.:
d
Each contour line : 0.0025 (a.u.)
3.0
2.2
E d
l a 9
7-
x 1.3
7
x 1.4
.:
j
W W
; 0.7
u
w
W
W
0.6
-0.2
-0.5 -1.5
-1.0 -0.75 q,
0.0 ( a.u. x
0.75
1.5
.I a.m.u.1
Figure 9. Contour plot of the approximate ground-state vibrational wave function and the NVP in the two-dimensional model system. The axis system is identical to that in Figure 3. The contours were drawn at equal spacing in arbitrary units of wave function amplitude.
TABLE V Tunneling Splitting, Zero-Point Energy, and the Reaction Barrier Height Calculated from the Adiabatic Pseudo Valley Approximation in the Two-Dimensional Model Svstem“
MPP NVP
tunneling splitting, cm-I
zero-point energy, cm-I
barrier ht, kcal/mol
0.92 (0.78) 0.63 (0.65)
276.4 (284.4) 230.2 (283.6)
8.9 (9.1) 11.9 (12.1)
“The results in parentheses are from the adiabatic calculations.
discrepancy actually found between “SLP” and “EXACT” is due to the use of a Taylor’s expansion to solve for the eigenstates of the reaction path Hamiltonian.) The incorporation of the diabatic effects improves the tunneling splittings from the adiabatic level of the predictions, except for the NVP. Even at the diabatic level, the MEP gives a poor prediction of the tunneling splitting. Pseudo Valley Approximation. An actual reaction is promoted by the delicate balance of a potential energy and a kinetic energy, and, in general, a MPP does deviate significantly from the floor valley on a potential energy surface. However, it always passes through the ridge of a nuclear wave function because of its definition. If we assume a virtual harmonic valley of a potential energy surface where a MPP lies on the floor, the curvature of the pseudo valley may be determined so as to reproduce the original nuclear wave function. The actual procedure we have employed is as follows: (1) Calculate the Hessian of the nuclear wave function, (8+(q)/d92)/+(q) at selected points of the MPP. Basically, this property corresponds to the force constant matrix of the virtual potential energy. (2). Project out the tangent component of the MPP from the Hessian. (3). Diagonalize the projected Hessian matrix and determine the curvature of the local normal modes bounded in the pseudo harmonic valley. Our pseudo valley approximation is a local harmonic approximation determined from a nuclear wave function. In this pseudo potential valley, the MPP always lies on the floor and the adiabatic separation or the reaction path Hamiltonian method7 can be justified as in the case of IRC. It is noted again that in our pseudo valley approximation, the first and second derivatives of a potential energy surface are no longer required, but the derivatives of the wave function are used, instead. Once a proper wave function is obtained, the procedure to calculated the pseudo valley is quite straightforward with no extra calculational effort, such as the derivative calculations of a potential energy surface. The results of our pseudo valley approximation are demonstrated in Table V. In these calculations, the pseudo valleys are
-2.0
-1.0
0.0
1 .o
2.0
9, ( a.u. x .IXiiK.) Figure 10. Contour plot of the potential energy surface and the four
reaction paths in the proton transfer of H2N-H---NH2 compound.
treated adiabatically. For convenience, the results of our adiabatic calculations are given in parentheses (these were already presented in Table 111). Somewhat surprisingly, the calculated properties are very close to the usual adiabatic results. Once a NVP is justified as a reasonable reference path, the pseudo valley can be easily obtained from the corresponding approximate wave function, with no extra calculation. 6. Proton-Transfer Reaction of H2N-H-*NH2-
The vibrational calculations discussed in the preceding section were for a simple two-dimensional model system. To assess the reliability of our concept in realistic systems, we have carried out the same type of calculations on the proton transfer reaction in the H2N-H-NH2- system:
For this purpose, we have first calculated the two-dimensional potential energy surface on a reaction surface. To define the reaction surface, we have employed the minimum energy surface criterion using symmetry-adapted internal coordinates ( p ) as follows:
(6.lb) where ml and m2 are the appropriate mass scaling factors. For the calculation of the potential energy surface, we used the ab initio SCF method with constrained geometry optimizations. We used the MIDI-4 basis set by Tatewaki et al.I9 augmented by p type polarization functions (orbital exponent = 0.68) on the moving proton. The vibrational calculation of this system was carried out within the reaction surface Hamiltonian framework at the adiabatic level of approximation. Figure 10 shows the contour plot of the potential energy surface of the H2N-H-NH2- compound on the reaction surface with the zero-point energy correction from the bath, together with the four reaction paths. The classical reaction barrier height on this potential energy surface is 5.0 kcal/mol, which is relatively small compared to other examples discussed so far. This will thus (19) Tatewaki, H.; Huzinaga, S. J . Comput. Chem. 1980, 3, 205.
Tunneling Paths in Intramolecular Proton Transfer
The Journal of Physical Chemistry, Vol, 95, No. 25, 1991 10463 3.0
TABLE VI: Tunneling Splitting, Zero-Point Energy, and the Reaction Barrier Height of the Various Calcrilational Schemes in H.N-H--.NH."
splitting, cm-I
zero-point energy, cm-'
barrier ht, kcal/mol
0.1 0.3 0.1 0.7 35.5
112.0 116.1 112.0 1418.3 1509.1
5.0 5.9 5.2 15.9 5.0
tunneling
MEP MPP NVP SLP exact
"The degrees of freedom transversal to the reaction path are ignored completely.
:1
d 2*2
5
x 1.4
i
# v
0.6
vw
-Os21 , 0 . 1( a.u. x
9,
GXii.)
Figure 12. Contour plot of an approximated ground vibrational wave
function obtained from the normal vibrational wave functions and the NVP in the proton transfer of H2N-H--NH2 compound. The axis system is identical to that in Figure 8. The contours were drawn at equal spacing in arbitrary units.
v
a" -0.2 O I
TABLE VII: Tunneling Splitting, Zero-Point Energy, and the Reaction Barrier Height of the Various Calculational Schemes in HZN-H--*.NH~"
-1.0
-1.0 0.0 1.0 2.0 9, ( a.u. x d-). Figure 11. Contour plot of the ground vibrational wave function and the MPP in the proton transfer of H2N-H---NH2 compound. The axis system is identical to that in Figure 8. The contours were drawn at equal spacing in arbitrary units. -2.0
MEP MPP NVP SLP exact
tunneling splitting, cm-'
zero-point energy, cm-I
barrier ht, kcal/mol
2.1 10.9 2.1 0.8 35.5
116.4 121.8 70.3 1405.1 1509.1
1.8 2.9 2.0 15.5
5.0
" The degrees of freedom transversal to the reaction path are treated adiabatically.
provide a critical test of our concept, since the tunneling effect is expected to be large in this case and it may be difficult to determine the appropriate tunneling path or to describe the dynamics on that reaction path. The calculated tunneling splitting, the zero-point energy correction, and the reaction barrier height are summarized in Table VI. For the MPP, thc NVP, and SLP, and the MEP calculations in Table VI, the bath modes are completely ignored. The tunneling splitting in the two-dimensional space ("EXACT") is relatively large, 35.5 cm-I. All the approximate treatments perform uniformly poorly. Figure 10 demonstrates that the MPP does not deviate strongly from the MEP, as compared with the previous examples (see Figures 3 , 6 , and 7). The NVP is approximately halfway between the MEP and the MPP. The exact and the approximate wave functions are plotted in Figures 11 and 12, together with the MPP and the NVP, respectively. Since the potential energy surface is so flat, the minima corresponding to reactants and products in Figure 11 are somewhat shifted from the equilibrium geometries in the direction of the geometrical saddle point. The wave function in Figure 12 does not closely resemble the exact wave function in Figure 11, presumably because the nuclear Hessian matrix at equilibrium does not contain sufficient information to construct the proper wave function. The adiabatic contribution from the bath mode is presented in Table VII. When the adiabatic correction is included, the tunneling splittings increase remarkably and are closer to the exact value of 35.5 cm-', except for the SLP. Especially for the MPP, the tunneling splitting increases by 2 orders of magnitude from the adiabatic level of the prediction. The results from the diabatic calculations are summarized in Table VIII. When the transversal
TABLE VIII: Tunneling Splitting and the Zero-Point Energy of the Various Calculational Schemes in H2N-H--.NH2"
tunneling splitting, cm-'
zero-point energy, cm-'
2.3 188.3 7.0 146.7 35.5
1814.7 1251.4 1560.8 1086.1 1509.1
MEP MPP NVP SLP exact
"The degrees of freedom transversal to the reaction path are treated diabatically. TABLE I X Tunneling Splitting, Zero-Point Energy, and the Reaction Barrier Height Calculated from the Adiabatic Pseudo Valley Approximation for H2N-H.-NH
