Quantum Instanton Evaluations of the Thermal Rate Constants for
Complex Systems

Zhao, Yi; Wang, Wenji

doi:https://doi.org/10.1155/2012/483504

Advances in Physical Chemistry

On this page

Abstract Introduction Results Conclusion Acknowledgments References Copyright Related Articles

Special Issue

Accurate Potential Energy Surfaces and Beyond: Chemical Reactivity, Binding, Long-Range Interactions, and Spectroscopy

View this Special Issue

Review Article | Open Access

Volume 2012 | Article ID 483504 | https://doi.org/10.1155/2012/483504

Quantum Instanton Evaluations of the Thermal Rate Constants for Complex Systems

Yi Zhao¹and Wenji Wang²

Academic Editor: Laimutis Bytautas

Received28 Jun 2011

Revised26 Aug 2011

Accepted10 Sept 2011

Published11 Dec 2011

Abstract

Quantum instanton (QI) approximation is recently proposed for the evaluations of the chemical reaction rate constants with use of full dimensional potential energy surfaces. Its strategy is to use the instanton mechanism and to approximate time-dependent quantum dynamics to the imaginary time propagation of the quantities of partition function. It thus incorporates the properties of the instanton idea and the quantum effect of partition function and can be applied to chemical reactions of complex systems. In this paper, we present the QI approach and its applications to several complex systems mainly done by us. The concrete systems include, (1) the reaction of , (2) the reaction of , (3) H diffusion on Ni(100) surface; and (4) surface-subsurface transport and interior migration for H/Ni. Available experimental and other theoretical data are also presented for the purpose of comparison.

1. Introduction

The accurate and efficient evaluation of chemical reaction rate constant is one of prime objectives of theoretical reaction dynamics. Since rigorous quantum mechanical approaches are limited to small molecular (several atoms) reactions, a variety of approximation approaches have been proposed. Benefited from the small recrossing dynamics at not-too-high temperatures, the transition state theories (TSTs), originally proposed by Eyring [1, 2] and Wigner [3],have become a possible and popular way to estimate rate constants. Due to their practical simplicity, they have been broadly applied to numerous reactions. The TST is inherently a classical theory and suitable at sufficiently high temperatures, where the classical description of nuclear motions may be adequate. At low temperatures, especially for the reactions involving the motions of light atoms (i.e., hydrogen), however, quantum effects become quite significant. To make the TST still valid for such low temperature reactions, many approaches have been proposed to quantize it [4–14]. However, there is no absolutely unambiguous way to do it.

To develop a more accurate and less ad hoc quantum version of TST, with a specific focus on the tunneling regime, Miller et al. [15–17] have proposed a quantum instanton (QI) approach recently. The QI is based on an earlier semiclassical (SC) TST [18] that became known as the instanton [19, 20]. The similarity between the QI and SC instanton lies in using the steepest descent approximation to evaluate relevant integrals in the quantum rate formula, while the crucial difference is that the Boltzmann operator is evaluated by the quantum mechanics and semiclassical approximation, respectively. The QI theory thus incorporates the tunneling, corner cutting [21–24], and anharmonicity correctly and is expected to overcome the quantitative deficiency of the SC instanton model. In particular, the QI theory considers all tunneling paths and automatically gives each path its naturally weight factor from the quantum Boltzmann operator, instead of choosing a single optimal tunneling path, which is taken into account in the SC instanton and TSTs with SC tunneling corrections. Indeed, it has been numerically demonstrated that the QI predicts accurate quantum rates for one-dimensional and two-dimensional models within 20% error over a wide temperature range, from the deep tunneling to overbarrier regimes.

A lot of developments and applications [25–39] have been made since the QI theory is proposed. The original QI [15] involves the second time derivative of the flux-flux correlation at time zero. It has been further improved [29] by taking into account the higher derivatives of the flux-flux function. For the 1D and collinear reactions, the improved model is considerable accuracy, giving the rates to within 5%–10% errors. For a practical purpose, a simple and general way for choosing dividing surfaces used in the QI is suggested [31], namely, using the family of (hyper)planes normal to the minimum energy path at various distances s. A “simplest” QI model [26] has also been suggested with one dividing surface, slightly less accurate than the original QI. To reveal the relationship with conventional TSTs, the classical limit of the QI has been derived [37]. It is found that the classical TST is just a special case of the QI in high temperature limit; moreover, the quantum correction of the prefactor is more important than that of the activation energy in the TST.

Since the QI solely involves the Boltzmann operator and its relevant quantities, it can be applied to quite complex molecular systems (from gas phase [17, 37], liquid [28], to surface [38, 39]) via well-established imaginary time path integral techniques. The first implementation of QI with path integral Monte Carlo and adaptive umbrella sampling techniques is applied to the three-dimensional hydrogen exchange reaction D + H₂ → HD + H [16]. Soon, the techniques are further extended to the reaction of [17]. The thermodynamic integration with respect to the mass of the isotope and the inverse temperature is also proposed to compute the kinetic isotope effects [30] and rate constants [34] directly. To improve the convergence of the Monte Carlo simulation, the efficient “virial” estimators [32] have been derived from the logarithmic derivatives of the partition function and the delta-delta correlation functions, and it is found that their statistical errors do not increase with the number of discrete time slices in the path integral. Most recently, the QI has been compared with other conventional approaches [36] for an intramolecular proton transfer on a full-dimensional potential energy surface that incorporates high-levels ab initio calculations along the reaction path. The obtained kinetic isotope effects from the QI are in reasonable agreement with those from the path-integral quantum TST.

In this paper, we firstly illustrate the QI formula and its path integral representation. Then, we display several applications, which are mainly done by ourselves. The systems include two gas phase reactions [17] and [37], H diffusion on Ni(100) surface [38], and surface-subsurface transport and interior migration for H/Ni [39].

2. Method

In this section, we summarize the rate formula for the QI evaluation. The detailed derivation can be found in [15–17]. The QI model proposes the following thermal rate constant: Here, is the reactant partition function per unit volume. is zero time value of the flux-flux correlation function where is the inverse temperature , is the Hamiltonian operator of the reaction system, and and are the flux operators given by with . In (3), is the step-side function, represents the Cartesian coordinates of the reaction system, and and define two separate dividing surfaces via the equations and , both and being positive (negative) on the product (reactant) side of the dividing surfaces.

in (1) is a specific type of energy variance given by In order to get the correct free particle (high temperature) limit (that would be 25% too large otherwise), an ad hoc term is added to , which has very little effect in the low temperature regime. and are zero time value and its second time derivative, respectively, of the “delta-delta” correlation function where the generalized delta-function operator is Here, is the total number of atoms, , denotes the Cartesian coordinates of the th atom, and is its atomic mass.

The dividing surfaces are determined by the stationary condition where is a collection of parameters that is involved in the location of the dividing surfaces. This condition originates from the SC instanton model, and the resulting dividing surfaces correspond qualitatively to the turning points of the periodic orbit that runs on an upside down PES in imaginary time (see Appendix A in [15, 18]).

Since all the relevant quantities in the QI expression (1) involved only the quantum Boltzmann operator, they can be readily evaluated using imaginary time path integral Monte Carlo (PIMC) [40] method.

We begin with the simplest quantity, , which can be discretized according to the standard procedure [41–43] to give where is a multiplicative constant, the number of imaginary time slices, and the Cartesian coordinates of the system associated with the sth time slice. is the discretized action given by where and represents . Path integral expressions for and are somewhat more complicated but can be obtained in a straightforward manner. The appropriate expressions are with with with being the total number of degrees of freedom (i.e., ).

In realistic calculations, we rewrite (1) as the product of several ratios

The terms of and are directly calculated as a constrained average over the same ensemble of paths [16, 17] with

The evaluation of , however, meets a challenge because the is the quantity associated with the transition state, while with the asymptotic reactant domain, we evaluate it using adaptive umbrella sampling techniques [44].

The QI approximation uses the short-time information of the flux-flux correlation function. Predescu and Miller [45] demonstrate that in the classical limit, Wigner’s variational principle and the quantum variational criterion based on the minimization of flux-flux correlation function produce the same optimal surface. Recently, Wang et al. [37] have shown that in the classical limit, the QI formula is the same as the classical transition state theory. These situations motivate us to write the QI formula (13) as the Arrhenius form Here, the free energy is defined by This two-dimensional free energy is related to (the reactant partition function) and , and it is corresponding to the quality of probability density at . has a similar property to the partition function at the transition state. The prefactor is given by The advantage of (16) allows us to investigate the respective quantum contributions to the rates from the quantum free energy and prefactor by comparing their quantum and classical values, since other factors such as the vibrational-rotational coupling and anharmonicity are automatically involved.

3. Results and Discussions

3.1. The Reaction of

The reaction is a prototype of polyatomic hydrogen abstraction reaction. Quantum dynamical studies of this reaction have become possible only recently, because it involves 12 internal degrees of freedom and thus poses difficulties to quantum dynamics calculations as well as construction of the potential surface. We apply the QI methodology to this reaction using the potential energy surface constructed by Espinosa-García [46]. All calculations are performed in terms of the Cartesian coordinates of all the atoms (i.e., 18 degrees of freedom).

In the path integral simulations, the number of imaginary time slices is chosen to be 20 and 100 at temperatures K and 200 K, respectively, while Monte Carlo cycles are run to achieve <10% statistical convergence.

3.1.1. Reaction Coordinate

For this reaction, we define a generalized reaction coordinate , where is an adjustable parameter that shifts the location of the dividing surface (defined by ). is defined by a linear interpolation between two constituent reaction coordinates and through the parameter , is a reaction coordinate whose dividing surface is designed to pass through the top of the classical potential barrier, which is defined here as with being the reaction coordinate that describes the abstraction process of one of the methane hydrogens by the incident one where denotes the interatomic distance between atoms and and is the value at the transition state geometry. , on the other hand, describes a dividing surface that is located far in the asymptotic reactant valley, which is given by Here, R is the scattering vector that connects the incident hydrogen and the center of mass of the methane. is an adjustable parameter which is chosen to be 9 Å in order to guarantee that the interaction potential energy between H and is negligible.

Now, the term (5) becomes a function of two parameters, and , as follows: and thus, one seeks a stationary point of in the two-dimensional space to obtain the corresponding “optimal” values.

3.1.2. Free Energy Surface

The quantity varies exponentially as a function of , it is convenient to define a quantum “free energy surface” as follows: and locate the saddle point of by visual inspection.

Figure 1 shows the free energy surfaces for the reaction. It exhibits a barrier-like profile along the direction , while it grows approximately quadratically with the increasing of the absolute value of . From this figure, it is seen that at a higher temperature (Figure 1(a)), there appears only a single saddle point at , while at a low temperature (Figure 1(b)), the saddle point bifurcates into a distinct pair at , which indicates the existence of nonnegligible tunneling effects in the rate constant [15, 16].

(a)

(b)

3.1.3. Rate Constants

Having obtained the “optimal” values of the at each temperature, one can now compute the quantum instanton rate by combining various quantities as in (8). The calculated QI rates are tabulated in Table 1 as well as other theoretical and experimental ones.

Comparing the quantum instanton rates with others, we find that is in good agreement with the experimental data more specifically, it is closer to the rates obtained by Baulch et al. [47] than those by Sutherland et al. [48], and agrees with within 10% for the temperature range T = 600–2000 K, but it becomes somewhat larger than the latter as the temperature is decreased (the deviation becomes 30% and 55% for and 300 K, resp.). It should be noted that the differences between and are much smaller than the uncertainty of the experimental data.

Figure 2 displays the Arrhenius plots of the rate constants. In Figure 2, we also plot the accurate quantum dynamics results of MCTDH (the multiconfigurational time-dependent Hartree approach) [49]. Compared to the MCTDH ones, our QI rate constants are larger by factors of about 2 to 3 over the temperature range 300–400 K. This difference may partly be due to the recrossing effect which is not considered in QI theory and partly arise from the use of the J-shifting approximation and the neglect of the vibrational angular momenta Hamiltonian in MCTDH method [49].

3.2. The Reaction of

The reaction is an important step in the radical mechanism of thermal decomposition of monosilane. We calculate the rates and kinetic isotope effects (KIEs) of this reaction with the quantum instanton method in full Cartesian space, on the basis of analytical potential energy surface constructed by Espinosa-García et al. [50].

The reaction coordinate of has the same form as that of (Section 3.1.1). In our QI calculations, the number of time slices, P, is set to be 20 (1000 K) and 120 (200 K) for the quantum evaluations. In our classical evaluations, the formula is the same as QI, but the number oftime slices in the path integral is set to be 1. The number of Monte Carlo is about for computing a single ensemble average, and it converges the relevant quantities within 10% statistical errors.

3.2.1. Free Energy and Prefactor

We have rewritten the QI formula (8) in the Arrhenius form (16), which consists of the free energy and prefactor. In this section, we calculate the corresponding quantum and classical quantities so as to investigate the respective quantum contributions to the rates from the quantum free energy and prefactor. Our calculated results are plotted in Figures 3 and 4.

Now, we look into the quantum effects from both the free energy and prefactor. Figure 3 plots the quantum and classical free energies as a function of temperature at the optimized stationary point. One immediately observes that the quantum effect becomes significant at T < 600 K, whereas both quantum and classical results nearly coincide at T > 600 K. As expected, the quantum effect always decreases the classical free energy. The difference of the contributions to the rates is about 40% at 200 K.

Figure 4 shows the temperature dependence of the quantum and classical prefactors. Compared with Figure 3, one easily finds that this quantum contribution is much larger than that for the free energy correction. Even at 1000 K, the difference is observable. At 200 K, the difference is about several orders of magnitude.

In the QI theory, it is not possible to explicitly distinguish the quantum contributions from the partition function and nuclear tunneling, but we can conclude that it is insufficient to estimate the accurate rate by only replacing the activation energy in the TST with its quantum analog, because the quantum prefactor plays more important rule in determining the quantum rate.

3.2.2. Rate Constants

In Table 2, we summarize the rate constants obtained from the QI theory and classical VTST (this kind of results are performed by setting the time slices in the path integral to be 1 in the QI calculations) as well as those from the canonical variational TST with the centrifugal-dominant small curvature SC adiabatic ground-state (CVT/CD-SCSAG) approach [50], the conventional TST with simple Wigner tunneling factor [51] and the experiment [52]. The corresponding Arrhenius plots are displayed in Figure 5.

Table 2 and Figure 5 display comparable results of the CVT/CD-SCSAG to the QI values, with 36% maximal errors at 200 K and a slightly different Arrhenius slope. Both results are in good agreement with experimental data in the tested temperature range. This manifests that the PES used is reasonable accurate. However, the values from the classical VTST are always smaller than the QI results especially in the deep tunneling regions. The conventional TST results have similar tendency to those from the classical VTST. Although the classical VTST can be much improved by evaluating the partition functions quantum mechanically, we do not focus on this improvement. However, it is noted that the anharmonicity, rotational-vibrational coupling are involved in the classical simulation. It is thus expected that these large errors come from the pure quantum effects.

3.2.3. Kinetic Isotope Effects

Kinetic isotope effect (KIE) is the characteristic of chemical reactions which may reveal the quantum effect. We consider the following isotopic reactions: The calculated values are tabulated in Table 3.

The KIEs of and in the temperature range of 200–1000 K are summarized in Table 3. Besides the QI values, this table also displays comparable results of the CVT/CD-SCSAG [50]. It is easily found that the KIEs predicted by the QI are smaller than 1, in agreement with the CVT/CD-SCSAG values. The detailed comparison reveals that the CVT/CD-SCSAG predicts smaller values than those from QI theory. Again, the maximal error occurs at 200 K and is about 23%. For KIEs, we find that although both approaches predict the values larger than 1 at 200 to 1000 K, the QI values are smaller than those from CVT/CD-SCSAG. Espinosa-García and coworkers [50] have pointed out that the CVT/CD-SCSAG rates may have been overestimated because of the high vibrational and tunneling contributions. This manifests that the QI approach indeed correctly accounts for the quantum effects.

3.3. H Diffusion on Ni(100) Surface

Diffusion plays a fundamental role in surface process. It reveals characteristics about the underlying surface potential and is intimately involved in determining the kinetics of surface catalyzed chemical reactions. The hydrogen atom and its isotopes are ideal candidates to exhibit quantum tunneling behavior due to their small masses. We explore the evaluation of the quantum instanton approximation to the process of H diffusion on Ni(100) surface using the EAM4 potential energy surface constructed by Truong and Truhlar [53].

In the path integral calculations, the numbers of time slices, and , for the degrees of freedom of H and quantized Ni atoms, respectively, are set to (24–40, 6–8) in the temperature range of 40–300 K. The number of Monte Carlo cycle employed is (1–10) × , which converges most of the relevant quantities within 10% statistical error (some of the statistical errors are within 20% below 100 K).

3.3.1. Model

Nickel crystallizes in a face-centered-cubic (fcc) lattice structure. The structural model used consists of 162 Ni atoms over four layers: 40 atoms are in each of the first and third layers, and 41 atoms are in each of the second and fourth layers.

Figure 6 shows the platform and profile chart of Ni(100) lattice. In the simulations, 89 Ni atoms (orange), in the four sides of the bulk metal and the bottom layer, are fixed. The 48 atoms (blue) closest to the fixed ones are considered to be movable and treated classically. The last 25 ones (red), surrounding the reactant and product sites and lying directly beneath the reactant and product sites, are treated quantum mechanically.

In our calculations, the lattice (see Figure 6) is rotated by 45 degrees counterclockwise, and we chose the coordinate of H atom (along the reaction path) to be reaction coordinate.

3.3.2. Probability Distribution of Paths

The path integral has an advantage that the paths of the particles can display the character of the diffusive process. We extract the probability distributions for the paths of H and Ni atoms (two reaction coordinate beads ( and ) of H path are fixed at the two dividing surfaces) for the purpose of guaranteeing the instanton property (the instanton is a periodic orbit between the two dividing surfaces). The dividing surfaces can be obtained from the free energy surface (24), their values at different temperatures are shown in Table 4.

Figure 7 displays the probability distributions for both H and Ni atoms at several temperatures. It is found that the probability distribution of H is localized at 300 K, and it becomes more and more delocalized as the temperature decreases. Below 80 K, the distribution rapidly becomes obvious in both sides of the transition state, and the path begins to continuously distribute between the two hollow sites at lower temperatures. This phenomenon indicates that H tunneling becomes remarkable at the temperatures lower than 80 K. The probability distributions of Ni atoms have small changes in the whole temperature range, and they seem to be frozen at very low temperature.

3.3.3. Diffusion Coefficients

According to the hopping model [54], the diffusion coefficient D is related to the rates through where b is the hop length and is equal to 2.489 Å for the Ni(100) surface. We have calculated the diffusion coefficients at temperatures ranging from 40 to 300 K, and the results are plotted in Figure 8.

We plot two kinds of QI diffusion coefficients in Figure 8, one is the result on a rigid lattice, the other is the result on a quantized lattice. At high temperatures, these two kinds of diffusion coefficients are nearly the same; however, at low temperatures, the ones on a rigid lattice are bigger than that on a quantized lattice, and this is mainly due to the fact that the free energy on a quantized lattice is higher than that on a rigid lattice [38]. Comparing the QI diffusion coefficients with the CVT/SCSAG ones [55] on a rigid lattice, we find that our QI results on a rigid lattice are very similar to those from the CVT/SCSAG. Now, we compare the present theoretical results with experimental ones. It is found that the QI results as well as the CVT/SCSAG results are in good agreement with the experiments of George et al. [56] and Mullins et al. [57] at 200 K and 300 K. However, at low temperatures, the theoretical results are significantly larger than Lee et al. [58] and Lin and Gomer [59]. The experimental transition temperature (100 K–160 K) is also different from the present calculations. The QI predicts it to be 70 K for the quantized lattice, while it is about 66 K for the rigid lattice obtained from both the QI and the CVT/SCSAG approaches. These discrepancies between theories and experiments may be attributed to the interaction potential. It is likely that the accuracy of the potential energy surfaces needs further improvement.

3.4. Surface-Subsurface Transport and Interior Migration for H/Ni

A full-dimensional potential energy surface (EAM6), constructed by Wonchoba and Truhlar [60], using the embedded diatomics-in-molecules (EDIMs) [61] potential energy function, is used in this QI calculations.

In path integral calculations, the numbers of time slices, and for the degrees of freedom of the H and quantized Ni atoms, respectively, are set to in the temperature range of 100–400 K. The number of Monte Carlo is about for computing a single ensemble average. It converges most of the values within 10% statistical errors (some of the statistical errors are within 20% at 100 K).

3.4.1. Model

Figures 9 and 10 show H diffusion processes in the nickel crystal with a face-centered-cubic (fcc) lattice structure. For a given rate process, we construct a lattice cell (each cell contains more than 200 Ni atoms) in which all atoms are treated to be movable to incorporate the effect of the crystal fluctuation on the rates. To be concrete, the Ni atoms in the four sides of the bulk metal and at the bottom layers are fixed, the Ni atoms surrounding the reactant and product sites and lying along the reaction path are treated quantum mechanically, and the others are treated classically. It should be noted that the lattice for the interior diffusion process has a structure of sphere and the outer layers are fixed.

In the QI calculations, we need to define the reaction coordinate operators (6). For the systems considered in this section, the hydrogen coordinates are essentially good choices. We thus adopt the following reaction coordinates for different rate processes. In the H diffusion on Ni(111) surface from a hcp site to a fcc site (), the coordinate of the H atom is chosen, whose direction is showed in Figure 9 as the black line connecting the A_hcp site to A_fcc site. In the H resurfacing from a subsurface octahedral vacancy to the fcc site () and from a subsurface tetrahedral vacancy to the hcp site (), the coordinate of the H atom (vertical to Ni(111) surface) is taken. In the H diffusion between the adjacent subsurface octahedral and tetrahedral vacancies () and between the adjacent interior octahedral vacancy and interior tetrahedral vacancy (), the reaction coordinates are along the directions that perpendicular to the planes of Ni₁–Ni₂–Ni₃ and Ni₄–Ni₅–Ni₆, respectively.

3.4.2. Free Energy and Prefactor

In order to investigate the quantized Ni lattice effect on the rates, we have recast the QI formula in (1) into Arrhenius form, which consists of the free energy and the prefactor (17) and (18).

In the free energy calculations, firstly, we calculate the free energy profile for each step of the H hopping paths, and , with the reaction coordinate defined in Section 3.4.1. Then, we connect these free energy profiles one by one, the final free energy profiles for the whole processes are displayed in Figures 11 and 12.

Figures 11 and 12 display the calculated free energy profiles for the surface-subsurface and interior processes, respectively, with the quantized lattice, the classical lattice and the rigid one at room temperature (300 K). The corresponding free energy barriers, prefactors, and rates are tabulated in Table 5.

Figures 11 and 12 clearly show that the classical lattices always reduce the free energy barriers compared with the rigid lattices, however, the differences between the free energy barriers with the quantized lattices and the ones with the classical lattices are very small. For different processes; however, the relaxation effect on the free energies is very different. For hydrogen diffusion on Ni(111), the classical lattice only slightly lowers the free energy barrier, while it decreases the barriers by more than one-half in subsurface and interior processes. More careful analysis from Figure 11 reveals that the two preferred Ni(111) surface binding sites, that is, hcp and fcc hollows have symmetric wells, manifesting that the motions of Ni atoms in the layer beneath the surface have little influence on the surface free energies despite the fact that the Ni atoms beneath the hcp and fcc hollows have different arrangements. It is also found that the site has a deeper well than the site and the well of the site nearly disappears as relaxed Ni atoms are considered. It manifests that the hydrogen at the site is very unstable and can easily move to the site or resurface to the hcp site. In the interior processes, although the hydrogen in the site is much less stable than in the site (see Figure 12), the well at the is obvious. One thus expects that a two-step reaction process can be used for the reaction from one to the other.

Next, we consider the prefactor. Table 5 shows that the prefactors of the classical lattices are much smaller than those of rigid ones except for the process of H diffusion on Ni(111), and the prefactors of quantized lattices are always smaller than that of classical ones, but their amplitudes have the same order. This may be explained by the fact that the quantum motions of the lattice atoms can induce the dissipative effect on the tunneling degrees of freedom [62, 63], because the prefactor essentially incorporates the dynamical effect. It is well known that pure dissipation in the overdamping regime always hinders the reaction rates for a given reaction barrier. The present results thus are consistent with above analysis.

The rate is determined by both the prefactor and the free energy barrier. In Table 5, the rate with the rigid lattice is smaller than the one on the classical lattice, which is because the prefactor changes a little and the free energy determines the rate. However, the rate with the quantized lattice is smaller than that of the classical lattice, which is due to the fact that quantized lattice has a smaller prefactor while the free energies are similar. Generally speaking, for the quantized lattices, the rates are lower by 20%–40% when compared to the ones on the classical lattices.

Another important feature of the free energy is its temperature dependence. Figure 13 displays the free energy profiles with both H- and Ni-treated quantum mechanically at several temperatures. Generally speaking, the free energies have a slight difference at 300 K and 400 K, whereas this difference becomes pronounced for 100 K and 200 K, and the barrier positions move to the directions of shallow well for asymmetric reactions. These properties can be explained by the hydrogen tunneling effect. At lower temperatures, the tunneling plays a more important role. Indeed, the barrier heights decrease with decreasing temperature except for among 200 K to 400 K. This special case may be due to the special structure of the lattice. The thermal average displacements of and (in Figure 9) vertical to the Ni(111) surface increase with increasing temperature. Thus, H goes through reaction bottleneck easier at a higher temperature, which makes the barrier decrease. Figure 13 also displays that the free energy barriers of change little in the temperature range of 200–400 K, while those of the become smaller and smaller with decreasing temperature. The corresponding free energy barriers are 0.47, 1.09, 1.47, and 1.66 kcal/mol at 100, 200, 300, and 400 K, respectively. It manifests that the diffusing H atom may not equilibrate in the interior tetrahedral vacancy at very low temperatures.

(a)

(b)

(c)

(d)

3.4.3. Surface-Subsurface Transport

In the resurfacing process, does not show an obvious barrier (see Figure 11), as the lattice atoms are treated quantum mechanically. This step thus can be thought as a barrierless process. Table 6 tabulates the rate constants for the other resurfacing and subsurface processes in the temperature range of 100 to 600 K. We also list the available CVT/SCT results for and its reverse reaction [60]. Again, the CVT/SCT rates are close to the QI rates except at 100 K. We think our much bigger rate constant at 100 K is due to the contribution of tunneling.

3.4.4. Interior Migration

In the interior of bulk Ni, the two most stable sites to cage H are symmetric octahedral vacancies (see Figure 12). H diffusion between them has been measured experimentally [64–67]. Several theoretical calculations have also been proposed to investigate this diffusion process. Wimmer et al. [68] calculate the diffusion coefficient via two-step reactions by using a transition state theory together with accurate ab initio energies, while Wonchoba and Truhlar [60] consider the kinetic step as a direct process with a double maximum barrier and calculate the diffusion coefficient using the CVT/SCT. In the present QI calculations, Figure 13 has explicitly shown that the free energies have a well at the tetrahedral site from 200 K to 400 K. It is thus reasonable to assume that the diffusing H atom temporarily equilibrates in the tetrahedral site before jumping forward or backward to a neighboring octahedral site. The free energy well, however, becomes very shallow at 100 K. In this case, the direct reaction from the octahedral site to the other one may be acceptable. Here, we only calculate the diffusion coefficients via the two kinetic steps at 200–400 K.

The temperature dependence of the diffusion coefficients are commonly fitted to the Arrhenius equation where is the gas constant and and are the pre-exponential factor and the activation energy, respectively. The QI calculations predict and . In the calculations, the rates of are used to obtain the diffusion coefficients for , because this process is much slower than that of and it determines the total reaction rates.

Table 7 tabulates the pre-exponential factors and activation energies coming from available experiments and theories. It is found that both and from the QI calculations are close to Ebisuzaki’s experimental data [65]. Further tracking down the comparisons with experiments is nontrivial, because the accuracy of the diffusion coefficients is much dependent of the potential energy surface. However, we may make a quantitative comparison for the QI and CVT/SCT results. With use of the diffusion coefficients from 295 K to 300 K obtained by the CVT/SCT, Wonchoba and Truhlar [60] predict 11.1 kcal/mol for , and for , respectively. These values are observably different from the QI calculations. and are 0.8 kcal/mol larger and 3 times smaller than those from the QI calculations, respectively. The origin of these differences can be explained by that Wonchoba and Truhlar have treated the processes and as a single kinetic step rather than as two kinetic steps, which results in a much longer tunneling path than that of the two steps. Compared with the results reported by Wimmer et al. [68], the activation energy is about 0.6 kcal/mol larger than the present one, and the pre-exponential factor is 10 times larger, which is also larger than all available experimental data. These differences may come from both the different potential energy surfaces and rate methods.

4. Conclusion

We have presented the basic principle of the quantum instanton (QI) approximation and its applications to chemical reactions from gas phase to surface systems. The applications demonstrate that the QI method makes it possible to treat hundreds of atoms because of the well-established techniques of imaginary time path integrals. For instance, more than 200 atoms have been incorporated in the present study of the H diffusion processes.

The QI approximation is a kind of ‘‘quantum transition state theory’’ in that there is no account of “recrossing” dynamics in the description. The recrossing effects on the quantum instanton rate constants have been quantified for several collinear reaction by Ceotto and Miller [25], and it is found especially evident for the collinear heavy-light-heavy reactions. Fortunately, the recrossing effects become generally less important in higher dimensions [69]. Therefore, the QI approach may become a suitable tool for the calculation of chemical reaction rate constants of complex systems.

Compared to conventional TST theories, the QI approximation involves two dividing surfaces, which are quantum analogs of the two turning point surfaces of the imaginary time trajectory in the semiclassical instanton theory. At high-temperature limit, these dividing surfaces coalesce into the one, the same as the dividing surface from Wigner’s variational principle. In this case, the QI approximation becomes exactly the same as the classical TST. As the tunneling corrections are incorporated, the CVT rates are consistent with the QI rates except at deep tunneling regime, where QI rates are generally greater than the rates from CVT with various tunneling corrections, since the CVT method uses an optimized tunneling path, while the QI method considers all tunneling paths and automatically gives each path its natural weight by the quantum Boltzmann factor. For the reaction of , accurate quantum dynamics rate constants are obtained with the MCTDH method. Compared to the MCTDH ones, the QI rate constants are larger by factors of about 2 to 3 over the temperature range 300–400 K. This difference may partly be due to the recrossing effect which is not considered in QI theory and partly arise from the use of the J-shifting approximation and the neglect of the vibrational angular momenta Hamiltonian in MCTDH method.

Acknowledgments

This work was financially supported by the National Science Foundation of China (Grant nos. 20833004 and 21073146), National Key Basic Research Foundation Program of China (2007CB815204), and the Scientific Research Foundation of Northwest A&F University (Grant no. Z109021103).

References

H. Eyring, “The activated complex in chemical reactions,” The Journal of Chemical Physics, vol. 3, no. 2, pp. 63–71, 1935.
View at: Google Scholar
H. Eyring, “The theory of absolute reaction rates,” Transactions of the Faraday Society, vol. 34, pp. 41–48, 1938.
View at: Publisher Site | Google Scholar
E. Wigner, “The transition state method,” Transactions of the Faraday Society, vol. 34, pp. 29–41, 1938.
View at: Publisher Site | Google Scholar
G. A. Voth, D. Chandler, and W. H. Miller, “Time correlation function and path integral analysis of quantum rate constants,” Journal of Physical Chemistry, vol. 93, no. 19, pp. 7009–7015, 1989.
View at: Google Scholar
G. A. Voth, D. Chandler, and W. H. Miller, “Rigorous formulation of quantum transition state theory and its dynamical corrections,” The Journal of Chemical Physics, vol. 91, no. 12, pp. 7749–7760, 1989.
View at: Google Scholar
S. Jang and G. A. Voth, “A relationship between centroid dynamics and path integral quantum transition state theory,” Journal of Chemical Physics, vol. 112, no. 20, pp. 8747–8757, 2000.
View at: Google Scholar
W. H. Miller, “Generalization of the linearized approximation to the semiclassical initial value representation for reactive flux correlation functions,” Journal of Physical Chemistry A, vol. 103, no. 47, pp. 9384–9387, 1999.
View at: Google Scholar
E. Pollak and J. L. Liao, “A new quantum transition state theory,” Journal of Chemical Physics, vol. 108, no. 7, pp. 2733–2743, 1998.
View at: Google Scholar
J. Shao, J. L. Liao, and E. Pollak, “Quantum transition state theory: perturbation expansion,” Journal of Chemical Physics, vol. 108, no. 23, pp. 9711–9725, 1998.
View at: Google Scholar
K. Yamashita and W. H. Miller, “‘Direct’ calculation of quantum mechanical rate constants via path integral methods: application to the reaction path Hamiltonian, with numerical test for the $H + H_{2}$ reaction in 3D,” The Journal of Chemical Physics, vol. 82, no. 12, pp. 5475–5484, 1984.
View at: Google Scholar
J. W. Tromp and W. H. Miller, “New approach to quantum mechanical transition-state theory,” Journal of Physical Chemistry, vol. 90, no. 16, pp. 3482–3485, 1986.
View at: Google Scholar
N. F. Hansen and H. C. Andersen, “A new formulation of quantum transition state theory for adiabatic rate constants,” The Journal of Chemical Physics, vol. 101, no. 7, pp. 6032–6037, 1994.
View at: Google Scholar
G. Krilov, E. Sim, and B. J. Berne, “Quantum time correlation functions from complex time Monte Carlo simulations: a maximum entropy approach,” Journal of Chemical Physics, vol. 114, no. 3, pp. 1075–1088, 2001.
View at: Publisher Site | Google Scholar
E. Sim, G. Krilov, and B. J. Berne, “Quantum rate constants from short-time dynamics: an analytic continuation approach,” Journal of Physical Chemistry A, vol. 105, no. 12, pp. 2824–2833, 2001.
View at: Google Scholar
W. H. Miller, Y. Zhao, M. Ceotto, and S. Yang, “Quantum instanton approximation for thermal rate constants of chemical reactions,” Journal of Chemical Physics, vol. 119, no. 3, pp. 1329–1342, 2003.
View at: Publisher Site | Google Scholar
T. Yamamoto and W. H. Miller, “On the efficient path integral evaluation of thermal rate constants within the quantum instanton approximation,” Journal of Chemical Physics, vol. 120, no. 7, pp. 3086–3099, 2004.
View at: Publisher Site | Google Scholar
Y. Zhao, T. Yamamoto, and W. H. Miller, “Path integral calculation of thermal rate constants within the quantum instanton approximation: application to the $H + {CH}_{4} \to H_{2} + {CH}_{3}$ hydrogen abstraction reaction in full Cartesian space,” Journal of Chemical Physics, vol. 120, no. 7, pp. 3100–3107, 2004.
View at: Publisher Site | Google Scholar
W. H. Miller, “Semiclassical limit of quantum mechanical transition state theory for nonseparable systems,” The Journal of Chemical Physics, vol. 62, no. 5, pp. 1899–1906, 1975.
View at: Google Scholar
S. Coleman, “Fate of the false vacuum: semiclassical theory,” Physical Review D, vol. 15, no. 10, pp. 2929–2936, 1977.
View at: Publisher Site | Google Scholar
A. J. Leggett, S. Chakravarty, A. T. Dorsey, M. P. A. Fisher, A. Garg, and W. Zwerger, “Dynamics of the dissipative two-state system,” Reviews of Modern Physics, vol. 59, no. 1, pp. 1–85, 1987.
View at: Publisher Site | Google Scholar
R. A. Marcus, “On the analytical mechanics of chemical reactions. quantum mechanics of linear collisions,” The Journal of Chemical Physics, vol. 45, no. 12, pp. 4493–4499, 1966.
View at: Google Scholar
E. A. Mccullough and R. E. Wyatt, “Dynamics of the collinear $H + H_{2}$ reaction. I. Probability density and flux,” The Journal of Chemical Physics, vol. 54, no. 8, pp. 3578–3591, 1971.
View at: Google Scholar
D. G. Truhlar and A. Koppermann, “Exact and approximate quantum mechanical reaction probabilities and rate constants for the collinear $H + H_{2}$ reaction,” The Journal of Chemical Physics, vol. 56, no. 5, pp. 2232–2252, 1972.
View at: Google Scholar
T. V. George and W. H. Miller, “Classical S-matrix theory of reactive tunneling: linear $H + H_{2}$ collisions,” The Journal of Chemical Physics, vol. 57, no. 6, pp. 2458–2467, 1972.
View at: Google Scholar
M. Ceotto and W. H. Miller, “Test of the quantum instanton approximation for thermal rate constants for some collinear reactions,” Journal of Chemical Physics, vol. 120, no. 14, pp. 6356–6362, 2004.
View at: Publisher Site | Google Scholar
C. Venkataraman and W. H. Miller, “The Quantum Instanton (QI) model for chemical reaction rates: the “simplest” QI with one dividing surface,” Journal of Physical Chemistry A, vol. 108, no. 15, pp. 3035–3039, 2004.
View at: Publisher Site | Google Scholar
Y. Li and W. H. Miller, “Different time slices for different degrees of freedom in Feynman path integration,” Molecular Physics, vol. 103, no. 2-3, pp. 203–208, 2005.
View at: Publisher Site | Google Scholar
T. Yamamoto and W. H. Miller, “Path integral evaluation of the quantum instanton rate constant for proton transfer in a polar solvent,” Journal of Chemical Physics, vol. 122, no. 4, Article ID 044106, pp. 1–13, 2005.
View at: Publisher Site | Google Scholar
M. Ceotto, S. Yang, and W. H. Miller, “Quantum reaction rate from higher derivatives of the thermal flux-flux autocorrelation function at time zero,” Journal of Chemical Physics, vol. 122, no. 4, Article ID 044109, pp. 1–7, 2005.
View at: Publisher Site | Google Scholar
J. Vaníček, W. H. Miller, J. F. Castillo, and F. J. Aoiz, “Quantum-instanton evaluation of the kinetic isotope effects,” Journal of Chemical Physics, vol. 123, no. 5, Article ID 054108, pp. 1–14, 2005.
View at: Publisher Site | Google Scholar
Y. Li and W. H. Miller, “Using a family of dividing surfaces normal to the minimum energy path for quantum instanton rate constants,” Journal of Chemical Physics, vol. 125, no. 6, Article ID 064104, 2006.
View at: Publisher Site | Google Scholar
S. Yang, T. Yamamoto, and W. H. Miller, “Path-integral virial estimator for reaction-rate calculation based on the quantum instanton approximation,” Journal of Chemical Physics, vol. 124, no. 8, Article ID 084102, 2006.
View at: Publisher Site | Google Scholar
J. Vaníček and W. H. Miller, “Efficient estimators for quantum instanton evaluation of the kinetic isotope effects: application to the intramolecular hydrogen transfer in pentadiene,” Journal of Chemical Physics, vol. 127, no. 11, Article ID 114309, 2007.
View at: Publisher Site | Google Scholar
M. Buchowiecki and J. Vaníček, “Direct evaluation of the temperature dependence of the rate constant based on the quantum instanton approximation,” Journal of Chemical Physics, vol. 132, no. 19, Article ID 194106, 2010.
View at: Publisher Site | Google Scholar
T. Zimmermann and J. Vaníček, “Three applications of path integrals: equilibrium and kinetic isotope effects, and the temperature dependence of the rate constant of the [1,5] sigmatropic hydrogen shift in (Z)-1,3-pentadiene,” Journal of Molecular Modeling, vol. 16, no. 11, pp. 1779–1787, 2010.
View at: Publisher Site | Google Scholar
K. F. Wong, J. L. Sonnenberg, F. Paesani et al., “Proton transfer studied using a combined ab initio reactive potential energy surface with quantum path integral methodology,” Journal of Chemical Theory and Computation, vol. 6, no. 9, pp. 2566–2580, 2010.
View at: Publisher Site | Google Scholar
W. Wang, S. Feng, and Y. Zhao, “Quantum instanton evaluation of the thermal rate constants and kinetic isotope effects for ${SiH}_{4} +H \to {SiH}_{3} + H_{2}$ reaction in full Cartesian space,” Journal of Chemical Physics, vol. 126, no. 11, Article ID 114307, 2007.
View at: Publisher Site | Google Scholar
W. Wang and Y. Zhao, “Path integral evaluation of H diffusion on Ni(100) surface based on the quantum instanton approximation,” Journal of Chemical Physics, vol. 130, no. 11, Article ID 114708, pp. 1–10, 2009.
View at: Publisher Site | Google Scholar
W. Wang and Y. Zhao, “Quantum instanton evaluations of surface diffusion, interior migration, and surface-subsurface transport for H/Ni,” Journal of Chemical Physics, vol. 132, no. 6, Article ID 064502, 2010.
View at: Publisher Site | Google Scholar
B. J. Berne and D. Thirumalai, “On the simulation of quantum systems: path integral methods,” Annual Review of Physical Chemistry, vol. 37, pp. 401–424, 1986.
View at: Google Scholar
R. P. Feynman and A. R. Hibbs, Quantum Mechanics and Path Integrals, McGraw-Hill, New York, NY, USA, 1965.
R. P. Feynman, Statistical Mechanics, Addison-Wesley, New York, NY, USA, 1972.
L. S. Schulman, Techniques and Applications of Path Integrals, John Wiley & Sons, New York, NY, USA, 1986.
M. Mezei, “Adaptive umbrella sampling: self-consistent determination of the non-Boltzmann bias,” Journal of Computational Physics, vol. 68, no. 1, pp. 237–248, 1987.
View at: Google Scholar
C. Predescu and W. H. Miller, “Optimal choice of dividing surface for the computation of quantum reaction rates,” Journal of Physical Chemistry B, vol. 109, no. 14, pp. 6491–6499, 2005.
View at: Publisher Site | Google Scholar
J. Espinosa-García, “New analytical potential energy surface for the CH₄+H hydrogen abstraction reaction: thermal rate constants and kinetic isotope effects,” Journal of Chemical Physics, vol. 116, no. 24, pp. 10664–10673, 2002.
View at: Publisher Site | Google Scholar
D. L. Baulch, C. J. Cobos, R. A. Cox et al., “Evaluated kinetic data for combustion modelling,” Journal of Physical and Chemical Reference Data, vol. 21, no. 3, pp. 411–734, 1992.
View at: Google Scholar
J. W. Sutherland, M. C. Su, and J. V. Michael, “Rate constants for H + CH₄, CH₃ + H₂, and CH₄ dissociation at high temperature,” International Journal of Chemical Kinetics, vol. 33, no. 11, pp. 669–684, 2001.
View at: Publisher Site | Google Scholar
S. Andersson, G. Nyman, A. Arnaldsson, U. Manthe, and H. Jonsson, “Comparison of quantum dynamics and quantum transition state theory estimates of the H + CH₄ reaction rate,” Journal of Physical Chemistry A, vol. 113, no. 16, pp. 4468–4478, 2009.
View at: Publisher Site | Google Scholar
J. Espinosa-García, J. Sansón, and J. C. Corchado, “The ${SiH}_{4} + H \to {SiH}_{3} + H_{2}$ reaction: potential energy surface, rate constants, and kinetic isotope effects,” Journal of Chemical Physics, vol. 109, no. 2, pp. 466–473, 1998.
View at: Publisher Site | Google Scholar
K. D. Dobbs and D. A. Dixon, “Ab initio prediction of the activation energies for the abstraction and exchange reactions of H with CH₄ and SiH₄,” Journal of Physical Chemistry, vol. 98, no. 20, pp. 5290–5297, 1994.
View at: Google Scholar
A. Goumri, W. J. Yuan, L. Ding, Y. Shi, and P. Marshall, “Experimental and theoretical studies of the reaction of atomic hydrogen with silane,” Chemical Physics, vol. 177, no. 1, pp. 233–241, 1993.
View at: Google Scholar
T. N. Truong and D. G. Truhlar, “Effect of steps and surface coverage on rates and kinetic isotope effects for reactions catalyzed by metallic surfaces: chemisorption of hydrogen on Ni,” Journal of Physical Chemistry, vol. 94, no. 21, pp. 8262–8279, 1990.
View at: Google Scholar
J. V. Barth, “Transport of adsorbates at metal surfaces: from thermal migration to hot precursors,” Surface Science Reports, vol. 40, no. 3, pp. 75–149, 2000.
View at: Publisher Site | Google Scholar
T. N. Truong and D. G. Truhlar, “The effects of steps, coupling to substrate vibrations, and surface coverage on surface diffusion rates and kinetic isotope effects: hydrogen diffusion on Ni,” The Journal of Chemical Physics, vol. 93, no. 3, pp. 2125–2138, 1990.
View at: Google Scholar
S. M. George, A. M. DeSantolo, and R. B. Hall, “Surface diffusion of hydrogen on Ni(100) studied using laser-induced thermal desorption,” Surface Science, vol. 159, no. 1, pp. L425–L432, 1985.
View at: Google Scholar
D. R. Mullins, B. Roop, S. A. Costello, and J. M. White, “Isotope effects in surface diffusion: hydrogen and deuterium on Ni(100),” Surface Science, vol. 186, no. 1-2, pp. 67–74, 1987.
View at: Google Scholar
A. Lee, X. D. Zhu, L. Deng, and U. Linke, “Observation of a transition from over-barrier hopping to activated tunneling diffusion: H and D on Ni(100),” Physical Review B, vol. 46, no. 23, pp. 15472–15476, 1992.
View at: Publisher Site | Google Scholar
T. S. Lin and R. Gomer, “Diffusion of ¹H and ²H on the Ni(111) and (100) planes,” Surface Science, vol. 255, no. 1-2, pp. 41–60, 1991.
View at: Google Scholar
S. E. Wonchoba and D. G. Truhlar, “General potential-energy function for H/Ni and dynamics calculations of surface diffusion, bulk diffusion, subsurface-to-surface transport, and absorption,” Physical Review B, vol. 53, no. 16, pp. 11222–11241, 1996.
View at: Google Scholar
T. N. Truong, D. G. Truhlar, and B. C. Garrett, “Embedded diatomics-in-molecules: a method to include delocalized electronic interactions in the treatment of covalent chemical reactions at metal surfaces,” Journal of Physical Chemistry, vol. 93, no. 25, pp. 8227–8239, 1989.
View at: Google Scholar
R. Baer and R. Kosloff, “Quantum dissipative dynamics of adsorbates near metal surfaces: a surrogate Hamiltonian theory applied to hydrogen on nickel,” Journal of Chemical Physics, vol. 106, no. 21, pp. 8862–8875, 1997.
View at: Google Scholar
R. Baer, Y. Zeiri, and R. Kosloff, “Hydrogen transport in nickel (111),” Physical Review B, vol. 55, no. 16, pp. 10952–10974, 1997.
View at: Google Scholar
K. Yamakawa, “Diffusion of deuterium and isotope effect in nickel,” Journal of the Physical Society of Japan, vol. 47, no. 1, pp. 114–121, 1979.
View at: Google Scholar
Y. Ebisuzaki, W. J. Kass, and M. O'Keeffe, “Diffusion and solubility of hydrogen in single crystals of nickel and nickel-vanadium alloy,” The Journal of Chemical Physics, vol. 46, no. 4, pp. 1378–1381, 1967.
View at: Google Scholar
W. Eichenauer, W. Loser, and H. Witte, “Solubility and rate of diffusion of hydrogen and deuterium in nickel and copper single crystals,” Z Metallk, vol. 56, pp. 287–293, 1965.
View at: Google Scholar
L. Katz, M. Guinan, and R. J. Borg, “Diffusion of H₂, D₂, and T₂ in single-crystal Ni and Cu,” Physical Review B, vol. 4, no. 2, pp. 330–341, 1971.
View at: Publisher Site | Google Scholar
E. Wimmer, W. Wolf, J. Sticht et al., “Temperature-dependent diffusion coefficients from ab initio computations: hydrogen, deuterium, and tritium in nickel,” Physical Review B, vol. 77, no. 13, Article ID 134305, 2008.
View at: Publisher Site | Google Scholar
S. Chapman, S. M. Hornstein, and W. H. Miller, “Accuracy of transition state theory for the threshold of chemical reactions with activation energy. Collinear and three-dimensional $H + H_{2}$ ,” Journal of the American Chemical Society, vol. 97, no. 4, pp. 892–894, 1975.
View at: Google Scholar

Copyright

Copyright © 2012 Yi Zhao and Wenji Wang. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

PDF Download Citation

Download other formats

Order printed copies

Views

1340

Downloads

785

Citations