4430
J. Phys. Chem. 1996, 100, 4430-4436
Path Integral Calculation of Quantum Nonadiabatic Rates in Model Condensed Phase Reactions Maria Topaler and Nancy Makri* School of Chemical Sciences, UniVersity of Illinois, 505 S. Mathews AVenue, Urbana, Illinois 61801 ReceiVed: June 16, 1995; In Final Form: September 30, 1995X
We present accurate path integral calculations of quantum rate constants for model nonadiabatic reactions in condensed matter. The model is described by two coupled diabatic potential surfaces interacting linearly with a bath of harmonic oscillators. The rate constant is obtained from the time integral of the flux-flux correlation function which is evaluated by the quasi-adiabatic propagator path integral method. We study the dependence of the reaction rate on friction, temperature, and exothermicity and compare with predictions of analytical theories. In particular, we observe a broad golden rule plateau as well as rate enhancement due to quantum resonances for low friction in agreement with the semiclassical analysis of Onuchic and Wolynes.
I. Introduction The rates of chemical reactions can be affected by a number of quantum phenomena such as tunneling, curve crossing, and quantum interference. Reactions in condensed matter such as solids, liquids, glasses, and biomolecules are also influenced by the coupling to the environment. Therefore, an accurate description of reactive processes must incorporate both quantum effects and dissipation. Quantum theories of dissipative rate processes have been significantly advanced during the past years.1-3 Most of them are based on transition state pictures, semiclassical techniques, or the path integral formalism. For common adiabatic dissipative processes the dependence of the quantum reaction rate on friction is generally nonmonotonic and qualitatively similar to that predicted from the classical Kramers theory.1-3,4 The reaction rate increases at small friction (energy-diffusion regime) and eventually falls off with friction when the latter becomes large enough to ensure rapid energy exchange (spatial-diffusion regime). However, the quantum transmission coefficient can be very large compared to its classical value due to tunneling, and the turnover from energy to spatial diffusion regime moves to very small values of friction at low temperature so that it may not be observable.1 In a recent paper5 we presented accurate calculations of quantum mechanical reaction rates for a model adiabatic condensed phase reaction using the quasiadiabatic propagator path integral methodology6-10 to evaluate the quantum reactive flux correlation function.11,12 We presented converged results for the rate constant over a wide range of temperature and friction and in general obtained very good agreement with analytical theories. We also observed an additional interesting quantum effect: in the low temperature-weak damping regime, vibrational coherences caused a significant quantum dynamical enhancement of the rate above the centroid density quantum transition state value.13,14 This effect manifested itself as an upward step structure of the integrated flux correlation function with the step length approximately equal to the vibrational period in the reactant well. Similar enhancement of the rate due to quantum resonances has been predicted earlier by Onuchic and Wolynes15 on the basis of a semiclassical calculation for a model nonadiabatic reaction comprising two intersecting harmonic potential curves. The present numerical study was motivated X
Abstract published in AdVance ACS Abstracts, February 1, 1996.
0022-3654/96/20100-4430$12.00/0
partly by the wealth of intriguing quantum mechanical effects revealed in that semiquantitative analysis of dissipative curvecrossing dynamics. The study of reaction rates for nonadiabatic condensed phase processes poses a very interesting and important problem. Some of the most important chemical reactions cannot be described in terms of the motion on a single Born-Oppenheimer potential surface. In particular, in condensed phase reactions such as electron transfer and sometimes proton transfer the nonadiabatic coupling of electronic surfaces plays a very significant role. A standard picture of such a process involves transitions between diabatic surfaces corresponding to the reactant and product sites induced by coupling between them. The Hamiltonian is
H0 )
(
ps2 V (s) Vd(s) 1+ 1 Vd(s) V2(s) 2m0
)
(1.1)
Here 1 is a unit matrix and s is a nuclear reaction coordinate that may correspond to some important intramolecular mode or may represent some collective coordinate. The potential curves are shown in Figure 1. The diabatic barrier V0 is equal to the energy difference between the crossing point and the reactant minimum, and the exothermicity ∆V is the energy difference between the reactant and the product sites. By diagonalizing the potential energy part of H ˆ 0 one obtains two adiabatic potential energy surfaces,
V((s) ) 1/2[V1(s) + V2(s)] ( /2{[V2(s) - V1(s)]2 + 4Vd2(s)}1/2 (1.2)
1
If the energy gap between the adiabatic surfaces is large compared to the thermal energy kBT then the adiabatic description of the reaction is adequate. When the electronic coupling is relatively small one needs to consider both potential surfaces and the curve crossing description is used. The most dramatic deviation from the Born-Oppenheimer single-surface dynamics occurs in the region of configuration space where two potential energy surfaces nearly cross. The dynamics in the vicinity of this avoided crossing point is strongly quantum mechanical. The pioneering work of Landau, Zener and Stueckelberg16-18 led to the following expression for the hopping probability between two diabatic surfaces: © 1996 American Chemical Society
Calculation of Quantum Nonadiabatic Rates
J. Phys. Chem., Vol. 100, No. 11, 1996 4431
PLZ ) 1 - e-δ
(1.3a)
with
δ)
2πVd2
(1.3b)
pV|V′2(0) - V′1(0)|
where we assume that the avoided crossing takes place at s ) 0 and V is the velocity of the particle at this point. Although eq 1.3 was obtained assuming constant diabatic coupling and linear diabatic surfaces, it often provides a good approximation because nonadiabaticity is only important in the vicinity of the crossing point where these assumptions usually hold. For undamped classical motion along the reaction coordinate the parameter δ determines whether the transition is adiabatic or nonadiabatic. In the absence of dissipation, the populations of the reactant and product sites for the system depicted in Figure 1 will remain almost constant or oscillate, depending on whether or not the corresponding vibrational levels are resonant. Coupling of the reaction coordinate to the environment generally causes irreversible transition from reactants to products described by a rate law. We follow many other authors (see, e.g. refs 19, 20, and 10) and assume that the two electronic states are directly coupled through one nuclear mode which is the reaction coordinate s introduced above. The reaction coordinate is in turn coupled to the remaining intramolecular and solvent degrees of freedom which will be referred to as the environment. In many cases, and in particular when electron transfer processes are considered, it is possible to represent the environment by a collection of harmonic oscillators coupled linearly to the reaction coordinate.21 Deviations from the harmonic picture are expected in processes that are not accurately described by linear response theory and have been observed in classical trajectory simulations.22 With these assumptions the total Hamiltonian takes the common system-bath form
(
Pj2
)
cjs 1 + mjωj Qj + H ) H0 + ∑ 2 j 2mj mω2 j
j
2
(1.4)
Equation 1.4 includes standard counterterms21 quadratic in s in order to ensure that the barrier height is independent of the coupling parameters cj. The reduced dynamics of H0 are influenced by the bath only through the spectral density function:21
J(ω) )
π
cj2
2
j
∑j m ω δ(ω - ωj)
(1.5)
j
A large amount of theoretical work has been devoted to the study of nonadiabatic reaction dynamics in dissipative environments. Zusman has used a stochastic Liouville equation to study the case of a classical reaction coordinate with strong dissipation.23 A qualitative description of nonadiabatic rate constants for classical nuclear motion and in particular their dependence on friction has been given by Fraunfelder and Wolynes and compared to the adiabatic case.24 The only quantum effect introduced there was surface hopping, i.e., electronic tunneling, and thus the approach of Fraunfelder and Wolynes is very close in spirit to numerical surface hopping molecular dynamics methods.25-27 Garg et al.19 and Wolynes20 used the path integral method to generalize the stochastic approach to the quantum mechanical case and obtained adiabaticity criteria similar to
Figure 1. Sketch of potential surfaces for the Hamiltonian of eq 1.1. Diabatic and adiabatic surfaces are shown by solid and dashed lines, respectively. Various parameters mentioned in the text are indicated.
those of Zusman. They also used instanton methods to advance toward the low-temperature regime. A major conclusion of all the above studies is that dissipation causes the reaction to become more adiabatic. Nonadiabatic rates have also been studied using imaginary time path integral methods by Wolynes,28 Chandler and co-workers,29-31 and many others. As we already mentioned at the beginning of this introduction, an additional quantum effect on nonadiabatic rates in condensed phases has been discussed in a semiquantitative semiclassical study by Onuchic and Wolynes15 which explored the interplay of friction, electronic nonadiabaticity, and quantum interference effects. This work led to qualitative arguments for how the transmission coefficient for a nonadiabatic reaction can be affected by resonances between reactant and product states and predicted that in the low damping regime, where dephasing rates are small, these resonance effects may lead to significant increase of the rate above the value predicted from the classical surface hopping model. In this paper we present accurate numerical calculations of quantum rate coefficients for the model nonadiabatic barrier crossing problem defined by eq 1.4. Our approach consists in evaluating the quantum flux correlation function using the quasiadiabatic propagator path integral (QUAPI) methodology developed earlier in our group.6-10 This work is the continuation of our earlier study of quantum mechanical reaction rates for the adiabatic case5 with the emphasis of the present work on those effects that are not observable in single-surface models. In section II we describe briefly the quantum flux correlation function formalism that we adopt for the nonadiabatic barrier crossing and review the quasiadiabatic propagator path integral method. In section III we present our numerical results on thermal rate constants and compare with some of the earlier analytical work. At high temperature where nuclear tunneling is insignificant, we observed in agreement with references15,24 a broad golden rule plateau where the rate is independent of friction and is given approximately by the Landau-Zener result. At lower temperatures this plateau is washed out due to nuclear tunneling effects and the transmission coefficient increases significantly above the classical surface hopping value. In the low damping regime we observe a step structure in the flux correlation function caused by constructive or destructive interference that leads to the resonance effects on the rate constant predicted by Onuchic and Wolynes.15 These effects are very similar to those that we found for adiabatic quantum reaction rates in the low damping-low temperature regime. Finally, some concluding remarks are presented in section IV. II. Methodology A. Flux Correlation Function. Following Miller, Schwartz, and Tromp12 the rate coefficient is given by the time integral
4432 J. Phys. Chem., Vol. 100, No. 11, 1996
Topaler and Makri
of the flux autocorrelation function
k ) Z-1∫0 Cf(t) dt tp
(2.1)
Here Z is the partition function of the reactants, tp is the conventional plateau time,32,33 and the correlation function is given by
Cf(t) ) Tr[FeiHtc*/pFe-iH c/p] t
(2.2)
In eq 2.2 tc ) t - ipβ/2 is a complex time that arises from combining the time evolution operator with the Boltzmann operator, β ) 1/kBT, and F is the symmetrized flux operator that measures the flux between reactants and products through a given dividing surface. This operator can be expressed12 as a commutator between the Hamiltonian of the reactive system and the projection operator Pprod that projects on the product part of configuration space:
F ) (i/p)[H,Pprod]
(2.3)
For an adiabatic reaction this projector is equal to a step function in position space that projects to the right of the surface separating reactants from products. Assuming the latter to be orthogonal to the reaction coordinate at the saddle point s ) 0 one obtains
Pprod )
{
0, s < 0 1, s > 0
any approximations. While this procedure reduces the dimensionality of the problem dramatically, it introduces interactions nonlocal in time into the dynamics of the system such that the latter is no longer solvable by conventional iterative onedimensional techniques. An iterative approach which took advantage of the finite memory length in the dissipative influence functional has been developed recently by Makarov and Makri,36,37 allowing fully quantum mechanical propagation over very long time. For the present purpose of the reaction rate calculations, which require evaluation of the flux correlation function for intermediate time length, we have found the direct global evaluation of the path integral more efficient. Below we review briefly the QUAPI methodology and the numerical techniques that we adopt for the calculation of quantum nonadiabatic reaction rates. The propagators in eq 2.7 are expressed in path integral form by dividing the complex time tc into N time slices of length ∆tc and approximating each short time evolution operator using the quasi-adiabatic propagator splitting6
exp(-iH∆t/p) ≈ exp(-iHenv∆t/2p) exp(-iH0∆t/p) × exp(-iHenv∆t/2p) (2.8) Here H0, given by eq 1.1, is the Hamiltonian for the system dynamics along the vibrationally adiabatic path from the avoided crossing point to the potential minima and defined by Qj ) cjs/ mjωj2, and
Henv ≡ H - H0 ≡ ∑Hj
(2.4)
Using this projector the expressions for the flux operator and for the flux correlation function can be evaluated in the coordinate representation.12,8 For the nonadiabatic reaction with the potential surfaces shown in Figure 1, another form of the projector seems more appropriate. Here we are interested not in the transition through a given dividing surface in nuclear configuration space but rather in transitions between the two electronic states. Denoting by |1〉 and |2〉 the electronic states corresponding to the diabatic surfaces V1(s) and V2(s), respectively, the projection operator Pprod that projects on the potential surface of the products has the form
Pprod ) |2〉〈2|
(2.5)
describe displaced harmonic oscillators of the bath which depend parametrically on the reaction coordinate s. Use of eq 2.8 leads to the following expression for the propagator products that appear in eq 2.7:
Trb[〈σ2N,s2N|eiHtc*/p|σN,sN〉〈σN,sN|e-iHtc/p|σ0,s0〉] )
∫-∞ds1 ... ∫-∞dsN-1∫-∞dsN+1 ... ∫-∞ds2N-1 ∞
2
∞
2
∞
2
〈1,s|F|2,s′〉 ) -〈2,s|Fˆ |1,s′〉 ) -(i/p)Vd(s)δ(s - s′)
(2.6)
so that only vertical transitions are permitted. Then the correlation function is given by the expression
∞
2
∑ ...
σ1)1
2N
∑ ∑ ... σ ∑)1 k)N+1 ∏ 〈σk,sk|eiH ∆t */p|σk-1,sk-1〉 × σ )1 σ )1 0
N-1
N+1
c
2N-1
N
〈σk,sk|e-iH ∆t /p|σk-1,sk-1〉I(s0,s1,...,s2N) ∏ k)2 0
We evaluate eq 2.2 for the two-surface problem in a product basis composed of diabatic electronic states {i〉, (i ) 1, 2)} and position eigenstates {s〉}. The matrix elements of the flux operator that are diagonal in the electronic degree of freedom vanish in this basis. The off-diagonal elements are given by
(2.9)
j
c
(2.10)
Here I(s0,s1,...,s2N) is the influence functional which arises from integrating out the harmonic bath. Notice that the influence functional does not depend on the electronic state index σ because the system-bath coupling has the same form for both electronic surfaces. For this reason the influence functional is formally identical to that of a forced harmonic oscillator38,39 and can be easily evaluated. The resulting expression is8,5 2N+12N+1
I ) ∏Ij0 exp{-(1/p) ∑
Cf(t) ) (2/p )Re∫-∞ds∫-∞ds′ Vd(s)Vd(s′) Trb[〈2,s|e 2
〈1,s′|e
∞
-iHtc/p
∞
iHtc*/p
|2,s〉 + 〈2,s|e
-iHtc/p
|1,s′〉〈1,s′|e
iHtc*/p
|1,s′〉
|1,s〉] (2.7)
where the symbol b stands for all bath degrees of freedom. B. Quasi-Adiabatic Propagator Path Integral (QUAPI) Methods. Feynman’s path integral method34,35 is ideally suited for calculations on problems of the type described by Hamiltonian of eq 1.4 because it allows the harmonic bath degrees of freedom to be integrated out analytically without introducing
j
∑
Bkk′sksk′}
(2.11)
k)0 k′)0
where
Bkk′ )
2 ∞ J(ω) ∫ dω ω2 Rkk′(ω;∆tc,β) π 0
(2.11a)
and the explicit expressions for the functions Rkk′(ω;∆tc,β) were given in refs 5 and 8. The propagators for H0 are evaluated exactly in terms of the numerically obtained eigenstates φm and
Calculation of Quantum Nonadiabatic Rates
J. Phys. Chem., Vol. 100, No. 11, 1996 4433
eigenvalues Em of the system Hamiltonian H06
〈σ′′,s′′|exp(-iH0∆tc/p)|σ′,s′〉 ) M0
∑ φm(s′′,σ′′)φm*(s′,σ′) exp(-iEm∆tc/p)
(2.12)
m)1
This guarantees that expression 2.10 is exact for any time step ∆tc in the absence of coupling between system and bath. For nonzero coupling the maximum time step allowed in eq 2.10 is often constrained by the highest frequency present in the problem. The quasi-adiabatic propagator provides a good reference because it treats accurately the high-frequency bath degrees of freedom which usually correspond to the shortest time scale in the problem. The quasi-adiabatic representation introduces no difficulty into the treatment of slow bath oscillators because the time scale of the latter is much longer than the typical time step. As a result, the QUAPI converges fast over a wide range of bath frequencies and allows large time steps so that the total number of integrations can be kept modest. C. Multidimensional Integration Schemes. Here, as in our previous work,5,7-9 we use both Monte Carlo and quadrature methods for the evaluation of the multidimensional integrals in eqs 2.7 and 2.10. Boltzmann factors in these expressions lead to damping of the integrand so that it becomes localized in a relatively small part of the multidimensional integration space. This damping of the integrand is especially pronounced when the coupling between reaction coordinate and bath is relatively strong. In that case Monte Carlo evaluation of the path integral is efficient for moderately long times even with relatively small time steps necessary for the strong friction regime. The Monte Carlo scheme employs a multidimensional random walk in the space of the path integral variables, which include the reaction coordinate and the electronic state index. The reaction coordinate is sampled on a grid.8 The calculation consists of three stages. In the first stage the absolute value of the integrand is used as the unnormalized sampling function. In the second stage the normalization integral is computed using the absolute value of one of the system propagator products as the sampling function whose normalization is evaluated in the third stage by matrix multiplication. However, as the system-bath coupling is decreased, the performance of Monte Carlo methods deteriorates significantly because the integrand is more delocalized and oscillatory. We have developed a very efficient quadrature scheme that allows evaluation of the path integral in this case.5 It employs a systemspecific discrete variable representation (DVR).40-48 These schemes aim at providing representations that combine the advantages of compactness characteristic of energy basis sets while retaining the attractive features of coordinate states. We have adopted a system-specific potential-optimized DVR47 to construct a highly efficient quadrature suitable for evaluation of the discretized QUAPI. In this context the DVR basis is obtained by performing a unitary transformation on the basis {φi} of M lowest eigenstates of H0 M
|ui〉 ) ∑ Lii′|φi′〉
(2.13)
i′)1
The transformation matrix {Lii′} is chosen to diagonalize the position operator s in the DVR basis:
〈ui|s|ui′〉 ) siDVRδii′
(2.14)
The DVR states obtained in this way are the discrete analog of the ordinary coordinate eigenstates, and the eigenvalues siDVR form the DVR grid. This grid typically consists of only 2-20 points and allows quadrature evaluation of the path integral for times sufficiently long for the reactive flux to reach a plateau.9,5 Equation 2.10 in this discrete variable representation takes the form
Trb[〈σ2N,s2N|eiHtc*/p|σN,sN〉〈σN,sN|e-iHtc/p|σ0,s0〉] ) M
M
M
M
... ∑ ∑ ... ∑ 〈σ2N,s2N|eiH ∆t */p|ui ∑ i )1 i )1 i )1 i )1 0
c
2N-1
1
N-1
N+1
〉×
2N-1
2N-1
∏
k)N+2
〈uik|eiH0∆tc*/p|uik-1〉〈uiN+1|eiH0∆tc*/p|σN,sN〉 × N
〈σN,sN|eiH0∆tc*/p|uiN-1〉 ∏ 〈uik|e-iH0∆tc/p|uik-1〉 × 〈ui1|e
k)2 |σ0,s0〉I(s0,siDVR ,...,siDVR ,sN,siDVR ,...,siDVR ,s ) 1 N-1 N+1 2N-1 2N
iH0∆tc*/p
(2.15) Here the propagators in the DVR basis are obtained from eqs 2.12 and 2.13 and the influence functional is still given by eq 2.11. We have found the DVR evaluation of the QUAPI extremely efficient; in typical situations the number M of grid points required for convergence is only slightly larger than the number M0 of system eigenstates surviving in the complextime propagator, eq 2.12. III. Results and Discussion In this section we apply the above methodology to study quantum reaction rates for transitions between two diabatic potential surfaces interacting with a dissipative environment according to eq 1.4. For simplicity in this study we chose the diabatic surfaces to be harmonic with frequency ω0 and the electronic coupling Vd(s) to be constant; clearly though, our methodology can be applied in the same way to the case of anharmonic potentials and coordinate-dependent coupling. The parameters that we use here are ω0 ) 500 cm-1, Vd ) 0.1 pω0, V0 ) 3 pω0, m0 ) 1836 au, and the exothermicity ∆V varies between zero and 0.5 pω0. The bath is described by an ohmic spectral density with exponential cutoff21
J(ω) ) ηωe-ω/ωc
(3.1)
with ωc ) ω0. The parameter η is the standard friction coefficient appearing in the (generalized) Langevin equation.49 We present our results in the form of transmission coefficients. Here we define the transmission coefficient as the ratio of the rate constant to its “classical” transition state theory value that we approximate by
kTST ≈
1 kBT -βV0 ω0 -βV0 ) e e 2πp Z 2π
(3.2)
It should be noted, though, that the classical rate for transitions between surfaces is strictly zero. The rate constants are computed according to eqs 2.1, 2.7, and 2.15, while the reactant partition functions Z are evaluated using similar (Monte Carlo or DVR quadrature) evaluation of the corresponding QUAPI expression. Figures 2-4 show the transmission coefficient as a function of the dimensionless damping parameter η/m0ω0 for different values of the temperature and the exothermicity ∆V. Figure 2
4434 J. Phys. Chem., Vol. 100, No. 11, 1996
Figure 2. Quantum transmission coefficient as a function of damping strength at relatively high temperature kBT ) 0.64pω0. (a, top) Golden rule, plateau. Circles and triangles correspond to results for ∆V ) 0 and 0.5pω0, respectively. (b, bottom) Blow-up of the small damping region, showing also intermediate values of exothermicity. Circles, ∆V ) 0; inverted triangles, ∆V ) 0.1pω0; squares, ∆V ) 0.2pω0; rhombi, ∆V ) 0.3pω0; filled triangles, ∆V ) 0.5pω0.
Topaler and Makri diabatic surfaces until it approaches the crossing region. During each passage through this region the particle can change its electronic state with probability given by the Landau-Zener expression, eq 1.3. The value of friction, η/m0ω0, determines the average number, Nc, of recrossings of the transition state region that the particle makes before it loses enough energy to the bath and settles in one of the wells, and the relevant adiabaticity parameter is PLZNc. For very low or very high friction the number of recrossings is very large, such that PLZNc > 1 and many electronic transitions are made. In this case the electronic state is randomized during each sequence of recrossings and the transmission coefficient κ ) 1/2Nc is independent of the electronic coupling, i.e., is equal to that for the adiabatic reaction, increasing with η for low friction and decreasing in the strong friction regime. However, when the number of recrossings is small or moderate, i.e., PLZNc , 1, then for each sequence of recrossings the probability of the reaction is proportional to the Landau-Zener probability times the number of attempts. This leads to κ ) 2PLZ which is independent of damping strength and smaller than the corresponding adiabatic value. These qualitative arguments have been confirmed by the stochastic surface hopping calculation of Cline and Wolynes.50 A simple interpolation formular between these two regimes has been proposed by Straub and Berne51,15
1 1 + 2Nc ≈ κ 2PLZ
(3.3)
In agreement with these results one observes in Figure 2 a golden rule plateau, i.e., a wide range of friction spanning several orders of magnitude where the transmission coefficient is essentially constant. This behavior is different from that for an adiabatic reaction with similar parameters where the rate has a maximum at η/m0ω0 ≈ 1/βV0. For the nonadiabatic reaction the transmission coefficient starts decreasing only at very large values of friction. Estimating the hopping probability using eq 1.3 and average thermal velocity
ν ) (2kBT/πm0)1/2 Figure 3. Quantum transmission coefficient as a function of damping strength at intermediate temperature kBT ) 0.32pω0. Circles and triangles correspond to results for ∆V ) 0 and 0.5pω0, respectively.
Figure 4. Quantum transmission coefficient as a function of damping strength at low temperature kBT ) 0.16pω0. Circles and triangles correspond to results for ∆V ) 0 and 0.5pω0, respectively.
corresponds to high temperature (pω0 ) 0.64kBT) where nuclear tunneling does not play a significant role and the only important quantum effect appears to be electronic curve crossing. The major contribution to the reaction rate comes from the vibrational levels that are close in energy to the crossing point. A simple qualitative description of nonadiabatic reaction rates in this regime has been given by Fraunfelder and Wolynes.24 Their argument is based on a quasiclassical surface hopping picture in which the particle moves classically on each of the
we obtain 2PLZ ≈ 0.04 which is in very good agreement with the plateau value κ ≈ 0.05 of the numerical results displayed in Figure 2. The small difference between numerical and analytic estimates of κ can be attributed partly to nuclear tunneling. When the temperature is lowered, kBT , pω0, nuclear tunneling along the reaction coordinate begins to play an important role. In this regime the main contribution to the reaction rate comes from vibrational levels that lie below the crossing point and the difference between nonadiabatic and adiabatic reactions becomes less pronounced. In the limit of very low temperature only the two lowest vibrational levels and their tunneling splitting will determined the reaction rate, and the same formal description will be valid both for adiabatic and nonadiabatic situation. In Figures 3 and 4 we show the transmission coefficient as a function of friction for kBT ) 0.32pω0 and kBT ) 0.16pω0, respectively. Indeed, we can observe in Figure 3 that the golden rule plateau is now narrower and in Figure 4 it essentially disappears. The transmission coefficient increases very rapidly when the temperature is lowered and becomes much larger than the classical hopping value of 2PLZ. Also shown in Figures 2-4 are rates for different values of exothermicity. At moderate to strong friction and high temperature the rate is essentially independent of ∆V, in agreement with eq 3.3. Indeed, all the dynamics here takes place in the
Calculation of Quantum Nonadiabatic Rates
J. Phys. Chem., Vol. 100, No. 11, 1996 4435
vicinity of the crossing point where the potential is very little influenced by the value of ∆V. For lower temperatures and strong friction one can observe a small dependence of the rate on the value of ∆V, the rate being lower for the symmetric ∆V ) 0 case. This has a simple explanation: the exothermicity is varied in our model, such that the barrier height V0 is kept constant but the width of the barrier at the same energies becomes smaller for the exothermic reaction compared to the symmetric case. Because at low temperatures the rate is largely determined by tunneling (and the tunneling rate is smaller for wider barrier), this difference in the barrier width leads to the difference in the reaction rates which increases with decreasing temperature because of the increasing role of tunneling. However, for small damping the rate in the symmetric potential is significantly larger than that in the asymmetric one for all the temperatures considered. For the resonant ∆V ) 0 case the transmission coefficient at low damping exceeds by a significant factor the surface hopping value of the golden rule plateau region. The above effect for nonadiabatic condensed phase reactions has been predicted by Onuchic and Wolynes15 on the basis of a qualitative semiclassical argument. At low damping dephasing processes are relatively slow and phase relationships are preserved on the time scale of several vibrational periods in the reactant well. These phase relationships lead to quantum interference effects that arise from the fact that quantum dynamics is described by probability amplitudes whose phases can constructively or destructively interfere. Onuchic and Wolynes generalized the surface hopping model to include these effects. They replaced the probability of electronic transition during each traversal of the crossing region, 2PLZ in eq 3.3, by the reaction probability during one coherent sequence of recrossings, Pcoh, computed from semiclassical probability amplitudes, and Nc by the average number of such coherent sequences. Their expressions give the dependence of the transmission coefficient on friction that is qualitatively the same as shown in Figure 2. The symmetric, resonant ∆V ) 0 reaction corresponds to the case of constructive interference while the antiresonant ∆V ) 0.5pω0 reaction to destructive interference. As can be seen from Figure 2a, intermediate values of exothermicity fall in between these two limiting situations. For the resonant case we have not been able to observe the transmission coefficient decrease below its resonant values when the friction is lowered further.52 The absence of that regime here may be due to the relatively small barrier height of our model. In order to illustrate these effects we examine the time dependence of the reactive flux. Figure 5 shows the timedependent transmission coefficient, defined as the (indefinite) time integral of the flux correlation function,
κ(t) ≡ (ZkTST)-1∫0 Cf(t′) dt′ t
(3.4)
In agreement with the above discussion, for strong friction κ(t) is essentially identical for the resonant and antiresonant reactions. It is an almost monotonic function of time and a welldefined plateau is reached within the first vibrational period in the reactant well. In the low-damping regime the transmission coefficient behaves differently for different values of ∆V. It exhibits a step structure with the step length equal to the vibrational period in the reactant well. For the resonant case we observe an upward step structure, each subsequent step being smaller than the previous one until the final plateau is reached. In the antiresonant case κ(t) makes steps in both directions but their total contribution is negative and leads to the decrease of
Figure 5. The time integral κ(t) of the flux correlation function normalized by the classical transition state theory value (see eq 3.4) for kBT ) 0.64pω0. Solid line, ∆V ) 0; dashed line, ∆V ) 0.2pω0; chain-dotted line, ∆V ) 0.5pω0. (a, top) η/m0ω0 ) 0.013, (b, middle) η/m0ω0 ) 0.026, (c, bottom) η/m0ω0 ) 0.75.
the rate compared to the “incoherent” short-time value. When the value of exothermicity is intermeidate between the resonant and antiresonant limits, the second and further steps in the correlation function lead to only a small correction to the rate coefficient. The step structure of the transmission coefficient at low damping is very similar to that we observed in the adiabatic case.5 There, at low temperatures when the probability of the reaction at each attempt to “cross” the transition state is small (as in the nonadiabatic case), we also observed an upward step structure of the integrated flux correlation function which increased the rate above the “incoherent” quantum transition state theory value. Therefore, it appears that these effects have similar origin for adiabatic and nonadiabatic reactions, arising from multiple coherent recrossings of the transition state region. IV. Conclusions In this paper we have presented accurate quantum mechanical calculations of rate constants for a model nonadiabatic reaction in dissipative environments. Our results exhibit the standard decrease of the rate from its adiabatic value by the LandauZener factor, nuclear tunneling at low temperature, and quantum interference effects. In particular, at high temperature we find a broad golden rule plateau where the rate coefficient is independent of friction and is determined by the Landau-Zener probability. At weak damping where the dephasing rate is slow, we observe in agreement with the semiclassical analysis of
4436 J. Phys. Chem., Vol. 100, No. 11, 1996 Onuchic and Wolynes15 a large enhancement of the reaction rate for the symmetric reaction which arises from constructive interference of quantum mechanical amplitudes. In the same regime, destructive interference leads to a decrease in the rate coefficient in the off-resonant case. These interference effects manifest themselves in a step structure of the integrated flux correlation function with step length equal to the vibrational period of the reactant potential well similar to those noted in our previous work on adiabatic reactions. While the calculations presented here employed a simple model of two harmonic potential surfaces, the effects of dissipation, nonadiabaticity, tunneling, and phase interference on the rate coefficient in nonadiabatic processes are quite general. The golden rule plateau and the crossover to the adiabatic regime at strong friction are relevant to proton and electron transfer processes in viscous solvents. The interference effects in the weak dissipation regime are very similar in nature to resonance structures observed in gas-phase experiments. Apart from providing quantitative measures for the effects discussed above, our results show that accurate quantum mechanical simulation of those condensed phase reactions that can be modeled in terms of effective harmonic bath modes is now feasible with standard computer resources. Using the quasi-adiabatic propagator path integral scheme we have been able to obtain converged results for the reaction rate over a wide range of temperature and for damping strength that spans four orders of magnitude. These calculations required up to several hours of CPU time on an IBM RISC 370 workstation for each set of parameters. At present, the basic limitation of our path integral approach is its restriction to reactions that can be described in terms of a system coupled to a harmonic bath. Within this class of models, the QUAPI methodology is simple as well as economical and thus readily applicable to the study of quantum rate processes in condensed phase or biological systems. Acknowledgment. This work was supported by the David and Lucile Packard Foundation through a Packard Fellowship for Science and Engineering. We are grateful to IBM for a one-time equipment award through the Shared University Research program. We also thank Peter Wolynes for interesting discussions. References and Notes (1) Hanggi, P.; Talkner, P.; Borcovec, M. ReV. Mod. Phys. 1990, 62, 251. (2) ActiVated barrier crossing; Fleming, G. R., Hanggi, P., Eds.; World Scientific: Singapore, 1993; p 82. (3) See the issue devoted to the fiftieth anniversary of Kramers’ reaction rate theory: Ber. Bunsenges. Phys. Chem. 1991, 95, No. 3. (4) Kramers, H. A. Physica (Utrecht) 1940, 7, 284. (5) Topaler, M.; Makri, N. J. Chem. Phys. 1994, 101, 7500. (6) Makri, N. Chem. Phys. Lett. 1992, 193, 435. Makri, N. In TimeDependent Quantum Molecular Dynamics; Broeckhove, J., Lathouwers, L., Eds.; Plenum Press: New York, 1992; p 209. (7) Topaler, M.; Makri, N. J. Chem. Phys. 1992, 97, 9001. (8) Topaler, M.; Makri, N. Chem. Phys. Lett. 1993, 210, 285. (9) Topaler, M.; Makri, N. Chem. Phys. Lett. 1993, 210, 448. (10) Makarov, D. E.; Makri, N. Phys. ReV. A 1993, 48, 3626. (11) Yamamoto, T. J. Chem. Phys. 1960, 33, 281.
Topaler and Makri (12) Miller, W. H.; Schwartz, S. D.; Tromp, J. W. J. Chem. Phys. 1983, 79, 4889. (13) Voth, G. A.; Chandler, D.; Miller, W. H. J. Chem. Phys. 1989, 91, 7749. (14) Voth, G. A. Ber. Bunsenges. Phys. Chem. 1991, 95, 393. (15) Onuchic, J. N.; Wolynes, P. G. J. Phys. Chem. 1988, 92, 6495. (16) Landau, L. D. Z. Sowjun. 1932, 2, 46. (17) Zener, C. Proc. R. Soc. A 1932, 137, 696. (18) Stueckelberg, E. C. G. HelV. Phys. Acta 1932, 5, 369. (19) Garg, A.; Onuchic, J. N.; Ambegaokar, V. J. Chem. Phys. 1985, 83, 4491. (20) Wolynes, P. G. J. Chem. Phys. 1987, 86, 1957. (21) Caldeira, A. O.; Leggett, A. J. Ann. Phys. (N.Y.) 1983, 149, 374. (22) Straub, J. E.; Borcovec, M.; Berne, B. J. J. Phys. Chem. 1987, 91, 4995. Berne, B. J.; Tuckerman, M. E.; Straub, J. E.; Bug, A. L. R. J. Chem. Phys. 1990, 93, 5084. Straub, J. E.; Berne, B. J.; Roux, B. J. Chem. Phys. 1990, 93, 6804. Tuckerman, M.; Berne, B. J. J. Chem. Phys. 1993, 98, 7301. (23) Zusman, L. D. Chem. Phys. 1980, 49, 295. (24) Fraunfelder, H.; Wolynes, P. G. Science 1985, 228, 337. (25) Tully, J. C.; Preston, R. K. J. Chem. Phys. 1971, 55, 562. (26) Webster, F. A.; Rossky, P. J.; Friesner, R. A. Comput. Phys. Commun. 1991, 63, 494. (27) Coker, D. F.; Xiao, L. J. Chem. Phys. 1995, 102, 496. (28) Wolynes, P. G. J. Chem. Phys. 1987, 87, 6559. (29) Bader, J. S.; Kuharski, R. A.; Chandler, D. J. Chem. Phys. 1990, 93, 230. (30) Marchi, M.; Chandler, D. J. Phys. Chem. 1991, 95, 889. 1992, 96, 1748. (31) Gehlen, J. N.; Chandler, D. J. Chem. Phys. 1992, 97, 4958. (32) Kapral, R.; Hudson, S.; Ross, J. J. Chem. Phys. 1970, 53, 4387. (33) Chandler, D. Introduction to modern statistical mechanics; Oxford University Press: Oxford, U.K., 1987. (34) Feynman, R. P. ReV. Mod. Phys. 1948, 20, 367. (35) Feynman, R. P.; Hibbs, A. R. Quantum Mechanics and Path Integrals; McGraw-Hill: New York, 1965. (36) Makarov, D. E.; Makri, N. Chem. Phys. Lett. 1994, 221, 482. (37) Makri, N.; Makarov, D. E. J. Chem. Phys. 1995, 102, 4600. Makri, N.; Makarov, D. E. J. Chem. Phys. 1995, 102, 4611. (38) Feynman, R. P. Statistical Mechanics; Addison-Wesley: Reading, MA, 1972. (39) Feynman, R. P.; Vrenon, Jr., F. L. Ann. Phys. 1963, 24, 118. (40) Harris, D. O.; Engerholm, G. G.; Gwinn, W. D. J. Chem. Phys. 1965, 43, 1515. (41) Dickinson, A. S.; Certain, P. R. J. Chem. Phys. 1968, 49, 4209. (42) Lill, J. V.; Parker, G. A.; Light, J. C. Chem. Phys. Lett. 1982, 89, 483. Light, J. C.; Hamilton, I. P.; Lill, J. V. J. Chem. Phys. 1985, 82, 1400. Lill, J. V.; Parker, G. A.; Light, J. C. J. Chem. Phys. 1986, 85, 900. (43) Choi, S. E.; Light, J. C. J. Chem. Phys. 1990, 92, 2129. Bacic, Z.; Light, J. C. J. Chem. Phys. 1986, 85, 4594. Bacic, Z.; Light, J. C. J. Chem. Phys. 1987, 86, 3065. Bacic, Z.; Light, J. C. Annu. ReV. Phys. Chem. 1989, 40, 469. (44) Colbert, D. T.; Miller, W. H. J. Chem. Phys. 1992, 96, 1982. Auerbach, S. M.; Miller, W. H. J. Chem. Phys. 1993, 98, 6917. (45) Hermann, M. R.; Fleck, Jr., J. A. Phys. ReV. A 1988, 38, 6000. Le Quere, F.; Leforestier, C. J. Chem. Phys. 1990, 92, 247. Bentley, J. A.; Wyatt, R. E.; Menou, M.; Leforestier, C. J. Chem. Phys. 1992, 97, 4255. (46) Das, S.; Tannor, D. J. J. Chem. Phys. 1990, 92, 3403. Williams, C. J.; Qian, J.; Tannor, D. J. J. Chem. Phys. 1991, 95, 1721. (47) Echave, J.; Clary, D. C. Chem. Phys. Lett. 1992, 190, 225. (48) Sim, E.; Makri, N. J. Chem. Phys. 1995, 102, 5616. (49) Zwanzig, R. J. Stat. Phys. 1973, 9, 215. (50) Cline, R. E.; Wolynes, P. G. J. Chem. Phys. 1987, 86, 3836. (51) Straub, J. E.; Berne, B. J. J. Chem. Phys. 1987, 87, 6111. (52) Because of the very large number of recrossings taking place at low friction, the correlation function needs to be computed for very long time. In this regime we employed the tensor multiplication scheme36,37 to calculate the quantum reactive flux for several hundred vibrational periods. The results obtained seem to indicate that at least for the parameters chosen here the transmission coefficient for the resonant case keeps increasing until the friction is lowered to values for which populations begin to exhibit coherent oscillations and a rate constant no longer exists.
JP951673K
