Software News and Updates

Software News and Updates DOIT: A Program to Calculate Thermal Rate Constants and Mode-Specific Tunneling Splittings Directly from Quantum-Chemical Ca...
Author: Fay Stanley
7 downloads 2 Views 243KB Size
Software News and Updates DOIT: A Program to Calculate Thermal Rate Constants and Mode-Specific Tunneling Splittings Directly from Quantum-Chemical Calculations

ABSTRACT: In this contribution we discuss computational aspects of a recently introduced method for the calculation of proton tunneling rate constants, and tunneling splittings, which has been applied to molecules and complexes, and should apply equally well to bulk materials. The method is based on instanton theory, adapted so as to permit a direct link to the output of quantum-chemical codes. It is implemented in the DOIT (dynamics of instanton tunneling) code, which calculates temperature-dependent tunneling rate constants and mode-specific tunneling splittings. As input, it uses the structure, energy, and vibrational force field of the stationary configurations along the reaction coordinate, computed by conventional quantum-chemical programs. The method avoids the difficult problem of calculating the exact least-action trajectory, known as the instanton path, and instead focusses on the corresponding instanton action, because it governs the dynamic properties. To approximate this action for a multidimensional system, the program starts from the one-dimensional instanton action along the reaction coordinate, which can be obtained without difficulty. It then applies correction terms for the coupling to the other vibrational degrees of freedom, which are treated as harmonic oscillators (transverse normal modes). The couplings are assumed linear in these modes. Depending on the frequency and the character of the transverse modes, they may either decrease or increase the action, i.e., help or hinder the transfer. A number of tests have shown that the program is at least as accurate as alternative programs based on transition-state theory with tunneling corrections, and is also much less demanding in computer time, thus allowing application to much larger systems. An outline of the instanton formalism is presented, some new developments are introduced, and special attention is paid to the connection with quantum-chemical codes. Possible sources of error are investigated. To show the program in action, calculations are presented of tunneling rates and splittings associated with triple proton transfer c 2001 John Wiley & Sons, Inc. J Comput Chem 22: in the chiral water trimer. 787–801, 2001 Keywords: thermal rate-constants; mode-specific tunneling; quantum-chemical calculations; DOIT

Correspondence to: Z. Smedarchina; e-mail: zorka. [email protected] Contract/grant sponsor: Xunta de Galicia; contract/grant number: XUGA20903A98

Journal of Computational Chemistry, Vol. 22, No. 7, 787–801 (2001) c 2001 John Wiley & Sons, Inc.

SOFTWARE NEWS AND UPDATES

Introduction

R

eactions in which a proton or a hydrogen atom is exchanged between two molecules or between two functional groups in the same molecule are among the most common in chemistry and biology. Because protons are light particles, their reactions may show strong quantum effects, such as tunneling through potential-energy barriers, and hence, cannot be treated by classical methods. This means that standard transition-state theory (TST), the method generally adopted to calculate reaction rates, is not directly applicable, but requires correction for tunneling.1 – 6 While in standard TST, the effect of the other degrees of freedom on the reaction coordinate is included simply through partition functions, in the quantum treatment appropriate for tunneling these effects enter through coupling terms, which need to be treated explicitly. As an alternative, classical trajectory calculations have been adapted to tunneling.7 – 12 The resulting multidimensional theories exist in various forms, depending on the chosen tunneling trajectory. The main drawback of these methods is that they become computationally very demanding when the tunneling trajectory deviates strongly from the minimumenergy path, because in that case the tunneling corrections are large. For this reason most applications are limited to small systems. The classicaltrajectories method, originally developed to calculate rate constants, has been extended to tunneling splittings8 but only the results for zero-point levels have proved to be reliable.9, 13 In a recent series of articles13 – 22 we have introduced an alternative method to deal with proton transfer processes, based on the instanton method. Although it has been known for some time that this is a very powerful method to describe tunneling,23, 24 quantitative applications to chemical reaction dynamics have been attempted only recently, possibly because at first sight it may appear that the method is ill suited to practical calculations due to its mathematical complexity. Benderskii and coworkers25 – 28 were the first to show, using exact analytical solutions of simple models, that the method offers a superior tool to deal with hydrogen tunneling in polyatomic molecules. They calculated “exact” instanton trajectories, i.e., least-action tunneling paths, for two-dimensional approximations to molecular potentials, and showed how these trajectories varied from nearly linear tunneling paths at low temperatures to the minimum-energy path

788

over the potential energy barrier at high temperatures. To make this method suitable to large molecules, complexes, and the condensed phase, it is necessary, however, to go beyond a two-dimensional approximation. Unfortunately, in practice, it is not possible to extend the calculations of least-action trajectories to more than two or three dimensions. From a chemical point of view, it is highly desirable to include all vibrational degrees of freedom, because the motion of the tunneling hydrogen is generally coupled to many other vibrations. This would allow the direct use of the output of quantum-chemical computations in the dynamics calculations. To achieve this generalization of the Benderskii approach we introduced a number of approximations, which allow us to avoid direct calculation of the instanton trajectory. Instead, we calculated the multidimensional instanton action in an approximate form, based on a generalization of the model used by Benderskii et al.27, 28 Assuming weak to moderate coupling between the tunneling mode and the transverse modes, we calculated the multidimensional instanton action by systematically correcting its one-dimensional (1D) variant for coupling to the transverse modes. The approximations used in our treatment made it imperative to test the method in detail. In a series of articles, we have used the method to calculate tunneling rate constants and tunneling splittings, including mode-specific tunneling splittings and isotope effects, in a number of molecules and complexes and compared the results with the available experimental data. We also carried out benchmark comparisons with tunneling methods based on transition-state theory with multidimensional tunneling corrections. The method is now available to other reseachers in the form of a computer program to be used in combination with a standard quantum-chemistry code. It is the purpose of the present article to describe this DOIT (dynamics of instanton tunneling) program.29 The program is written as a direct dynamics program, which uses as input the output of the quantum-chemistry code calculated at the stationary points of the chosen adiabatic surface to generate the full-dimensional potential-energy surface. The mode with imaginary frequency in the transition state is taken as the reaction coordinate, and the potential along this coordinate is derived from the positions, energies, and curvatures of the stationary points. The remaining 3N − 7 transverse modes of the N-atomic molecule or complex are taken to be harmonic, and their linear couplings to the

VOL. 22, NO. 7

SOFTWARE NEWS AND UPDATES tunneling coordinate are calculated from their displacements between the equilibrium configurations and the transition state. It follows that in its present form the program deals with transitions involving two potential-energy minima separated by a barrier, i.e., transitions inside a molecule or complex or in a condensed phase. In general, this includes all hydrogen-bonded systems, but not bimolecular reactions with short-lived collision complexes or reactions subject to strong variational effects. The latter limitation is usually not serious for hydrogenbonded systems, which tend to have “tight” transition states. In the following sections we describe the basic instanton method without going into mathematical details; for these we refer to the literature (e.g., see the reviews in refs. 25, 26, and 30, and the original papers cited therein). We discuss the approximations introduced to make the method applicable to multimode systems but, to avoid undue repetition, we frequently refer to our earlier articles in which these approximations are derived in more detail and illustrated by practical applications. In this contribution we focus on the techniques we use to compute the desired results from the input parameters. To demonstrate the operation of the program, we treat intermolecular hydrogen transfer in the chiral water trimer as an example, and attach the output of the program as supplementary material.

The Approximate Instanton Method The instanton approach is a quasiclassical method to calculate the dynamics of tunneling through a potential-energy barrier that is part of a multidimensional adiabatic potential-energy surface.24 – 26 It represents the overall result of tunneling along numerous paths of varying length and energy by that of tunneling along a single trajectory, called the instanton path, defined by the condition that at any temperature T it minimizes the Euclidian action. It is, thus, an extremal trajectory. The Euclidian action concept arises in the formulation of the classical action in the upside-down barrier, where the tunneling particle moves in imaginary time. Accordingly, it is evaluated with the Hamiltonian, instead of the Lagrangean. The tunneling rate constant is expressed via the “instanton action” SI along this trajectory by the relation  k(T) = A(T)e−SI (T) = P(T) i0 /2π e−SI (T) , (1) where A(T) is a parameter reflecting the contribution of a Gaussian distribution of trajectories close to

JOURNAL OF COMPUTATIONAL CHEMISTRY

the instanton path. This parameter is conveniently expressed in units i0 /2π, which represent the frequency with which the particle attacks the barrier, so that P(T) is a dimensionless parameter. Here and in the following the action is expressed in units h¯ . This rate constant has a non-Arrhenius temperature behavior, with a finite zero-temperature limit related to tunneling from the lowest energy level in the initial well. In general, direct evaluation of the instanton trajectory for a molecule with N atoms requires the solution (for every temperature) of 3N − 6 coupled equations of motion, run in imaginary time. This is not feasible, not even for relatively small systems; so far analytical instanton results in the presence of coupling that enhances tunneling have been reported only for two model double-well potentials in which the reaction coordinate was coupled to a single “effective” mode.27, 28 However, in typical proton transfer reactions, the heavy atoms tend to undergo substantial reorganization during the process, so that it will be necessary to include all the vibrational degrees of freedom of the system. This is possible only through the introduction of approximations. Some of these approximations are standard, such as the assumption that the transverse modes are harmonic and that their coupling to the tunneling mode is linear in the transverse mode coordinates. Coupling between transverse modes along the reaction coordinate will be included only in calculations of state-dependent properties such as mode-specific tunneling splittings. To treat this standard model, additional, less transparant approximations will be introduced that require careful testing. In our approximate instanton method14, 21 we concentrate on the evaluation of the instanton action without a direct search for the instanton trajectory; below, we describe the corresponding computational scheme, as implemented in the computer code DOIT 1.2.29 The instanton action in eq. (1) is found from the corresponding 1D action S0I (T), which is evaluated numerically for the potential along the reaction coordinate, through analytical corrections that represent its coupling to the 3N − 7 transverse modes. The dynamics is formulated in terms of the normal modes of the transition state, which is the configuration of highest symmetry, and the mode with imaginary frequency is chosen as the reaction coordinate. The linear coupling constants of the reaction coordinate to the transverse modes, treated as independent harmonic oscillators, are obtained directly from standard electronic structure and force-field calculations at the stationary points.

789

SOFTWARE NEWS AND UPDATES The double-minimum potential along the reaction coordinate is derived from the positions and curvatures of these stationary points; the curvature of this potential at the stable configurations is obtained by a unitary coordinate transformation. Thus, to run tunneling dynamics with the DOIT code, the minimum set of quantum-chemical input parameters needed are the energies, geometries, and force fields at the stationary points along the reaction coordinate. Because the actual dynamics calculations are computationally undemanding, the small number of configurations required make it possible to apply the method to reasonably large systems, including biopolymers and solids. The method allows a full-dimensional calculation of proton tunneling rate constants in a wide range of temperatures, and of tunneling splittings of zeropoint and vibrationally excited levels. It provides a clear physical picture of the contributions made to the tunneling by helping and hindering transverse vibrations. To test its accuracy, calculations have been carried out on several molecules for which reliable experimental data are available.21 These tests depend on the accuracy with which the quantumchemical methods used can reproduce the potentialenergy surfaces, and thus have no absolute validity unless the system is so small that the potentials can be assumed “exact.” However, such systems usually concern bimolecular reactions or couplings that are either very strong or involve anharmonic largeamplitude motions. The most successful application thus far is the reproduction of mode-specific tunneling splittings relative to zero-point splittings and deuterium isotope effects in tropolone and 9-hydroxyphenalenone.13, 17 This achievement, which has not yet been matched by any other available dynamics method, shows that the program deals effectively with couplings between the tunneling mode and transverse modes, at least for low vibrational quantum numbers. Results reproducing observed rates of inversion of small molecules have been calculated that are at least as accurate as those obtained with transition-state theory with tunneling corrections based on the same potential-energy surfaces. However, the effective dimensionality of these moleculess is low, so that they do not provide a stringent test for the program. Further testing will be necessary before the accuracy of the method can be reliably assessed. THE POTENTIAL-ENERGY SURFACE In our approach, which was recently reviewed,21 the reaction Hamiltonian is formulated in terms of

790

the (mass-weighted) normal coordinates of the transition state {x, y}, with frequencies {ω}, the mode x with imaginary frequency ω∗ being the reaction coordinate. In this formulation the 3N − 7 transverse vibrations y can be divided into two groups, depending on whether or not they undergo a displacement between the equilibrium configurations and the transition state. Only those that are displaced are linearly coupled to the reaction coordinate and will be considered further, because modes that are not displaced, for example, out-of-plane modes in planar molecules, do not contribute to the linear coupling. If the double well is symmetric, the point group of the transition state will have a corresponding symmetry element with respect to which the displaced modes are either symmetric (s) or antisymmetric (a). In the simple case of a point group without three- or higher fold axes, this symmetry element is a mirror plane. The program implements this classification by a procedure explained below. The multidimensional potential is formulated in terms of these two classes of displaced modes14 2 1X 2 U(x, y) = UA (x) + ωa ya ± xCa /ωa2 2 a X 2 1 + ωs2 ys − x2 Cs /ωs2 , (2) 2 s where UA (x) is the 1D double-well adiabatic potential along the reaction coordinate and +/− refers to the left/right well, the transition state being at zero. The parameters Ca,s are the coupling constants expressed in terms of the frequencies ωa,s and the displacements 1x, 1ya,s, derived, respectively, from the force fields and the geometries of the transition state and the stable configurations: Ca = ωa2 1ya /1x,

Cs = ωs2 1ys /(1x)2 .

(3)

Only antisymmetric modes can be displaced between the initial and final configuration, and the corresponding reorganization introduces a Franck– Condon factor into the expressions for the tunneling rate and tunneling splitting. To find the coupling constants in eq. (3), one needs to evaluate the 3N − 6 displacements between the transition state and the minima. Their relation to the Cartesian coordinates of the initial configuration i }, the transition state X = {Xα,j }, (reactant) Xi = {Xα,j and/or the final configuration (product), if differf ent, Xf = {Xα,j }, is given by 1xi,f =

3N 3 X X

 i,f M1/2 α,jj Xα,j − Xα,j Lα,j1 ,

(4)

α=1j=1

VOL. 22, NO. 7

SOFTWARE NEWS AND UPDATES i,f

1yk =

3N 3N−6 3 X X X

 i,f M1/2 α,jj Xα,j − Xα,j Lα,jk ,

(5)

α=1j=1 k=2

where α = x, y, z; M is the 3N × 3N diagonal matrix of the masses; and L the 3N × (3N − 6) matrix of the eigenvectors in the transition state, in which the first column corresponds to the mode with imaginary frequency. The total displacement |1xi | + |1x f | along the reaction coordinate is used to characterize the width of the potential-energy barrier, and the adiabatic barrier height U0 = UA (0) is used to characterize its height. The average distance travelled by the tunneling particle, measured along actual trajectories, will depend on the temperature and the (effective) mass; because it is not an observable, it is not calculated explicitly. To classify the displaced modes in a symmetric double well as either symmetric or antisymmetric, the program starts with the mode with imaginary frequency, which is obviously antisymmetric, and forms the vector r1 = (rα,1 ) whose components are the sums of the components of the elements of the eigenvector: rα,1 =

N X

Lα,i1 .

(6)

i=2

This defines the mirror plane perpendicular to r1 through the transition state. The program then forms similar vectors rj for the 3N − 7 remaining normal modes of the transition state and projects these vectors onto this plane through the normalized scalar product  (7) pj = (r1 , rj ) |r1 ||rj |. If the mode is symmetric, pj = 0; if it is antisymmetric, pj = 1. With this procedure there is no need to calculate the final state to classify the modes in the symmetric double well. In an asymmetric double-well potential the normal modes have partial symmetric and antisymmetric character. Clearly, the Franck–Condon factor is determined by the total reorganization, so that the antisymmetric component must be of the form f f (8) 21ya = 1yk − 1yik = 1yk + 1yik . From the analogy with the symmetric case it follows that the symmetric component must be the alternative combination, which implies  f (9) 1ys,a = 12 1yk ± 1yik . Similarly, we obtain for the reaction coordinate  (10) 1x = 12 1x f − 1xi .

JOURNAL OF COMPUTATIONAL CHEMISTRY

The potential of eq. (2) near the barrier top, i.e., for x ' 0, assumes the form UA (x) ' U0 − 12 (ω∗ x)2 ,

(11)

where U0 is the adiabatic barrier height. Close to the minima, where x ' xi,f on the reactant and product side, respectively, the potential can be approximated by 2 2 (12) UA (x) ' 12 i0 x − 1xi ,   2 f 2 (13) UA (x) ' 12 0 x − 1x f , i,f

0 being the corresponding effective frequency associated with the reaction coordinate. The potentialenergy surface of eq. (2) can be rewritten in terms of the normal modes of the equilibrium configurations, zi,f . The coordinates zi are related to the coordinates {x, y} by the (3N − 6) × (3N − 6) unitary matrix Gi = L† Li ,

(14)

where L† is the transposed matrix of the transitionstate eigenvectors and Li the matrix of the reactant eigenvectors. The columns of Gi represent the coefficients that relate the modes of the initial configuration to a given transition-state mode. If the first column is chosen as the mode with imaginary frequency, the effective frequency of the transfer mode in the initial state is represented by the combination "3N−6 #1/2 X  i i i 2 . (15) G1k ωk 0 = k

A similar relation can be written for the final state. These equations are approximate, because they ignore the fact that in the transformed coordinates the kinetic energy is not diagonal. For the calculation of temperature-dependent rate constants and zero-point tunneling splittings, we include off-diagonal elements only of the form Gi1k , i.e., we neglect mixing between the transverse modes. However, for mode-specific tunneling splittings (and mode-specific rate constants, when implemented), we include all off-diagonal matrix elements of the modes populated in the state whose splitting (or rate constant) is being considered. All displacements and transformations are carried out in such a way that the Eckart conditions31 are obeyed. This is necessary to conserve the linear and angular momenta, and thus prevent mixing of the vibrations with translations or rotations. These conditions are implemented in the program by the procedure described in ref. (32).

791

SOFTWARE NEWS AND UPDATES valid in the limit T = 0 and is directly related to the zero-point tunneling splitting.

EVALUATION OF THE MULTIDIMENSIONAL INSTANTON ACTION To evaluate the multidimensional instanton action, we apply an approximation scheme that amounts to a multidimensional generalization of an exact solution for a two-dimensional model system obtained by Benderskii et al.27 These authors considered tunneling in a 1D symmetric potential formed by two crossing parabolas coupled to a symmetric harmonic mode. If the mode is slower than the tunneling motion, the instanton action can be transformed so as to separate the terms related to the 1D tunneling potential from those related to the coupling. In the opposite case the mode is treated in the adiabatic limit, which results in massrenormalization of the tunneling particle.27 We generalize this result to the multidimensional case by assuming that both coupling effects are additive, because the transverse modes are independent. This scheme yields an expression for the instanton action, in which the numerator is represented by the instanton action for the 1D potential along the reaction coordinate, evaluated for a particle with renormalized mass, representing coupling to highfrequency modes; the coupling to low-frequency symmetric modes enters the denominator via correction terms. Because we have shown that these terms are independent of the shape of the 1D potential, at least in the zero-temperature limit, we reasonably assume that this approach need not be limited to the potential formed by two crossing parabolas, but should apply equally well to any well-behaved 1D potential coupled to a system of independent harmonic modes. Low-frequency antisymmetric vibrations are included through Franck– Condon factors of the type F(T) = e−δa (T)

(16)

and thus simply add a term to the instanton action. This leads to a final approximate expression for the multidimensional instanton action of the form X S0I (T) P + αs δa (T), (17) SI (T) = 1 + s δs (T) a where the parameters δs,a represent the contributions of the transverse vibrations that are slower than the tunneling motion. The contributions of the faster vibrations, which can be treated adiabatically, are included in the 1D action S0I (T) in the form of an effective mass meff . All these contributions are functions of the couplings Ca and Cs . The factor αs represents the modulation of the antisymmetric coupling by the symmetric coupling.27 Equation (16) remains

792

THE ONE-DIMENSIONAL POTENTIAL-ENERGY BARRIER By definition, the 1D instanton path minimizes the Euclidian action evaluated in imaginary time in the upside-down 1D barrier. The details of its evaluation in the present approach are given in the Thermal Rate Constants section. Here we proceed with the definition of the one-dimensional adiabatic potential UA (x) for which the 1D instanton action is evaluated. This potential should obey the conditions ∂U(x, ya,s)/∂ya,s = 0 for all modes ya,s, so that the gradient is directed along x. It is evident that points along the minimum-energy path do not obey these conditions, because our reaction coordinate is the mode with imaginary frequency in the transition state and, except for this point, does not coincide with that path; however, it can be obtained by projecting the minimum-energy path onto the mode with imaginary frequency. Because the curvatures at the stationary points have been evaluated, the direction of this gradient is known near these points. Therefore, a simple way to approximate the potential UA (x) for intermediate points is to connect the regions of the stationary points by their common tangents, i.e., the tangent to the parabolas corresponding to frequency ω∗ at the transition state i,f and to 0 at the equilibrium configurations. This procedure proved sufficient in all cases studied so far. If a more accurate potential is needed, it is necessary to find intermediate configurations that obey the above conditions and to determine their energy; this problem will be discussed elsewhere. In simple cases it may be sufficient to use an analytical interpolation scheme; for instance, to calculate tunneling splittings, one may only need an accurate representation of the barrier shape deep inside the well. In that case a two-parameter quartic expression of the form  2 (18) UA (x) = U0 1 − (x/1x)2 , may be satisfactory. A possible generalization of this potential for an asymmetric well is given by   2 2 , x ≤ 0, (19) UA (x) = U0 1 − x/1xi    2 2 , UA (x) = 1E + (U0 − 1E) 1 − x/1x f x ≥ 0,

(20)

where 1E is the difference in energy between the two stable configurations (i.e., the exothermicity of the reaction). For thermal rate-constant calculations,

VOL. 22, NO. 7

SOFTWARE NEWS AND UPDATES the shape of the 1D potential becomes more important, and the potential with common tangents should be generated, using the curvatures of the stationary points in addition to the barrier width, barrier heigh, and the exothermicity. This potential is generated from the parabolas represented by eqs. (11), (12), and (13), and their first derivatives. i,f The linkage points x± are given by the equations i,f

x± = 1xi,f where Ai,f =

Bi,f

Ai,f ± Bi,f i,f

1 + (0 /ω∗ )2

,

 i,f  1xi,f 0 2 , |1xi,f | ω∗

   2Ei,f 1 2 1/2 1 = 1− + ∗2 (ω ) (1xi,f )2 (i,f )2

(21)

(22)

0

with E = U0 and E = U0 + 1E. In addition to the quartic and common-tangent potentials, the program allows the input of any 1D potential in analytical or numerical form. For some applications (see Thermal Rate Constants section), it turns out to be useful to include the zero-point energy difference in the definition of the barrier. This vibrationally adiabatic barrier height is given by i

UVA,0

f

M−1 M 1X 1X i 1 h¯ ωm − h¯ ωm + h¯ i0 = U0 + 2 m 2 m 2

(23)

where m runs over all M modes linearly coupled to the reaction coordinate. THE TRANSVERSE MODES Depending on their frequency, the transverse modes are divided into two groups. Modes that are faster than the tunneling motion contribute to the mass renormalization27 X X 2 2 meff = 1 + Ca /ωa2 + 4x2 Cs /ωs2 , (24) a

s

and modes that are slower than the tunneling motion contribute to the correction terms δa,s in eq. (17). The symmetric correction term assumes the form

terms of the coupling constant:16, 17  δs (0) = 14 i0 /ωs γs2 , where

(26)

 γs = ωs2 1ys 2 1x

(27)

is a dimensionless coupling parameter. In this expression  is the scaling frequency, defined as  = q

UC0 /1x, where UC0 is the crude-adiabatic barrier height, obtained from eq. (2) by setting x = 0, y = 1y. If it is lower than the imaginary frequency ω∗ , as is usually the case,  rather than ω∗ is used to divide the transfer modes into fast modes contributing via eq. (24) and slow modes contributing via eq. (25). In other words, the “cutting” frequency is given by c = min[ω∗ , ]. For strongly displaced symmetric modes we introduce a new coupling scheme based on the expansion in normal-mode displacements of the instanton action for the model of two crossing parabolas 27 referred to above. In this scheme eq. (27) is replaced by    g κ2 2 4 − 15 + 10κ − κ δs (T) = −3g 1 − 3 6   h¯ ωs 2g 1 − 8g 1 + 5κ 2 coth , (28) + κ 2ηkB T where g=

γs2 ; 8(1 − κ 2 )2

κ=

ωs . i0

(29)

The antisymmetric correction term assumes the form . h¯ ωa , ωa ≤ c , (30) δa (T) = η(1ya )2 a20a coth 2ηkBT a0a = (h¯ /ωa )1/2 being the zero-point amplitude of mode ya . At high temperatures this expression turns into the usual activation term δa (T) = Ua /kB T, where Ua = ωa2 (1ya )2 /2 represents the activation energy for the final reorganization of this mode. THE ZERO-TEMPERATURE LIMIT

(25)

Within the instanton formalism the rate constant of eq. (1) has a non-zero limit at T = 0:  (31) k(0) = P(0) i0 /2π e−SI (0) ,

where η = 1 + 1xe /1x f , 1xe being the value of x at energy E = 0 at the product side of the barrier, i.e., the point where the proton exits the barrier. The zero-point correction term δs (0) is expressed in

where the instanton action of eq. (17) at zero temperature equals X 2SC (E0 ) P + αs δa (0), (32) SI (0) = 1 + s δs (0) a

δs (T) = δs (0) coth

h¯ ωs , 2ηkB T

JOURNAL OF COMPUTATIONAL CHEMISTRY

793

SOFTWARE NEWS AND UPDATES 2SC (E0 ) being the classical action evaluated for the lowest energy level E0 in the well Z x2    1/2 dx, (33) 2meff UA (x) − E0 SC (E0 ) = 2 x1

where x1 and x2 are the classical turning points for the zero-point energy level. TUNNELING SPLITTINGS The zero-level tunneling splitting is calculated as  1/2 i  −S (0)/2 . (34) 0 /π e I 1(0) = P(0) Our approach also allows evaluation of modespecific tunneling splittings. To relate the excitation of a given mode ωki of the initial equilibrium configuration to the level vik , to the quantum numbers vs (symmetric) and va (antisymmetric) of the corresponding transition-state modes, which are the reference modes in our approach, we use the G-matrix vs =

3N−6 X

2

Gisk vik ,

va =

k=1

3N−6 X

2

Giak vik .

(35)

k=1

For symmetric modes, the effect of the coupling is proportional to the squared vibrational amplitude; for low vibrational quantum numbers where anharmonic effects may be neglected, the excitation of such a mode increases linearly with the vibrational quantum number, and therefore, the correction term increases accordingly: δs (vs) = (2vs + 1)δs (0).

(36)

The excitation of an antisymmetric mode depends on the magnitude of the displacement 1ya relative to the zero-level amplitude a0a of the mode: if 1ya ≥ a0a , the Franck–Condon effect for an excited level va is weaker that that of the ground-state level as follows from the expression −αs δa (va )

Fa (va ) = e where

,

 δa (va ) = 2δa (0) (2va + 1).

(37) (38)

If 1ya ≤ a0a , the excitation of the mode leads to a stronger suppressing effect given by the expression   (39) Fa (va ) = Lv αs δa (0) e−αs δa (0) , where Lv is the Laguerre polynomial of order v = vik . In addition to the effects described above, the excitation of a given mode in the equilibrium configuration generally results in (partial) excitation of the reaction coordinate, in which case the zero-level energy E0 should be replaced by the effective energy 2 (40) Eeff = 12 + vik Gi1,k h¯ i0 .

794

THERMAL RATE CONSTANTS To achieve maximum clarity we start the explanation of our computational scheme with an elementary application of the instanton method to tunneling in a 1D potential UA (x); for more detailed treatments, we refer to several recent reviews21, 25, 26, 30 and the original articles cited therein. The transfer rate constant is a Boltzmann average over the partial rate constants for every bound level in the well: ∞ X −1 0 k(Ei )e−βEi (41) k (T) = Z0 i=1

where Z0 is the one-dimensional partition function, given approximately by  i (42) Z−1 ¯ 0 /2 . 0 ' 2 sinh β h The rate constant consists of a tunneling component, k0tun , and an over-barrier (“classical”) component, k0cl , which comprise the effect of levels below and above the barrier, respectively. To illustrate the temperature dependence of the rate constant, we replace the sum in eq. (41) by the integral Z ∞ ρ(E)W(E)e−βE dE, (43) k0 (T) = Z−1 0 0

where the quasiclassical tunneling probability is evaluated as 1 , (44) ρ(E)W(E) = (2π h¯ )[1 + e2SC (E) ] SC (E) being the classical action of eq. (33) evaluated at energy E. At sufficiently high temperatures k0 (T) shows Arrhenius-type behavior and can be treated satisfactorily by classical 1D transition-state theory −βU0 . k0cl (T) ' k0TST (T) = (1/2π h¯ β)Z−1 0 e

(45)

At low temperatures tunneling prevails and the temperature dependence of the rate constant deviates from the Arrhenius equation. The integral in eq. (43) represents competition between long paths with small tunneling probabilities but favorable Boltzmann factors and short paths with large tunneling probabilities but less favorable Boltzmann factors. Given the exponential dependence of the rate on these parameters, this means that the tunneling is dominated by a narrow band of levels around a single energy level Ea (T), characterized by a period of motion in the upside-down barrier 2 = h¯ β.23 The method of steepest descent allows direct integration of eq. (43) which yields [cf. eq. (33)] 1/2 −2S (Ea )−βEa  e C . (46) k0tun (T) = Z−1 ¯ /|∂2/∂Ea | 0 2π h

VOL. 22, NO. 7

SOFTWARE NEWS AND UPDATES This expression can be written in an Arrhenius-like form k(T) = C(T)e−βEa , Ea (T) being the apparent activation energy, which decreases with temperature. When T → 0, both Ea and C(T) tend to zero yielding a finite rate constant k0 (0) of the order28 k0 (0) ' (i0 /2π)e−2SC (E0 ) ,

(47)

E0 being the lowest energy level in the well. As the temperature increases, the vibrational period 2 decreases but cannot become smaller than 2π/ω∗ . Thus, T∗ = h¯ ω∗ /2πkB

(48)

the so-called crossover temperature, defines the threshold above which the periodic orbit ceases to exist. Equation (46) is equivalent to the standard instanton formula which relates the Gamov-type rate of decay of a metastable state to the imaginary part of the free energy k0tun (T)

=

Z−1 0

∞ X

0(Ei )e−βEi

i=1

= 2(Im Z)/βZ0 = (2/h¯ ) Im F.

(49)

where 0(Ei ) is the “width” of the level Ei . Through the path-integral formulation of the partition function Z, this expression relates the rate constant to the instanton trajectory, which is a saddle point of the Euclidian action. Thus, the tunneling rate constant takes the form  0 (50) k0tun (T) = P(T) i0 /2π eSI (T) , where S0I (T) is the 1D Euclidian action corresponding to the instanton path, and the prefactor P(T) represents the contribution of a narrow Gaussian channel of trajectories around it. The equivalence of eq. (50) to the periodic-orbit solution of eq. (46) allows evaluation of the 1D instanton action S0I (T) in eq. (17). However, this approach can be applied only at temperatures below T∗ , where such a periodic orbit exists. As we demonstrate below, above T∗ tunneling still prevails and a different method has to be found for the evaluation of the tunneling rate constant. For this purpose we return to eq. (41), and instead of approximating it by eq. (43), we represent the 1D rate constant by Z ∞ ρ(E)W(E)e−βE dE, (51) k0 (T) ' k0 (0) + Z−1 0 E0

which yields correct temperature behavior in the limits of low and high temperature and thus offers a unified expression for the whole temperature range.

JOURNAL OF COMPUTATIONAL CHEMISTRY

The validity of this representation has been tested for several model potentials and a range of parameter values. Below T∗ it generally yields good agreement when compared with the results of direct diagonalization and the periodic-orbit calculations of eq. (46), and at high temperatures it is found to link smoothly to the classical rate expression of eq. (45). This is illustrated in Figure 1 for the rate constant of thermally activated escape from the well of a cubic parabola UA (x) = 12 20 x2 − 13 kx3 , which is an example of an asymmetric potential.33 Similar behavior was found for symmetric potentials, for example, the quartic potential of eq. (18). Figure 1 indicates that the term “crossover temperature” used for the temperature T∗ , above which the instanton path ceases to exist, is a misnomer, because tunneling tends to remain the dominant transfer mechanism up to temperatures of about 2T∗ . This is of great significance for proton transfer across hydrogen bonds, where T∗ is often found to be close to room temperature, while isotope effects and non-Arrhenius temperature dependence clearly show the importance of tunneling at these temperatures. Unlike the periodic-orbit solution of eq. (46), the tunneling part of the rate constant in our representation can be evaluated well above T∗ , up to tem-

FIGURE 1. Rate constants of thermal escape from the metastable state of a cubic parabola, U0 /¯h0 = 3,

T ∗ = 343 K. Open circles represent the rate constant of eq. (41) for i = 0, 1, 2 obtained via direct diagonalization.33 The dot-dash line represents the rate constant obtained by classical transition-state theory. The dashed line corresponds to the periodic-orbit solution of eq. (46) in the temperature region where it exists. The thick solid line represents eq. (51) adopted in the DOIT code. The two vertical lines indicate the region between T ∗ and 2T ∗ .

795

SOFTWARE NEWS AND UPDATES peratures where overbarrier transitions take over. This allows us to define a 1D “instanton-like” action above T∗ by rewriting the tunneling component of eq. (51) in the form Z U0 ρ(E)W(E)e−βE dE k0tun (T) ' k0 (0) + Z−1 0 E0

 0 = P(T) i0 /2π e−SI (T) .

(52)

At high temperatures, typically above 2T∗ , overbarrier transfer will become dominant and the tunneling contribution becomes negligible. Thus, we can turn this equation into an expression valid for arbitrary temperature by adding a term that represents overbarrier 1D transfer defined by eq. (45):  0 (53) k0 (T) = P(T) i0 /2π e−SI (T) + k0TST (T). Although for simplicity this derivation was carried out for a 1D potential, this partitioning of the rate constant can be readily extended to multidimensional potentials, because the basic arguments do not depend on the dimensionality. The extrema of the Euclidian action include two saddle points, one corresponding to classical overbarrier transfer governed by the Boltzmann factor, and one corresponding to through-barrier tunneling, governed by the instanton action. This is an exact result, independent of the dimensionality. We, therefore, can replace S0I (T) by the multidimensional instanton action of eq. (17), and k0TST (T) by its multidimensional equivalent kTST (T). Next we consider the evaluation of the preexponential factor A(T) = P(T)(i0/2π). It has been shown28 that in the adiabatic limit for the transverse modes P(T) tends to unity, provided that the 1D instanton action S0I (T) in eq. (17) is evaluated with the vibrationally adiabatic potential UVA (x), obtained from UA (x) by the application of zero-point energy corrections of eq. (23). Then the rate constant assumes the form  (54) k(T) = i0 /2π e−SI (T) + kTST (T) with a zero-temperature limit defined by eq. (31). Note that the approximation scheme leading to eq. (54) uses the adiabatic approximation for the transverse modes only to evaluate the preexponential factor; for the evaluation of the most important parameter, the instanton action, the effect of coupling of the transverse modes with the tunneling mode is fully taken into account. The instanton approach reduces the (3N − 6)-dimensional problem effectively to a one-dimensional problem, because the approximation scheme that transforms one-dimensional into multidimensional

796

instanton action is carried out analytically. This greatly simplifies the dynamical calculations and reduces the laborious part of the computations to that of finding the proper potential-energy surface. This means that the method can be directly applied to any system for which the structure and force field can be evaluated at the stationary points of the reaction coordinate. SOLVENT EFFECTS Many reactions of biological interest take place in solution or in a polar environment, such as a cell or a protein molecule. The code DOIT 1.2 includes the option of incorporating dynamic solvent effects. Although long-range solvent effects can be treated by a continuum model, short-range effects generally require the explicit evaluation of the structure of the complex of the molecule associated with a number of discrete solvent molecules. This leads to a model of a specific supermolecule (molecule + solvent shell) in a cavity surrounded by a dielectric continuum, treated by a self-consistent reaction field model. The simplest of the continuum models is the Onsager model.34 The supermolecule is treated as an isolated species, but it is convenient to combine its lowest frequency modes, namely those with frequencies h¯ ωc ≤ 4kB T, with the slow motions responsible for the reorganization of the continuum and treat this combination classically. Their thermal excitation can be described by a Marcus-type activation energy Ur0 = (λ + 1G0 )2 /4λ,

(55)

where 1G0 is the free-energy of the reaction, λ = λd + λb is the sum of the reorganization energy of the supermolecule X ωc2 1y2c , (56) λd = 2 c

and that of the bulk. In the special case of a spherical cavity, the bulk term assumes the form  2  (µ Ei − µ E f )2 ε − 1 n2 − 1 e − 2 . (57) λb = 3 R 2ε + 1 2n + 1 4πε0 E f are the dipole moment In this expression µ E i and µ vectors of the initial and final states, respectively, of the supermolecule, R is the radius of the effective spherical cavity containing it (calculated, for instance, with the Onsager model), ε and n the dielectric constant and the refractive index of the solvent, e the electron charge, and ε0 the permittivity of the vacuum; the high-frequency modes of the supermolecule remain unactivated. The Franck– Condon factor in the presence of a strong solvent

VOL. 22, NO. 7

SOFTWARE NEWS AND UPDATES effect thus assumes the form

P −αs [ a-c δa-c (0)+Ur0 /kB T]

F(T) = e

.

(58)

In the program this scheme is used if the reorganization energy of the bulk is comparable to that of the supermolecule, i.e., if the relation λb ≥ λd /4 holds; in the opposite case the contribution due to the bulk is neglected and the program uses the same scheme as for the isolated complex.

Example: Proton Transfer and Tunneling Splitting in the Chiral Water Trimer To demonstrate the properties of the method and the code, we consider the hydrogen-bonded ring of three water molecules illustrated in Figure 2. This water trimer is of general interest as an example of a hydrogen-bonded network and has received much attention in the recent literature. 35 – 38 Our own interest in these systems is in their relation to water wires proposed as proton conductors in biological systems. For this reason we focus on the coherent proton exchange of the three hydrogenbonding protons in the ring. Although a spectroscopic study of the chiral water trimer by Pugliano and Saykally36 had provided experimental evidence for this interconversion of the two enantiomers, Monte Carlo simulations led Gregory and Gary 37 to predict that such a process would not take place. Recently, Loerting and coworkers35 addressed this problem by means of variational transition-state theory (variational effects turned out to the negligible) with small-curvature and large-curvature tunneling corrections and reported values for the zero-point tunneling splittings and the corresponding coherent transfer rate constant. They concluded that large-curvature tunneling was the better approach.

FIGURE 2. The two enantiomers of the chiral water trimer (left and right) and the transition state (center), calculated at the B3LYP/cc-pVTZ level. All distances are in Å, angles in degrees. The absolute energies in atomic unints are −229.410435 for the equilibrium configurations (symmetry C1 ) and −229.380092 for the transition state (Cs ).

JOURNAL OF COMPUTATIONAL CHEMISTRY

However, there are good reasons to repeat the calculations with a different dynamics method, because of an apparent discrepancy between the values reported for the low-temperature tunneling rates k(T → 0) and the zero-point tunneling splittings 1(0). The two quantities are approximately related through the simple expression  2 (59) k(0) ' (π/2ω0 ) 1(0) , where ω0 is the frequency along the reaction coordinate in the initial state. Within the formalism used in ref. 34, it is the lowest frequency along the minimum-energy path. Substitution of this frequency ω0 = 231 cm−1 and the zero-point splitting of 1.9 × 10−5 cm−1 calculated by Loerting et al.35 into eq. (59) yields a rate constant in the zerotemperature limit k0 = 0.5 s−1 . This value is inconsistent with the value k(100) ' 10−8 s−1 calculated by these authors in the large-curvature approximation. Because the program they used does not converge to the correct zero-temperature limit, we believe that the latter value is incorrect. We, therefore, repeat the calculations with the instanton method, as implemented in the DOIT 1.2 code, at the B3LYP/cc-pVTZ level.39, 40 The quantum-chemical calculations were performed with the GAUSSIAN 98 program.41 The resulting structural and energetic parameter values, displayed in Figure 2 and Table I agree well with those reported by Loerting et al.,34 except for the slightly lower value of the barrier height (25.3 against 27.0 kcal/mol) and the corresponding reduction in the imaginary frequency (1807i against 1833i cm−1 ). The reaction coordinate corresponding to the mode with imaginary frequency is illustrated in Figure 3a. The instanton calculations based on our parameter values yield a zero-point splitting of 3.0 × 10−4 cm−1 , which is to be compared with the value 1.9 × 10−5 cm−1 reported by Loerting et al. If we adjust our barrier to their higher value, our splitting is reduced to 1.2 × 10−4 cm−1 , leaving an unexplained discrepancy of a factor of 6. This is not an exceptionally large value given the sensitivity of the splitting to the structural and energetic parameters and the approximations inherent in both calculations. In our treatment the effective frequency of the mode representing the reaction coordinate in the equilibrium configuration, 0 , equals 3154 cm−1 , as follows from eq. (15). As expected, it coincides almost completely with the antisymmetric OH-stretching mode, the corresponding diagonal element of the G matrix being 0.953, while the offdiagonal elements are about 0.2 or smaller. The diag-

797

SOFTWARE NEWS AND UPDATES TABLE I. Vibrational Frequencies, Assignments, and Displacements of the Normal Modes of the Chiral Water/Heavy Water Trimer in the Transition State Obtained at the B3LYP/cc-pVTZ level; ωi Is the Frequency and Gia,s;k the G-Matrix Element of the Major Component in the Stable Configuration.

ω, cm−1

Assignment

1y, Å amu1/2

ωi , cm−1

Gia,s;k

a0 3852/2805 3847/2799 2010/1425 1690/1229 1624/1184 1411/1007 1253/906 731/704 671/536 553/504 502/362 441/319

OH(D) stretch CC + CO stretch OH(D) stretch OH(D)O in-plane def out-of-plane OH(D)O bend OH(D) stretch OH(D) stretch OH(D)O breathing wag OH(D)O def twist + rock rock

0.020/0.041 0.010/0.011 0.157/0.211 0.067/0.088 0.008/0.010 0.074/0.093 0.050/0.061 1.560/1.714 0.599/0.074 0.020/0.223 0.104/0.175 0.012/0.076

3869/2822 3864/2817 966/696 1652/1205 1659/1211 642/463 3675/2592 231/221 392/183 196/288 263/183 219/156

0.852/0.816 0.968/0.957 0.764/0.770 0.598/0.567 0.603/0.554 0.820/0.817 0.842/0.819 0.799/0.941 0.674/0.523 0.708/0.738 0.806/0.699 0.661/0.757

3850/2801 1581/1143 1539/1115 1259/909 780/572 688/529 548/511 447/327 1807i/1291i

OH(D) stretch out-of-plane OH(D)O bend OH(D)O def OH(D) stretch wag wag OH(D)O def + H-bond stretch rock OH(D) stretch

0.003/0.002 0.007/0.022 0.009/0.024 0.003/0.009 0.088/0.155 0.099/0.066 0.040/0.121 0.127/0.199 0.811/1.136

3866/2819 1659/1211 711/516 3682/2596 490/366 388/198 198/280 198/151 3495/2536

0.843/0.807 0.526/0.619 0.739/0.673 0.882/0.846 0.816/0.806 0.725/0.568 0.618/0.767 0.643/0.761 0.920/0.898

a00

onal elements of G are collected in Table I, together with the normal-mode displacements between the initial state and the transition state for both the light and heavy water trimer. The largest displacement is that of the symmetric breathing mode of 731 cm−1 , illustrated in Figure 3b, which helps the tunneling by shortening the O—O distance. The corresponding correction factor δs (0) = 0.300, derived from eq. (28), reduces the instanton action by 30%, which represents a very large effect. In heavy water this factor assumes a value of 0.465, which is close to the limit where the method ceases to be reliable. The coupling effect of antisymmetric modes, such as the deformation mode illustrated in Figure 3c is very small for the trimer, so that the overall effect of all couplings amounts to a strong enhancement of the tunneling. As follows from Table I, the major component of the strongly coupled symmetric breathing mode has a frequency of 231 cm−1 in the equilibrium configuration, which is the reaction-coordinate mode of the minimum-energy path. Because G2sk = 0.64 for this mode, the corresponding correction factor

798

for its first excited level vik = 1 will be δs (1) = 0.684, leading to a mode-specific tunneling splitting of 1.1 × 10−2 cm−1 , a factor of 40 larger than that calculated for the zero-point level, which should make it much easier to detect experimentally. Excitation of the mode of 490 cm−1 , which represents the major component of the antisymmetric mode illustrated in Figure 3c, has little effect on the splitting. A summary of the theoretical predictions is collected in Table II. The calculated tunneling rates as a function of temperature are displayed in Figure 4. Specifically, the graph shows that k(0) ' 8.6 s−1 . For the tunneling-mode frequency 0 = 3154 cm−1 , this result agrees with the value derived from eq. (59) and disagrees strongly with the rate constant reported in ref. 35. Hence, although we agree qualitatively with the zero-point splittings obtained by these authors, we conclude that their calculated rate constants must be in error by many orders of magnitude. According to our calculations, which yield results for the tunneling rate constants and tunneling splittings that are mutually compatible, the

VOL. 22, NO. 7

SOFTWARE NEWS AND UPDATES

FIGURE 4. Arrhenius-type plot of the triple proton tunneling rate constant governing tautomerization of the chiral water trimer. The thick solid line depicts the calculated rate constant for a barrier height of 25.3 kcal/mol; the thin solid line depicts the rate constant without coupling to transverse modes, and the dot-dash line the rate constant obtained by standard transition-state theory.

FIGURE 3. The water trimer transition state with three normal modes involved in the transfer between the two enantiomers: (a) The mode with imaginary frequency ω∗ = 1807 cm−1 ; (b) a totally symmetric breathing mode with frequency ωs = 731 cm−1 , which helps tunneling; and (c) an antisymmetric mode with frequency ωa = 780 cm−1 , which hinders tunneling.

cause the actual dynamics calculations take only a small fraction of the computer time necessary for the quantum-chemical calculations and the quantumchemical calculations required represent a minimum in the sense that no dynamics method can produce accurate results without at least this information. The accuracy of the results will of course depend on the quality of the potential, but the aim is for the dynamics method not to introduce errors larger than those resulting from inaccuracies in the potential. Whether this has been achieved remains open to question, although the tests carried out thus far indicate that inaccuracies in the calcu-

lifetime of an enantiomer of the chiral water trimer should be of the order of 30 µs at room temperature.

Discussion The basic aim of the method presented here is to extent our ability to calculate and interpret dynamic properties of protons and hydrogen atoms to all systems for which the appropriate potentials can be calculated. This has certainly been achieved, beTABLE II.

Zero-Level and Mode-Specific Tunneling Splittings 1, in cm−1 , and Thermal Rate Constants kH,D (T), in s−1 , for Hydrogen (H) and Deuterium (D) Exchange in the Chiral Water and Heavy Water Trimers. Band

Freq.

1H a

1D a

1H b

0(0,0) 5(0,1) 9(0,1)

— 231 490

3.0 · 10−4 [1.2 · 10−4 ] 1.1 · 10−2 [5.1 · 10−3 ] 1.9 · 10−4 [1.1 · 10−4 ]

6.6 · 10−7 [1.9 · 10−7 ] — —

1.9 · 10−5 — —

kD a (100) 2.5 · 10−2 [7.1 · 10−6 ]

kH a (300) 3.1 · 104 [5.0 · 103 ]

kH b (100) '10−8

kH b (300) 3.4 · 10−4

kH a (100) 10.2 [2.9]

Numbers in square brackets refer to values with the barrier height adjusted to 27.0 kcal/mol, the value of ref. 35. a This study. b Large-curvature tunneling, ref. 35.

JOURNAL OF COMPUTATIONAL CHEMISTRY

799

SOFTWARE NEWS AND UPDATES lated barrier heights are the main source of error. The most accurate tests are those concerned with tunneling splittings. The method has been very successful in calculating isotope effects and the effect of vibrational excitations on these splitting,13 – 15, 17 – 19 which are the effects most strongly influenced by the dynamics. In most cases reproduction of the absolute value of the splittings required adjustment of the quantum-chemically evaluated barrier heights. Comparison with the results of the dynamics methods and variation of the level of quantum chemistry used indicates that this is a quantum-chemical rather than a dynamics problem, because the program accounts satisfactorily for properties such as the kinetic isotope effect and the relative splitting of vibrationally excited levels, even when the absolute values are inaccurate. The program, therefore, offers the option of determining the barrier height that yields the best fit to the observed dynamics. Unreliable barrier heights are a problem especially for hydrogen-bonded systems. This diagnosis calls for a greater effort to test calculated barrier heights through the evaluation of tunneling splittings and rate constants. The harmonic-oscillator approximation for transverse modes is essentially unavoidable for directdynamics methods requiring the output of standard quantum-chemical programs. The linear coupling approximation for couplings with the tunneling mode is also common and is, in fact, equally unavoidable, because without it, the integral over the modes in the instanton formalism cannot be solved analytically, with the result that the treatment becomes unwieldy. The model does not include mixing of transverse modes, except for modespecific properties; this is not considered serious for temperature-dependent properties, because this mixing would only redistribute couplings between similar modes without much effect on the overall coupling. The coupling of low-frequency antisymmetric modes enters through a Franck–Condon factor, the accuracy of which should be good unless the factor is so small that anharmonic effects become important. This has not been the case in the examples studied. The coupling to symmetric modes enters in a form derived from an exact solution of a simple two-dimensional model; it has successfully predicted the tunneling splittings of vibrationally excited levels of these modes,13 – 15, 17 – 19 although thus far the applications are limited to low vibrational quantum numbers. The one-dimensional instanton action to which these couplings apply is calculated very accurately. The pre-exponential factor in the instanton rate expression is also an approxi-

800

mation, but the error is expected to be negligible.28 The method is not suitable for reactions where variational effects are significant. For many purposes, for instance biological applications, the absolute accuracy of the calculated rate constants is not an important issue. It is more significant to obtain semiquantitative insight into processes and mechanisms. The instanton approach is particularly well suited for this purpose because it gives a physically accurate picture of the tunneling event and the role of various coupled vibrations in helping or hindering the transfer. It has also been successful in describing transfers involving two and three protons,20, 22 an example of which is presented earlier. The observation that these protons can transfer either coherently or incoherently is of special interest for biological processes, where strings of water molecules may function as proton conduits. Our calculations have shown that such proton transfer at biological temperatures is likely to be dominated by quantum-mechanical tunneling, and may involve the coherent transfer of several protons.42 The program permits calculation of these transfer rates at arbitrary temperatures. We, therefore, conclude that the method and program described here provide a promising new approach to tunneling rates and tunneling splittings in molecules43 and molecular assemblies.44 – 46 The results are at least as accurate as those of any other program presently available for this purpose, and allow a complementary range of practical applications, partly because there are additional options, such as the calculation of vibronic tunneling splittings, and partly because there is much less demand on computer time, which allows applications to larger systems.

Supplementary Material A standard output of DOIT 1.2 for the rate constant and tunneling splitting calculations of the water trimer is available as a test result in Ref. 29.

References 1. Truhlar, D. G.; Garrett, B. C.; Klippenstein, S. J. J Phys Chem 1996, 100, 12771. 2. Truhlar, D. G.; Isaacson, A. D.; Garrett, B. C. Theory of Chemical Reaction Dynamics; Baer, M., Ed.; CRC Press: Boca Raton, FL, 1985; vol. 4. 3. Garrett, B. C.; Truhlar, D. G. J Chem Phys 1983, 79, 4931. 4. Miller, W. H.; Handy, N. C.; Adams, J. E. J Chem Phys 1980, 72, 99.

VOL. 22, NO. 7

SOFTWARE NEWS AND UPDATES 5. Skodje, R. T.; Truhlar, D. G.; Garrett, B. C. J Phys Chem 1981, 85, 3019.

32. Bunker, P. R.; Jensen, P. Molecular Symmetry and Spectroscopy; NRC Research Press: Ottawa, 1998; 2nd ed.

6. Garrett, B. C.; Joseph, T.; Truong, T. N.; Truhlar, D. G. Chem Phys 1989, 136, 271.

33. Hontscha, W.; Hänggi, P.; Polak, E. Phys Rev B 1990, 41, 2210.

7. Waite, B. A.; Miller, W. H. J Chem Phys 1980, 73, 3813.

34. Onsager, L. J Am Chem Soc 1936, 58, 1486.

8. Makri, N.; Miller, W. H. J Chem Phys 1989, 91, 4026.

35. Loerting, T.; Liedl, K. R.; Rode, B. M. J Chem Phys 1998, 109, 2672. 36. Pugliano, N.; Saykally, R. J. Science 1992, 257, 1938.

9. Sewell, T. D.; Guo, Y.; Thompson, D. L. J Chem Phys 1996, 103, 8557. 10. Guo, Y.; Thompson, D. L. J Chem Phys 1996, 105, 1070.

37. Gregory, J. K.; Clary, D. C. J Chem Phys 1995, 102, 7817.

11. Guo, Y.; Thompson, D. L. J Chem Phys 1996, 105, 7480.

38. Liedl, K. R.; Sekusak, S.; Kroemer, R. T.; Rode, B. M. J Phys Chem A 1997, 101, 4707.

12. Guo, Y.; Li, S.; Thompson, D. L. J Chem Phys 1997, 107, 2853. 13. Fernández-Ramos, A.; Smedarchina, Z.; Zgierski, M. Z.; Siebrand, W. J Chem Phys 1998, 109, 1004. 14. Smedarchina, Z.; Siebrand, W.; Zgierski, M. Z.; Zerbetto, F. J Chem Phys 1995, 102, 7024. 15. Smedarchina, Z.; Caminati, W.; Zerbetto, F. Chem Phys Lett 1995, 238, 279. 16. Smedarchina, Z.; Siebrand, W.; Zgierski, M. Z. J Chem Phys 1995, 103, 5326. 17. Smedarchina, Z.; Siebrand, W.; Zgierski, M. Z. J Chem Phys 1996, 104, 1203. 18. Smedarchina, Z.; Fernández-Ramos, A.; Rios, M. A. J Chem Phys 1997, 106, 3956. 19. Smedarchina, Z.; Zerbetto, F. Chem Phys Lett 1997, 271, 189. 20. Smedarchina, Z.; Zgierski, M. Z.; Siebrand, W.; Kozlowski, P. M. J Chem Phys 1998, 109, 1014. 21. Siebrand, W.; Smedarchina, Z.; Zgierski, M. Z.; FernándezRamos, A. Int Rev Phys Chem 1999, 18, 5. 22. Fernández-Ramos, A.; Smedarchina, Z.; Siebrand, W.; Zgierski, M. Z.; Rios, M. A. J Am Chem Soc 1999, 121, 6280.

39. Becke, A. D. J Chem Phys 1993, 98, 5648. 40. Kendall, R. A.; Dunning, T. H., Jr.; Harrison, R. J. J Chem Phys 1992, 96, 6796. 41. Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Zakrzewski, A. G.; Montgomery, J. A., Jr.; Stratmann, R. E.; Burant, J. C.; Dapprich, S.; Millam, J. M.; Daniels, A. D.; Kudin, K. N.; Strain, M. C.; Farkas, O.; Tomasi, J.; Barone, V.; Cossi, M.; Cammi, R.; Mennucci, B.; Pomelli, C.; Adamo, C.; Clifford, S.; Ochterski, J.; Petersson, G. A.; Ayala, P. Y.; Cui, Q.; Morokuma, K.; Malick, D. K.; Rabuck, A. D.; Raghavachari, K.; Foresman, J. B.; Cioslowski, J.; Ortiz, J. V.; Stefanov, B. B.; Liu, G.; Liashenko, A.; Piskorz, P.; Komaromi, I.; Gomperts, R.; Martin, R. L.; Fox, D. J.; Keith, T.; Al-Laham, M. A.; Peng, C. Y.; Nanayakkara, A.; Gonzalez, C.; Challacombe, M.; Gill, P. M. W.; Johnson, B.; Chen, W.; Wong, M. W.; Andres, J. L.; Gonzalez, C.; Head-Gordon, M.; Replogle, E. S.; Pople, J. A. Gaussian 98, Revision A.3; Gaussian, Inc.: Pittsburgh, PA, 1998. 42. Smedarchina, Z.; Siebrand, W.; Fernández-Ramos, A.; Gorb, L.; Leszczynski, J. J Chem Phys 2000, 112, 566.

23. Miller, W. H. J Chem Phys 1975, 62, 1899.

43. Zgierski, M. Z.; Grabowska, A. J Chem Phys 2000, 112, 6329; 113, 6845.

24. Callan, C. G.; Coleman, S. Phys Rev D 1977, 16, 1762; Langer, G. S. Ann Phys 1969, 54, 258; Polyakov, A. M. Nucl Phys 1977, 121, 121, 429.

44. Fernández-Ramos, A.; Smedarchina, Z.; Zgierski, M. Z. J Chem Phys 2000, 113, 2662.

25. Benderskii, V. A.; Goldanskii, V. I.; Makarov, D. E. Phys Rep 1993, 233, 195. 26. Benderskii, V. A.; Makarov, D. E.; Wight, C. H. Adv Chem Phys 1994, 88, 1. 27. Benderskii, V. A.; Goldanskii, V. I.; Makarov, D. E. Chem Phys 1991, 154, 407. 28. Benderskii, V. A.; Goldanskii, V. I.; Grinevich, P. G. Chem Phys 1993, 170, 275. 29. Smedarchina, Z.; Fernández-Ramos, A.; Zgierski, M. Z.; Siebrand, W. DOIT 1.2, a Computer Program to Calculate Hydrogen Tunneling Rate Constants and Splittings, National Research Council of Canada; http://gold.sao.nrc. ca/sims/software/doit/index.html. 30. Grabert, H.; Schramm, D. N.; Ingold, G. L. Phys Rep 1988, 168, 115. 31. Eckart, C. Phys Rev 19XX, 47, 552.

JOURNAL OF COMPUTATIONAL CHEMISTRY

45. Fernández-Ramos, A.; Smedarchina, Z.; Siebrand, W.; Zgierski, M. Z. J Chem Phys 2000, 113, 9714. 46. Fernández-Ramos, A.; Smedarchina, Z.; Rodríguez-Otero, J. J Chem Phys 2001, 114, 1567.

Z ORKA S MEDARCHINA1 A NTONIO F ERNÁNDEZ -R AMOS1,2 W ILLEM S IEBRAND1 1 Steacie Institute for Molecular Sciences National Research Council of Canada Ottawa, K1A 0R6 Canada 2 Department of Physical Chemistry Faculty of Chemistry University of Santiago de Compostela 15706 Santiago de Compostela, Spain

801