Tunneling Dynamics Using Classical-like Trajectories with an Effective

Letter pubs.acs.org/JPCL

Tunneling Dynamics Using Classical-like Trajectories with an Effective Quantum Force Li Huaqing, Jens Poulsen, and Gunnar Nyman* Department of Chemistry and Molecular Biology, University of Gothenburg, 41296 Gothenburg, Sweden ABSTRACT: We develop an effective quantum force to propagate trajectories in Wigner phase space. In this way, the quantum entanglement between trajectories is represented by the effective quantum force, which is related to the variance of the momentum distribution. Reaction probabilities are obtained for three model potentials. For a cubic potential, the reaction probabilities are found to continuously grow in time long after the classical reaction probabilities have become stationary. This growth is due to tunneling and is in good agreement with numerically exact quantum results. The variance of the momentum distribution is updated at each time step of the trajectories, which is the key to the increase of the reaction probabilities at longer times and allows a nice interpretation of tunneling.

SECTION: Molecular Structure, Quantum Chemistry, and General Theory

E

suffer from issues such as unstable prefactors, cumbersome Taylor expansions of the potential energy surface, or multidimensional curve-fitting procedures that have plagued many of the other methods. Of course, in comparison, CW, CMD, and RPMD can be argued to be inferior, since they do not account for interference effects. The CW method is in fact so simple that it utilizes only “ordinary” classical trajectories and all quantum effects that it treats are included by simply adopting a Wigner distribution for the sampling of the initial conditions. The purpose of this paper is to “upgrade” the CW model by adding dynamical quantum effects to it. This is accomplished by replacing the classical force by an effective time-dependent one. A key parameter of the model is the momentum dispersion of the particles, which is being used as an inverse measure of quantum delocalization and goes into the expression for the effective force. The momentum dispersion is computed from the spread of the momentum, which is obtained on-the-fly from nearby trajectories. The Wigner function corresponding to the density operator ρ̂(t) is22

ssentially all quantum methods aimed for describing quantum effects in anharmonic many-dimensional adiabatic systems rely on defining appropriate “classical trajectories” for the system of interest. This is so since the exponential scaling of basis set methods must be avoided. Among the trajectory-based methods we find ring polymer molecular dynamics (RPMD),1,2 centroid molecular dynamics (CMD),3 the semiclassical initial value representation (SCIVR),4 entangled trajectory molecular dynamics (ETMD),5−9 Bohmian mechanics (BM),10,11 quantized Hamilton dynamics,12 forward−backward IVR,13 the classical Wigner method (CW),14−18 coupled coherent states (CCS),19 and frozen Gaussian cellular dynamics (FGCD).20 Utilizing trajectories for all degrees of freedom is advantageous for several reasons. First, the problem of how to mix quantum and classical degrees of freedom is avoided. This is a central issue in surface hopping schemes,21 where the primary quantum degree of freedom is not described by trajectories. Another benefit is the possibility of sampling trajectories by Monte Carlo; see, e.g., applications of RPMD, CMD, and CW. The above methods differ in how the trajectories account for quantum effects. For instance, the nonlocality of quantum mechanics is in BM, ETMD and CCS described by introducing interactions between the trajectories. On the other hand, semiclassical methods such as the SC-IVR and the FGCD methods use uncoupled trajectories, but instead augment each trajectory with a phase and a preexponential factor. Of the mentioned methods, only RPMD, CMD, and CW have been applied to large quantum systems composed of hundreds of degrees of freedom with chemically realistic potentials. The main reason is that these methods do not © 2013 American Chemical Society

ρ W (q , p , t ) =

1 2π ℏ

∫

q + y/2|ρ ̂(t )|q − y/2 e−ipy / ℏ dy (1)

with q and p being the position and momentum coordinates. y is related to the importance of the off-diagonal elements ⟨q + y/2|ρ̂|q − y/2⟩. Received: June 18, 2013 Accepted: August 20, 2013 Published: August 20, 2013 3013

dx.doi.org/10.1021/jz4012564 | J. Phys. Chem. Lett. 2013, 4, 3013−3018

The Journal of Physical Chemistry Letters

Letter

Figure 1. Reaction probability versus time for the cubic potential at three different average energies. The blue curves represent the accurate quantum mechanical results.7 The red curves represent the results obtained using our effective force. The green curves represent the ETMD results.7 The pink curves represent results obtained using the classical force.

as shown by Martens and collaborators.7 In these equations ṗ is the force, i.e., the time derivative of the momentum and q̇ is the velocity. From eq 2 we can see that the force depends on both the momentum and the position of “other trajectories”. This is obviously different from the classical force. The trajectories are said to be entangled with each other and thus the name entangled trajectory molecular dynamics (ETMD) is used to describe dynamics based on eq 2.8 The ETMD method has been applied to several model potentials and found to yield reaction probabilities in good agreement with numerically exact results.5,7 These examples involve tunneling so they illustrate situations where this phenomenon is handled quite well via the entanglement in the ETMD method. To calculate the ETMD force for a point in phase space, there are two integrations to carry out, which might be numerically expensive. Also, eq 2 for the force diverges when ρW(q, p, t) → 0 and there is a difficulty in regions where the Wigner function is negative. These last two issues can be overcome by introducing the Husimi representation whereby no negative regions remain as a result of averaging over regions too small to be distinguishable when accounting for the uncertainty principle.9,23 In the classical Wigner (CW) model14,16−18 classical trajectories are run but with the initial conditions sampled from a Wigner distribution. Thus quantum effects enter through the initial conditions but the dynamics is classical. Therefore it is not capable of describing deep tunneling nor interference effects. The classical trajectories may for example be used to evaluate a flux−flux correlation function in the short time limit, which in turn can be used to give thermal rate constants. We have in previous work24,25 modified the force used in the CW model to become an effective quantum force. In ordinary CW calculations, the force is evaluated at

Figure 2. The explanation of the tunneling effect for the cubic potential.

Requiring that the time evolution of the Wigner function obeys the continuity equation, the quantum mechanical equations of motion for trajectories may be formally written ṗ = −

1 ρ (q , p) W

∫ dy

⎡ py ⎤ 1 exp⎢ −i ⎥ ⎣ ℏ ⎦ 2π ℏ

∫

V (q + y/2) − V (q − y/2) y ⎡ yξ ⎤ W dξ exp⎢i ⎥ρ (q , ξ) ⎣ ℏ⎦

(2)

and q̇ =

p m

(3) 3014



Letter

Figure 3. The classical cubic potential and the corresponding effective potentials at three different times.

the middle of the two Feynman paths that enter the expression for the flux−flux correlation function. We instead evaluated the force from the difference of the potential energy at two points. The two points are separated by a characteristic distance that is chosen to correspond to a separation between the Feynman paths, which contributes most to the thermal rate constant. In the present work we begin from eq 2 and we are interested in propagating the Wigner distribution corresponding to the density operator, rather than to the flux operator, forward in time. We wish to retain an approximate treatment of dynamical quantum effects but without the force being dependent on the momenta of the trajectories, i.e., we wish to substitute the entanglement of the trajectories by a simple delocalization parameter. In eq 2, the integration over ξ is the inverse Fourier transform of the Wigner function which is defined in eq 1. Thus 1 ṗ = W ρ (q , p)

∫

In the limit where the off diagonal elements of the density matrix are zero, classical dynamics becomes exact. This is an assumption of the CW method,14,16−18 which uses the classical force to propagate trajectories in Wigner space. Classical dynamics does not treat tunneling. In line with our earlier work,24 we wish to find a characteristic value y0(q) to approximately treat quantum effects. As before, we set V (q + y/2) − V (q − y/2) y =

ṗ = −

−V ′(q) W

ρ (q , p)

= −V ′(q)

×

∫

(6)

V (q + y0 (q)/2) − V (q − y0 (q)/2) y0 (q)

(7)

rather than the classical expression in eq 5. As already mentioned, in earlier work we found the characteristic distance y0 as the value of y corresponding to the largest element of the matrix ⟨q + y/2|F̂(β/2)|q − y/2⟩. Since the present work concerns the density rather than the flux and the matrix elements of the density operator peak along the diagonal, we need another way to find the characteristic distance, or else we will obtain ordinary classical dynamics. We choose to relate the characteristic distance to the width of the Wigner distribution function corresponding to the density operator. For this purpose we locally in q approximate the Wigner function to be a Gaussian. This distribution fulfills the minimum uncertainty principle and the variance of the momentum σp(q) is related to the to the variance of the position σq(q) by σp(q)σq(q) = ℏ/2. The smallest phase space area compatible with the uncertainty principle is called a quantum blob.26 For a

(4)

The elements of the density matrix peak along the diagonal, thereby amplifying the importance of small values of y in the integration in eq 4. We can therefore Taylor expand the potential difference to obtain V(q + y/2) − V(q − y/2) ≈ V′(q) × y. The force then becomes ṗ =

y0 (q)

so that the force in eq 4 becomes

V (q + y/2) − V (q − y/2) dy y

⎡ py ⎤ × exp⎢ −i ⎥ q + y/2|ρ ̂|q − y/2 ⎣ ℏ⎦

V (q + y0 (q)/2) − V (q − y0 (q)/2)

⎡ py ⎤ dy exp⎢ −i ⎥ q + y/2|ρ ̂|q − y/2 ⎣ ℏ⎦ (5) 3015



Letter

Figure 4. Reaction probability versus time for four different initial energies. (a) Symmetric Eckart barrier. (b) Asymmetric Eckart barrier. Red curves are results obtained using the effective quantum force. Green curves are results from the ETMD method.7−9 The blue curves represent accurate results.7−9

we define y0(q) ≡ ℏ/2σp(q) = σq(q). We can update σq(q) by first calculating the variance of the momentum at the position q. The only property that we have used from the Gaussian distribution is that it fulfills the minimum uncertainty principle. We use the above discussion as a motivation for our calculation of y0 for our actual distribution, for which we obtain the variance of the momentum, Δp, as

Gaussian distribution function, the quantum blob has an area of ℏ/2. If we for a position q find the variance of the momentum to be σp(q), then the resolution of this position is limited to ℏ/2σp(q). To evaluate the force in eq 7 requires two points in phase space that are separated by a distance y0(q). We choose to determine this distance by considering the area of the quantum blob. Let us first see what happens if we assume that we have a Gaussian distribution. With arguments from above, 3016


The Journal of Physical Chemistry Letters Δp = p̅ = p2 =

Letter

p2 − (p ̅ )2

Ψ(q , 0) =

∫ dpρ W (q , p , t )p/∫ dpρ W (q , p , t ) ∫

dpρ W (q , p , t )p2 /

∫

dpρ W (q , p , t )

(11)

where the q0 corresponds to the initial center of the position and p0 to an initial momentum. The initial Wigner distribution function is then

(8)

⎛ (q − q )2 (p − p0 )2 ⎞ 1 0 ⎜ ⎟ − ρ (q , p , 0) = exp⎜ − πℏ 2σq2 2σp2 ⎟⎠ ⎝

We then simply set

W

y0 (q) = ℏ/2Δp

(9)

y0 (q) = 2 p2

ℏ ⎛ ⎜1 − ⎝

7(t ) = ∞

=

⎛ (p )2 ⎞ ⎜⎜1 + ̅ ⎟⎟ 2 p2 ⎠ 2 p2 ⎝

∫q

1/2

∞

P(q , t )h(q) dq div

∞

∫−∞ dq ∫−∞ dpρW (q , p , t )h(q − qdiv)

(13)

where qdiv is the dividing surface that separates reactants from products. P(q,t) is the probability distribution for position q at time t. A. Cubic Potential. First we calculate reaction probabilities for a particle moving in a cubic potential given by V (q) =

1 1 mω02q2 − bq3 q < qd 2 3 V (q) = Vd

q > qd

(14)

Using atomic units throughout this paper we choose m = 2000, ω0 = 0.01, b = 0.2981, Vd = −0.015 and qd = 1.12556.7 This gives a minimum V = 0 at q = 0 and a maximum V0 = V(qdiv) = 0.015 at qdiv = 0.6709. The initial momentum is set to be zero (p0 = 0) and three initial positions are chosen in order to produce three different initial energies. The initial positions are q0 = −0.2, q0 = −0.3, and q0 = −0.4 corresponding to the average energies E0 = 0.75V0, E0 = 1.26V0, and E0 = 2.0V0, respectively.27 Ten thousand trajectories were run with the time step δt = 5. We use 300 boxes to discretize the region (−1.5,1.5) in order to calculate the y0 values. The resulting reaction probabilities are presented in Figure 1. The present calculations using the classical Wigner effective quantum force (CWEQF) in eq 7 reproduce the accurate quantum mechanical results quite well, as do the ETMD results. The difference to the classical results is due to tunneling, for which the effective quantum force expression offers an interesting explanation. The dominant part of y0(q) is obtained from ℏ/(2(p2)1/2), as seen in eq 10. The trajectories with the highest momenta tend to pass over the barrier to become products early on and thereafter they are not included when calculating p2. Thus, the average kinetic energy, or p2, of the reactants decreases with time. The decrease in kinetic energy results in an increase of y0(q), which in turn reduces the effective potential barrier height. This allows more trajectories to leave the well, thus increasing the accumulated reaction probability and also further lowering the average kinetic energy of the remaining reactants. This procedure continues, thereby allowing trajectories that were initially trapped in the well to later pass over the barrier as its height effectively is being reduced with time, considering how the effective quantum force is

(p ̅ )2 ⎞ p2

(12)

where σq = (ℏ/(2mω)) and σp = (ℏmω/2) which specifies the width of the distribution along q and p. The reaction probability is expressed as 1/2

Thereby we are able to calculate the force in eq 7. To evaluate this force does not require the double integration as in eqs 2 and 4, thereby being numerically much faster. For each spatial position in phase space, we need to integrate over the momentum space only once to get the variance of the momentum. The force generated in the ETMD method describes “nonlocal” contributions via entanglement. In the present work, the entanglement between the trajectories with different momenta is absorbed into the parameter y0(q), which is used to obtain the quantum force. Classical trajectories are completely local in nature. y0(q) has the effect of bringing in the “nonlocality” or delocalization of quantum mechanics in an approximate way. For illustration we will calculate reaction probabilities for three simple one-dimensional model problems. We will initiate the trajectories by Monte Carlo sampling the initial conditions. To obtain the value y0(q), we need the variance of the momentum at the position q. To do this we divide q into a number of boxes, which each contains a number of trajectories. Some boxes may, however, contain only a few or even zero trajectories at a particular time. In order to avoid numerical trouble in calculating the variance of the momentum, we adopt the following Taylor expansion procedure for obtaining the value of y0. Writing Δp = (p2 − (p)̅ 2)1/2 = (p2)1/2(1 − ((p)̅ 2/p2))1/2 and remembering that (p̅)2 ≤ p2, eq 9 can be approximated as

≈

⎛ mω ⎞1/4 ⎞ ⎛ mω ⎜ ⎟ exp(ip0 q) exp⎜ − (q − q0)2 ⎟ ⎠ ⎝ 2ℏ ⎝ πℏ ⎠

⎟

⎠

ℏ

(10)

This expression is numerically stable even if there is only one trajectory in a box (zero trajectories in a box is a noninteresting case). In our earlier work,24 we used a short time approximation where we obtained the characteristic distance from the initial distribution. In the present work we, however, update y0(q) at each time step, which is a key aspect as it is imperative for handling the tunneling at long times in the present application to the cubic potential. We illustrate the effect of using eq 7 by appling this effective quantum force in calculating reaction probabilities for a cubic potential and a symmetric Eckart and an asymmetric Eckart potential. In these applications, we use a Gaussian minimum uncertainty wave packet 7,8 to initiate the trajectories. This Gaussian corresponds to a harmonic oscillator with frequency ω and mass m. The initial wave function is 3017



■

evaluated. The mechanism is depicted in Figure 2. The decrease of the effective potential with time is illustrated in Figure 3. The reduction of the effective potential is by no means a new explanation for the quantum tunneling effect. The clear connection to the decrease in momentum of the remaining reactants we have however not encountered before. We think that this gives an interesting picture of the tunneling process. B. Symmetric Eckart Barrier. The symmetric Eckart barrier problem that we study here7 is defined by the potential V(q) = Vbsech2(aq), where the parameter a = 0.5 and V(qdiv) = Vb = 0.0364, where Vb is the height of the barrier and qdiv = 0. 600 boxes have been used in the range (−8,8). The average initial position of the wave packets is q0 = −7.0, just as in ref 7. We choose four different average energies, viz., E0 = 0.4Vb, E0 = 0.6Vb, E0 = 0.8Vb, and E0 = 1.0Vb, to allow comparison with the ETMD results of ref 7. We run 10 000 Monte Carlo sampled trajectories using a time step δt = 5. The resulting reaction probabilities are shown in Figure 4a. We find satisfactory agreement between the effective quantum force results, the ETMD results and the accurate quantum dynamical results. C. Asymmetric Eckart Barrier. As a third illustration we apply the effective quantum force to an asymmetric Eckart barrier V(q) = (A/(1 + exp(−aq)) + (B/(4cosh2(aq/2)), where A = V1 − V2, B = ((V1)1/2 + (V2)1/2)2, a = 1, V1 = 0.0364 and V2 = 0.0146.7 Ten thousand Monte Carlo sampled trajectories were run using a time step δt = 5. 600 boxes have been used in the range (−8,8). In Figure 4b, we compare the results obtained using the effective quantum force that we defined in the present work with results obtained using ETMD7 and accurate results.7 We see that the effective quantum force captures the entanglement between the trajectories quite well and produces quite satisfactory results. In conclusion, the effective quantum force that we constructed in this paper captures the essential “nonlocal” contributions to the dynamics. It can thus be used to approximately treat quantum dynamics. Further, the effective potential that we have constructed and the expression for the characteristic distance y0 offer nice explanations for the tunneling effect, particularly for the cubic potential employed here. We end by discussing the computational scaling of the method used here with respect to the number of degrees of freedom. A main cause for concern is the computation of the position-dependent momentum variance, Δp. It could be feared that an exponential increase of the number of required grid-boxes would occur. However, we argue that this is not so. To see this, first remember that the initial phase-space points are sampled by Monte Carlo. Hence, every sampled point has an appreciable density of neighbor points close to it. Since the Liouville theorem holds (due to a momentum-independent force), it follows that also later, when a point has been propagated, it will still have the same density of neighbors. If the Wigner density is not “squeeze-stretched” in an extreme fashion, the conservation of density should indeed make the computation of p2 possible also for larger dimensions. In this connection, a neighbor list approach similar to that employed in conventional molecular dynamics simulations can be applied to keep track of which particles are close to each other and this allows for the calculation of p2 at each position.

Letter

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■

REFERENCES

(1) Habershon, S.; Manolopoulos, D. E.; Markland, T. E.; Miller, T. F. Annu. Rev. Phys. Chem. 2013, 64, 387−413. (2) Craig, I. R.; Manolopoulos, D. E. J. Chem. Phys. 2004, 121, 3368−3373. (3) Voth, G. A. Adv. Chem. Phys. 2007, 93. (4) Miller, W. H. J. Chem. Phys. 2006, 125, 132305. (5) Wang, L.; Martens, C. C.; Zheng, Y. J. Chem. Phys. 2012, 137, 034113. (6) Donoso, A.; Zheng, Y.; Martens, C. C. J. Chem. Phys. 2003, 119, 5010. (7) Wang, A.; Zheng, Y.; Martens, C. C.; Weiyi, R. Phys. Chem. Chem. Phys. 2009, 11, 1588. (8) Donoso, A.; Martens, C. C. Phys. Rev. Lett. 2001, 87, 223202. (9) Lopez, H.; Martens, C. C.; Donoso, A. J. Chem. Phys. 2006, 125, 154111. (10) Hughes, K. H.; Wyatt, R. E. J. Chem. Phys. 2004, 120, 4089− 4097. (11) Garashchuk, S.; Mazzuca, J.; Volkov, M. V. J. Chem. Phys. 2012, 137, 074115. (12) Prezhdo, O. Theor. Chem. Acc. 2006, 116, 206−218. (13) Sun, X.; Miller, W. H. J. Chem. Phys. 1999, 110, 6635. (14) Poulsen, J. A.; Nyman, G. N.; Rossky, P. J. J. Chem. Phys. 2003, 119, 12179. (15) Heller, E. J. J. Chem. Phys. 1975, 62, 1544. (16) Wang, H.; Sun, X.; Miller, W. H. J. Chem. Phys. 1998, 108, 9726. (17) Wang, H.; Sun, X.; Miller, W. H. J. Chem. Phys. 1999, 110, 4828. (18) Sun, X.; Wang, H.; Miller, W. H. J. Chem. Phys. 1998, 109, 4190. (19) Shalashilin, D. V.; Child, M. S. J. Chem. Phys. 2008, 128, 054102. (20) Heller, E. J. J. Chem. Phys. 1991, 94, 2723−2729. (21) Tully, J. C. J. Chem. Phys. 1990, 93, 1061. (22) Wigner, E. Phys. Rev. 1932, 40, 749. (23) Husimi, K. Proc. Phys.-Math. Soc. Jpn. 1940, 22, 264. (24) Poulsen, J. A.; Li, H.; Nyman, G. J. Chem. Phys. 2009, 131, 024117. (25) Li, H.; Poulsen, J. A.; Nyman, G. J. Phys. Chem. A 2011, 115, 7338. (26) Dragoman, D. Phys. Scr. 2005, 72, 290−296. (27) The values of q0 have been chosen as in ref 5. For q0 = −0.3 this gives E0 = 1.26V0 rather than E0 = 1.25V0 as stated in ref 5.

3018


Tunneling Dynamics Using Classical-like Trajectories with an Effective

Recommend Documents