Semiclassical Surface Hopping HâK Propagator: Application to Two

6562

J. Phys. Chem. B 2001, 105, 6562-6569

Semiclassical Surface Hopping H-K Propagator: Application to Two-Dimensional, Two-Surface Problems† Guangcan Yang and Michael F. Herman* Department of Chemistry, Tulane UniVersity, New Orleans, Louisiana 70118 ReceiVed: December 15, 2000; In Final Form: April 2, 2001

A semiclassical surface hopping Herman-Kluk (HK) propagator is developed for multisurface, multidimensional nonadiabatic problems. Its application for two-surface model problems is presented for two-dimensional collinear systems. Both noncrossing and crossing cases are investigated. The calculation on the basis of Monte Carlo method shows that the general multisurface semiclassical HK propagator method is efficient and accurate.

I. Introduction It is quite common to separate the degrees of freedom of a system into a fast subsystem and a slow subsystem in many reactive or nonreactive scattering problems. Within the adiabatic approximation, the dynamics are effectively confined to one adiabatic state for the fast degrees of freedom. In these cases, semiclassical methods often provide an efficient method for the calculation of physically important quantities of the system such as wave functions, energies, transition probabilities, etc, for the slow degrees of freedom.1-3 In other cases, such as electrontransfer reactions, the slow degrees of freedom are not confined on a single adiabatic energy surface, and transitions between different adiabatic quantum states for the fast electronic degrees of freedom occur. Effective path methods, complex path methods, and surface hopping procedures are among the semiclassical methods for many surface problems.4-34 In the present paper, we focus on the surface hopping methods. Many surface hopping approaches have been developed with different criteria for hopping between potential surfaces.11-34 If the nonadiabatic transitions are localized about avoided crossing seams, then simple surface hopping procedures generally work well. However, there are cases in which the nonadiabatic interactions are not localized, such as vibrational energy relaxation in condensed matter.35,36 A general formulation of the surface hopping propagator has been previously presented and analyzed as a semiclassical solution of the time-dependent Schro¨dinger equation (TDSE).18 The semiclassical propagator was derived by the dual expansion in powers of p and in powers of the nonadiabatic coupling. This surface hopping propagator satisfies the Schro¨dinger equation through first order in p, which is the same order as the familiar single surface propagator, and through all orders of the nonadiabatic coupling. However, the theoretical analysis presented previously focused on the formal analysis in order to provide a basis for understanding surface hopping techniques. The present paper develops the theory further to provide a numerical method that is applicable and practical. The main numerical difficulty with the previous multisurface propagator comes the use of the van Vleck expression for the single surface semiclassical propagator in the surface hopping propagator. In a numerical calculation, it would be necessary †

Part of the special issue “Bruce Berne Festschrift”. * To whom correspondence should be addressed.

to find all trajectories originating from position x0 at t ) 0 and ending at positions xt at time t. An additional difficulty arises from the singularities in the van Vleck formula when the trajectories encounter caustics. In order overcome the numerical difficulties arising from the use of the van Vleck propagator, the HK initial value representation (IVR) form of the single surface propagator37-45 is used in this paper to develop a surface hopping theory. The paper is organized as follows: In section II.A, we briefly review the semiclassical surface hopping theory by dual expansions in p and nonadiabatic coupling. The HK propagator theory is briefly outlined in section II.B. Then, we develop the semiclassical HK surface hopping propagator by combining the dual expansion surface hopping theory with the HK propagator method in section II.C. To test the accuracy and effeciency of the method, numerical calculations are presented for two level collinear systems in section III. Both noncrossing and crossing cases are considered. Finally, in section IV, we present a summary of our results and some further discussion. II. Theory A. Semiclassical Surface Hoppping Propagator. Consider a system with the fast variables r and the slow variables R. For example, in many cases, the electronic coordinates are the fast variables and the positions of the nuclei are the slow variables. In the adiabatic approximation, the wave function for the slow variables is confined on a single potential surface Vi(R). The propagator for the slow variables Kijs(R1,R0,t) is given by the standard WKB expression1

K0ij(R0,R1,t) ) δijAii exp[iSii/p]

(1)

where Sii is the action for a classical path from R0 to R1 in time t under the influence of the potential Vi(R). The prefactor Aii is given by the well-known van Vleck expression

[|

|

∂2Sii Aii ) /(-2πip)d ∂R1 ∂R0

]

1/2

(2)

where d is the dimensionality of R. In some cases, the action of the gradient operator ∇s of slow variables on the wave function ψif of fast variables cannot be neglected; that is, the influence of the nonadiabatic coupling

10.1021/jp004509s CCC: $20.00 © 2001 American Chemical Society Published on Web 05/15/2001

Semiclassical Surface Hopping H-K Propagator

ηij )

∫ dr1 ψfj(r1|R1)*∇sψfi(r1|R1)

J. Phys. Chem. B, Vol. 105, No. 28, 2001 6563

(3)

is important. When this is the case, we should modify the semiclassical propagator K(R1,R0,t) for slow variables, eq 1, to ensure that the full propagator K(r1,R1,r0,R0,t) still satisfies the Schro¨dinger equation to order p. The nonadiabatic semiclassical propagator for slow variables K(R1,R0,t) can be written as the summation of the contribution of hopping trajectories18

Kij(R1,R0,t) ) δijAii exp(iSii/p) +

∫RR

∫RR

1

0

du τijATij exp(iSTij /p) + ∞

1

0

du FijARij exp(iSRij /p) +

K(n) ∑ ij (R1,R0,t) n)2

(4)

The second and third term on the rhs of eq 4 are the contribution from the trajectores undergoing a single hop. The second term gives the contributions from trajectories for which the sign of P‚ηij is the same before and after the hop, where P is the momentum for the slow variable coordinates. These are denoted as transmission or T-type trajectories. The third term provides the contribution from trajectories that have a change in the sign of P‚ηij at the hop. These are referred to as reflection or R-type trajectories. All of these trajectories travel on surface Vi from R0 to a hopping point Rh, hop to Vj, and then travel from Rh to R1. The u integration in eq 4 is over all possible hopping points for the given type of trajectory. The set of hopping points for all T-type or all R-type single hop trajectories that start at R0 on Vi and end on Vj at R1 at time t form a curve that goes from R0 to R1.18 The u integration is over this set of hopping points, and du is defined as the change in the component of the hopping point, Rh, parallel to ηij along this curve.18 The phase factor Sijγ, where γ ) T or R, is the classical action for the trajectory, and the prefactor Aijγ obeys the continuity equation and is continuous at the hopping point. Because the phase factor, S ) ∫(T - V) dt′, includes the contribution from the kinetic energy of slow degrees of freedom, phase coherence is handled exactly within the semiclassical approximation in the method. The prefactor is given by eq 2 multiplied by |Piη/Pjη|1/2, where Piη and Pjη are the magnitudes of component of the momentum on the initial potential surface Vi and final potential surface Vj in the direction of nonadiabatic coupling, respectively, at the hopping point.18 The amplitudes for the T-type and R-type transitions are given by18

τij ) Fij ) -

(Piη + Pjη) ηij 2Pjη (Piη - Pjη) ηij 2Pjη

∑ ∑β ∫R j ,j ...j

R1 0

1 2

(5)

du1 du2... dun Aijβ1j2..jn-1j

n-1

exp(iSijβ1j2..jn-1j/p)σij1σij2...σjn-1j

K(R,R0,t) )

(2πp1 ) ∫ da C(a,t)e

(6)

The β summation is over all sequences of T-type (no sign change in P‚ηij) and R-type hops (sign change in P‚ηij). The

d

iS(a,t)/p

g(R,at) g*(R0,a) (7)

where d is the number of degrees of freedom in the system, a ) (Ra,Pa) denotes a set of initial coordinates and momenta for a classical trajectory, da equals dRa dPa, and at ) (Rt,Pt) represents the phase space point obtained by classically evolving these variables to time t. S(a,t) is the classical action along this trajectory. The function

g(R,a) ) (γ/2π)d/2 exp[-γ(R - Ra)2 + (i/p)Pa‚(R - Ra)] (8) is a Gaussian wave packet with average position Ra and average momentum Pa with width parameter γ. The prefactor C(a,t) is definded by

{ [(

C(a,t) ) det

where ηij is the magnitude of ηij multiplied by the sign of Pi‚ ηij and Pi is the momentum just before the hop. K(n) ij (R1,R0,t) in eq 4 is the contribution to the propagator from trajectories undergoing n hops. Its general form can be written as18

K(n) ij (R1,R0,t) )

index jk designates the potential surface after the kth hop, and σab (where a ) jk-1 and b ) jk with the condition i ) j0 and j ) jn) is either τ or F, depending on whether the kth hop is a T-type or R-type event. S is the classical action of the trajectory undergoing n hops and A is the corresponding amplitude. They should be continuous at each hop. In order for A to be continuous, it is given by eq 2 multiplied a factor |Plη/Pkη|1/2 for each hop, where k and l denote the initial and final potential surfaces for the hop, respectively.18 Usually, the contribution of the transmission trajectories to the propagator is much bigger than that from reflection trajectories. In some surface hopping problems, we omit the reflection trajectories to simplify the calculation. B. HK Propagator. The surface hopping propagator described in the previous section includes the van Vleck propagator expression on a single potential surface. In practical calculations, the evaluation of the van Vleck propagator involves searches for the trajectories satisfying the desired boundary conditions at times 0 and t. Such searches can be very cumbersome. An additional difficulty comes from the singularities in the van Vleck formula, when the trajectories encounter caustics. To overcome the numerical difficulties, the semiclassical IVR techniques can be adopted, which generally express the propagation of the system as an integration over the initial positions, momenta, or phase space points for a set of classical trajectories.37-50 Among the IVR methods, the HK approach has been shown to yield high accuracy and relatively quick convergence.39-45 The HK propagator can be expressed as

∂Rt 1 ∂Pt ∂Rt 1 ∂Pt + - 2ipγ 2 ∂Pa ∂Ra ∂Pa 2ipγ ∂Ra

)]}

1/2

(9)

where the branch of the square root is chosen so that C(a,0) ) 1 and C(a,t) is a continuous function of time. Formally, γ must have a positive real part. In practical calculations, there is a certain degree of latitude in the choice of its value, but its optimization is both possible and desirable. The HK formula for the propagator can be used to determine the wave function that evolves semiclassically at time t from the initial state ψ0. The wave function at time t is

ψ(R,t) ) 1 d 2πp

( ) ∫ da dR C(a,t)e 0

iS(a,t)/p

g(R,at) g*(R0,a) ψ0(R0,0) (10)

At t ) 0, we have (Rt,Pt) ) (Ra,Pa), S(a,0) ) 0 and C(a,0) ) 1. In this case, eq 10 is simply the expansions of ψ0 in terms of

6564 J. Phys. Chem. B, Vol. 105, No. 28, 2001

Yang and Herman

the overcomplete Gaussian basis g(R,a). As t increases, the average position and momentum of each Gaussian wave packet travel along a classical trajectory (Rt,Pt) with the width parameter γ held fixed, and this travelling Gaussian is multiplied by the factor C(a,t) exp[iS(a,t)/p]〈g(R0,a)|ψ0(R0)〉, where 〈‚‚‚〉 indicates integration over R0, to obtain its contribution to ψ(R,t). C. HK Surface Hopping Propagator. For a problem involving multiple quantum states, the component of wave function for state j at time t can be written as

ψj(R,t) )

∑i ∫ dR0 Kij(R,R0,t) ψi(R0,0)

(11)

where

Kij(R,R0,t) )

∫ dr dr0 ψfj(r|R)*K(r,R,r0,R0,t) ψfi(r0,R0)

(12)

and K(r,R,r0,R0,t) is the Feynman propagator for the full system. The summation in eq 11 is over all initial states. For simplicity, we assume the initial wave function corresponds to a specific state i, and the summation drops out. A HK-IVR version of the semiclassical surface hopping propagator can be developed by following the same steps taken in the original derivation of the HK propagator,39 except that the multistate surface hopping generalization of the van Vleck semiclassical propagator is employed, rather than the single surface form. The nonhopping trajectories are treated as zeroth order hopping trajectories, and the zeroth order term is just the familiar single surface HK propagator. For the convenience of further derivation, we initially consider only the term from the first-order T-type trajectories. The first-order contribution of wave function in eq 4 is

∫ dR0 K1ij(R,R0,t) ψi(R0,0) ) ∫dR0duATij exp(iSTij/p)τijψi(R0,0)

ψ1j (R,t) )

(13)

where STij is the classical action of the T-type hopping trajectory and ATij is still given by eq 2 but multiplied by |Pjη/Piη|1/2 at each hopping point, where Piη and Pjη are the magnitudes of component of momenta along the direction of nonadiabatic coupling vector η before and after the hop. As is well-known, the basis of Gaussians, g(R,a), form an over complete set for any real γ > 0. The resolution of the identity in terms of Gaussians is

δ(R - R′) ) (2πp)-d

∫ da1 g(R,a1) g*(R′,a1)

ψ1j (R,t) )

∫ dR0 dRm dRn du da1 da2 g(R,a2)

g*(Rn,a2)ATij exp(iSTij /p)τijg(Rm,a1) g*(R0,a1) ψi(R0,0) (15) The integrand has the form of A exp[φ]ψi(R0,0), where the exponent φ contains iSij/p and the exponents from the Gaussian functions. We change the integration variables from P1 and P2 to new variables q1 and q2 defined by

0

∂Pjt+ ∂zia

h

)

∂zjt+ ∂Rit∂Rit- ∂zia

+

∂zjt+ ∂Pit∂Pit- ∂zia

(18)

where z denotes R either P, whereas t- ) th - δ and t+ ) th + δ as δ f 0+ and th is the time of the hop. The subscripts i and j are attached to quantities to denote that the quantity corresponds to a time before or after the hop, respectively. The derivatives ∂zjt/∂zjt+ and ∂zit-/∂zia are evaluated as in the single surface case. The only remaining derivative in eq 19 is ∂zjt+/ ∂zit-, which can be obtained as follows. Assume the hop from Vi to Vj occurs at time th in the interval between t and t + ∆t. Pi and Pj are the momenta on potential surfaces Vi and Vj at the hopping point. Then we have

Rt+∆t ) Rt +

(16)

then perform the R1, R2, Rm, and Rn integrations by the

d

∂zjt ∂Rjt+ ∂zjt ∂Pjt+ ∂zjt ) + ∂zia ∂Rjt+ ∂zia ∂Pjt+ ∂zia

P1 ) q1 - 2iγpR1 P2 ) q2 - 2iγpR2

Piη τ′ g(R,at) Cij(a,t) µ ij exp[iSa(t)/p]g*(R0,a) ψi(R0,0) (17)

(2πp1 ) ∫ da dR dt

where τ′ij ) (Pjη/Piη)1/2τij and τij is given by eq 5. The factor (Pjη/Piη)1/2 in τ′ij comes from the prefactor ATij in eq 13. In the surface hopping propagator discussed in section II.A, this factor multiplies the van Vleck form of the prefactor, eq 2. In this section, we move the (Pjη/Piη)1/2 factor to the transition amplitude, whereas the van Vleck expression for the prefactor is needed in order to obtain the HK prefactor Cij(a,t). The prefactor Cij(a,t) has the same form as eq 9, where the initial and final position and momenta are those for the hopping trajectories. This expression for Cij(a,t) can be obtained from the expression for the van Vleck prefactor, eq 2, the second derivative of the φ factor from the stationary phase integrations, and the Jacobian for the coordinates transformation from (q1,q2) to (Ra,Pa).40,41 The integration over hopping points has been replaced in eq 17 with an integration over all hopping times, th, between zero and t with the change of variable du ) dt Piη/µ. This change is simple in the IVR, because du is the change in the component of the hopping point Rh parallel to ηij, and the integration over hopping points is performed for a fixed value of the initial phase space point a. (Here we have assumed a single mass for simplicity. If this is not the case, mass weighted coordinates can be used.) This simplification of the integration over hopping trajectories is an advantage of the IVR form over the previously given form of the surface hopping propagator.18 We can evaluate the derivatives in Cij(a,0) at a hopping point by a limit procedure. We assume the trajectory starts on surface i at time ti and travels on this surface until hopping to surface j at time th, and then it travels on surface j until time tj. According to the chain rule, we have

(14)

where d is the dimensionality of the system and a1 ) (R1,P1). Following the derivation of the single surface HK-IVR39-41 and inserting (14) into (13) twice, we have

ψ1j (R,t) )

stationary phase method holding R0, q1, and q2 fixed. After completing the stationary phase integration and transforming the integration variables q1 and q2 to the initial position and momentum a ) (Ra,Pa), we have39

and

Pi Pj ∆th + (∆t - ∆th) + O[(∆t)2] µ µ

(19)



Pt+∆t ) Pt + (Pfη - Piη)eη - ∇Vi∆th - ∇Vf(∆t - ∆th) + 2

O[(∆t) ] (20) where ∆th ) th - t, eη is the unit vector parallel to the nonadiabatic coupling vector ηij with the direction defined such that Pi‚eη > 0, and Piη and Pjη are the components of the momentum in this direction before and after the hop, respectively. The value of Pjη is chosen so that energy is conserved during the hop. The second term of the rhs in eq 20 comes from the momentum change at the hop. Using eq 19, we have the derivative of the final position with respect to the initial position of the interval

Pi - Pf ∂th ∂Rt+∆t )I+ + O(∆t) ∂Rt µ ∂Rt

(21)

where the third term of the rhs will tend to zero as ∆t f 0. The changes in the hopping point Rh accompanying a change in the initial position or momentum for the hopping trajectory are always in a direction perpendicular to ηij in the semiclassical surface hopping expansion. (This condition is imposed so that the action satisfies the classical derivative expressions ∂S/∂Ri ) - Pi, ∂S/∂Rj ) -Pj, and ∂S/∂t ) -E.)18 Therefore, the time that it takes for the trajectory to go from Rt to the hopping point is given by

∆th )

µ e ‚(R - Rt) Piη η h

(22)

derivation can also be applied for the first-order R-type surface hopping term in eq 4, and this results in a similiar replacement. This argument is also applicable to the higher terms in eq 4. Collecting the terms of all orders, we obtain a full expression for the surface hopping HK propagator similiar to eq 4 but with the replacement of the van Vleck expression with the HK formula for hopping trajectories. In a practical application, the th integration in eq 17 is evaluated numerically by considering hops in successive time intervals. If the nonadiabatic coupling is large in time intervals, as could be the case near an avoided crossing, then multiple hop trajectories can made significant contributions, even in a small time interval. However, the integrand in eq 17 can be approximated further to account for multiple hops in the small time intervals. We develop this approximation for a two states nonadiabatic coupling case. Consider a short time interval from ti to tj. The first-order contribution to the propagator from trajectories that start at ti and phase space point a and hop in this time interval is

k1ij(ti,tf) )

(23)

Inserting the expression into eq 21, we arrive at

∂Rt+ ∂Rt-

) lim

∆tf0

eη ∂Rt+∆t ) I - (Pi - Pf) ∂Rt Piη

(24)

By a similiar derivation, we can obtain the other derivatives in eq 18

∂Rt+ ∂Pt∂Pt+ ∂Rt-

i

Piη C(a,t)τijeiSij/p µ

(26)

k1ij(ti,tf) ) C(a,t)τijeiSij(m)/p

∫t t dt f

i

Piη i∆Sij/p τ′ e µ ij

(27)

where ∆Sij is the difference between the classical actions of a trajectory that hops at th and one that hops at tm and Sij(m) is the action for the trajectory that hops at tm. Only the trajectories with an odd number of hops contribute to the propagator with different initial and final potential surfaces, if only two surfaces make important contributions in the small interval. The next term of the transition propagator comes from trajectories with three hops. For the three hops trajectory, we have

k3ij(ti,tf) ) )0

∫t t dt1 f

i

) -[∇Vi - ∇Vj] eη

{

∂th ∂eη + (Pjη - Piη) + ∂Rt ∂Rt

[

Piη - Pjη ∂eη ∂th Pi (e ‚∇Vi) Pjη ∂Rt ∂Rt η

]}

eη ) I + (Piη - Pjη) eη ∂PtPfη

Piη Pjη Piη dt2 dt3 Cij(a,t)τ′ij(t1)τ′ji(t2)τ′ij(t3)eiSij/p (28) µ µ µ

with the ordered time t1< t2< t3. We can again pull out the prefactor and exponent at midpoint from the integrand. Then eq 28 can be approximated as

-

∂Rh µ ∇(Vj - Vi) Pjη ∂Rt ∂Pt+

f

The trajectory can hop at any time in the interval. We approximate the prefactor as constant with its value taken to be that for the trajectory that hops at the time tm at the middle of the small time intervel. In this approximation the integration becomes

where only the lowest order terms in ∆th have been kept in obtaining eq 22. This gives

∂th µ ) - eη ∂Rt Piη

∫t t dt

(25)

The prefactor in eq 4 for a hopping trajectory can be calculated from eq 18 and these relations. Comparing the first-order T-type surface hopping HK propagator in eq 17 with its counterpart in eq 4, we can see that the above derivation just results in the replacement of the van Vleck expression with the corresponding HK formula. The same

k3ij(ti,tf) ≈ -

(∫

1 C(a,t)eiSij(m) 3!

tf

ti

dt


)

3

(29)

where we have used the relation τ′ji ) -τ′ij. In eq 29, the factor 1/3! is introduced to compensate for the change from the time ordered integrations in eq 28 to integrations over all time orderings. The main approximation in eq 30 involves the action difference ∆Sij, which is defined as in the single hop calculation in eq 27. The correct ∆S for a three hop trajectory is not the same as the sum of ∆S for three single hop trajectory. However, when the interval is small, we expect the error from the phase factor to be quite small and the approximation to be reasonable in practical calculations.


Yang and Herman

Defining

Tij )

∫ dth


(30)

the contribution to the full i f j propagator involving any number of hops in the short time interval is approximated as

(

kij(ti,tf) ≈ C(a,t)eiSij(m) Tij -

1 3 1 5 T + T + ... 3! ij 5! ij

)

) C(a,t)eiSij(m) sin(Tij)

(31)

On the other hand, the corresponding contribution to the propagator on the same surface includes only even terms and can be approximated as

kij(ti,tf) ≈ C(a,t)e

(

iSii(m)

1 1 1 - T2ij + T4ij + ... 2 4!

) (32)

This numerical approximation is particularly useful near avoided crossings where the nonadiabatic couplings can be become very large and multiple hop terms can make significant contributions within a time interval of fairly modest width. The total contribution from each trajectory has the g(R,at) and 〈g(R,a0)|ψ0〉 factors as well as the prefactor C(a,t) and phase factor eiS/p, and there is a cos(Tij) or sin(Tij) for each small time interval. The trajectory, for which the contribution is evaluated, either stays on the same surface throughout a time interval or hops at the time in the middle of that interval for each small time interval between 0 and t. In each interval, a normalized hopping probability is defined as sin|Tij|/(sin|Tij | + cos|Tij |), and the probability of staying on the same surface is cos|Tij |/(sin|Tij | + cos|Tij |). A hop-or-not decision is made by comparing a random number between 0 and 1 with the hopping probability. An extra factor of ((sin|Tij | + cos|Tij |) is added from each interval to the weight for the trajectory to compensate for the use of these normalized branching probabilities. III. Numerical Calculations In this section, we present numerical calculations for a model collinear atom-diatom system to illustrate and test some of the issues discussed in the previous sections. The particles in the system have the masses of a helium and a hydrogen molecule with the atomic masses given by MH ) 1836.2 a.u. and MHe ) 4MH (although the model system does not correspond to an actual H2/He system). The model is defined by the following diabatic Hamiltonian matrix elements for the fast variable subsystem

1 H f22(R1,R2) ) k2(R1 - R2e)2 + A2e-B2R2 + ∆ 2

{

[(

)

R1e + R2e 2 + Be(R2 - Ce) 2

where the “-” sign is used for j ) 1 and “+” sign for j ) 2. The adiabatic wave functions are expressed in terms of the diabatic wave functions as

()(

)( )

(35)

]

(36)

where the angle θ is given by

[

2H f12 1 -1 θ ) tan 2 H f11 - H f22

The nonadiabatic coupling can be obtained from eqs 3 and 35 as η12 ) ∇θ. The initial wave packet is assumed to be a Gaussian wave packet with average position R0 and average momentum P0. The initial wave packet is entirely on the lower potential surface V1. Thus, we have

ψ1(0) ) (R/2π)d/4 exp[-R(R - R0)2 + (i/p)P0‚(R - R0)] ψ2(0) ) 0

(37)

where the number of degrees of freedom d is 2 and the width parameter R ) 0.02. In present case, we take R0 ) (1.0, 5.0), and the initial wave packet has an average vibrational kinetic energy corresponding to the vibrational ground state and an average incoming kinetic energy Tin ) 1.5, so that P0 ) (6.9, -22.1). Fully quantum time dependent wave packet calculations are compared with semiclassical calculations. The quantum calculations are performed in the diabatic representation, and wave functions are transformed to the adiabatic representation before evaluating the time dependent transition probability. The quantum calculations are performed by a generalized splitoperator method.51-53 The semiclassical time-dependent evolution of the Gaussian wave packet can be evaluated in terms of the initial wave function and the semiclassical nonadiabatic propagator

∫K11(R1,R2,t) ψ0(R1,R2) dR ψ2(R1,R2,t) ) ∫K12(R1,R2,t) ψ0(R1,R2) dR ψ1(R1,R2,t) )

1 H f11(R1,R2) ) k1(R1 - R1e)2 + A1e-B1R2 2

Aae 1 - tanh Ae R1 -

1 Vaj ) [(H f11 + H f22)2 - x(H f11 - H f22)2 + 4(H f12)2] (34) 2

d ψf1 cos θ sin θ φ1 ) -sin θ cos θ φd2 ψf2

) C(a,t)eiSii(m) cos(Tij)

H f12(R1,R2) )

between the diabatic surfaces and is set to ∆ ) 0.5. With this choice, the two diabatic surfaces H f11 and H f22 do not cross. The coordinate R1 is the diatomic bond length, and R2 is the distance between the incoming atom and the atom in the diatomic closest to it. The equilibrium position parameters R1e and R2e for coordinate one are set to 1.1 and 1.2, respectively. The adiabatic surfaces obtained by diagonalization of the diabatic Hamiltonian matrix are

]}

(38)

where the initial wave packet is given in eq 37 and nonadiabatic HK propagator is presented in section II.C. The integrals over a ) (Ra,Pa) are evaluated by the Monte Carlo method, sampling the initial phase space points from the distribution

(33)

where the parameters are set to k1 ) 10.0, k2 ) 30.0, A1 ) 1.0, B1 ) 1.5, A2 ) 2.0, B2 ) 2.0, Ae ) 1.0, Aae ) 0.1, Be ) 2.0, and Ce ) 2.0. The parameter ∆ is the asymptotic separation

F(a) ) (2πσRσP)-2 exp[-(Ra - R0)2/2σR2 (Pa - P0)2/2σP2] (39) where (R0,P0) is the average phase space position of the initial



Figure 1. Noncrossing two surface collinear system. The transition probability P(t) from the lower potential surface to the upper potential surface is plotted versus time. The quantum result (solid line) is based on the generalized split-operator method. The semiclassical results are presented for calculations employing 50 000 Gaussians. The dash line corresponds to the calculation including multihop contributions in each time step, whereas the short dash line corresponds to the single hop case with the same other parameters.

wave packet, σR ) (γ + R)/2γR, and σP ) 2(γ + The sampling distribution is proportional to |〈g(R,a)|ψ0〉|. The calculation of classical trajectory is completed by a RungeKutta numerical integration procedure, and each Gaussian wave packet travels along a classical trajectory with a multiplying factor of C(a,t) exp[iS(a,t)/p]. By the Monte Carlo procedure, the final wave function is given by the following summation 2

ψ1(R,t) ) ψ2(R,t) )

R)/p2.

2

A

∑i gi(R,a)eiφCi(a,t)eiS (a,t)/p i

N A

∑j gj(R,a)eiφCj(a,t)eiS (a,t)/p

N

j

(40)

where subscripts i and j denote the trajectories ending on the surface 1 and 2, respectively. The factor A and eiφ in eq 40 arise because the overlap between the initial wave packet and the Gaussian g(R,a) is proportional to, but not the same, as the Monte Carlo sampling function, i.e., 〈g|ψ0〉 ) AF(a)eiφ. N is the total number of sampled trajectories. The final transition probabilities can be obtained by the double summation over the trajectories ending on a particular surface using the fact that the integration ∫g*(R,a1) g(R,a2) dR is easily performed analytically. A. Numerical Results. Given the diabatic Hamiltonian matrix elements in eq 33, H f11(R) and H f22(R) do not cross. Furthermore, the adiabatic coupling in this case is not highly localized. As we mentioned in the Introduction, some surface hopping methods have difficulties with this kind of delocalized nonadiabatic coupling problem. In the other hand, the application of the HK surface hopping propagator presented in this work is not constrained by the localized or delocalized nature of nonadiabatic coupling. The semiclassical transition probabilities calculated by the HK surface hopping propagator are compared with the quantum transition probabilities in Figure 1. Only T-type transitions are included in these calculations. The probability of the transition is evaluated in each time step, and a decision whether the trajectory hops or not is made along the trajectories, as described in section II.C. This hop-or-not-hop time step is chosen to be 0.6 a.u. (A smaller time step of 0.06 a.u. is used for the Runge-Kutta integration of the trajectories.) To estimate the importance of multiple hops within an interval,

Figure 2. Classical trajectory of the center of initial Gaussian wave packet ψ0 on the lower surface in the noncrossing surface problem. One of its branches (dash line) undergoes a hop to the upper surface at t ) 35 a.u., and the other one (short dash line) hops to the upper surface at t ) 140 a.u. The time interval between two adjacent × marks is 20 a.u.

results are also shown in Figure 1 for calculations which account for only single hops in each interval by replacing sin(Tij) in eq 31 with Tij and cos(Tij) in eq 32 with unity. We can see that the single hop calculation is satisfactory, but it can be improved further by including multiple hops in each interval. In the calculation, we have taken the width parameter of Gaussian functions to be R ) γ ) 0.02. The relative sampling error of semiclassical calculation is 5% as estimated from 10 runs using different random number sequences. The classical trajectory that starts with the average position and momentum of the initial wave function on the lower potential surface is plotted in Figure 2, whereas two branching trajectories which hop to the upper surface at t ) 35 and t ) 140 respectively are plotted in the same figure. We can see that the hopping trajectories deviate their original trajectory dramatically. The wave packet travels inward to the interaction area and oscillates in the R1 direction. At the early times, t < 25, the wave packet has largely not yet reached the interaction region and almost no transition occurs. From t ) 25 to t ) 45, the wave packet travels through the interaction area and part of the wave packet hops to the upper surface. After passing through the interaction area, the wave packets on the lower and upper surfaces continue their evolution until they are reflected back by the potential barriers. They continue to propagate, largely on a single surface, until they enter the interaction region again. In this small R2 propagation, almost no interaction occurs between the two surfaces and the transition probabilities remain almost constant. At about t ) 130, the wave packet enters the interaction region again. Much more of the wave packet on the lower surface jumps to the upper surface on this second pass through the interaction region, and the transition probability to the upper surface increases rapidly from less than 0.01 to about 0.05. Then the wave packet leaves the interaction region, and the transition probabilities remain almost unchanged. In this numerical example, the semiclassical calculation employing 50 000 Gaussians takes about 12 h 30 min CPU time on a IBM RS6000 43P-260 workstation, and the corresponding quantum calculation using 1024 × 1024 FFT grid takes about 10 h 40 min on the same machine. Because computational effort required by the quantum method is strongly dimension-dependent, we can expect the present semiclassical method will be the numerically more preferable procedure in higher dimensional calculations.


Yang and Herman

Figure 3. Surface crossing collinear system. The transition probability P(t) from the lower potential surface to the upper potential surface is plotted versus time. The quantum result (solid line) is based on the generalized split-operator method. The semiclassical results are presented for calculation employing 20 000 Gaussians. The dash line corresponds to the calculation including multihop contributions in each time step, whereas the short dash line corresponds to the single hop case with the same other parameters.

If we modify H f11 in eq 33 to the following form

1 H f11(R1,R2) ) k1(R1 - R1e)2 + A1e-B1R2+3.0 2

(41)

then the surfaces H f11 and H f22 cross near R2 ) 1.0. In this case, the nonadiabatic coupling is localized near the avoided crossing curve. The other parameters are the same as those in the noncrossing case. The semiclassical transition probabilities are compared with the result of the full quantum calculation for this crossing problem in Figure 3. The relative sampling error of semiclassical calculation is 8%, as estimated from eight runs using different random number sequences. The transition probability displays pronounced oscillations in this surface crossing example because of the interference between the contributions to the transition amplitudes from different paths. IV. Discussion The semiclassical surface hopping HK propagator developed in this paper generalizes our previous work.18 The present surface hopping HK propagator can be more easily used in practical calculations than the previous based on the van Vleck expression.18 The present formulation is an IVR technique, which yields converged results from comparatively small samples of trajectories. Because it is an IVR, it avoids the difficult root searches for hopping trajectories obeying double ended boundary conditions. This feature results in a considerable numerical simplification. In the previous van Vleck type approach, the integration over hopping points, even in the firstorder term, is a difficult task. For fixed values of the initial point Ri, final point Rf, and time t for the trajectory, all hopping trajectories which hop at a time between zero and t have to be included. The hopping points for this set of trajectories form a curve running from Ri to Rf. In the IVR form presented here, the hopping point integration is over all hopping trajectories with initial phase space points a. Because these trajectories all have the same initial position and momentum, the integration over hopping points, Rh is easily transformed into an integration over hopping times, th, with the simple conversion of the integration variable du ) dt Piη/µ, resulting in an equation which can be implemented numerically in a straightforward manner.

Our model calculations consider a two-dimensional collinear system, but the method can be applied in the same mannner for higher dimensional problems. We summarize our numerical Monte Carlo implementation of the HK surface hopping propagator as follows: 1. The initial phase points are sampled by a Monte Carlo algorithm54 so as to obtain a distribution which is proportional to the overlap of the initial wave packet and Gaussian functions |〈g(R,a)|ψ0〉|. The classical trajectories are then run by integrating Hamiltion’s motion equations, making a hop-or-not decision in each time step. 2. The prefactor C(a,t) is calculated numerically according to eq 9. Because the trajectories involve hops between different potential surfaces, more care must be taken in the calculation of the derivatives at hopping points. This is discussed in section II.C. 3. The transition probability can be obtained by the integration on the squared amplitude of wave functions. Because the wave function consists of a summation of Gaussian functions, the integration can be achieved analytically. In summary, an analysis of the semiclassical surface hopping HK propagator for multisurface nonadiabatic problems in many dimensions has been presented. We emphasize that it provides not only a useful formal framework for the analysis of surface hopping problems but also a feasible computational approach for many dimensional problems. Acknowledgment. This research is supported by NSF Grant CHE-9816040. References and Notes (1) Miller, W. M. AdV. Chem. Phys. 1974, 25, 69; 1975, 30, 77. (2) Heller, E. J. Chem. Phys. 1975, 62, 1544; 1981, 75, 2923; 1991, 94, 2723. (3) Herman, M. F. Annu. ReV. Phys. Chem. 1994, 45, 83. (4) Delos, J. B.; Thorson, W. R. Phys. ReV. A 1972, 6, 720; 728. (5) Miller, W. M.; McCurdy, C. W. J. Chem. Phys. 1978, 69, 5163. (6) McCurdy, C. W.; Meyer, H. D.; Miller, W. M. J. Chem. Phys. 1979, 70, 3177. (7) Meyer, H. D.; Miller, W. M. J. Chem. Phys. 1979, 70, 3214; 71, 2156; 1980, 72, 2272. (8) McCurdy, C. W.; Miller, W. M. J. Chem. Phys. 1980, 73, 3191. (9) Micha, D. A. J. Chem. Phys. 1983, 78, 7138. (10) Miller, W. M.; George, T. F. J. Chem. Phys. 1972, 56, 5637. (11) Preston, R. K.; Tully, J. C. J. Chem. Phys. 1971, 54, 4297. (12) Tully, J. C.; Preston, R. K. J. Chem. Phys. 1971, 55, 562. (13) Pechukas, P. Phys. ReV. 1969, 181, 166; 174. (14) Laing, J. R.; Freed, K. F. Chem. Phys. 1977, 19, 91. (15) Herman, M. F. J. Chem. Phys. 1982, 76, 2949. (16) Herman, M. F. J. Chem. Phys. 1983, 79, 2771; 1984, 81, 754, 764. (17) Herman, M. F. J. Chem. Phys. 1985, 82, 3666. (18) Herman, M. F. J. Chem. Phys. 1995, 103, 8081. (19) Tully, J. C. J. Chem. Phys. 1990, 93, 1061. (20) Hammes-Schiffer, S.; Tully, J. C. J. Chem. Phys. 1994, 101, 4657. (21) Yang, J. Y.; Hammes-Schiffer, S. J. Chem. Phys. 1997, 107, 5727. (22) Webster, F.; Rossky, P. J.; Friesner, R. A. Comput. Phys. Commun. 1991, 63, 494. (23) Webster, F.; Schniker, J.; Friedrich, M. S.; Friesner, R. A.; Rossky, P. J. Phys. ReV. Lett. 1992, 66, 3172. (24) Webster, F.; Tang, E. T.; Rossky, P. J.; Friesner, R. A. J. Chem. Phys. 1994, 100, 4835. (25) Bittner, E.; Rossky, P. J. J. Chem. Phys. 1995, 103, 8130. (26) Prezhdo, O. V.; Rossky, P. J. J. Chem. Phys. 1997, 107, 5863. (27) Space, B.; Coker, D. F. J. Chem. Phys. 1991, 94, 1976. (28) Batista, V. S.; Coker, D. F. J. Chem. Phys. 1997, 106, 6923. (29) Blais, N. C.; Truhlar, D. G. J. Chem. Phys. 1983, 79, 1334. (30) Blais, N. C.; Truhlar, D. G.; Mead, C. A. J. Chem. Phys. 1988, 89, 6204. (31) Topaler, M. S.; Hack, M. D.; Allison, T. C.; Lui, Y. P.; Mielke, S.; Schwenke, D. W.; Truhlar, D. G. J. Chem. Phys. 1997, 106, 8699. (32) Sun, X.; Miller, W. H. J. Chem. Phys. 1997, 106, 6346. (33) Muller, U.; Stock, G. J. Chem. Phys. 1997, 107, 6230.

Semiclassical Surface Hopping H-K Propagator (34) 4713. (35) (36) (37) (38) (39) (40) 326. (41) (42) (43) (44) (45) (46)

Kohen, D.; Stillinger, F. H.; Tully, J. C. J. Chem. Phys. 1998, 109, Herman, M. F.; Arce, J. C. Chem. Phys. 1994, 183, 335. Arce, J. C.; Herman, M. F. J. Chem. Phys. 1994, 101, 7520. Miller, W. M. J. Chem. Phys. 1970, 53, 3578. Miller, W. M. J. Chem. Phys. 1991, 95, 9428. Herman, M. F.; Kluk, E. Chem. Phys. 1984, 91, 27. Kluk, E.; Herman, M. F.; Davis, H. L. J. Chem. Phys. 1986, 84, Herman, M. F. J. Chem. Phys. 1986, 85, 2049. Kay, K. G. J. Chem. Phys. 1994, 100, 4377, 4432. Kay, K. G. J. Chem. Phys. 1994, 101, 2250. Kay, K. G. J. Chem. Phys. 1997, 107, 2313. Elran, Y.; Kay, K. G. J. Chem. Phys. 1999, 110, 3653. Sepulveda, M. A.; Heller, E. J. J. Chem. Phys. 1994, 101, 8004.

J. Phys. Chem. B, Vol. 105, No. 28, 2001 6569 (47) Spath, B. W.; Miller, W. H. J. Chem. Phys. 1996, 104, 96; Chem. Phys. Lett. 1996, 262, 486. (48) Keshavamurthy, S.; Miller, W. H. Chem. Phys. Lett. 1994, 218, 189. (49) Compolieti, G.; Brumer, P. Phys. ReV. A 1994, 50, 997; 1996, 53, 2958; J. Chem. Phys. 1997, 107, 791. (50) Provost, D.; Brumer, P. Phys. ReV. Lett. 1995, 74, 250. (51) Feit, M. D.; Fleck, J. A.; Steige, A. J. Comput. Phys. 1982, 47, 412. (52) Alvarellos, J.; Metiu, H. J. Chem. Phys. 1988, 88, 4957. (53) Balakrishnan, N.; Kalyanaraman, C.; Sathyamurthy, N. Phys. Rep. 1997, 280, 79. (54) Kalos, M. H.; Whitlock, P. A. Monte Carlo Methods, Volume 1: Basics; Wiley-Interscience: New York, 1986.

Semiclassical Surface Hopping HâK Propagator: Application to Two

Recommend Documents

Semiclassical Surface Hopping HâK Propagator: Application to Two