New Research In
Physical Sciences
Social Sciences
Featured Portals
Articles by Topic
Biological Sciences
Featured Portals
Articles by Topic
 Agricultural Sciences
 Anthropology
 Applied Biological Sciences
 Biochemistry
 Biophysics and Computational Biology
 Cell Biology
 Developmental Biology
 Ecology
 Environmental Sciences
 Evolution
 Genetics
 Immunology and Inflammation
 Medical Sciences
 Microbiology
 Neuroscience
 Pharmacology
 Physiology
 Plant Biology
 Population Biology
 Psychological and Cognitive Sciences
 Sustainability Science
 Systems Biology
Quantum mechanical reaction rate constants by vibrational configuration interaction: The OH + H_{2}→H_{2}O + H reaction as a function of temperature

Edited by Bruce J. Berne, Columbia University, New York, NY, and approved February 10, 2005 (received for review October 28, 2004)
Abstract
The thermal rate constant of the 3D OH + H_{2}→H_{2}O + H reaction was computed by using the flux autocorrelation function, with a timeindependent squareintegrable basis set. Two modes that actively participate in bond making and bond breaking were treated by using 2D distributed Gaussian functions, and the remaining (nonreactive) modes were treated by using harmonic oscillator functions. The finitebasis eigenvalues and eigenvectors of the Hamiltonian were obtained by solving the resulting generalized eigenvalue equation, and the flux autocorrelation function for a dividing surface optimized in reduceddimensionality calculations was represented in the basis formed by the eigenvectors of the Hamiltonian. The rate constant was obtained by integrating the flux autocorrelation function. The choice of the final time to which the integration is carried was determined by a plateau criterion. The potential energy surface was from Wu, Schatz, Lendvay, Fang, and Harding (WSLFH). We also studied the collinear H + H_{2} reaction by using the Liu–Siegbahn–Truhlar–Horowitz (LSTH) potential energy surface. The calculated thermal rate constant results were compared with reported values on the same surfaces. The success of these calculations demonstrates that timeindependent vibrational configuration interaction can be a very convenient way to calculate converged quantum mechanical rate constants, and it opens the possibility of calculating converged rate constants for much larger reactions than have been treated until now.
The evaluation of the thermal reaction rate constant k(T) from the quantum mechanical flux provides an efficient alternative to the computation of rate constants by means of the scattering matrix. The quantum mechanical formulation of k(T) in terms of flux autocorrelation functions C _{f}(t) was presented by Yamamoto (1) and Miller et al. (2, 3), and there have been several applications to calculate thermal rate constants for specific systems. The flux operator can be used to compute the cumulative reaction probability N(E) or the flux autocorrelation function, and either of these functions can be used to compute the thermal rate constant. Various approaches (4–24) involving basis functions, path integrals, and wave packet propagation methods have been used. Recently, Manthe and coworkers (19, 22) have calculated the thermal rate constants for the CH_{4} + H→CH_{3} + H_{2} and CH_{4} + O→CH_{3} + OH reactions by calculating N(E) as a function of energy E by using the multiconfiguration timedependent Hartree (MCTDH) method. Earlier work on triatomic reactions showed that accurate results can be obtained with an approach based on diagonalizing the timeindependent Hamiltonian (4, 8, 10). One advantage of this formulation is that the variational principle is used to identify the relevant subspace of the basis set. This approach is appealing in terms of its generality and straightforward extension to larger systems, and it is extended to polyatomic reactions in this article.
In this work, we have used flux autocorrelation functions to compute the thermal rate constants of two benchmark reactions: collinear H + H_{2}→H_{2} + H (which is used as a test of our computer program) and fulldimensional OH + H_{2}→H_{2}O + H. Both of these reactions have been studied extensively in the past by using various potential energy surfaces (8, 11, 13, 16, 20, 22–37), although most calculations involve approximations in the dynamics. In this work, we have used a timeindependent squareintegrable (L ^{2}) basis set to represent the Hamiltonian and the flux operator, and we formulated the method in a way that should be applicable to general polyatomic reactions. The basis functions are expressed in terms of massscaled normal mode coordinates defined at the saddle point or at a variational transition state. We have used 2D distributed Gaussian functions to represent the two modes that actively participate in the bondforming and bondbreaking process. For the collinear H + H_{2}→H_{2} + H reaction, there are only two modes, and both of these modes were treated by using 2D distributed Gaussian functions. The use of distributed Gaussian functions allows us to saturate the basis space in the strong interactions region (i.e., on and around the transition state). A similar strategy has been used for calculating scattering matrices (38, 39).
The OH + H_{2}→H_{2}O + H reaction has become a benchmark reaction for fouratom systems. Recently, two potential energy surfaces (33, 34) have been developed for this reaction. We have used the Wu–Schatz–Lendvay–Fang–Harding (WSLFH) potential energy WSLFH surface (33) for our work, and the resulting thermal rate constants are compared with earlier wave packet calculations by Goldfield and Gray (35). This application will be of special interest because there are few (if any) reactions with four or more atoms for which results of two different methods for calculating converged thermal rate constants have been shown to agree. The demonstration that converged results can be obtained by timeindependentbasisfunction methods is of even more interest because the advantages of this approach have been largely overlooked.
There are six normal mode coordinates at the saddlepoint geometry, and the Hamiltonian is represented as a function of the six massscaled normal mode coordinates. Two stretching modes that represent the bondmaking and bondbreaking process are treated by using 2D distributed Gaussian functions, and the remaining four modes are treated by using harmonic oscillator (HO) basis functions. We form 6D basis functions by taking a direct product of the 2D distributed Gaussian functions with the HO functions, and the matrix elements of the Hamiltonian operator are evaluated in this basis. Because the 2D Gaussian functions are not orthogonal to each other, the overlap matrix is computed, and the generalized eigenvalue problem is solved to obtain the eigenvalues and eigenvectors. The eigenvalues and eigenvectors are used to compute the flux autocorrelation function and the thermal rate constant.
Methods
Quantum Mechanical Theory. The thermal rate constant can be expressed in terms of the quantum mechanical flux operator via the flux autocorrelation function (1–3). The derivation leading to this result for a bimolecular reaction has been presented earlier (2), and thus, we will simply summarize the important relations here. For further details, see refs. 3 and 8.
The expression for the symmetric flux operator F is where h̄ is Plank's constant divided by 2π, H is the Hamiltonian operator, θ is the Heaviside unit step function, s is the reaction coordinate, and s = 0 defines a dividing surface separating reactants from products. The flux autocorrelation function C _{f} at time t for a given temperature T is where Tr{} represents a quantum mechanical trace, and k _{B} is Boltzmann's constant. The thermal rate constant k(T) is related to the flux autocorrelation function by where σ is the symmetry number of the reaction, is the electronic degeneracy of the potential energy surface on which the reaction occurs (in this case, is 2), Φ^{R} is the distinguishableparticle reactant partition function per unit volume, and L(T)isthe Laplace transform of the cumulative reaction probability given by One method of computing the flux correlation function is to evaluate the trace in Eq. 2, in the basis formed from the eigenvectors of the Hamiltonian, and the resulting expression is as follows: where E_{i} and Ψ _{i} are the eigenvalues and eigenvectors of the Hamiltonian operator, respectively. Inserting the expression of the flux operator from Eq. 1, we can rewrite Eq. 6 as follows: where θ _{ij} is the matrix element of the Heaviside step function in the basis formed from the eigenvectors of the Hamiltonian operator.
In Eq. 3, the presence of σ indicates that we assume that L and Φ^{R} are calculated without considering identical particle symmetry, and the presence of d ^{TS} _{el} indicates that we assume reaction occurs on a single Born–Oppenheimer potential energy surface of degeneracy d ^{TS} _{el}. Carrying out the integral of Eq. 5 analytically for a finite upper limit t, we obtain the following: where The upper limit t should be chosen to be large enough that the correlation function has practically decayed to zero and the area under the C _{f} curve remains unchanged with time. Because total angular momentum J is a good quantum number, we calculate the contributions of each J value to L or I(t) separately. Therefore, we write
To evaluate the rate constant by using Eq. 9, we need the eigenvectors E_{i} and eigenvectors Ψ _{i} of the Hamiltonian operator. We expand the eigenvectors in a nonorthogonal basis as follows: where the overlap matrix S has the following elements: The eigenvalues were obtained by solving the following generalized eigenvalue equation: where H is the Hamiltonian matrix in the nonorthogonal basis, and c _{i} is the eigenvector with elements c_{ki} .
Collinear H + H_{2} System. The collinear H + H_{2} reaction has been studied extensively (6, 25, 28), and it is used here to validate the method. The details of the basis functions and the calculations are given in Supporting Appendix A and Supporting Appendix B, respectively, which are published as supporting information on the PNAS web site. The results agree with those calculated by scattering theory (28) within 1%.
OH + H_{2} Calculation. The three hydrogen atoms in the OH + H_{2} system were labeled as follows: The OH + H_{2} system has six vibrational degrees of freedom. Normal mode analysis was performed at the saddlepoint geometry, and the resulting frequencies are 3,675, 2,439, 1,191, 690, and 573 cm^{–1} for modes 1–5, respectively, and 1,210 i cm^{–1} for mode 6. The modes Q _{1}, Q _{3}, Q _{4}, and Q _{5} represent the spectator O—H stretch, the outofplane bend, inplane bend, and torsion, respectively. The two stretch modes that actively participate in the bondmaking and bondbreaking process are Q _{2} and Q _{6}. All Q_{m} (m = 1, 2,..., 6) are zero at the saddle point. The two reactive modes were represented by using 2D distributed Gaussian functions. The remaining four spectator modes were represented by HO functions. The 6D basis was formed by taking a direct product between 2D distributed Gaussian functions and the HO functions. Earlier work (31) has shown that the contribution of the vibrational angular momentum term is very small for this system and can be dropped from the expression of the Hamiltonian. Therefore, the Hamiltonian in massscaled (27) normal coordinates is where μ is the scaling mass. In these calculations, the value of μ was set at 1 atomic mass unit. The matrix elements of the Hamiltonian operator were evaluated by using 6D basis functions, and the resulting generalized eigenvalue equation was solved.
Dividing surface. All rate calculations (except for those used in the check calculations described in the last paragraph of Results for OH + H_{2} →H_{2}O + H) were carried out with a dividing surface optimized in reduced dimensionality calculations. The dividing surface was defined in terms of the reactive modes by s = 0, where where η and c are parameters. Variations in η and c allow rotation and translation of the dividing surface, respectively. We generated 15 such dividing surfaces by taking various combination of η and c parameters, and they were tested in 2D calculations; the details are provided in Supporting Appendix C, which is published as supporting information on the PNAS web site. In this section, we discuss only the optimal dividing surface, which is shown in Fig. 1. This surface is labeled D_{1} , and it has η = π/4 and c = 0. The figure also shows the dividing surface of conventional transition state theory; it separates the reactant from the product region near the saddlepoint region but does not separate the reactants from products in the region of large Q _{2} and small Q _{6}. As shown in Fig. 1, this problem is corrected by rotating the dividing surface by π/4.
Basis functions. A 2D cut of the 6D internal coordinate space was obtained by varying the R _{OB} and R _{BC} distances and keeping all of the other degrees of freedom fixed at the saddlepoint geometry. A square grid with spacing d was formed by placing M points, from a minimum R _{Min} to a maximum R _{Max}, along each axis. Of the total of M ^{2} points, all points with potential energy of ≤E _{Cut} were selected. The selected points were then transformed to normal coordinates and were used as centers for 2D distributed Gaussian functions with x = Q _{6} and y = Q _{2}. The number of 2D Gaussian functions formed is called N _{g}. All N _{g} Gaussian functions were assigned the same width parameter A such that the maximum overlap between any two Gaussian functions is ≤S _{Max}.
The number N of basis functions used for the 6D calculations is directly proportional to the number N _{g} of 2D Gaussian functions. It is desirable to have a small value of N _{g} so that the number of basis functions for the 6D calculations is affordable. Therefore, we carried out reduceddimensional calculations whose purpose was to optimize the 2D distributed Gaussian functions.
The reduceddimensional calculations were performed in the Q _{2} Q _{6} subspace, with the other vibrational degrees of freedom frozen at their saddle point values. The calculations were repeated with various sets of 2D Gaussian functions, and the parameters were varied to find small basis sets that yield converged results in 2D. Table 1 lists parameters for two different sets of 2D Gaussian basis functions, labeled as G_{1} and G_{2} , that yield converged flux autocorrelation functions in 2D for the temperature range of 300–1,000 K. The G_{2} set has a larger value for the cutoff energy E _{Cut} and a smaller value of grid spacing d and it contains more basis functions than the G_{1} set. The centers of the 243 distributed Gaussian functions in set G_{1} are shown in Fig. 1. This set was selected for use in 6D calculations.
The basis functions are functions of six normal coordinates and are given by where and ϕ _{nmk} is an HO function, with orders n_{mk} = 0,1,....Thefunction χ _{k} is a 2D Gaussian function centered at , with width parameter A; the details of the 2D Gaussian functions used in the present calculations are provided in Supporting Appendix A. Note that Eq. 17 does not include rotation. One could carry out calculations for nonzero total angular momentum J by multiplying Eq. 17 by a rotational function (40–43). In this article, we use basis functions only for J = 0. Because the J = 0 rotational eigenfunction is a constant, including rotation in Eq. 17 affects only the normalization. Note that is an eigenfunction of the 4D HO Hamiltonian defined at the saddlepoint geometry, with the following eigenvalue: where ω _{m} is the frequency of mode m at the saddle point. Note that the zero point energy of the 4D HO is included in the zero of energy. The 6D basis functions were formed by taking a direct product between N _{HO} 4D HO functions and N _{g} 2D Gaussian functions, resulting in N = N _{g} N _{HO} 6D basis functions. The procedure used to select the HO functions is discussed in Results for OH + H_{2} →H_{2}O + H.
Thermal rate constant calculations. The matrix elements needed in Eqs. 9 and 13 were computed by standard methods as explained in Supporting Appendix D, which is published as supporting information on the PNAS web site. The eigenvalues and eigenvectors of the Hamiltonian were obtained from Eq. 13, and the J = 0 contribution to the L^{J=0} (T) Laplace transform was obtained from Eqs. 8 and 9. The integral depends on the upper limit t of the integration, and we must choose the upper limit large enough that the results are converged (3, 8). However, because of the finite size of the basis, the integral does not really converge but only reaches a plateau. Fig. 2 shows C _{f}(t) and its time integral at 300 K for the time range of 0–60 fs. We see that in the 0–20 fs and 40–60 fs regions the value I(t) fluctuates with t, but it is very stable in the time range of 30–40 fs. One can use this plateau value of I(t) to compute thermal rate constants. Factors that affect the range of time over which I(t) is stable at its plateau value include the number of basis functions used for the calculation, the locations of the basis functions, and the position of the dividing surface. The type of basis functions is also an important factor; for example, it was shown in ref. 5 that 1D distributed Gaussian functions and sine functions are more efficient basis functions than HO functions for evaluating C _{f} for a 1D Eckart potential.
The value of the upper limit of time integration was determined in this work by finding the widest time interval over which I(t) is constant to within 1%. The center of this interval is called t _{f} and the width is called Δt. (Details of the algorithm for finding t _{f} are given in Supporting Appendix E, which is published as supporting information on the PNAS web site.) We take I(t _{f}) as L^{J=0} .
Rather than compute L by the exact relation of Eq. 10, weuse the separablerotation approximation (44), which yields the following: where A ^{TS}, B ^{TS}, and C ^{TS} are the rotational constants of the transition state evaluated at the saddlepoint geometry. The values of A ^{TS}/hc, B ^{TS}/hc, C ^{TS}/hc were 18.2 cm^{–1}, 2.81 cm^{–1}, and 2.44 cm^{–1}, respectively. These values are identical to the values used by Goldfield and Gray (35). It is known from refs. 20 and 35 that the separablerotation approximation overestimates the rate constant by ≈40% for this reaction.
Last, we consider three other factors in Eq. 3. The symmetry number σ is 2, the electronic degeneracy d ^{TS} _{el} is 2, and where Φ^{Rel} is the relative translational partition function per unit volume, and is the reactant electronic partition function given by where Δ is the spinorbit splitting of the OH radical and equals 140 cm^{–1} (45). is the product of two diatomic vibrational–rotational partition functions computed without considering identicalparticle symmetry or ^{–}nuclear spin (those effects are in σ). For consistency with Eq. 21, these partition functions are computed from the diatomic rotational constant and vibrational energy levels and diatomic rotational constant B ^{diat} by with Q ^{diat} _{vib} evaluated from the accurate diatomic energy levels for the given potential energy surface. The converged vibrational partition function for the reactants was computed by summing Boltzmann factors for the vibrational states whose energies were obtained by solving the 1D vibrational Hamiltonians for the reactant molecules by using HO basis functions.
Results
Results for OH + H_{2}→H_{2}O + H. In this section, we discuss the convergence with respect to the HO basis functions in the nonreactive degrees of freedom and compare the results to previous calculations.
The parameters for the basis functions used for the computation of thermal rate constants are shown in Table 2. The 4D spaces of basis sets B_{1} –B_{3} were obtained by a twostep procedure. the first In step, a cutoff parameter was defined, and all with were selected for the next step. Among these selected , only functions with three of the four n_{m} , with m = 1, 3, 4, and 5, equal to zero were selected. As an example, consider a basis set formed by using . There are 91 4D HO functions with excitations energies of ≤5,000 cm^{–1}. Of these 91 functions, there are only 21 functions with three or more n_{m} = 0. Basis sets B_{4} –B_{8} are variations on basis B_{1} . First, the maximum number of quanta in each of the four modes for the B_{1} set was labeled as , , , and , respectively. For B_{4} , three HO functions of the form with k = 1,..., 3 were formed and were combined with the 18 existing HO function in the B_{1} set to give a set of 21 HO functions. A similar procedure was carried out with , , and to yield B_{5} –B_{7} . Note that B_{4} –B_{7} all contain 5,103 6D functions, but they have different distributions of quanta in the four nonreactive modes. Table 3 shows good convergence of the thermal rate constants computed by using the B_{1} –B_{7} basis sets. Comparison of computed rate constant with experimental results (46, 47) and previous theoretical calculations (20, 35, 48) are shown in Fig. 3.
Convergence of the B_{1} set was also tested by using twomode excited HO functions in the basis set. There are the following six possible combinations for exciting any two modes at the same time: (n _{1}, n _{3}, 0, 0), (n _{1}, 0, n _{4}, 0), (n _{1}, 0, 0, n _{5}), (0, n _{3}, n _{4}, 0), (0, n _{3}, 0, n _{5}), and (0, 0, n _{4}, n _{5}). The three lowest energy states of the form (n _{1}, n _{3}, 0, 0) were combined with the 18 HO functions from the B_{1} set to obtain a new set of 21 HO functions. The direct product between these 21 HO functions and 243 distributed Gaussian functions was formed, and the resulting set of 5,103 functions was labeled as B_{8} . Similar calculations were also performed for the remaining four combinations, and the results for B_{9} –B_{13} were found to agree well with those for B_{1} –B_{7} . The largest deviation for any of the results with bases B_{4} –B_{13} from those obtained with our largest basis, B_{3} , is 1%.
If the wave function is represented adequately by using a large enough basis set, then the flux autocorrelation function decays to zero after a time of order (10, 17). As temperature goes down, one needs a larger basis set to converge the thermal rate constant. The basis sets used in the present calculations were found to be adequate in the temperature range of 290–1,100 K. For temperatures <290 K, we found that the time integral of flux autocorrelation function was oscillatory and incapable of providing converged rate constants. However, this phenomenon is due to finite size of the basis set and can be corrected by extending the boundaries of the basis set more deeply into the reactant and product regions.
We have also calculated the quantum mechanical rate constants at 300 and 500 K by using dividing surfaces based on variational transition state theory (VTST). The aim of the VTST calculations was to obtain the variational dividing surface to be used in the quantum mechanical rate constant calculations. The results obtained with the VTST dividing surfaces agree with those presented here within 1% for both temperatures. Fig. 4 shows plots of the flux autocorrelation function obtained by using the dividing surface at the saddle point and the variational dividing surface at 300 K. This comparison provides a numerical verification of the fact that the results are independent of the location of the dividing surface. Details of these check calculations are provided in Supporting Appendix F, which is published as supporting information on the PNAS web site.
The flux autocorrelation functions for the B_{1} set at 1,000 K is shown in Fig. 5, which is published as supporting information on the PNAS web site.
Not only are the results well converged, but they agree well with the wave packet results of ref. 35. Over the temperature range of 300–1,000 K, the largest deviation and the mean unsigned deviation for any of the B_{1} –B_{13} basis sets from the results obtained in ref. 35 were 7% and 2%, respectively. (In comparing with ref. 35, we compare to their separablerotation results because the emphasis in this article is on calculation of the Laplace transform and not on improving on the separablerotation approximation.)
A side product of the calculations presented here is that we can test VTST. VTST calculations with microcanonically optimized multidimensional tunneling (49) were presented by Troya et al. (48). Their rate constants are larger than the present converged rate constants by 54%, 40%, and 75% at 300, 500, and 1,000 K, respectively.
Although the calculations reported here are based on an analytic potential energy surface, the method can be used equally well with direct dynamics (50–53).
Conclusions
The thermal rate constant was computed for a fulldimensional fourbody reaction by using timeindependent squareintegrable basis functions, and the Laplace transform was found to be well converged and in good agreement with previous calculations based on timedependent wave packets. Rate constants for a system with more than three atoms have been calculated here by using timeindependent basis functions. It was shown that an efficient basis set can be formed by using Gaussian basis functions for the two active stretch modes and HO functions for the remaining modes. We found that the number of HO functions needed in the nonreactive degrees of freedom to get converged results is very small as compared with the number of 2D distributed Gaussian functions in the bondforming and bondbreaking modes.
Although much of the effort in quantum dynamics is currently focused on wave packet methods, the success of the present calculation opens the possibility of treating general multidimensional reactions by convenient timeindependent basis set methods, and we anticipate that further reduction in the computational cost could be achieved by combining the present methods with new schemes (54, 55) for reducing the size of Gaussian basis sets for nuclear motion and with hierarchical representations of the potential (40–43). This comparison of timeindependent to timedependent approaches bears elaboration. It is a general question of working in the time domain or the energy domain, which are complementary in a Fourier sense. In the early days of quantum mechanical collision theory, timeindependent methods, especially the closecoupling method, received almost all of the attention and effort (56–58). Later, attention turned more to timedependent quantum dynamics (59, 60), which is now considered by most workers to be the method of choice for accurate polyatomic dynamics (61, 62). Timedependent quantum mechanics is especially efficient for generating results at a series of total energies, as required for thermally averaged rate constants, because a single wave packet carries information about a wide range of energies (60). However, the present application shows that timeindependent quantum mechanics can also be used for accurate polyatomic reaction dynamics. It is particularly important to note that the method used here is not special to fourbody reactions, and in fact, it is based on straightforward use of vibrational configuration interaction. Thus, by taking advantage of recent advances in vibrational configuration interaction calculations (41  43) and the fast convergence of the calculations with respect to completing the basis in nonreactive degrees of freedom, it should be possible to extend the method used here to much larger systems. In this regard, we note that use of the flux autocorrelation method (1  3) is a general method for taking advantage of the fact that reaction rates are often dominated by shorttime dynamics in a localized region around a transition state. This makes timedependent methods like the timedependent Hartree method (18–20, 22, 61–63) more affordable because the wave packet needs to be represented over only a short period, and timeindependent methods like vibrational configuration also benefit greatly by the need to represent the wave function over only a localized region of space. We anticipate that vibrational–rotational configuration interaction calculations will provide a powerful general tool for calculating flux autocorrelation functions for many other polyatomic reactions.
Acknowledgments
We thank Stephen Gray for sending the reaction rate constants of ref. 35 in tabular form. This work was supported in part by National Science Foundation Grant CHE0349122.
Footnotes

↵ * To whom correspondence should be addressed. Email: chakra{at}comp.chem.umn.edu.

Author contributions: D.G.T. designed research; A.C. and D.G.T. performed research; A.C. and D.G.T. analyzed data; and A.C. and D.G.T. wrote the paper.

This paper was submitted directly (Track II) to the PNAS office.

Abbreviations: HO, harmonic oscillator; VTST, variational transition state theory.
 Copyright © 2005, The National Academy of Sciences
References
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵

Mielke, S. L., Peterson, K. A., Schwenke, D. W., Garrett, B. C., Truhlar, D. G., Michael, J. V., Su, M.C. & Sitherland, J. W. (2003) Phys. Rev. Lett. 91 , 1–4.

↵
Zhang, J. Z. H., Li, Y. M., Wang, L. & Xiang, Y. (2004) in Modern Trends in Chemical Reaction Dynamics, eds. Yang, X. & Liu, K. (World Scientific, Singapore), pp. 209–248.
 ↵
 ↵
 ↵
 ↵

Carter, S., Bowman, J. M. & Handy, N. C. (1998) Theor. Chem. Acc. 100 , 191–198.
 ↵
 ↵

↵
Herzberg, G. (1950) Molecular Spectra and Molecular Structure. I. Spectra of Diatomic Molecules (Van Nostrand, Princeton).
 ↵

↵
Baulch, D. L., Cobos, C. J., Cox, R. A., Esser, C., Frank, P., Just, Th., Kerr, J. A., Pilling, M. J., Troe, J., Walker, R. W. & Warnatz, J. (1992) J. Phys. Chem. Ref. Data 21 , 411–734.
 ↵
 ↵
 ↵

Truhlar, D. G. & Gordon, M. S. (1990) Science 249 , 491–498.

↵
Truhlar, D. G. (1995) Understanding Chem. React. 16 , 229–255.
 ↵
 ↵

↵
Bernstein, R. B., ed. (1979) AtomMolecule Collision Theory (Plenum, New York).

Clary, D. C., ed. (1986) The Theory of Chemical Reaction Dynamics (Reidel, Dordrecht, The Netherlands).

↵
Truhlar, D. G. (1994) Comput. Phys. Commun. 84 , 79 90.

↵
Kulander, K. C., ed. (1991) TimeDependent Methods for Quantum Dynamics (Elsevier, Amsterdam).

↵
Wyatt, R. E. & Zhang, J. Z. H., eds. (1996) Dynamics of Molecules and Chemical Reactions (Dekker, New York).
 ↵

↵
Meyer, H.D. & Worth, G. A. (2003) Theor. Chem. Acc. 109 , 251–267.
 ↵