Abstract
The quantum-mechanical tunneling is often important in low-energy reactions, which involve motion of light nuclei, occurring in condensed phase. The potential energy profile for such processes is typically represented as a double-well potential along the reaction coordinate. In a potential of this type defining reaction probabilities, rigorously formulated only for unbound potentials in terms of the scattering states with incoming/outgoing scattering boundary conditions, becomes ambiguous. Based on the analysis of a rectangular double-well potential, a modified expression for the reaction probabilities and rate constants suitable for arbitrary double- (or multiple-) well potentials is developed with the goal of quantifying tunneling. The proposed definition involves energy eigenstates of the bound potential and exact quantum-mechanical transmission probability through the barrier region of the corresponding scattering potential. Applications are given for several model systems, including proton transfer in a HO–H–CH3 model, and the differences between the quantum-mechanical and quasiclassical tunneling probabilities are examined.
1. Introduction
There is great interest in understanding the role of quantum-mechanical (QM) effects on reactivity at low temperature in complex systems, such as liquids [1], proteins [2], at surfaces [3], and so forth. One example, which motivated many theoretical investigations [4], is the unusually high experimental kinetic isotope effect () of the transferring hydrogen/deuterium substitution on the reaction rate constant in soybean lipoxygenase-1 [5], attributed to QM tunneling as the dominant reaction mechanism. A few other examples, where tunneling defines the reaction rate are isomerization of methylhydroxycarbene [6] and reaction of hydroxyl radical with methanol [7]. To understand the role of tunneling for a specific system a typical approach is to construct an electronic potential energy profile of a double-well character, along one-dimensional reaction coordinate, and to compute tunneling probabilities using the quasiclassical Wentzel-Kramers-Brillouin (WKB) approximation. This was done, for example, in recent studies of the barrier shape effect on the classical and QM contributions to reactivity in enzyme catalysis [8–10] and in a study of the tunneling control of reactivity in carbenes [6]. Our goal here is to define and analyze QM reaction probabilities and thermal reaction rate constants in the bound potentials in a more rigorous way than the WKB theory. There are two QM effects that should be considered: (i) decay of the wavefunction inside a barrier for a range of energies, treated approximately in the WKB approach, and (ii) quantization of the energy levels due to the bound character of a potential. (The latter is not included in the WKB except via the harmonic zero-point energy of the reactant well.)
The ultimate information about a reactive process in gas phase is contained in the scattering matrix defined for a given collision energy . Its elements , subscripts and label specific internal states of reactants and products, are defined in terms of the energy eigenstates with incoming/outgoing boundary conditions in the reaction coordinate for a potential which vanishes at large reactant-product separation [11]. From the scattering matrix the cumulative reaction probability, , may be computed as Furthermore, determines the thermal reaction rate constant as the Boltzmann average over the collision energy [12]: where is the inverse temperature, . (The atomic units are used throughout.) The cumulative reaction probability can also be equivalently determined from expressions involving correlation functions of the flux operator eigenstates evolved in time [12]. These flux-based approaches naturally connect to quantum and classical transition state theories (beyond the scope of this paper).
Reactions in condensed phase are typically described by the potential energy profile along the reaction coordinate of a double-well character. Therefore, definitions of probabilities or rate constants related to scattering eigenstates are inapplicable. The flux-based approaches involving time-correlation functions are not formally guaranteed to reach a plateau (and do not reach it in one-dimensional models) at long times without inclusion of other degrees of freedom. In the full-dimensional treatment the other “bath” degrees of freedom are responsible for the wavefunction decoherence and for the energy dissipation from the reactive degree of freedom. The plateau behavior of correlation functions is required for a rigorous definition of thermal reaction rate constants. It has been argued, though not formally proven, that already in two degrees of freedom oscillatory behavior of correlation functions is sufficiently damped with time to make the flux-based approaches work [13]. Of course, high-dimensional models give more realistic representation of reactions in condensed phase but are often impractical. If the coupling of the reaction coordinate to the bath degrees of freedom is weak, or the bath degrees of freedom are not involved in rearrangement of atoms, then a reliable estimate of quantum tunneling probability as a characteristic of the barrier and of the potential energy profile along the reaction coordinate is important in itself. Given that significance of QM tunneling is the highest in the regime of just a few reactant/product eigenstates contributing to reactivity, a simple expression for the QM tunneling probabilities, which can be implemented by solving the time-independent Schrodinger equation, is useful. In the time domain, treating a chemical reaction as an event proceeding from reactants to products, rather than the transition state flux, can be useful as well, especially if the potential depends on time, for example, mimicking the effect of substrate motion in reactions on surfaces [14].
For bound potentials the QM reaction probabilities and rate constants are usually estimated using the transition state theory [15] or using the quasiclassical WKB expression for the tunneling probability evaluated at the zero-point energy of the reactant well [16]: The transition state estimate is accurate for a single barrier at energies comparable to its height, while the WKB expression is accurate at low energies, when transmission from products to reactants is justifiably neglected. However, the product region and a possibility of a reverse reaction should be incorporated for a more rigorous treatment; the energy level splittings, measured spectroscopically, are due to tunneling and depend on all regions of the potential.
Even though the approximate treatments of tunneling can be surprisingly accurate, the search for more conceptually rigorous expressions led to the development of less traditional approaches, such as the stationary state decomposition for quadratic potentials [17] via the quantum Hamilton-Jacobi equation adapted for energy eigenstates [18, 19], the exterior complex scaling method [20], and, more recently, the quantum instanton theory [21], ring polymer dynamics [22], and perturbative “beyond instantons” correction to the WKB approximation [23]. The quantum tunneling is known to be very sensitive to the barrier shape and energy, as shown, for example, in the context of enzyme catalysis by Hay and coworkers [9]. Thus, to assess the importance of QM tunneling we propose an adaptation of definitions of the QM tunneling probability and reaction rate constant based on exact energy levels of the bound system and exact QM tunneling through the barrier of the corresponding scattering system. The expressions are obtained from the analysis of a rectangular double-well potential (Section 2) and its unambiguous mapping onto a scattering problem. (For general potentials there is some ambiguity associated with the definition of the reactant and product regions.) After model applications to rectangular and piecewise quadratic potential, we use the approach to compute KIE on the reaction rate constant of proton/deuteron transfer in a model HO–H–CH3 along a linear reaction path with realistic potential (Section 3). Section 4 concludes.
2. Formalism
The modified expression for the reaction rate constants is based on the common functional form of the energy eigenstates for the open scattering system and the double-well system (or any sequence of wells and barriers) consisting of rectangular potential steps. Let us consider a one-dimensional open scattering system with a potential being constant in the asymptotic regions as and arbitrary in the interaction region, referred to as the “barrier” region: Describing scattering in the energy domain, the energy eigenstates with appropriate boundary conditions, such as a wavefunction incoming from the left in the direction from reactants to products, have the following asymptotic form for ( is the particle mass): The wavefunction in the barrier region is an unspecified solution to the one-dimensional Schrodinger equation which smoothly connects solutions in the two asymptotic regions and does not explicitly enter the expression for the transmission probability : For the double-well potential the energy eigenstates with scattering boundary conditions [11] are undefined making (5) and (7) inapplicable. Weiner [17, 24] extended the concept of scattering states to a piecewise quadratic double-well potential, defining the incoming and outgoing boundary conditions in the parabolic wells. These asymptotically diverging solutions were determined from the quantum Hamilton-Jacobi equation [25], resulting in impractical implementation besides being limited by the functional form of the potential.
Dynamics: the probability of finding a particle, initially localized in the reactant well, on the product side oscillates in time without converging to a well-defined value, unless energy dissipation and interaction of the reactive systems with the environment are fully included.
We will consider the rectangular double-well potential, which is easy to analyze in terms of the plane-wave solutions. For such a potential , shown in Figure 1, the energy eigenstates and the corresponding transmission probabilities for the open scattering system defined by (4) and for the full double-well potential of (8) are the same at the energy levels of the full potential in the region of finite . Another feature of this bound potential is that the reactant and product regions are uniquely defined. Thus, the expression for the thermal reaction rates of (2) can be modified (i) using the transmission probability in the open system having the same barrier region as the double-well potential and (ii) using the energy eigenstates of the double-well potential: The summation goes over the energy eigenstates of the bound system: and is the probability of finding a particle on the reactant side for each state : Operator projects a wavefunction onto the reactant region defined as .
In case of the rectangular double well of (8) the projections of the eigenstates on the reactant well are unambiguous (): Generally, a definition of requires specification of the reactant region. Some possibilities are as follows: (i) is a position of the barrier top for any and (ii) is a classical turning point; that is, is a root of in the reactant well. The first option is reasonable for simple barriers with a single maximum, while the second one is more general and consistent with the WKB tunneling expression.
To summarize the formalism, for general potentials with reactant and product well regions the transmission probability is defined by the properties of the open scattering system, that is by the barrier, with the set to constants beyond the reactant/product well minima, and , respectively. This is shown in Figure 2 for a piecewise quadratic continuous potential. The double-well potential simply specifies which energies out of the continuum of the scattering energy states to include into the Boltzmann averaging in (9). For an asymmetric potential, the reactant well may have eigenstates below the asymptote of the product well. Those eigenstates should be included into the Boltzmann averaging in (9) with the zero transmission probability. The same arguments and expressions are applicable to more complicated bound potentials, as long as the reactant and product regions are specified.
3. Examples and Discussion
In this section we examine reaction probabilities and rate constants computed from (9) for several model systems with the emphasis on comparison of the QM and WKB tunneling probabilities. We also consider better-than-WKB approximations relevant to each model: a “minimalistic” two-state description for a symmetric rectangular well (Section 3.1), a mixed QC/QM calculation for an asymmetric well (Section 3.2), and a parabolic barrier approximation for a piecewise quadratic potential (Section 3.3). Accuracy of transmission probabilities and the energy level splittings following from the quasiclassical WKB expression is examined. Example of the KIE calculation for a model proton transfer system is described in Section 3.4.
3.1. Comparison of the QM, WKB, and Two-State Model Probabilities for a Symmetric Double Well
Let us apply (9) to a symmetric rectangular double well and compare reaction probabilities computed using the exact QM and approximate energy levels and transmission probabilities. For a general rectangular barrier, given by (8) without the infinitely high walls, the transmission probability is In the expression above , and is the barrier width; and are given in (5). Analytical expression for QM transmission probability allows high precision calculations performed with Maple [26]. The particle mass is a. u.
We are considering a symmetric barrier of fixed height and variable width. The width of the wells , , is fixed at . The barrier height is set to . The barrier width is in the range of . In the limit of infinitely wide barrier the reactant and product wells have two bound states. The transmission probabilities, , and the energy level splitting, , associated with both asymptotic energy levels are computed by considering two lowest pairs of energy eigenstates of the finite-width potentials. Exact QM transmission probabilities are compared with those obtained using the quasiclassical WKB approach (labeled QC) and the two-state description (labeled ) of this system. The WKB tunneling probabilities are used in comparison, even though the WKB method is inaccurate for potentials with discontinuous derivatives, because of its routine use in tunneling calculations for chemical systems.
In a symmetric potential the ground and first excited states have even/odd symmetry, so the half-sum/half-difference of the two states represents a particle located on the reactant/product side. Therefore, for consistent comparison of exact and approximate expressions the QM transmission is defined as to account for localization of the “initial” wavefunction in the reactant well. The difference between (9) and (14) is that the latter neglects the effect of temperature on the occupation of and states. This omission is accurate in the limit of small energy level splittings. Equation (14) is also consistent with the two-state representation of a wavefunction in the basis formed by the ground eigenstates of reactant and product wells, and , respectively. These functions are the lowest energy eigenfunctions of the reactant and product wells in the limit of an infinitely wide barrier. (The eigenfunctions of the isolated wells with infinitely high walls cannot be used for a rectangular potential, because such eigenfunctions will have zero overlap.) In the two-state representation the energies of the ground and excited states are the generalized eigenvalues of the 2-by-2 Hamiltonian matrix with the overlap matrix . For a symmetric potential the matrix elements are The energy level splitting is .
Following [24], we use the quasiclassical relation between the energy level splitting and the tunneling probability: with the ground state energy used instead of the frequency of the harmonic reactant well. In the two-state representation, the energy level splitting, , following from (17) yields an estimate for the tunneling . The same (18) gives estimates of the quasiclassical energy level splitting, , once the quasiclassical WKB tunneling probability, , is computed from (3). Two pairs of states of the full potential , and , are analyzed; the higher energy states, and , are described in the basis of the first excited states of the reactant and product wells.
Table 1 lists the ground and first excited energy levels obtained exactly, using quasiclassical WKB approximation and within the two-state approximation labeled , , and respectively, are listed in the upper half of Table 1. The same quantities for the next highest pair of levels () are given in the lower half of Table 1. In the quasiclassical WKB treatment the energy levels are defined by the width of the reactant well and, therefore, do not depend on the barrier width. The ground state energy in the two-state approximation is remarkably accurate even for narrow barriers, while the accuracy of level deteriorates for smaller barrier widths. The ratios of the quasiclassical and two-state estimates of tunneling probability to the exact QM values as a function of the barrier width are plotted in Figure 3 for the lowest pair of states. The agreement between and is excellent for the barrier width . The quasiclassical WKB results underestimate the tunneling probability by a factor of . This discrepancy can be understood by comparing quasiclassical WKB and QM tunneling at the energies of the asymptotic eigenstates, shown in Figure 4 for the barrier width . The transmission probabilities at these energies differ by a factor of . For a wider well of the width (not shown here) the ground state energy happens to be near the intersection of the QM and quasiclassical WKB tunneling curves at , which coincidentally yields much better agreement than the results shown in Figure 3. Therefore, the accuracy of the quasiclassical results depends on the reactant well width, even if the accuracy of the quasiclassical tunneling probability does not. At very low energies the discrepancy between the QM and QC tunneling is in orders of magnitude, but fortunately the energy regime below the zero-point energy of reactants does not contribute to the tunneling in a double-well system.
To summarize, the two-state tunneling probabilities for the lowest pair of eigenstates are remarkably accurate except for very narrow barriers and may provide useful estimates if the rest of the states are of much higher energy. The barrier can be considered sufficiently wide when the overlap of the reactant and product well eigenstates, , is less than a few percent (5% percent in the current example). The QC expression, (18), relating the energy level splitting and tunneling is inaccurate for the higher energy pairs of “split” levels. The accuracy of the QC approximation depends on both the accuracy of QC tunneling probability and discrepancy between the energy levels of the full potential and of the asymptotic reactant well. For one-dimensional systems the use of exact energy eigenstates of the full potential rather than those of an isolated reactant well in (9) is more appropriate even with approximate calculation of the tunneling probability.
3.2. Asymmetric Rectangular Barrier
Now let us consider an asymmetric rectangular double-well potential, given by (8) and sketched in Figure 1. There are 7 states with energies below the barrier top. The ground and first excited states are localized on the reactant and product sides, respectively. The eigenstates of the full potential with even quantum numbers correlate with the eigenstates of the reactant well; the difference in energies of the full and asymptotic potentials increases with the quantum number. The parameter values for the potential are , , and . The well width is . The particle mass is a. u.
Exact QM energy levels and eigenfunctions, needed to compute projections on the reactant well, are obtained as eigenvalues and eigenvectors of the Hamiltonian matrix in the Colbert-Miller discrete variable representation [27].
Exact QM transmission probability is given by (13). Quasiclassical transmission probability is defined by (3) and set to (1) for energies above the barrier top. The exact and approximate rate constants computed from (9) are listed in Table 2 for temperatures measured in units of the barrier top . Results of the fully QM, fully QC, and mixed QC/QM (QC transmission evaluated at exact QM energy levels) descriptions are shown in Figure 5 and in Table 2. At the lowest temperature, equivalent to 1% of the barrier height, the QC rate constant is 93 times lower than the QM result, while the mixed QC/QM estimate is only 2 times lower. At high temperatures, when the high energy eigenstates significantly contribute to the reaction rate constant, the discrepancy between exact and approximate rate constants becomes smaller, with the mixed QC/QM estimate being closer to the exact result at all energies. One concludes that, generally, the accuracy of depends both on the accuracy of QC transmission probability and on the difference between exact and asymptotic-well energy levels. Using exact QM energy levels and functions in conjunction with quasiclassical transmission probabilities is more accurate than the fully QC description.
Note, that for the rectangular potential considered here, for is inaccurate because QC approximation does not describe above-the-barrier reflection leading to oscillations of . For general potentials the reaction probability may be computed using, for example, the time-dependent QM wave packet approach [28] used in the next example, yielding reaction probabilities for a range of energies from a single calculation, or with a time-independent method for a few energy values.
3.3. Piecewise Quadratic Potential
Next, we examine a piecewise quadratic potential sketched in Figure 2. The reactant well and barrier are quadratic functions in . The reaction rate constants are obtained using the QM and quasiclassical WKB expressions and using the analytical QM transmission: through a parabolic barrier. Equation (19) is more accurate than the QC expression near the barrier top [16]. A continuous piecewise quadratic potential is defined as follows: with the proper choice of parameters—, , , and —the potential, and its first derivative is continuous function at and . Here a symmetric double well, , is considered; the remaining parameter values are , , , and . The particle mass is a. u. (If the Hamiltonian were scaled by the mass of the proton, would be equivalent to 0.3 eV or 7 kcal/mol.) The barrier parameters are chosen such that there are only few energy levels under the barrier top, which is typical for a chemical reaction at low temperature.
The time-dependent wave packet correlation approach [28] has been used to calculate QM transmission probability, entering (9), since we examine dependence of the reaction rate on the accuracy of the energy levels used in (9). The scattering matrix element describes transmission from reactants to products: The subscripts refer to reactant/product reaction channels and refers to the asymptotically incoming/outgoing wave relative to the barrier. An incoming wave packet , taken as a Gaussian function in the left asymptotic region of , is evolved in time in the “unfolded” potential shown in Figure 2 with a dash. Another Gaussian function, is placed in the product region of the potential. The time-dependent overlap of evolving with the stationary , or correlation function , is computed and Fourier-transformed into the energy domain. The denominator in (21) accounts for the distributions of energy eigenstates in the reactant and product wave packets at time ; these distributions depend on the wave packet localization and kinetic energy: A single wave packet propagation, accomplished on a grid using the split-operator method [29, 30], gives transmission probability for a range of energies represented in the reactant/product wave packets, which is convenient if multiple eigenstates contribute to the rate constant. At low energies it is hard to converge transmission probability using time-dependent dynamics methods, but the energies below the zero-point energy of the bound potential do not contribute to . The exact QM energy levels and projections of eigenfunctions on the reactant well are performed as in Section 3.2. There are five energy levels below the barrier top, but many more states contribute to at higher temperatures.
The energy-resolved exact QM, parabolic barrier approximation, and QC transmission probabilities are shown in Figure 6(a). Tunneling probabilities are underestimated in the QC approximation, while the parabolic barrier expression gives accurate results for energies above (the barrier is in fact parabolic in this energy range). The discrepancy grows at energy below . This energy range, however, is below the zero-point energy and does not contribute to the thermal rate constants shown in Figure 6(b). The approximate rate constants differ from exact QM results by at most 25–30%.
(a)
(b)
3.4. Proton Transfer
As a chemically relevant model, we consider the proton transfer in the HO–H–CH3 system for a constrained, collinear O–H–C geometry. In this model system, the proton is transferred from a donor carbon to acceptor oxygen. The potential energy surface is obtained from density function theory (DFT) electronic structure calculations, in particular at the B3LYP/6-31G(d,p) level of theory [31]. The energies were compared to those obtained at the CCSD(T)/aug-cc-pVDZ theory level for the same geometries and were found in excellent agreement. In this system, collinear donor-proton-acceptor arrangement is a good approximation to the fully optimized reaction path as determined by a set of constrained geometry optimizations. In these optimizations was fixed to a value in the range 0.8–2.0 Å: the average deviation of ∠CHO from linearity was less than 2° and the maximum deviation was 10°.
One-dimensional potential energy surfaces of 40 points were generated as functions of for a fixed distance. Proton transfer on three surfaces for is analyzed below. The surfaces were parametrized as sixth degree polynomials in . The potential energy curves, shown in Figure 7(a), have characteristic asymmetric double-well shape.
(a)
(b)
We examine the effect of donor-acceptor distance, , on the tunneling rate constants, using QM and QC transmission probabilities and exact and approximate energy levels in (9). QM transmission probability is calculated using the wave packet correlation approach outlined in Section 3.3. QC transmission is defined by (3). The calculation of reaction rate constants using either of these approaches requires eigenstate projections on the reactant region. The eigenstates are computed exactly, as outlined in Section 3.2 defining the reactant region to the right of the barrier top.
The results are shown in Figure 7(b) as a function of temperature. At low temperatures, proton rate constants are several orders of magnitude higher than those for the deuteron, and this gap decreases as temperature rises. As a result, the kinetic isotope effect should be largest at low temperatures, and this trend is shown in Figure 7(b). The KIE calculated for each surface remains nearly constant in the low-temperature (0–300 K) region and then begins to approach 1 as the temperatures rises. These very large KIE values point to a reaction dominated by quantum tunneling at low temperatures. The largest discrepancy in the KIE between different surfaces is seen at very low temperatures as well, and by increasing by only 0.1 Å, the KIE is enhanced by an order of magnitude.
When calculating QC rates, the energy levels of the isolated donor well are traditionally used as rather than the energy corresponding to as we have defined it in (11). In the three double-well potentials considered here, the QC ground state energy is lower than that of the QM calculations. Thus, the QC calculations yield lower rate constants by up to a factor of 5 as shown in Table 3. The table compares the ground state contributions to and . Despite the discrepancy in rate constants, the KIE predicted by the QC method is within a factor of 2 of the QM results due to cancellation of error.
4. Conclusions
Reliable estimates of the QM tunneling probabilities through a barrier along the reaction path are often used in studies of reactions proceeding in condensed phase. The formal definition of the QM reaction probability and rate constant based on asymptotic scattering states [11] cannot be used for bound potentials representing such processes. In dynamics of a wave packet representing reactants, this aspect manifests itself as persisting-in-time oscillations of the reaction probabilities. A proposed modification of the QM expression for bound potentials, (9), which addresses this problem, is based on the analysis of a rectangular double-well potential. For this potential the rate constant expression separates into (i) the reactant/product transmission probability through a barrier of a scattering system with the same barrier region as the full potential and (ii) the eigenstate energies and eigenstate projections on the reactant region of the full bound potential. Exact QM or quasiclassical (or other approximate) methods may be used to estimate the tunneling probability and energy levels. For example, for the rectangular double-well potential the two-state representation gave fairly accurate estimates of tunneling probabilities derived from the quasiclassical relationship between the energy level splitting and tunneling (18). The relationship worked for the lowest pair of eigenstates even for narrow barriers but did not hold for higher energy pairs of eigenstates.
The quasiclassical WKB estimates of rate constants are shown to depend on both the accuracy of the transmission probabilities and positions of the energy levels: use of exact QM energy eigenstates is preferred. The QC rate constants are lower at low energy and higher at energies comparable to the barrier top in comparison to QM results. Performing QM scattering calculations to obtain transmission probabilities is more expensive than QC estimates, but it makes the approach generalizable to more than one dimension. For a smooth potential (piecewise quadratic potential) the approximate transmission probabilities were quite accurate. The accuracy of the reaction rate constant in the parabolic approximation to the barrier was better than 4% and better than 40% for the quasiclassical WKB approximation. In all cases we find that the accuracy of rate constants is improved when exact eigenstates are used with approximate probabilities in (9). For the proton transfer model for the HO–H–CH3 system with constrained donor-acceptor distance, the QC approximation gave reasonable estimates of the tunneling: the QC reaction rate constants were approximately 4 times smaller than the exact QM counterparts, and their ratio, the KIE, was within 50% of the exact QM value due to cancellation of errors. While multidimensional dynamics is preferable for a rigorous theoretical study of a reaction in condensed phase, this simple approach of computing reaction probabilities and thermal rate constants in bound potentials may be used to analyze barriers and to assess importance of QM tunneling for a given system.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
This material is based upon work supported by the National Science Foundation under Grant no. CHE-1056188. The authors thank Vitaly Rassolov for active discussions.