The Classical Wigner Method with an Effective Quantum Force

ARTICLE pubs.acs.org/JPCA

The Classical Wigner Method with an Effective Quantum Force: Application to the Collinear H þ H2 Reaction Huaqing Li, Jens Aage Poulsen, and Gunnar Nyman* Department of Chemistry, University of Gothenburg, SE-412-96, Gothenburg, Sweden ABSTRACT: To improve the classical Wigner (CW) model, we recently proposed the classical Wigner model with an effective quantum force (CWEQF).1 The results of the CWEQF model are more accurate than those of the CW model. Still the simplicity of the CW model is retained. The quantum force was obtained by defining a characteristic distance η0 between two Feynman paths that enter the expression for the fluxflux correlation function. η0 was considered independent of the position along the reaction path. The CWEQF leads to a lowering of the effective potential barrier. Here we develop the method to use position dependent η0 values. The method is also generalized to two dimensions. Applications are carried out on one-dimensional model problems and the two-dimensional H þ H2 collinear reaction.

I. INTRODUCTION To obtain rate constants is one of the fundamental tasks in chemical reaction dynamics. Methods developed for computing reaction rates may rely on quantum mechanics, classical mechanics or a mix of these. Full-dimensional quantum methods are exact but can so far in general only be applied to reactions with a small number of atoms. Classical molecular dynamics methods can be applied to large systems but are unreliable if quantum effects are significant. Ways to bridge the capabilities of the quantum and classical methods are thus desirable. Wyatt et al.25 have implemented Bohmian mechanics,6,7 which bridges classical and quantum mechanics through a quantum potential in configuration space. The quantum potential is constructed from the propagated wave function, which is a solution to the Schr€odinger equation. However, this method experiences numerical difficulties due to the nodes of the wave function which causes the quantum potential to diverge. As an alternative, dynamics is also studied in phase space, particularly Wigner space.8 The Wigner function is a quasi-probability function in phase space which allows the study of both pure and mixed quantum states. Martens et al.9 proposed the so-called entangled trajectory molecular dynamics (ETMD) scheme. For each trajectory in phase space, the propagation is here entangled with other trajectories. The ETMD method captures the essential quantum effects due to the entanglement in the dynamics but is numerically heavy and thus difficult to implement in multidimensional cases. Also, it is, strictly speaking, only applicable to positive Wigner functions. Liu and Miller10 proposed ways to simplify the ETMD method by using an effective potential that depends on density, position, and momentum. To simplify the calculation, they used the equilibrium distribution approximation (EDA) and thermal Gaussian approximation (TGA).11,12 r 2011 American Chemical Society

Another option that introduces quantum effects into classical simulations and also reduces the numerical cost is the classical Wigner (CW) model.1317 The time propagation of each phasespace point is in the CW model performed using classical dynamics, which is independent of the other trajectories so that dynamical quantum effects are not captured. The CW model is exact for harmonic potentials and is otherwise most accurate for short time problems and at high temperature.16,18 When dynamical quantum effects are significant, which can be expected at low temperature, classical trajectory-based methods are generally not accurate. To improve the CW model without increasing the numerical cost, we proposed the classical Wigner model with an effective quantum force (CWEQF).1 It provides an effective potential to substitute for the classical potential in order to approximately include quantum effects in the dynamics. The results of the CWEQF model are more accurate than those of the CW model. Still the simplicity of the CW model is retained. The question then arises: can we extend this CWEQF model from onedimensional model systems to larger systems? This is the subject of the present paper where we apply the CWEQF model to simple reactions rate problems and compare with results from other methods, such as ring polymer molecular dynamics (RPMD)19,20 and the quantum instanton (QI) method.21 This paper is structured as follows: Section II describes a generalized theory for the one-dimensional CWEQF model. In Section III, extension to the two-dimensional H þ H2 collinear Special Issue: J. Peter Toennies Festschrift Received: January 27, 2011 Revised: May 6, 2011 Published: May 26, 2011 7338

dx.doi.org/10.1021/jp200886v | J. Phys. Chem. A 2011, 115, 7338–7345

The Journal of Physical Chemistry A

ARTICLE

Figure 1. Two extrema, (2.55, 2.55) and ( 2.55, 2.55), are seen in this picture of the matrix element Æx|F^(β/2|yæ, where x = q η/2 and y = q þ η/2, for β = 12 for the symmetric Eckart barrier. The η0 value can be calculated as η0 = 2.55 (2.55) = 5.1.

reaction is described. Section IV contains results, and finally we summarize our findings in Section V.

II. THEORY A. The CW Model. The thermal rate constant k(T) can be expressed as22 Z 1 ¥ kðTÞ ¼ dtCf f ðtÞ ð1Þ Qr 0

where Qr is the reactant partition function, and Cff(t) is the fluxflux autocorrelation function, which is given by ^ ^ F^ðβ=2Þ expð iHt=pÞ ð2Þ Cf f ðtÞ ¼ Tr½F^ðβ=2Þ expðiHt=pÞ with β = 1/kBT. The half-Boltzmannized flux operator F^(β/2) is β^ ^ β^ ^ ð3Þ F ðβ=2Þ ¼ exp H F exp H 4 4 where the flux operator F^ is given by 1 F^ðqds Þ ¼ ½^pδð^q qds Þ þ δð^q qds Þ^p 2m

ð4Þ

Here ^p and ^q correspond to the momentum and position operators respectively. qds specifies the dividing surface between the reactant and product. From refs 15 and 16, for any two ^ and B ^, we have operators A ^ expðiHt=pÞ ^ ^ expð iHt=pÞg ^ TrfA B Z Z 1 ¼ dp0 dq0 Aw ðq0 , p0 ÞBw ðq0 , p0 , tÞ 2πp Z Z 1 dp0 dq0 Aw ðq0 , p0 ÞBw ðqt , pt Þ 2πp

ð5Þ

where we have utilized the CW approximation in the last equality. ^ is defined as The Wigner transformation of an operator C Z ¥ ^ jq þ η=2 dη exp½ipη=p q η=2jC ð6Þ Cw ðq, pÞ ¼ ¥

path integral, which involves a potential difference between the two real-time paths that result from the two time evolution operators. η is the distance between these Feynman paths.16 The CW model approximates the propagation of the Wigner function as Bw(q0,p0,t) = Bw(qt,pt), where the propagation from an initial point (q0,p0) to a final point (qt,pt) is performed using a classical force f(q), which originates from linearizing the aforementioned potential difference as V(q η/2) V(q þ η/2) ≈ f(q)*η (see eqs 7 and 13 in ref 16), which holds when η ≈ 0. The CW model is accurate when the characteristic value of η is small, i.e., when the temperature is high.1 Using the rate constant expression in eq 1 together with eqs 2 and 5, we get Z ¥ Z 1 kðTÞ dt dp0 Qr ðTÞh 0 Z dq0 Fw ððβ=2, qds Þ; qt , pt ÞFw ððβ=2, qds Þ; q0 , p0 Þ ð7Þ where Fw is the Wigner transform of F^(β/2). Standard classical trajectories are used to propagate Fw in time. Since Cff(t) depends on the dividing surface, we will follow the same procedure to choose the dividing surface as we did in ref 1, i.e., we minimize Cff(0). B. The CWEQF Model Using a Position Independent η0. The CW model is accurate in the high temperature limit, but becomes poorer as the temperature is lowered. In a previous paper,1 we managed to find a distance η0 that is characteristic for the distance between the two Feynman paths in the expression for the fluxflux correlation function, which appears in the rate constant expression in eq 1. The quantum effects in the dynamics are reflected by this parameter, which can be thought of as estimating a delocalization effect. The CW model assumes η0 to be infinitesimally small, thereby removing quantum effects from the dynamics. In ref 1 we defined η0 to be the value of η that extremizes the function g(q,η) = Æq η/2|F^(β/2)|q þ η/2æ. This is illustrated in Figure 1. We note that q η/2 and q þ η/2 are the initial positions of the two Feynman paths. We also used the short time approximation that η0 is representative for the Feynman path separation also for later times. The force at a position q is expressed as f ef f ðqÞ ¼

The trace in eq 5 can be expressed in terms of Feynman path integrals. The two time evolution operators then yield a Feynman 7339

V ðq η0 =2Þ V ðq þ η0 =2Þ η0

ð8Þ

dx.doi.org/10.1021/jp200886v |J. Phys. Chem. A 2011, 115, 7338–7345


ARTICLE

Figure 2. The blue and green curves represent the linearized potential difference using classical and quantum forces respectively, the red curve represents the true potential difference and the pink curve represents g(q,η); all as a function of η for β = 12. The position is chosen as q = 0.5, i.e. rather close to the saddle point. The linearized and true potential differences match for η = η0 = 5.1, which is the η-value for which the matrix elements have their two extrema.

Figure 3. Effective potentials and the original symmetric Eckart potential for different values of β.

Table 1. Transmission Coefficients j = k/kclTST for a Symmetric Eckart Barrier Potential V(q) = b/cosh2(q*c) in Natural Units,30 a k

β=2

4

6

8

10

12

accurate

1.22 2.1 5.2 21.8 162 1970

CWEQF using position-dependent η0

1.16 1.8 4.4

19

149 1846

CWEQF using position-independent η0 1.16 1.8 4.6

20

150 1890

CW

1.15 1.7 3.3 9.9

QTSTb

1.23 2.1 5.7 29.3 248 3058

49

398

m = p = ω# = 1, b = 6/π, and c = (π/12)1/2, ω# is the barrier frequency. k is the relevant rate constant, and kclTST is the rate constant obtained by ‡ classical transition state theory with kclTST = eβV0 /hβQr, where V‡0 is the potential at the dividing surface.14 The QTST results in Tables 1 and 2 are from Pollak31. b From ref 31. a

The two points that we use to linearize the potential difference are separated by the distance η0. We can plot the true potential difference V(q þ η/2) V(q η/2) and the linearized potential difference feff(q) η for a fixed value of q for different η values. We can also plot the function of g(q,η) for different η values for the same q value. This is seen in Figure 2, for the symmetric Eckart barrier defined in Table 1. If we fix q at a value close to its saddle point value, the η value where the linearized potential difference agrees with the true potential difference, agrees with the characteristic distance η0. The quantum effects in the dynamics are captured quite well by using η0 to obtain an effective quantum force. The effective potentials and the original potential for the symmetric Eckart barrier at different temperatures are plotted in Figure (3). The barrier height of the effective potential becomes lower as the temperature is lowered. This approximately accounts for the quantum effect of tunneling through the barrier. C. The CWEQF Model Using Position Dependent η0 Values. The CWEQF model with a position-independent η0 worked well for a symmetric Eckart barrier at all temperatures. Still, if we choose a position some distance away from the dividing surface,

Figure 4. The linearized potential difference using the quantum force (with a position-independent value of η0), the true potential difference, and g(q,η) are shown as a function of η for β = 12. The position is chosen as q = 1.8, whereby the corresponding agreement between the linearized and true potential differences occurs for η = 5.1. The matrix element is, however, at its peak when η = 3.6, so there is a mismatch.

for example, q = 1.8, and make the corresponding plot to Figure 2, we arrive at Figure 4. It is seen that there is a mismatch between where the curve for the true potential difference crosses the linearized potential difference between the Feynman paths and where the extrema of g(q,η) are found. This implies that the CWEQF model with a position independent η0 can be improved. In our previous paper,1 using a position-independent η0, we encountered a problem for an asymmetric Eckart barrier at low temperature (the potential is defined below in Table 2). The matrix elements g(q,η) (see Section II.B) sometimes contained two pairs of extrema instead of only one pair (see Figure 5). We handled this by splitting the flux distribution function into two different parts with different η0 values. If a trajectory moved between the two parts, we defined the η0 value to be the average of the two individual η0 values. However, this procedure is not easy to apply to general potentials. Here we improve the CWEQF model from this point of view by calculating η0 values as a function of q. For a certain value of q we search for the value 7340



ARTICLE

of η that corresponds to the extrema of g(q,η). This value is denoted η0(q). By using the short time approximation, this value is not changed for later time. The characteristic distances for a symmetric and an asymmetric Eckart barrier, for different temperatures, are shown in Figures 6 and 7, respectively. The effective force is then η ðqÞ η ðqÞ V q 0 V qþ 0 2 2 f ef f ðqÞ ¼ η0 ðqÞ

the following section, we apply the CWEQF model using position-dependent η0 values to a two-dimensional case.

ð9Þ

In the following, we omit the dependence of η0 on q in the notation. In Figure 8 we compare the linearized potential difference obtained from the quantum force in eq 9 with the true potential difference for q = 1.8. The function g(q,η) is also presented. It is seen that the η value where the curve for the true potential difference crosses the linearized potential difference agrees with where we find the extrema of g(q,η). An advantage of the CWEQF model using position-dependent η0 values is that it can be applied to potentials resulting in more than one pair of extrema. A drawback is that we have to specify the characteristic distances for different positions compared to finding only one such value in the CWEQF model using a position-independent η0 value. However, the computational effort to obtain all the η0 values for different positions is trivial compared to the effort required to obtain the Wigner function. In

Figure 6. η0 values as a function of q for a symmetric Eckart barrier for different values of β.

Table 2. Transmission Coefficients for an Asymmetric Eckart Potential V(q) = A/(1 þ exp(2q/C)) þ b/cosh2(q/c), with30 a = 18/π, b = 13.5/π, c = 8/(3π)1/2, and m = p = ω# = 1 k accurate CWEQF using position-dependent η

a

β=2 4

6

8

10

1.2 2.0 5.3 26 250 1.17 1.9 5.1 23.5 172

12 4100 3646

CWEQF using position-independent η 1.17 1.9 5.4 24

250

CW

1.16 1.8 4.1 21

374 29300

QTSTa

1.2

2.0 5.6 44.0 1100 28000

RPMDb

1.2

2.0 5.3 28

From ref 31. b From ref 19.

310

6800

5900

Figure 7. η0 values as a function of q for an asymmetric Eckart barrier for different values of β.

Figure 5. Two pairs of extrema for the matrix element Æx|F^(β/2|yæ for an asymmetric Eckart barrier at β = 12. 7341



ARTICLE

Figure 9. (ηx0,ηy0) for the collinear LSTH H þ H2 potential at 200 K. Figure 8. The linearized potential difference using the quantum force (with position dependent η0 values), the true potential difference and g(q,η) as a function of η for β = 12. The position is chosen as q = 1.8, for which the matching between the true and linearized potential difference occurs for η = 3.6. This η value also gives the extrema of the matrix element.

III. TWO-DIMENSIONAL H þ H2 RATE CONSTANT CALCULATION A. Theoretical Construction. Here we apply the CWEQF method to the two-dimensional benchmark HA þ HBHC f HAHB þ HC collinear reaction. To proceed with this, we specify the reaction coordinate x and the corresponding orthogonal coordinate y as23 rffiffiffi 1 ‡ ‡ ðRAB RAB y¼ þ RBC RBC Þ, 2 rffiffiffi 1 ðRBC RAB Þ x¼ ð10Þ 2

where RAB and RBC are the bond distances between different hydrogen atoms. (R‡AB, R‡BC) specifies the saddle point. The product side is defined as x > 0, and the reactant side has x < 0. For the collinear reaction we always have RAC = RAB þ RBC. The Hamiltonian for this system is23 H ¼

1 2 1 2 px þ p þ V ðx, yÞ, mx ¼ mH =3, 2mx 2my y my ¼ mH

ð11Þ

Following Yamamoto et al.,22 using the fact that the Boltzmann operator commutes with the time evolution operator, we can get a symmetrized fluxflux correlation function for the twodimensional case C2d f f ðtÞ ¼

1 ð2πpÞ2

Z

Z dpx

Z dx

dpy

Z

dyFw ððβ=2, xds Þ; xðtÞ, yðtÞ, px ðtÞ, py ðtÞÞFw ððβ=2, xds Þ; x, y, px , py Þ

ð12Þ The rate constant for the two-dimensional calculation is then Z ¥ 1 k2d ðTÞ ¼ dtC2d ð13Þ ff ðtÞ Qr ðTÞ 0

At initial time (t = 0), the integrand of eq 12 is always positive. The symmetrized fluxflux correlation function is therefore easier to converge than the fluxside correlation function.22 For a reaction in two dimensions, η0 becomes a two-dimensional vector, ηB0= (ηx0,ηy0), which is illustrated in Figure 9. It can be seen that the ηB0 vectors on the minimum energy path (MEP) point along it. For a point (x,y), we may search for the momentum vector (px,py) which corresponds to the extrema of Fw((β/2,xds);x,y,px,py). This vector is denoted B p 0. We have visually noticed that the directions of p 0 are the same. the ηB0 vector and B Let (x,y) and (x0 ,y0 ) define the points on the two real-time Feynman paths. To construct the effective quantum force is simple in the one-dimensional case. However, for the twodimensional case, where η0 is a vector, it is not obvious how to linearize the potential difference V(x,y) V(x0 ,y0 ) in order to construct the effective quantum force. We now define B r = (rx,ry) and ηB0 = (ηx,ηy) and we let rx = (x þ x0 )/2, ry = (y þ y0 )/2, ηx = x x0 , and ηy = y y0 . The potential difference becomes V( B r þ ηB/2) V( B r ηB/2) ≈ fBeff 3 ηB. The force vector contains two components. To determine them, we need two equations. The first equation is obtained by linearizing the potential difference between the two extrema of the elements, which are connected by the blue arrow in Figure (10). This gives one equation: ! ! ! ! ef f ! η 0ð B rÞ η 0ð B rÞ rÞ ¼ V B r V B r þ B f 3 η 0ð B 2 2 ð14Þ For each fixed B r , there is only one pair of extrema of the matrix elements when (ηx, ηy) is varied. However, eq 14 does not make the linearized surface unique. It can be rotated around the axis of ηB0. We resolve this problem by finding a pair of points along the orthogonal direction to ηB0. The characteristic distance in the orthogonal direction to ηB0 should be chosen infinitesimally small, as illustrated in Figure 10. This results in ! fxef f ð B rÞ ð! η ^0 ε þ ! η 0Þ rÞ 3 fyef f ð B ¼ Vð B r þ! η 0 =2 þ ! η ^0 ε=2Þ V ð B η ^0 ε=2Þ r ! η 0 =2 ! ð15Þ Here we have ε f 0 and ηB^0 = {ηy0/[(ηx0)2 þ (ηy0)2]1/2, ηx0/[(ηx0)2 þ (ηy0)2]1/2} is the unit vector in the direction 7342



ARTICLE

Figure 10. Plot of the matrix element Ærx ηx/2, ry ηy/2|F^(β/2)|rx þ ηx/2, ry þ ηy/2æ as a function of ηx/2 and ηy/2. The beginning and ending points of the blue arrow represent the extrema of the matrix elements, thereby defining ηx0 and ηy0. The yellow vectors specify the infinitesimal movement along the orthogonal direction to ηB0. T = 1000K, rx = 0, and ry = 0.5.

orthogonal to ηB0. Combining eqs 14 and 15, we obtain fx ðx, yÞ ðV ðx þ ηx0 =2, y þ ηy0 =2Þ V ðx ηx0 =2, y ηy0 =2ÞÞηx0

¼

ðηx0 Þ2 þ ðηy0 Þ2 0:5 ðVy ðx ηx0 =2, y ηy0 =2Þ þ Vy ðx þ ηx0 =2, y þ ηy0 =2ÞÞηx0 ηy0 0

þ

0

2

ðηx0 Þ þ ðηy0 Þ2 0:5 ðVx ðx ηx0 =2, y ηy0 =2Þ þ Vx ðx þ ηx0 =2, y þ ηy0 =2ÞÞðηy0 Þ2 0

0

2

ðηx0 Þ þ ðηy0 Þ2

fy ðx, yÞ ¼

ðV ðx þ ηx0 =2, y þ ηy0 =2Þ V ðx ηx0 =2, y ηy0 =2ÞÞηy0 ðηx0 Þ2 þ ðηy0 Þ2 0:5 ðVy ðx ηx0 =2, y ηy0 =2Þ þ Vy ðx þ ηx0 =2, y þ ηy0 =2ÞÞðηx0 Þ2 0

0

x, yjF^ðβ=2Þjx0 , y0

Z

ðηx0 Þ2 þ ðηy0 Þ2 0:5 ðVx ðx ηx0 =2, y ηy0 =2Þ þ Vx ðx þ ηx0 =2, y þ ηy0 =2ÞÞηy0 ηx0 0

þ

the remaining n 1 orthogonal directions and match the potential difference as in eq 15, which will result in other n 1 equations. The n-dimensional force will be obtained by solving these n equations. B. Numerical Scheme. In this section we give a recipe for the numerical application of the CWEQF model to the collinear H þ H2 reaction in four steps in order to obtain the thermal rate constant (eq 13) from the fluxflux correlation function (eq 12). The first three steps concern initial time t = 0, while the fourth concerns the dynamical aspects. This is followed by a brief description of an accurate calculation of the thermal rate constant. We use the LiuSiegbahnTruhlarHorowitz (LSTH) potential24,25 for this calculation. Step 1: Matrix Elements of ^F(β/2). The first step is to determine the matrix elements Æx,y|F^(β/2)|x0 ,y0 æ. We partition these matrix elements into a set of submatrix elements as follows:

0

¼

2

ðηx0 Þ þ ðηy0 Þ2

ð16Þ Here V0x means the derivative of the potential with respect to x with the corresponding notation for the y coordinate. This way of specifying eq 16 is not guaranteed to be optimal, and there may be other forms that give more accurate rate constants. It is also straightforward to generalize this method to n-dimensional (n > 2) problems. An n-dimensional ηB0-vector connects the two extrema of the matrix elements Æx0, x1, ...xn|F^(β/2)|x00 , x10 , ...xn0 æ and matches the potential difference as in eq 14. Then we can take infinitesimal steps along

D E 0 0 0 0 dx1 dy1 dx2 dy2 :::: dxn dyn dx1 dy1 ::::: dxn dyn x, yjeβF H=n jx1 , y1 D E D E x1 , y1 jeβF H=n jx2 , y2 ::::::: xn1 , yn1 jeβF H=n jxn , yn D ED 0 0 E 0 0 0 0 xn , yn jF^ðβflux Þjx1 , y1 x1 , y1 jeβF H=n jx2 , y2 ::::: D 0 ED 0 0 E 0 0 0 xn 1 , yn 1 jeβF H=n jxn , yn xn , yn jeβF H=n jx0 , y0 ð17Þ

We define βF and βflux as 2βF þ βflux ¼ β=2, βF =n ¼ βflux =2

ð18Þ

We choose n = 15 for 1000 K and n = 20 for 200 K and 300 K in order to keep the temperature in each submatrix high, which means that βF/n and βflux will be small. We transform the coordinates as x = (x þ x0 )/2, Δx = x x0 , y = (x þ x0 )/2, and 7343



ARTICLE

position by looping over ηB for each B r to find the ηB = ηB0 which gives the extrema for the matrix elements ÆrB ηB/2|F^(β/2)|rB þ ηB/2æ. Then we use these ηB0 vectors to construct the effective force according to eq 16. Step 4: Time Propagation. The effective quantum force is used to propagate the initial Wigner distribution in time so as to obtain Fw at time t and insert into eq 12. We use the velocity Verlet method26 for the time propagation. Using eq 12 in eq 13 completes the numerical procedure. Accurate Thermal Rate Constant. In two dimensions, it is feasible to obtain the accurate thermal rate constant for the collinear H þ H2 reaction, which can be compared with our CWEQF rate constant. An expression for the accurate rate constant is Z t 1 kðTÞ ¼ lim dtTrfF^ðβ=2ÞF^ðβi =2Þg Figure 11. Comparison of the classical and effective MEPs at T = 300 K.

Δy = y y0 . Then we have by standard manipulations Δxmx mx ðx xdiv Þ2 =2βflux p2 mx ðΔxÞ2 =2βflux p2 x, yjF^ðβflux Þjx0 , y0 e e πβ2flux p3 rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi my my ðΔyÞ2 =2βflux p2 βflux V ðx, yÞ e 2e 2πβflux p

¼

emy ðΔyÞ n=2βF p eβF V ðx, yÞ=n 2

2

ð20Þ

Here, xdiv is the position of the dividing surface. We define an effective MEP to be the path that follows the largest contribution to the matrix element of F^(β/2) as a function of ^r = (x,y). This effective MEP coincides with the classical MEP at high temperature, but the two curves deviate from each other at lower temperature. It is only necessary to calculate the matrix elements around the effective MEP. The grids for the x and y coordinates are chosen as (Δx,Δy) = (0.05,0.06) with the number of grid points (Nx,Ny) = (80,50) for 1000 K and (Δx,Δy) = (0.075,0.075) with the number of grid points (Nx,Ny) = (80,80) for 300 K and 200 K. To obtain the effective MEP with a limited CPU usage, we first run a small calculation to get the matrix elements on a sparse grid in the y coordinate, keeping its range large. In this way we generate a coarse effective MEP. Next we use this coarse effective MEP in another small calculation using the same number of points for the y coordinate but with reduced range. Finally we obtain the matrix elements Æx,y|F^(β/ 2)|x0 ,y0 æ on a dense grid around the effective MEP. The effective MEP and the classical MEP are plotted in Figure 11. Step 2: Wigner Transform of ^F(β/2). In this step, we numerically Wigner transform the initial (t = 0) matrix elements in eq 17 Z Z according to dηx dηy Fw ðx, y, px , py Þ ¼ expði B p 3! η =pÞ B r ! η =2jF^ðβ=2Þj B r þ! η =2 ð21Þ Step 3: The Effective Quantum Force. In this step, we find the effective quantum force. First we determine ηB0 as a function of

dt 0

Qr t f ¥ 0

Z

¥

Z dx

Z

dy

dx0

Z

dy0 x, yjF^ðβ=2Þjx0 , y0 x0 , y0 jF^ðβi =2Þjx, y

ð22Þ βi ¼ β þ 4it with

ð19Þ sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi D E sffiffiffiffiffiffiffiffiffiffiffiffiffiffi my n mx n mx ðΔxÞ2 n=2βF p2 βF H=n 0 0 x, yje jx , y 2e 2πβF p2 2πβF p

1 Qr

Z

! β i ^ F^ exp βi H ^ F^ðβi =2Þ ¼ exp H 4 4

ð23Þ

The matrix elements are calculated as in eq 17 in step (1) but changing β/2 to β/2 þ 2it. There is no approximation in this calculation. Using eq 22 to obtain the rate constant we reproduce the results from the accurate quantum dynamical calculation of Bondi et al.27

IV. RESULTS We first present one-dimensional transmission coefficients, which are directly related to the rate constants, for a symmetric Eckart barrier in Table 1 and for an asymmetric Eckart barrier in Table 2. From these tables it is seen that the CWEQF model with position-dependent η0 values and the model with positionindependent η0 values perform essentially equally well. Both of them give results that are substantially closer to the accurate ones than the CW and quantum transition state theory (QTST) models at low temperature. For the asymmetric Eckart barrier, there are also RPMD results available which are of comparable accuracy to the CWEQF results. Next we turn to the two-dimensional calculations for the collinear HþH2 reaction. We plot the autocorrelation functions at T = 1000 K and T = 200 K in Figures 12 and 13, respectively. At the lower of these temperatures the improvement of the autocorrelation function obtained for the CWEQF model compared to the CW model is apparent. The corresponding rates are shown in Table 3. It is seen that the CW model underestimates the rate constant, especially at the lower temperature 200 K, where the rate constant from the CW model is only one-fifth of the exact one. The CWEQF model works consistently better than the CW model, which is particularly noticeable at low temperature, which follows from the behavior of the autocorrelation functions. In Table 3 we also show results from previous work using the CW method based on the flux-side correlation function28 and QTST.23 These methods overestimate the accurate rate constant 7344



ARTICLE

Figure 12. Comparison of an accurate Cff with those for the CW and CWEQF models at T = 1000 K.

and QTST calculations. The results are of comparable accuracy as those of the RPMD model. An advantage of the RPMD model is its low demand for computational power. As for the two-dimensional H þ H2 calculation, we have shown that the CWEQF model consistently produces better results than the ordinary CW method and QTST calculations, while QI calculations are quite accurate. We note that the rate constants obtained from the CW model depends strongly on whether the fluxflux or the fluxside autocorrelation function is used as seen in Table 3. There may be better ways to construct the effective force for multidimensional calculations than what we have suggested here. A way to update the ηB0 vector that goes beyond the present short time approximation might further improve the accuracy of the results. However, the presented procedure is reasonably accurate and simple. In particular, given the ηB0 vector, the evaluation of the effective quantum force requires only the classical potential and its derivative. The main limitation for application in high dimensional calculations is the difficulty to obtain the converged multidimensional thermal flux Wigner distribution function. There has, however, been recent progress on this subject using a variational method.29 The dynamics can then be run using the effective quantum force in the manner presented here.

’ REFERENCES

Figure 13. Comparison of an accurate Cff with those for the CW and CWEQF models at T = 200 K.

Table 3. Rate Constants (in cm molecule1 s1) for the TwoDimensional H þ H2 Collinear Reaction k

200 K

300 K

1000 K

accuratea

6.2 102

4.8

6.7 103

CWEQF:

3.6((0.05) 102

2.7((0.03)

5.46((0.01) 103

CW

1.3((0.05) 102

1.6((0.03)

5.0((0.01) 103

1

CW using Cfs (1.74 ( 0.19) 10 b

1

2.3((0.4) 10

c

QTST

2

5.3 10

d

QI

6.31 ( 0.21 (6.92 ( 0.08) 103 9.2 ( 0.7

(7.1 ( 0.1) 103

5.19

7.24 103

From ref 27. CW model based on flux-side correlation functions from ref 28. c QTST from ref 23. d QI theory from ref 21. a

b

by a factor of 3 and 4, respectively at 200 K. There are also results from the QI method available. These are close to the accurate results for the temperature shown. However, contrary to the CW method, the QI model does not describe recrossing, whereby it does not produce the exact results in the high temperature limit.21

V. SUMMARY AND DISCUSSION It has been shown that the CWEQF model works well both for one and two-dimensional rate constant calculations. In particular, it retains the simple mechanism of the CW model and approximately accounts for quantum effects through the quantum force. In one dimension, the results we have obtained are highly satisfactory and more accurate than results from the CW model

(1) Poulsen, J. A.; Li, H. Q.; Nyman, G. J. Chem. Phys. 2009, 131, 024117. (2) Wyatt., R. E. J. Chem. Phys. 1999, 111, 4406. (3) Eric., R.; Bittner. J. Chem. Phys. 2000, 112, 9703. (4) Kendrick., B. K. J. Chem. Phys. 2004, 121, 2471. (5) Goldfarb, Y.; Tannor, D. J. J. Chem. Phys. 2007, 127, 161101. (6) Bohm, D. Phys. Rev. 1952, 85, 166. (7) Bohm, D. Phys. Rev. 1952, 85, 180. (8) Wigner, E. Phys. Rev. 1932, 40, 749. (9) Donoso, A.; Marten., C. C. Phys. Rev. Lett. 2001, 87, 223202. (10) Liu, J.; Miller, W. H. J. Chem. Phys. 2007, 126, 234110. (11) Frantsuzov, P.; Neumaier, A.; Mandelshtam, V. A. Chem. Phys. Lett. 2003, 381, 117. (12) Frantsuzov, P.; Mandelshtam, V. A. J. Chem. Phys. 2004, 121, 9247. (13) Heller, E. J. J. Chem. Phys. 1976, 65, 1289. (14) Wang, H.; Sun, X.; Miller, W. H. J. Chem. Phys. 1998, 108, 9726. (15) Qiang, S.; Eitan, G. J. Chem. Phys. 2003, 118, 8173. (16) Poulsen, J. A.; Nyman, G.; Rossky, P. J. J. Chem. Phys. 2003, 119, 12179. (17) Hernandez, R.; Voth, G. A. Chem. Phys. 1998, 233, 243. (18) Sun, X.; Wang, H.; Miller, W. H. J. Chem. Phys. 1998, 109, 4190. (19) Craig, I. R.; Manolopoulos, D. E. J. Chem. Phys. 2005, 123, 034102. (20) Collepardo-Guevara, R.; Suleimanov, Y. V.; Manolopoulos, D. E. J. Chem. Phys. 2009, 130, 174713. (21) Ceotto, M.; Miller, W. H. J. Chem. Phys. 2004, 120, 6356. (22) Yamamoto, T.; Wang, H.; Miller, W H. J. Chem. Phys. 2002, 116, 7335. (23) Liao, J.-L.; Pollak, E. J. Phys. Chem. A 2000, 104, 1799. (24) Liu, B. J. Chem. Phys. 1973, 58, 1925. (25) Liu, B.; Siegbahn, P. J. Chem. Phys. 1978, 68, 2457. (26) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Oxford University Press: Oxford, U.K., 1987. (27) Bondi, D. K.; Clary, D. C.; Connor, J. N. L.; Garrett, B. C.; Truhlar, D. G. J. Chem. Phys. 1982, 76, 4986. (28) Zheng, Y.; Pollak, E. J. Chem. Phys. 2001, 114, 9741. (29) Poulsen, J. J. Chem. Phys. 2011, 134, 034118. (30) Jang, S.; Voth, G. A. J. Chem. Phys. 2000, 112, 8747. (31) Pollak, E.; Liao, J.-L. J. Chem. Phys. 1998, 108, 2733. 7345


The Classical Wigner Method with an Effective Quantum Force

Recommend Documents