Materials Chemistry Research Dept. BL0111780 Research Paper

Frank H. Stillinger

Model Studies of Non­Adiabatic Dynamics

Daniela Kohen (a)*, Frank H. Stillinger (a);(b), and John C. Tully (c)

(a) Lucent Technologies, Bell Laboratories, 700 Mountain Avenue, Murray Hill, NJ 07974,
(b) Princeton Materials Institute, Princeton University, NJ 08544,
(c) Chemistry Department, Yale University, New Haven, CT 06511.
(*) Current address: Department of Chemistry, University of California, Irvine, Irvine, CA, 92697

January 1998

Abstract: Mixed quantum­classical methods are applied to an increasingly challenging series of model problems, and their accuracy examined. The models involve one light and one heavy degree of freedom, and exhibit substantial nonadiabatic behavior. In all of the models the coupling between the light and heavy particles is linear (harmonic). In addition, different external potentials are applied to the heavy particle only. The energies of the light particle quantum states, as a function of the position of the heavy particle, define a sequence of ground and excited Born­Oppenheimer potential energy curves. Because the light particle experiences a purely harmonic potential, the potential energy curves are parallel and equally spaced for all of the models. In addition, the nonadiabatic couplings among potential energy curves persists for all times due to the non­vanishing linear coupling between light and heavy particles. The model problems were used to test two strategies for carrying out mixed quantum­classical dynamics in systems involving nonadiabatic transitions: mean field and surface hopping. The model calculations reported here suggest that, in cases where linear couplings dominate, the mean field mixed quantum­classical method displays useful accuracy and is robust to the introduction of anharmonic heavy­particle interactions. The model calculations also reveal special situations in which the surface hopping approximation is inappropriate.

• Paper [Text Approximation]
• Paper [Postscript (235KB)]
• Figures
  • Fig. 1a [.jpg (19KB)]
  • Fig. 1b [.jpg (17KB)]
  • Fig. 1c [.jpg (31KB)]
  • Fig. 1d [.jpg (22KB)]
  • Fig. 1e [.jpg (32KB)]
  • Fig. 2a [.jpg (20KB)]
  • Fig. 2b [.jpg (16KB)]
  • Fig. 2c [.jpg (26KB)]
  • Fig. 3a [.jpg (22KB)]
  • Fig. 3b [.jpg (19KB)]
  • Fig. 3c [.jpg (19KB)]
  • Fig. 4a [.jpg (20KB)]
  • Fig. 4b [.jpg (13KB)]
  • Fig. 4c [.jpg (20KB)]
  • Fig. 4d [.jpg (5KB)]
  • Fig. 4e [.jpg (5KB)]
  • Fig. 5a [.jpg (19KB)]
  • Fig. 5b [.jpg (18KB)]
  • Fig. 5ca [.jpg (28KB)]
  • Fig. 5cb [.jpg (37KB)]
  • Fig. 5d [.jpg (6KB)]
  • Fig. 5e [.jpg (5KB)]

 
Daniela Kohen (a)*, Frank H. Stillinger (a);(b), and John C. Tully (c) 
(a) Lucent Technologies, Bell Laboratories, 
700 Mountain Avenue, Murray Hill, NJ 07974, 
(b) Princeton Materials Institute, Princeton University, NJ 08544, 
(c) Chemistry Department, Yale University, New Haven, CT 06511.
(*)Current address: Department of Chemistry, University of California,
   Irvine, Irvine, CA, 92697
January 1998
  
Abstract: 
Mixed quantum­classical methods are applied to an increasingly 
challenging series of model problems, and their accuracy examined. 
The models involve one light and one heavy degree of freedom, and 
exhibit substantial nonadiabatic behavior. In all of the models the 
coupling between the light and heavy particles is linear (harmonic). 
In addition, different external potentials are applied to the heavy par­ 
ticle only. The energies of the light particle quantum states, as a func­ 
tion of the position of the heavy particle, define a sequence of ground 
and excited Born­Oppenheimer potential energy curves. Because the 
light particle experiences a purely harmonic potential, the potential 
energy curves are parallel and equally spaced for all of the models. In 
addition, the nonadiabatic couplings among potential energy curves 
persists for all times due to the non­vanishing linear coupling be­ 
tween light and heavy particles. The model problems were used to 
test two strategies for carrying out mixed quantum­classical dynamics 
in systems involving nonadiabatic transitions: mean field and sur­ 
face hopping. The model calculations reported here suggest that, in 
cases where linear couplings dominate, the mean field mixed quantum­ 
classical method displays useful accuracy and is robust to the intro­ 
duction of anharmonic heavy­particle interactions. The model cal­ 
culations also reveal special situations in which the surface hopping 
approximation is inappropriate. 
2 

1 Introduction 
Many dynamical phenomena in physics and chemistry involve coupled mo­ 
tions of particles with widely differing masses. Electronic and nuclear degrees 
of freedom respectively, in single atoms and molecules, provide a vivid illus­ 
tration, which the Born­Oppenheimer separation scheme [1] exploits to good 
advantage. The transfer of light protons along hydrogen bonds between pairs 
of heavy electronegative atoms (e. g. N, O, F, Cl) constitutes another class of 
examples [2]. The interaction of liquid He 3 or He 4 with high­atomic­weight 
solids, both crystalline and amorphous, present further illustrations [3]. 
The presence of extreme mass ratios invites the use of a mixed quantum 
and classical description [4]. This is particularly appropriate for describing 
the (virtually classical) motions of electron­bearing heavy atoms or molecules 
that form nonreactive insulating media. That is precisely the basis of the 
molecular dynamics simulation method, for which the Born­Oppenheimer 
ground­electronic­state energy surface provides the potential energy function 
for numerical integration of Newton's equations [5, 6]. Assuming that elec­ 
tronic excitation remains absent, methods are available when required to 
generate quantum corrections to the strict classical dynamical description of 
3 

the heavier particles [7]­[10]. On the other hand, there are methods that can 
accurately treat the electron dynamics even if electronic excitations occur but 
provided that their impact on the motion of the heavier atoms is negligible 
[11]­[14]. 
Difficulties begin to arise in those cases where electronic (or other light 
particle) excitations out of the ground state become significant and these, in 
turn, influence the behavior of the heavier particles. One encounters this 
situation, for example, in photochemical reactions involving nonradiative 
transitions, and in oxidation­reduction or proton transfer reactions in liquid 
solutions. Accurate description of structure and dynamics in liquid metals 
and reactions at metal surfaces are also problematic, since a continuum of 
excited conduction electrons must be involved. 
The intrinsic power of molecular dynamics computer simulations to il­ 
luminate physical and chemical phenomena has lent considerable weight to 
extending its traditional classical version as described above, to cover cases 
with electronic excitation. Indeed several approximation schemes for real­ 
izing this goal treating the quantum and classical degrees of freedom self­ 
consistently have been proposed (for reviews, see [15]­[16]). These include 
both `mean field' [17]­[33] and `surface hopping' [34]­[42] approaches. But in 
4 

spite of extended effort no universally acceptable approximation strategy has 
yet emerged. The present paper reports results of a modest study intended 
to contribute some further understanding to this active research area. In par­ 
ticular, our efforts could be regarded as an attempt to understand to what 
degree the use of some of these methods is appropriate when the evolution 
of a heavy plus light particle system is studied using the particle coordinates 
directly (rather than some collective coordinates). For simple model prob­ 
lems it would be more appropriate to base a quantum­classical separation 
on normal mode coordinates in order to minimize the effect of the quantum 
coordinates on the classical ones, and vice versa. In realistic many particle 
simulations a normal mode separation is not practical, however, and the bare 
particle must generally be used. This introduces coupling between classical 
and quantum coordinates, even in harmonic systems, and this coupling can 
be problematical for mixed quantum­classical dynamics. The simplest model 
potentials that represent this situation are those of linearly coupled harmonic 
and anharmonic oscillators. In this context nonadiabatic transitions may oc­ 
cur between the adiabatic (Born­Oppenheimer) light­mass states simply be­ 
cause the bare particle coordinates do not correspond to the `normal modes' 
of the system. The following section 2 argues in favor of studying simple 
5 

tractable models of nonadiabatic effects, and introduces three such models, 
each involving a coupled pair of particles moving in one dimension under the 
influence of an external potential. Section 3 provides computational details 
about the respective dynamical investigations for each of these three models. 
Section 4 presents our numerical results for the three. The final section 5 dis­ 
cusses the implications of our results for future theoretical and simulational 
activity. 
2 Simple Models 
Each of the three cases to be examined involves a pair of spinless particles 
confined to move in one dimension only. The lighter particle of the pair 
(mass m, position x, momentum p) interacts with the heavier particle (mass 
M , position X, momentum P ) via potential v(x; X). At the same time the 
heavier particle experiences an external potential V (X). The two­particle 
Hamiltonian consequently has the form 
H = p 2 
2m + P 2 
2M + v(x; X) + V (X) (1) 
The three cases to be analyzed are distinguished by the corresponding as­ 
signments of v and V , to be specified below. 
6 

Reliance on one­dimensional models to illuminate three­dimensional phe­ 
nomena obviously must be viewed with caution. Nevertheless this strategy 
offers some advantages. In particular, some one­dimensional models can be 
solved analytically in closed form (see model I below), while others that are 
not fully solvable may still permit a more thorough analysis than do their 
more complicated three­dimensions counterparts. Furthermore, various con­ 
tributing effects in one dimension can usually be easily isolated from model 
to model, thereby simplifying the task of interpretation. We view the study 
of nonadiabatic effects in the one­dimensional context as a helpful precursor 
to the development of theoretical techniques that are broadly applicable to 
three­dimensional systems. 
Let / n (x; X) and E n be the exact eigenfunctions and eigenvalues for 
the time­independent wave equation corresponding to the Hamiltonian H in 
operator form: 
H/ n (x; X) = E n / n (x; X) (2) 
subject to given boundary conditions. The general time­dependent solution 
to this quantum­mechanical problem has the form 
\Psi(x; X; t) = 
X 
n 
A n exp[ 
with constant A n selected to satisfy initial conditions. One must keep in mind 
that this last equation must be modified to include an appropriate integral 
over continuum states when they exist for a case under consideration. 
The Born­Oppenheimer (BO) approximation replaces the full quantum 
problem eq. 2 by a sequential pair of single­particle problems. The first ex­ 
amines the light particle motion in the presence of a stationary heavy particle 
(i. e. fixed at position X) 
H (0) (xjX)OE j (xjX) = ffl (0) 
j (X)OE j (xjX) (4) 
where H (0) (xjX) simply contains X as a parameter: 
H (0) (xjX) = p 2 
2m + v(x; X) (5) 
The second utilizes ffl (0) 
j (X) as a supplement to V (X) for motion of the heavy 
particle: 
H (1) 
j (X)' l;j (X) = ffl (1) 
l;j (X)' l;j (X) (6) 
where now 
H (1) 
j (X) = P 2 
2M + V (X) + ffl (0) 
j (X) (7) 
The outcomes of this sequential approximation are estimates for the exact 
eigenfunctions 
/ n (x; X) ß OE j (xjX)' l;j (X) (8) 
8 

and for the exact eigenvalues 
E n ß ffl (1) 
l;j (9) 
for the full two­particle Hamiltonian H. In principle these estimates could 
be inserted into eq. 3 to approximate time dependence of the two­particle 
dynamics. 
For all three models to be examined, the interaction v(x; X) was assumed 
to depend only on relative separation x 
extra simplification, 
the first stage of the BO sequential approximation, eqs. 4­5, leads to the 
explicit results [43] 
OE j (xjX) = N j exp[ 
1 4 (x 
and ! (0) is the harmonic frequency for the light particle oscillating about 
the fixed heavy particle, ! (0) = (k=m) 1=2 . Note that m is larger than the 
reduced mass for the pair of particles, ¯ = mM=(m+M ); consequently ! (0) 
is smaller than ! = (k=¯) 1=2 , the natural oscillation frequency for the pair in 
free space. 
Our three models are now distinguished by the following V (X) choices: 
V (X) = 1 
2 
KX 2 model I 
V (X) = ffl[exp( 
es a Morse oscillator potential 
for the heavy particle, and can be viewed as a one dimensional analog of 
an atom­substrate binding potential of the type often encountered in surface 
physics and chemistry [44]. The third introduces quartic anharmonicity, and 
if A ? 0, B ! 0 it can represent bistable motion. 
Because the ffl (0) , eq. 12, are equally spaced and independent of X for the 
harmonic v choice, eq. 10, the adiabatic potential surfaces for motion of the 
heavy particle are parallel and equally spaced for all three models. 
10 

3 Solution Procedure 
Our objective for each of the three models is to determine the precise quan­ 
tum evolution of an initial wavepacket state, and to compare these results 
with predictions of two specific approximation schemes (Ehrenfest mean field 
[17]­[33] and fewest switches surface hopping [35]). 
The initial state is selected to conform to the BO `picture' reflecting 
the intention to use the heavy and light particle coordinates as the relevant 
coordinates. For the model potentials that are explored below the use of 
approximately separable linear combinations of the bare coordinates might 
simplify our efforts. But despite this possibility this alternative was rejected 
in order to retain generality, given that as the complexity of the system in­ 
creases it becomes more difficult to identify such combinations. In particular, 
the light particle is placed in one of the H (0) eigenstates, OE j , and the heavy 
particle in a translating Gaussian wave­packet. Consequently we set, 
\Psi(x; X; t = 0) = N OE j (xjX) exp[ 
values typically around 20 times the inverse of the initial momentum, 
11 

unless otherwise noted. The BO functions OE j ' l;j form a complete set for all 
three models, and could be used for propagation in time. Indeed one of our 
objectives is to monitor time dependence of the light particle quantum­state 
probabilities 
a i (t) = 
X 
j 
c \Lambda 
j;i (t)c j;i (t) (16) 
where 
c j;i (t) = 
Z 
dX ' \Lambda 
j;i (X) 
Z 
dx OE \Lambda 
i (xjX )\Psi(x; X; t) (17) 
that indicate transitions between the adiabatic potential surfaces. In addition 
to following the a i (t) we have found it useful to examine the time dependence 
of the heavy particle mean position and momentum: 
! X ? (t) = 
Z Z 
dx dX \Psi \Lambda (t) X \Psi(t); (18) 
! P ? (t) = 
Z Z 
dx dX \Psi \Lambda (t) P \Psi(t) (19) 
We present these results below as phase space diagrams, i. e. as curves in the 
X;P plane. 
Model I has the advantage that the fully quantum problem is solvable 
in terms of independent normal modes (see Appendix I). This is used when 
examining the system time evolution. First, the trivial time evolution in 
the normal modes basis set is calculated and then transformed into the BO 
12 

basis set to obtain the a i mentioned above. Note that these normal mode 
quantum states differ from the BO states by an amount dependent on the 
mass ratio m=M ; as the ratio goes to zero the two representations converge 
to one another. However when the ratio is non­zero the chosen initial state 
(eq. 15) cannot be a proper normal mode eigenstate, so that the subsequent 
evolution inevitably involves transitions between the BO surfaces. Therefore 
this model is a very simple case to study the extent to which the mixed 
quantum­classical approximation schemes under consideration are suitable 
to represent the system evolution using the particle coordinates explicitly. 
To calculate the full quantum evolution for Models II and III fast Fourier 
transform techniques are used. The procedure used to solve the 2­dimensional 
time dependent Schr¨odinger equation is that of Kosloff and Kosloff [45] gen­ 
eralized to n light particle states (the n accessible adiabatic states). To do 
13 

so \Psi(x; X; t) is expanded in the OE j (xjX) basis; it becomes 
\Psi(x; X; t) ! 
0 
B B B B B B B B B B B B B B B B B B @ 
\Theta 0 (X; t) 
\Theta 1 (X; t) 
\Theta 2 (X; t) 
::: 
\Theta n (X; t) 
1 
C C C C C C C C C C C C C C C C C C A 
(20) 
where 
\Theta i (X; t) = 
Z 
dx OE \Lambda 
i (xjX) \Psi(x; X; t) (21) 
is a wavefunction evolving on V i (X) = V (X)+ ffl (0) 
i (X), and the Hamiltonian 
becomes a non­diagonal matrix, 
H i;j (X) = 
Z 
dx OE \Lambda 
i (xjX) H(x; X) OE j (xjX): (22) 
Thus the two­dimensional time evolution problem has been transformed into 
n­coupled one­dimensional time dependent Schr¨odinger equation. Note that 
the coupling occurs because the BO basis functions, OE j (xjX), are not eigen­ 
vectors of the P operator included in H (see eqn. 1). As the m=M ratio goes 
to zero the coupling vanishes, and as in the case of Model I transitions cease 
to occur. 
14 

The precise full quantum results will be compared with results generated 
by the two approximate methods mentioned above, Ehrenfest mean field 
and the fewest switches surface hopping. These are two well­established 
approximation schemes that are designed to take advantage of the simplicity 
of classical dynamics for the heavier particle, while presuming to account 
for nonadiabatic quantum effects of the lighter particle in a self­consistent 
manner. These methods are described below. 
The Ehrenfest mean field approximation postulates classical dynamics for 
the heavier particle, subjected to an effective force F eff defined as 
F eff = \Lambda H (0) \Phi 
oe 
(23) 
where \Phi(xjX[t]), the light particle wavefunction in the presence of the clas­ 
sically moving heavy particle, evolves according to the following time depen­ 
dent Schr¨odinger equation: 
i¯h 
@ 
@t 
\Phi(xjX[t]) = H (0) (x; X[t])\Phi(xjX[t]): (24) 
These equations define the method and can be used directly. Alternatively 
the time­dependent solution for the motion of the light particle can formally 
15 

be expressed 
\Phi(xjX[t]) = 
X 
j 
c j (t)OE j (xjX); (25) 
using the previously introduced light particle BO basis functions, OE j (xjX). 
These c j 's are the coefficients that will be use to compute the light particle 
quantum­state probabilities that will be compared with the results of the 
fully quantum calculations. Clearly, the equations for the c j 's time evolution 
can be obtained using eqn. 24, 
i¯h — 
c j = ffl (0) 
j c j 
is then given by the following expression: 
F eff = 
c \Lambda 
j c l + c \Lambda 
l c j ][ffl (0) 
j 
contribution to F eff , while the third term comprises nonadiabatic effects 
(light particle quantum number changes). Observe that the effective force 
interpolates among the forces defined separately on each of the accessible 
states. 
Fewest switches surface hopping is an intrinsically stochastic method in 
which the heavier particle executes finite­time­interval classical motion on 
distinct surfaces. Eqn. 26, the time dependent Schr¨odinger equation for the 
light particle \Phi(xjX[t]) in terms of the light particle adiabatic basis set, is 
also used here, but the motion of the classical heavy particle evolves on only 
one adiabatic potential surface at each instant of time, with instantaneous 
`hopping' between states according to the fewest switches algorithm. Thus, 
in contrast to the method described above, the heavier particle is subjected 
to a force that is simply 
F = 
total en­ 
ergy when hops occur. Note that this method --unlike the mean­field method-- 
is defined upon a basis set which in this case is chosen to be the light par­ 
ticle adiabatic basis set, OE j (xjX). The fewest switches algorithm minimizes 
17 

the number of state switches while maintaining a statistical distribution in 
an ensemble of trajectories that closely reproduces the light particle state 
populations. To do so, instantaneous switching probabilities from the occu­ 
pied level i to all the others states j during the time interval t to t + \Delta are 
calculated as 
g k;j = \Deltab j;k 
jc k j 2 
(31) 
where b is defined by (see eqn. 26) 
@jc j j 2 
@t 
= 
X 
k 6=j 
b j;k (32) 
= 
h to state 2 will occur 
if ¸ ! g 1;2 . A switch to state 3 will occur if g 1;2 ! ¸ ! g 1;2 + g 1;3 , etc. 
In addition, for a hop to occur the kinetic energy of the classical particle 
must be large enough to compensate for the loss in potential energy that 
might be involved in the transition. To understand how the algorithm works 
note that if the system has only two levels the switching probability is the 
rate at which the quantum light particle population probability is increasing 
on the surface where the classical particle might hop. In the absence of 
18 

`forbidden' hops, i.e., the appearance of non­zero amplitudes in states that 
are energy forbidden, the fewest switches algorithm statistically partitions 
trajectories correctly among different potential energy surfaces according to 
the probabilities jc j j 2 . The treatment of forbidden hops is discussed below. 
The approximate mixed quantum­classical methods involve solving both 
the time dependent Schr¨odinger equation eqn. 26, and the Newton equations. 
Both were solved numerically using the Runga­Kutta­Gill method [46]. 
4 Numerical Results. 
4.1 Model I: Harmonic Oscillator. 
This model consists of harmonic potentials in both the heavier and lighter 
particle coordinates, with a coupling that could be viewed as bilinear (ex­ 
panding the square of (x 
e the evolution of a system under this Hamiltonian 
can be determined trivially by transforming the problem into its separable 
19 

(normal) modes, see Appendix I. This will be used to calculate the exact 
(fully quantum) evolution of the system. But we reiterate that the motions 
of the heavy and light particles will be viewed as such, drawing upon the 
Born­Oppenheimer representation. Figure 1a displays the potential surfaces, 
the BO eigenvalues plus the heavy particle external potential (see eqn. 12), 
V n (X) = 1 
2 KX 2 + (n + 1 
2 )¯h 
q 
k=m (35) 
V n (X) = 1 
2 15X 2 + (n + 1 
2 ) 
p 
5: (36) 
These eqns. indicate choices for some of the parameters used in our numerical 
calculations. In addition the masses m = 1 and M = 10 were assigned, 
and for convenience we use atomic units. The width parameter was chosen 
to be oe = 0:15 in eqn. 15 and is such that the results do not significantly 
change when oe is smaller. Initially the ground state is the only populated 
adiabatic state, and X o and 
falls clearly in a regime 
were the BO approximation should not hold and transitions between the 
levels must occur. The separation of the parallel adiabatic states assures 
20 

that when a transition occurs the amount of energy transferred from the 
quantum degree of freedom to the classical one affects the dynamics of the 
heavy particle. This problem then may represent a challenge to any of the 
approximation methods. 
First we concentrate on the mean field approximation. Figure 1b shows 
the phase space diagram. The solid line corresponds to the solution of the 
fully quantum problem (! X ? is plotted vs ! P ?). The mean field results 
(X is plotted vs P ) are essentially identical to the fully quantum results. 
The dotted line is the quantum adiabatic approximation solution, i. e. no 
transitions (! X ? is plotted vs ! P ?). Thus exact quantum and mean 
field versions are identical, even though as shown by the contrast with the 
adiabatic approximation, transitions occur to a great extent. This result can 
be understood with the aid of the Ehrenfest theorem [47] which states that 
for any Hamiltonian of the form 
H(q) = 
p 2 
q 
2m + V (q); (37) 
the mean values ! q ? and ! p q ? are given by 
d 
dt ! q ?= ! p q ? 
m ; d 
dt ! p q ?= 
When applied to the present completely harmonic model, eq. 34, we obtain 
d 
dt ! x ? = ! p ? 
m 
; d 
dt ! p ?= 
heavy parti­ 
cle motion, and then taking the classical limit, produces the following Newton 
equations of motion: 
d 
dt 
X = P 
M 
; d 
dt 
P = @ 
@X 
ae 1 
2 KX 2 + 
Z 
dx\Phi \Lambda H (0) \Phi 
oe 
= KX 
the second to the last 
equation. The dynamics of the quantum degree of freedom can be examined 
applying the Ehrenfest theorem to H (0) (see eqn. 4), this leads to 
d 
dt 
! x ?= ! p ? 
m 
; d 
dt 
! p ?= 
and 42, are identical to the solution of the fully quantum problem, eqns. 39 
and 40 (with X $! X ? and P $! P ? ), explaining why the mean field 
approximation is exact while calculating the evolution of the heavy particle. 
This is true only because the coupling is bilinear and the potentials quadratic 
in both coordinates, otherwise deviations occur as will be seen below. As a 
curiosity --it is hard to envision a practical direct use of what follows-- note 
that if the role of x and X are interchanged while applying the mean field 
method (the heavier particle is treated by quantum mechanics and the lighter 
by classical mechanics) their evolution in phase space would still be described 
exactly. 
Even for harmonic models, the mean­field mixed quantum­classical method 
is exact only in reproducing the first moments, ! x ? and ! p ?. Higher 
moments are not computed exactly. This is apparent when comparing the 
time dependence of the light particle quantum­state probabilities in figure 
1c. While the curves resemble each other, they are not identical. 
Overall the performance of the mean field method is remarkably good 
for model I. The harmonic oscillator problem is peculiar and therefore it is 
interesting to study to what degree this result survives when the harmonicity 
is somehow relaxed. Clearly, there are circumstances where the method will 
23 

eventually stop performing well. In particular, as has been described in 
detail elsewhere [15, 35, 42, 48] this method cannot describe accurately the 
evolution of minority reaction channels if these channels involve potential 
energy surfaces that differ considerably from the one experienced by the 
majority channel. It would be of little use to modify the potential just in 
order to see the mean field approximation fail; it seems more informative 
instead to concentrate on a specific class of potentials. Here, as mentioned 
earlier, we will focus on the class of potentials that are harmonic in the 
separation between x and X but let the V (X) change (the potential energy 
surfaces are then constrained to be parallel) and monitor the degree to which 
the given approximation holds. 
Now we focus on the performance of the second approximation method, 
the fewest switches surface hopping algorithm. The results shown correspond 
to averaging 1000 trajectories, all starting at (X 0 ; 
Figure 1d shows the resulting phase space 
diagrams. The different lines have a similar meaning to those in figure 1b 
although here the classical quantities represent an average over all the trajec­ 
tories. Note how the first two curves differ significantly, especially at longer 
times. Figure 1e shows the population evolutions. The solid line is the exact 
24 

result, the broken line the population calculated as the fraction of trajectories 
on each surface and the dotted line the population calculated as the average 
(over all the different trajectories) of c \Lambda 
i c i . At early times, the populations 
computed by surface­hopping are very similar to the mean field results of Fig. 
1c. At longer times, the results of surface hopping show damped oscillations 
in the populations and an inward spiraling in the phase space diagram. This 
can be attributed to the fact that different trajectories exhibit a different 
history of random hops. Over a period of time this causes the ensemble to 
lose coherence and therefore spread in phase space. A second factor that 
can introduce error in surface hopping is associated with `forbidden hops'. 
With the algorithm employed here, hops induced by the evolving quantum 
populations are tested for energetic feasibility. In the present calculation the 
classical particles were started at a turning point and immediately quantum 
population began to develop on the upper surface. On the other hand some 
time elapses before the classical particle reaches a position where it is ener­ 
getically possible for it to begin to evolve on the upper energy surface. As a 
consequence the fewest switches algorithm fails to statistically partition tra­ 
jectories correctly among different potential energy surfaces. This not appear 
to be the major source of error in this case, however. 
25 

A related issue is that of the treatment of a forbidden surface hop. It has 
been argued [49] that the velocity of the classical particle should be reversed 
after such an occurrence. More recently, M¨uller and Stock [48] have shown 
that it may be more accurate to continue motion in the same direction. 
The reason behind the unsatisfactory performance of the surface hop­ 
ping algorithm for longer times is not related to the precise nature of the 
potential in Model I. Whenever the potential involves a coupling that does 
not vanish outside a localized scattering region, surface hopping will intro­ 
duce unphysical loss of coherence. In addition, this class of interactions also 
encourages the presence of regions where the hops are classically forbidden 
but encouraged by quantum mechanics. This situation describes all of the 
models studied in this work; they involve heavy and light particles connected 
by a spring and therefore coupled for all times. The only difference between 
them is the profile of the scattering region. Therefore we will focus on the 
performance of the mean field approximation. 
26 

4.2 Model II: Morse Oscillator A. 
In relaxing the harmonic restriction in the form of V (X) an obvious choice 
is that of the Morse oscillator, 
H(x; X) = p 2 
2m + P 2 
2M + 
1 
2 k(x 
n (X) = ffl[exp( =m (44) 
V n (X) = 60[exp( 
and M = 10, and initially 
the ground adiabatic state is the only one populated in the lighter particle 
degree of freedom. 
The anharmonicity of the potential causes the quantum wavepacket to 
deform, losing its Gaussian shape as time progresses, unlike model I. There­ 
fore, when using this mixed quantum­classical method one trajectory is not 
sufficient. It is necessary to perform dynamics for a collection of trajectories 
differing in the initial position and momentum of the heavy particle. To 
27 

choose the initial conditions we have used a Wigner transform [50]. The 
Wigner transform carries a density operator in the coordinate representation 
ae(X; X 0 ) = j\Theta(X ) ?! \Theta(X 0 )j to phase space; 
ae(X; P ) = 1 
߯h 
Z 1 
ber of different initial conditions is such that 
the results are substantially independent of it, and this number depends not 
only on the potential surface but also on the initial conditions and length of 
the dynamics. 
Figure 2b shows phase space diagrams. The solid line corresponds to the 
solution of the fully quantum problem (solved numerically as described in 
section 3), the broken line to the ensemble averaged mean field results and 
the dotted line to the quantum adiabatic approximation solution. Figure 2c 
shows the time dependence of the light particle quantum­state probabilities. 
Even though the results are not identical, the method performance is still 
quite good. We emphasize that this result, although representative, was 
obtained keeping in mind the goal of finding the limitations of the mean field 
method for this class of potentials. Its success motivated the selection of 
28 

model III described in the next subsection. 
4.3 Model III: Double Well. 
The mean field method cannot properly describe situations for which the 
wave packet bifurcates into two or more distinct paths governed by different 
interactions. 
Model III introduces quadratic anharmonicity that can represent bistable 
motion (if A ? 0 and B ! 0). 
H(x; X) = p 2 
2m + P 2 
2M + 1 
2 k(x 
the potential surfaces, 
V n (X) = AX 4 +BX 2 + (n + 
1 
2 )¯h 
q 
k=m (48) 
V n (X) = X 4 
freedom is a Gaussian 
29 

wavepacket in the fully quantum calculation, while in the mixed classical­ 
quantum calculations it is described by a collection of X o and 
in 
section 3), the broken line to the ensemble averaged mean field results and 
the dotted line to the quantum adiabatic approximation solution. Figure 3c 
shows the time dependence of the light particle quantum­state probabilities. 
Note that the phase space trajectory, Fig. 3b, spirals inwards approaching 
an unsymmetric position of x about 0:5. Quantum mechanically, this is due 
to trapping of the fraction of the wave packet that has been excited to the 
n = 1 level. The initial conditions of the wave packet have been chosen 
so that there is sufficient energy for most of the wave packet to traverse 
the barrier in the ground (n = 0) state, but not in the first excited (n = 1). 
Tunneling through the excited state is sufficiently slow that the asymmetry in 
position persists over the time scale of the calculation. It is interesting to note 
that the mean field approximation reproduces this effect quite well, at least 
in an average way. There is a spread of initial energies of the independent 
30 

mean field trajectories, determined from the Wigner distribution. For the 
conditions of Fig. 3, trajectories in the low energy tail of this distribution 
get trapped by the potential barrier (slightly higher that the ground state 
barrier due to mixing with the excited state). The high energy majority 
traverse the barrier. The phase space asymmetry in the mean­field result 
of Fig 3 is thus largely determined by the width of the initial momentum 
distributions, not by the probability of excitation to the n = 1 state as in 
the fully quantum case. Nevertheless, on average, the mean field method 
quite accurately reproduces the asymmetry of the full quantum calculation, 
at least for the initial conditions chosen here. 
Even though the results are not perfect the performance of the mean field 
method is still quite good, despite the fact that during the evolution a small 
but significant portion of the population is in the excited state. Note that at 
the energy we selected, the potential does not permit exploring large regions 
of space. This may imply that even when the detailed evolution is not being 
described properly, the evolution of the averaged quantities that are being 
measured is accurate. Therefore, in the next section the model II will be 
reexamined but in a case were the evolution is not confined to a restricted 
region of space. 
31 

4.4 Model II revisited: Morse Oscillator B and C. 
The potential surfaces involved in the first example in this subsection are 
displayed in figure 4a, 
V n (X) = 10[exp( 
opulation 
might be transiently trapped in the well of the first excited potential curve. 
The masses are m = 1 and M = 10, and the initial state is the ground 
adiabatic state in the lighter particle degree of freedom, and in the heavier 
particle degree of freedom it is a gaussian wavepacket in the fully quantum 
calculation, and in the mixed classical­quantum calculations it is described 
by a set of X o and 
), is also shown in figure 4a, its baseline 
representing total energy of the system in this calculation. Note that outside 
the well region the total energy is less than V 1 (X). 
Figure 4b shows phase space diagrams. The solid line corresponds to the 
solution of the fully quantum problem (! X ? vs ! P ?), the broken line to 
the ensemble averaged mean field results and the dotted line to the quantum 
32 

adiabatic approximation solution. Figure 4c shows the time dependence of 
the light particle quantum­state probabilities. Note the persistent oscillations 
(clearly shown in the inserts displaying the 190­200th time units) in both the 
populations and in the phase space diagram, resulting from the fact that 
the coupling between the heavy and light particle does not vanish when the 
compound particle leaves the scattering region. The oscillations correspond 
to the periodic transfer of energy between the heavy and light particles, and 
far from the scattering region simply monitor the behavior of the heavier 
particle in the unperturbed composite particle. 
The amount of population transfer is rather small, despite our efforts 
to find parameters where the mean field approximation was likely to fail. 
This surprising situation is caused by the fact that if the initial momen­ 
tum is increased the size of the oscillations just described increase as well, 
thus masking any other effects. This adds to the fact that increasing the 
magnitude of the coupling does not have a significant net effect on the tran­ 
sition rates since k is also proportional to the separation of the levels and 
therefore inversely proportional to the transition probabilities. Despite this, 
the method performance is still better than expected even though the pop­ 
ulations in the different channels evolve in regions that are quite separated 
33 

spatially. 
This `puzzle' can be solved by recognizing that population can be evolving 
in the upper surface but not be trapped (and therefore not in a minority 
channel). Figure 4d and 4e show \Theta 0 (X; t) and \Theta 1 (X; t) respectively. The 
region enclosed in the quadrilateral figure in 4e is the only portion of the 
population effectively trapped, and this portion is almost non­existent. Most 
of the population in the upper surface behaves as a `ghost' of the population 
on the ground state, the behavior of the former merely mirroring that of the 
latter. This is a consequence of the persistent coupling between the light and 
heavy particles. 
Thus, to observe a significant failure the transition rate must be somehow 
increased. This may be accomplished by increasing the slope of the inner wall 
and therefore increasing the match between the frequency of the bounce and 
that of the oscillations. This creates a small fraction of the population with 
a significantly different amount of kinetic energy than that of the majority, 
consequently evolving under inappropriate forces. That is the case if the 
potential surfaces are 
V n (X) = 0:001[exp( 
as shown in figures 5a­e (the figures are equivalent to those in figure 4). 
In this example the well is extremely shallow and the inner wall is quite 
steep, therefore it might represent a collision with a hard wall. To encourage 
transitions given the absence of a substantial well, the population is initially 
in the first excited adiabatic state. An alternative would have been to give 
the heavier particle more initial kinetic energy, but this would have resulted 
in the appearance of masking oscillations. 
This final choice of parameters produces a bifurcation of the initial wave 
packet into two scattered portions of quite different velocity. An average 
path cannot be expected to properly describe the individual pathways, so the 
mean field method is expected to break down. The surface hopping method 
was developed in order to handle such bifurcation. However, this specific 
example, in addition to producing a bifurcation, involves a coupling between 
light and heavy that never vanishes. It is likely, therefore, that surface­ 
hopping would also have difficulty with this model, at least reproducing the 
long time oscillatory behavior. 
35 

5 Discussion and Conclusions 
We have examined the accuracy of mixed quantum­classical dynamics as ap­ 
plied to an increasingly challenging series of model problems. The models 
involve one light and one heavy degree of freedom, and exhibit substantial 
nonadiabatic behavior. In all of the models the coupling between the light 
and heavy particles is linear (harmonic). In addition, an external potential 
is applied to the heavy particle only. The models differ only in the mag­ 
nitude of the coupling force constant, and the form and magnitude of the 
external potential. The energies of the light particle quantum states, as a 
function of the position of the heavy particle, define a sequence of ground 
and excited Born­Oppenheimer potential energy curves. Because the light 
particle experiences a purely harmonic potential, the potential energy curves 
are parallel and equally spaced for all of the models. In addition, the nona­ 
diabatic couplings among potential energy curves persists for all times due 
to the non­vanishing linear coupling between light and heavy particles. 
The model problems were used to test two strategies for carrying out 
mixed quantum­classical dynamics in systems involving nonadiabatic tran­ 
sitions, mean­field and surface­hopping. With the mean­field approach the 
36 

classical motion of the heavy particle is governed by an effective potential 
energy surface resulting from an average over light particle quantum states, 
weighted by their time­varying populations. It is well known that this method 
is inadequate in cases where the quantum wave packet splits into parts that 
are governed by substantially different forces. The surface­hopping procedure 
for stochastically splitting a trajectory into multiple branches was developed 
in order to overcome this difficulty. In the present series of model studies, 
however, the forces associated with each quantum state are identical since 
the potential energy curves are parallel. The necessity for splitting trajec­ 
tories into branches is therefore not apparent. Furthermore, the fact that 
nonadiabatic transitions continue to occur for arbitrarily long times presents 
a significant obstacle for surface­hopping methods; the different stochastic 
histories of individual trajectories result in an unphysical loss of quantum 
coherence at long times. This was demonstrated by application to model I, 
where surface hopping proved to be an adequate approximation to the exact 
quantum results at short times, but was inadequate at long times and clearly 
inferior to the mean field method. 
The fact that the forces are identical on each potential energy curve does 
not ensure the success of the mean­field method. First, the classical approx­ 
37 

imation is certainly not exact even when motion is confined to a single po­ 
tential energy surface. Since our interest here is in mixed quantum­classical 
treatment of systems involving multiple potential surfaces, we selected model 
parameters so that classical mechanics provided an accurate description of 
single surface dynamics, as judged by comparison with exact quantum cal­ 
culations within the adiabatic approximation. Second, even with parallel 
potential energy surfaces, the mean­field approximation is not exact. The 
results reported above help delineate the range of validity of the method. 
For a completely harmonic system, Model I, we have shown in analogy 
with the standard Ehrenfest theorem that the mean­field classical heavy par­ 
ticle position X and momentumP reproduce exactly the corresponding quan­ 
tum mechanical expectation values ! X ? and ! P ?. Higher moments of 
X and P and quantum state probabilities are not reproduced exactly by 
the mean­field method, but they are quite accurate. With the addition of 
anharmonicity into the external potential that acts on the heavy particle, 
Models II and III, the Ehrenfest correspondence is no longer exact, and the 
mean­field results deviate somewhat more from the exact quantum results. 
Nevertheless, the mean­field method was of acceptable accuracy for all but 
one of the model systems studied. 
38 

For the double well potential, Model III, we expected that the mean­ 
field method would be incapable of correctly describing the bifurcation of 
the quantum wave packet into two parts localized respectively in the left and 
right potential wells. However, the fraction of mean­field trajectories trapped 
in each well agreed quite accurately with the quantum probabilities as did 
the mean values of X and P. The only case for which the mean­field method 
did not provide an acceptable description was the final parameterization of 
Model II, for which the quantum wave packet split into a scattered wave 
packet on the ground state potential curve and a scattered wave packet on 
the first excited potential. 
The class of models considered here, restricted to linear light­heavy cou­ 
pling, is admittedly artificial. The models were chosen to elucidate the ap­ 
plication of mixed quantum­classical dynamics to condensed phase systems 
involving many heavy (classical) degrees of freedom which are coupled to a 
few light (quantum) degrees of freedom. In such systems it is impractical to 
carry out dynamical simulations using normal mode coordinates. Rather, the 
bare coordinates of the heavy and light coordinates are generally used. In the 
bare particle coordinate system, the coupling between light and heavy par­ 
ticles will have a linear term. It is reasonable to suppose that in many cases 
39 

this linear term will dominate the coupling, particularly when the masses of 
the light and heavy particles do not differ greatly, e.g., hydrogen vs. car­ 
bon. The model calculations reported here suggest that, in cases where 
linear couplings dominate, the mean­field mixed quantum­classical method 
can be of useful accuracy and is robust to the introduction of anharmonic 
heavy­particle interactions. 
Acknowledgement 
We thank Dr. David Sholl for valuable discussions. 
Appendix I. Normal modes of Model I. 
The total Hamiltonian is, 
H(x; X) = p 2 
2m + P 2 
2M + 1 
2 k(x 
0 
B 
B 
B 
@ 
ff fl 
fl fi 
1 
C 
C 
C 
A 
0 
B 
B 
B 
@ 
q 
Q 
1 
C 
C 
C 
A 
(53) 
40 

where ff = k=m, fi = (k + K)=M and fl = 
fi 
fi 
fi 
fi 
fi 
fi 
fi 
fi 
fi 
ff 
\Sigma = ff + fi 
2 \Sigma 
s 
(ff 
(–+ 
modes coordinates i and j, 
H(x; X) = 
p 2 
j 
2 + 
p 2 
i 
2 + –+ j 2 
2 + – 
š j (j) = N j exp[ 2 ]H j (j 0 ); ffl j 
j = (j + 1 
2 
) 
q 
–+ ; (61) 
¸ i (i) = N i exp[ 
, H n are the Hermite polynomials, 
and N n are normalization constants. 
Figure Captions 
Figure 1. Numerical example for model I: Harmonic Oscillator. 
a. Potential surfaces. The small circle indicates the initial position of the 
heavy particle (in the mixed quantum­classical calculations) or the center of 
the wavepacket (in the fully quantum calculations). 
b. Mean Field results: Phase space diagram. The solid line corresponds 
to the solution of the fully quantum problem (! X ? is plotted vs ! P ?). 
For this case the mean field results (X vs P ) are identical to the fully quantum 
results. The dotted line corresponds to the quantum adiabatic approximation 
solution (! X ? is plotted vs ! P ?). The figure shows the first twenty 
time units of the evolution. 
c. Mean Field results: Evolution of populations. The solid line is the 
numerically exact ground state population and the broken line is the result 
42 

of the mean field approximation, both are plotted as a function of time. Only 
the ground state population is shown for clarity although there is population 
transfer to the first and second excited states. 
d. Surface Hopping results: Phase space diagram. The solid line corre­ 
sponds to the solution of the fully quantum problem, the broken line to the 
surface hopping results (average X is plotted vs average P ) and the dotted 
line to the quantum adiabatic approximation solution. The first twenty time 
units of the evolution are shown. 
e. Surface Hopping results: Evolution of populations. The solid line is the 
numerically exact ground state population, the broken line the population 
calculated as the fraction of trajectories on each surface and the dotted line 
the population calculated as the average (over all the different trajectories) 
of jc i j 2 . The last two are results of the surface hopping approximation. As in 
figure 1c only the ground state population is shown for clarity. 
Figure 2. Numerical example for model II: Morse Potential A. 
a. Potential surfaces. The initial wavepacket, \Theta 0 (X; t = 0) (see text) is 
also shown. 
b. Phase space diagram. The solid line corresponds to the solution of 
43 

the fully quantum problem (! X ? is plotted vs ! P ?), the broken line to 
the mean field results (ensemble averaged X vs ensemble averaged P ) and 
the dotted line to the quantum adiabatic approximation solution (! X ? 
is plotted vs ! P ?). The figure shows the first twenty time units of the 
evolution. 
c. Evolution of populations. The solid lines are the numerically exact 
populations and the broken lines are the result of the ensemble averaged 
mean field approximation, both are plotted as a function of time. Only the 
ground state population is shown for clarity although there is population 
transfer to the first and second excited states. 
Figure 3. Numerical example for model III: Double Well. 
a. Potential surfaces. The initial wavepacket, \Theta 0 (X; t = 0) (see text) is 
also shown. 
b. Phase space diagram. Same as Fig. 2. 
c. Evolution of populations. Same as Fig. 2. 
Figure 4. Numerical example for model II: Morse Potential B. 
a. Potential surfaces. The initial wavepacket, \Theta 0 (X; t = 0) (see text) is 
44 

also shown. 
b. Phase space diagram. Same as Fig. 2. 
c. Evolution of populations. Same as Fig. 2. 
The inserts in figures b and c display the evolution between the 190 and 
200 time units. 
d. Evolution of \Theta 0 (X; t). 
e. Evolution of \Theta 1 (X; t). The region enclosed in the quadrilateral is the 
only portion of the population effectively trapped. 
Figure 5. Numerical example for model II: Morse Oscillator C. 
a. Potential surfaces. The initial wavepacket, \Theta 1 (X; t = 0) (see text) is 
also shown. 
b. Phase space diagram. Same as in Fig. 2. 
c. Evolution of population on the ground state. 
d. Evolution of population on the second excited state. Note how when 
the mean field approximation is used the state to which most of the popula­ 
tion is transferred is not the correct one. 
e. Evolution of \Theta 0 (X; t). The region enclosed in the quadrilateral is the 
portion of the population evolving in the minority channel. 
45 

f. Evolution of \Theta 1 (X; t). 
46 

References 
[1] M. Born and J. R. Oppenheimer, Ann. Phys, 84, 457 (1927). 
[2] L. Pauling, The Nature of the Chemical Bond, Third edition (Cornell U. 
P., Ithaca, 1960), chap. 12. 
[3] K. R. Atkins, Liquid Helium (Cambridge U. P., Cambridge, 1959). 
[4] P. J. Rossky and J. Schnitker, J. Phys. Chem. 92, 4277 (1988) 
[5] M. P. Allen and D. J. Tildesley, Computer Simulation of Liquids (Claren­ 
don, Oxford, 1987). 
[6] L. M. Raff and D. L. Thompson, Theory of Chemical Reaction Dynam­ 
ics, Vol II, ed. M. Baer (CRC Press, Boca Raton, Fl, 1985), p 1. 
[7] E. J. Heller, Acc. Chem. Res. 14, 368 (1981). 
[8] W. H. Miller, J. Chem. Phys. 53, 1949 (1970). 
[9] J. H. Irving and R. W. Zwanzig, J. Chem. Phys. 19, 1173 (1951) 
[10] B. J. Berne and D. Thirumalai, Ann. Rev. Phys. Chem. 37, 401 (1986) 
[11] N. F. Mott, Proc. Cambridge Phil. Soc. 27, 553 (1931). 
47 

[12] E. E. Nikitin, Theory of Elementary Atomic and Molecular Processes in 
Gases, (Clarendon, Oxford, 1974). 
[13] A. G. Redfield, Adv. Magn. Reson. 1, 1 (1965). 
[14] J. M. Jean, R. A. Friesner and G. R. Fleming, J. Chem. Phys. 96, 5827 
(1992). 
[15] J. C. Tully, Modern Methods for Multidimensional Dynamics Com­ 
putation in Chemistry edited by D. L. Thompson, (World Scientific, 
Singapore) in press. 
[16] M. F. Herman, Annu. Rev. Phys. Chem. 45, 83 (1994) 
[17] A. D. McLachlan, Mol. Phys. 8, 39 (1964). 
[18] W. H. Miller and C. W. McCurdy, J. Chem. Phys. 69, 5163 (1978). 
[19] H. D. Meyer and W. H. Miller, J. Chem. Phys. 72, 2272 (1980). 
[20] D. A. Micha, J. Chem. Phys. 78, 7138 (1983). 
[21] Z. Kirson, R. B. Gerber, A. Nitzan and M. A. Ratner, Surf. Sci. 137, 
527 (1984). 
48 

[22] S. I. Sawada, A. Nitzan and H. Metiu, Phys. Rev. B. 32, 851 (1985). 
[23] S. Mukamel, Y. J. Yan and J. Grad, Stochasticity and Intramolecular 
Redistribution of Energy, ed. R. Lefebvre and S. Mukamel (D. Reidel, 
Dordrecht, Holland, 1986), pag. 109. 
[24] G. Wahnstrom, B. Carmeli and H. Metiu, J. Chem. Phys. 88, 2478 
(1988). 
[25] M. Amarouche, F. X. Gadea and J. Durup, Chem. Phys. 130, 145 
(1989). 
[26] A. T. Amos, K. W. Sulston and S. G. Davidson, Adv. Chem. Phys. 76, 
335 (1989). 
[27] H. J. C. Berendsen and J. Mavri, J. Phys. Chem. 97, 13464 (1993). 
[28] P. Bala, B. Lesyng and J. A. McCammon, Chem. Phys. Lett. 219, 259 
(1994). 
[29] G. D. Billing, Int. Rev. Phys. Chem. 13, 309 (1994). 
[30] G. Stock, J. Chem. Phys. 103, 2888 (1995). 
[31] M. Head­Gordon and J. C. Tully, J. Chem. Phys. 103, 10137 (1995). 
49 

[32] D. Antoniu and S. D. Schwartz, J. Chem. Phys. 104, 3526 (1996). 
[33] M. Hintenender, F. Rebentrost, R. Kosloff and R. B. Gerber, J. Chem. 
Phys. 105, 11347 (1996). 
[34] J. C. Tully and R. K. Preston, J. Chem. Phys. 55, 562 (1971). 
[35] J. C. Tully, J. Chem. Phys. 93, 1061 (1990). 
[36] J. C. Tully, Int. J. Quantum Chem. Symp. 25, 299 (1991). 
[37] F. Webster, P. J. Rossky and R. A. Friesner, Comput. Phys. Commun, 
63, 494 (1991). 
[38] F. J. Webster, J. Schnitker, M. S. Friedrichs, R. A. Friesner and P. J. 
Rossky, Phys. Rev. Lett. 66, 3172 (1991). 
[39] D. F. Coker and L. Xiao, J. Chem. Phys. 102, 496 (1995). 
[40] H. S. Mei and D. F. Coker, J. Chem. Phys. 104, 4755 (1996). 
[41] O. V. Prezhdo and P. J. Rossky, J. Chem. Phys. 107, 825 (1997) 
[42] C. C. Martens and J­Y. Fang, J. Chem. Phys. 106, 4918 (1997). 
50 

[43] C. Cohen­Tannoudji, B. Diu and F. Lalo¨e, Quantum Mechanics (Wiley, 
New York, 1977) vol 1, pag. 500. 
[44] S. Ross and J. P. Oliver, Solid Surfaces and the Gas­Solid Interface, Ad­ 
vances in Chemistry Series No. 33 (American Chemical Society, Wash­ 
ington, D. C., 1961), p. 303. 
[45] D. Kosloff and R. Kosloff, J. Comp. Phys. 52, 35 (1983). 
[46] M. J. Romanelli, Mathematical Methods for Digital Computers edited 
by A. Ralston and H. S. Wilf (Wiley, New York, 1960), Chap 9. 
[47] C. Cohen­Tannoudji, B. Diu and F. Lalo¨e, Quantum Mechanics (Wiley, 
New York, 1977) vol 1, pag. 242. 
[48] U. M¨uller and G. Stock, J. Chem. Phys. 107, 6230 (1997) 
[49] S. Hammes­Schiffer and J. C. Tully, J. Chem. Phys. 101, 4657 (1994). 
[50] E. P. Wigner, Phys. Rev., 40, 747 (1932); M. Hillery, R. F. O'Connell, 
M. O. Scully and E. P. Wigner, Phys. Rep., 106, 121 (1984) 
51