Quantum Instanton Evaluations of the Thermal Rate Constants for Complex Systems

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 H + CH4 → H2 + CH3, (2) the reaction of H + SiH4 → H2 + SiH3, (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.


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-toohigh 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][5][6][7][8][9][10][11][12][13][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][16][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][22][23][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 Advances in Physical Chemistry 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][26][27][28][29][30][31][32][33][34][35][36][37][38][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 2 → HD + H [16]. Soon, the techniques are further extended to the reaction of H+CH 4 → H 2 + CH 3 [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 H + CH 4 → H 2 + CH 3 [17] and H + SiH 4 → H 2 + SiH 3 [37], H diffusion on Ni(100) surface [38], and surface-subsurface transport and interior migration for H/Ni [39].

Method
In this section, we summarize the rate formula for the QI evaluation.The detailed derivation can be found in [15][16][17].The QI model proposes the following thermal rate constant: Here, Q r is the reactant partition function per unit volume.C f f (0) is zero time value of the flux-flux correlation function where β is the inverse temperature 1/(k B T), H is the Hamiltonian operator of the reaction system, and F a and F b are the flux operators given by with γ = a, b.In (3), h is the step-side function, r represents the Cartesian coordinates of the reaction system, and s a (r) and s b (r) define two separate dividing surfaces via the equations s a (r) = 0 and s b (r) = 0, both s a (r) and s b (r) being positive (negative) on the product (reactant) side of the dividing surfaces.ΔH(β) 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 2)/β, which has very little effect in the low temperature regime.C dd (0) and Cdd (0) are zero time value and its second time derivative, respectively, of the "delta-delta" correlation function where the generalized delta-function operator is Here, N is the total number of atoms, ∇ i = ∂/∂r i , r i denotes the Cartesian coordinates of the ith atom, and m i is its atomic mass.
The dividing surfaces are determined by the stationary condition where {c k } 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.
Path integral expressions for C f f (0) and Cdd (0) are somewhat more complicated but can be obtained in a straightforward manner.The appropriate expressions are × Δ s b r (P/2) exp −βΦ r (s) f v r (s) , (10) with with with f being the total number of degrees of freedom (i.e., f = 3N).
In realistic calculations, we rewrite (1) as the product of several ratios The terms of C f f (0)/C dd (0) and ΔH are directly calculated as a constrained average over the same ensemble of paths [16,17] with The evaluation of C dd (0)/Q r , however, meets a challenge because the C dd (0) is the quantity associated with the transition state, while Q r 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 Advances in Physical Chemistry Here, the free energy ΔF is defined by This two-dimensional free energy is related to Q r (the reactant partition function) and C dd (0; ξ a , ξ b ), and it is corresponding to the quality of probability density at (ξ a , ξ b ).C dd (0; ξ a , ξ b ) 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 vibrationalrotational coupling and anharmonicity are automatically involved.

The Reaction of H
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 P is chosen to be 20 and 100 at temperatures T = 1000 K and 200 K, respectively, while 3 × 10 7 Monte Carlo cycles are run to achieve <10% statistical convergence.

Reaction Coordinate.
For this reaction, we define a generalized reaction coordinate s(r; ξ), where ξ is an adjustable parameter that shifts the location of the dividing surface (defined by s(r; ξ) = 0).s(r; ξ) is defined by a linear interpolation between two constituent reaction coordinates s 0 (r) and s 1 (r) through the parameter ξ, s 1 (r) 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 s x (r)(x = α, β, γ, δ) being the reaction coordinate that describes the abstraction process of one of the methane hydrogens H x by the incident one H where r(X − Y ) denotes the interatomic distance between atoms X and Y and r † (X − Y ) is the value at the transition state geometry.s 0 (r), 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.R ∞ is an adjustable parameter which is chosen to be 9 Å in order to guarantee that the interaction potential energy between H and CH 4 is negligible.Now, the term C dd (0) (5) becomes a function of two parameters, ξ a and ξ b , as follows: and thus, one seeks a stationary point of C dd (0; ξ a , ξ b ) in the two-dimensional (ξ a , ξ b ) space to obtain the corresponding "optimal" values.

Free Energy Surface.
The quantity C dd (0; ξ a , ξ b ) varies exponentially as a function of (ξ a , ξ b ), it is convenient to define a quantum "free energy surface" as follows: and locate the saddle point of F(ξ a , ξ b ) by visual inspection.
Figure 1 shows the free energy surfaces for the H + CH 4 reaction.It exhibits a barrier-like profile along the direction ξ = (ξ a + ξ b )/2, while it grows approximately quadratically with the increasing of the absolute value of Δξ = ξ a − ξ b .From this figure, it is seen that at a higher temperature T = 1000 K (Figure 1(a)), there appears only a single saddle point at (ξ, Δξ) = (1.02,0.0), while at a low temperature T = 200 K (Figure 1(b)), the saddle point bifurcates into a distinct pair at (ξ, Δξ) = (1.1,±0.35), which indicates the existence of nonnegligible tunneling effects in the rate constant [15,16].

Rate Constants.
Having obtained the "optimal" values of the (ξ a , ξ b ) 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 k QI 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 k QI agrees with k CVT/μOMT 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 T = 500 and 300 K, resp.).It should be noted that the differences between k QI and k CVT/μOMT 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  ) is the saddle point of the quantum free energy surface.k QI is the quantum instanton rate constant.k CVT/μOMT is the rate constant of the canonical variational theory with microcanonical optimized multidimensional tunneling (CVT/μOMT) [46].k a exp t and k b exp t are the experimental Arrhenius fits, k(T) = 2.18 × 10 −20 T 3.0 exp(−4045/T) [47] and k(T) = 6.78 × 10 −21 T 3.156 exp(−4406/T) [48], respectively.Unit: cm 3 s −1 for rates.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].

The Reaction of H
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 H + SiH 4 has the same form as that of H + CH 4 (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 of time slices in the path integral is set to be 1.The number of Monte Carlo is about (6 − 10) × 10 6 for computing a single ensemble average, and it converges the relevant quantities within 10% statistical errors.

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 [47,48], respectively.The dot-dashed line is the result of the multiconfigurational time-dependent Hartree approach (MCTDH) [49].The open squares are the results of the canonical variational theory with microcanonical optimized multidimensional tunneling (CVT/μOMT) [46].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.

T (K)
Present 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 groundstate (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.

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 k QI (R1)/k QI (R2) and k QI (R1)/k QI (R3) 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 k(R1)/k(R2) 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 k(R1)/k(R3) 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.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].

H Diffusion on Ni(
In the path integral calculations, the numbers of time slices, P and P bath , for the degrees of freedom of H and quantized Ni atoms, respectively, are set to (P, P bath ) = (24)(25)(26)(27)(28)(29)(30)(31)(32)(33)(34)(35)(36)(37)(38)(39)(40)(6)(7)(8) in the temperature range of 40-300 K.The number of Monte Carlo cycle employed is (1-10) ×10 6 , which converges most of the relevant quantities within 10% statistical error (some of the statistical errors are within 20% below 100 K). 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 x coordinate of H atom (along the reaction path) to be reaction coordinate.

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 (x 0 and x P/2 ) 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.

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.constructed by Wonchoba and Truhlar [60], using the embedded diatomics-in-molecules (EDIMs) [61] potential energy function, is used in this QI calculations.

Surface-Subsurface Transport and Interior Migration
In path integral calculations, the numbers of time slices, P and P bath for the degrees of freedom of the H and quantized Ni atoms, respectively, are set to (P, P bath ) = (30−40, 6−8) in the temperature range of 100-400 K.The number of Monte Carlo is about (2 − 6) × 10 6 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).

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 s (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 (A hcp → A fcc ), the x 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 (Oc sub → A fcc ) and from a subsurface tetrahedral vacancy to the hcp site (Te sub → A hcp ), the z coordinate of the H atom (vertical to Ni(111) surface) is taken.In the H diffusion between the adjacent subsurface octahedral and tetrahedral vacancies (Oc sub → Te sub ) and between the adjacent interior octahedral vacancy and interior tetrahedral vacancy (Oc inte → Te inte ), the reaction coordinates are along the directions that perpendicular to the planes of Ni 1 -Ni 2 -Ni 3 and Ni 4 -Ni 5 -Ni , respectively.

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, A hcp → A fcc → Oc sub → Te sub → A hcp and Oc inte → Te inte → Oc inte , 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;

S ( Å)
Figure 12: Free energy profiles with respect to the processes in Figure 10 at 300 K.The green, the red, and the blue lines correspond to the rigid, classical, and quantized lattices, respectively.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 Oc sub site has a deeper well than the Te sub site and the well of the Te sub site nearly disappears as relaxed Ni atoms are considered.It manifests that the hydrogen at the Te sub site is very unstable and can easily move to the Oc sub site or resurface to the hcp site.In the interior processes, although the hydrogen in the Te inte site is much less stable than in the Oc inte site (see Figure 12), the well at the Te inte is obvious.One thus expects that a two-step reaction process can be used for the reaction from one Oc inte 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 Oc sub → Te sub among 200 K to 400 K.This special case may be due to the special structure of the lattice.The thermal average displacements of Ni 1 and Ni 3 (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 Oc inte → Te inte change little in the temperature range of 200-400 K, while those of the Te inte → Oc inte 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.

Surface-Subsurface Transport.
In the resurfacing process, Te sub → A hcp 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 Oc sub → A fcc 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.

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][65][66][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 Oc inte → Te inte → Oc inte 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 R is the gas constant and D 0 and E a are the preexponential factor and the activation energy, respectively.The QI calculations predict D 0 = 3.93 × 10 −3 cm 2 s −1 and E a = 10.26 kcal/mol.In the calculations, the rates of Oc inte → Te inte are used to obtain the diffusion coefficients for Oc inte → Oc inte , because this process is much slower than that of Te inte → Oc inte 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 D 0 and E a 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 E a , and 1.3 × 10 −3 cm 2 s −1 for D 0 , respectively.These values are observably different from the QI calculations.E a and D 0 are 0.8 kcal/mol larger and 3 times smaller than those from the QI calculations, respectively.The origin of these differences can be   Figure 13: The temperature dependence of the free energies.The blue, the black, the red, and the green lines correspond to free energy profiles at 100 K, 200 K, 300 K, and 400 K, respectively.explained by that Wonchoba and Truhlar have treated the processes Oc inte → Te inte and Te inte → Oc inte 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.

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 Advances in Physical Chemistry 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 H + CH 4 → H 2 + CH 3 , 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.

Figure 1 :
Figure 1: Local topography of the quantum free energy surface defined by (24) near the top of the barrier.(a) T = 1000 K; (b) T = 200 K.The cross symbols show the location of the saddle points.The values of ξ a and ξ b at the saddle points are used as input for computing the quantum instanton rate.

Figure 3 :
Figure 3: The temperature dependence of the free energy.The solid and dotted lines correspond to the quantum and classical calculations, respectively.

Figure 4 :Figure 5 :
Figure 4: The temperature dependence of the prefactor.The solid and dotted lines correspond to the quantum and classical calculations, respectively.

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 :
Figure 6: The platform and profile chart of H diffusion on Ni(100) lattice.The gray circles represent H atom on different surface sites, and the orange circles represent the fixed Ni atoms.Blue circles are Ni atoms treated classically, while the red ones are treated quantum mechanically.

Figure 7 :
Figure 7: Distribution of H and Ni quantum paths on the surface with two reaction coordinate beads (x 0 and x P/2 ) of H path fixed at the two dividing surfaces.Labels of <010> and <001> denote the crystal directions of nickel.The probabilities are normalized for both H and Ni.

Table 4 :Figure 8 :
Figure 8: Arrhenius plots of the diffusion coefficients in the range 40-300 K. Solid line: the QI results for the quantized lattice; dotted line: the QI results for the rigid lattice; dashed line: the CVT/SCSAG results for the rigid lattice [55].Pluses, circles, squares, and crosses are experimental results of Lee et al. [58], George et al. [56], Mullins et al. [57], and Lin and Gomer [59], respectively.

Figure 9 :
Figure 9: A lattice model with a few Ni atoms for the circle reaction processes of hydrogen.A hcp , A fcc , Oc sub , and Te sub are the abbreviations for a hcp hollow site, a fcc hollow site, a subsurface octahedral vacancy, and a subsurface tetrahedral vacancy, respectively.The black lines stand for the directions of the reaction paths.The atoms of Ni 1 , Ni 2 , and Ni 3 colored in orange are specially used to determine the reaction coordinate of the process Oc sub → Te sub .

Figure 10 :
Figure 10: A lattice model with a few Ni atoms for the process of H diffusion in interior of bulk Ni.Oc inte and Te inte are the abbreviations for an interior octahedral vacancy and an interior tetrahedral vacancy, respectively.The black lines stand for the general directions of the reaction paths.The atoms of Ni 4 , Ni 5 , and Ni 6 colored in orange are specially used to determine the reaction coordinate of the process Oc inte → Te inte .

Figure 11 :
Figure11: Free energy profiles with respect to the processes in Figure9at 300 K.The green, the red, and the blue lines correspond to the rigid, classical, and quantized lattices, respectively.

Table 1 :
Theoretical and experimental rate constants obtained for the H + CH 4 → H 2 + CH 3 reaction.(ξ a , ξ b plots the quantum and classical free energies as a function of temperature Figure 2: Arrhenius plots of the thermal rate constant for the H + CH 4 → H 2 + CH 3 reaction.The solid line is the quantum instanton rate.Dotted and dashed lines are the Arrhenius fits of the experimental data from
a: The results for a rigid lattice.b: The results for a classical lattice.c: The results for a quantized lattice.