combined quantum/classical dynamics scheme. Any average energy increase

(decrease) in the quantum degrees of freedom is compensated by a decrease

(increase) in energy of the classical degrees of freedom.7 Note that the

force on the classical particles (5.81) is the expectation of the force matrix,

ˆ

which is the expectation of the gradient of H and not the gradient of the

ˆ

expectation of H. The latter would also contain a contribution due to the

gradients of the basis functions in case these are a function of R. The force

calculated as the expectation of the negative gradient of the Hamiltonian is

called the Hellmann“Feynman force; the matrix elements Fnm in (5.81) are

the Hellmann“Feynman force matrix elements. The use of forces averaged

over the wave function is in accordance with Ehrenfest™s principle (see (3.15)

on page 43). Note that the energy conservation is exact, even when the basis

7 The DME method with average back reaction has been applied to a heavy atom colliding with

a quantum harmonic oscillator (Berendsen and Mavri, 1993) and to a quantum harmonic oscil-

lator in a dense argon gas by Mavri and Berendsen (1994). There is perfect energy conservation

in these cases. In the latter case thermal equilibration occurs between the harmonic oscillator

and the gas.

132 Dynamics of mixed quantum/classical systems

set is incomplete and the wave function evolution is therefore not exact.

The forces resulting from gradients of the basis functions are called the

Pulay forces; as shown above, they should not be included in the forces on

the classical particles when the perturbation term due to the non-adiabatic

coupling is included in the dynamics of the quantum subsystem. In the next

section we return to the Pulay forces in the adiabatic limit.

Proof First we split (5.84) into two parts:

dEtot ™

™

= tr (ρH) + tr (ρH), (5.85)

dt

and consider the ¬rst part, using (5.37):

i ™ ™

tr ([ρ, H] H) + R · tr ([ρ, d] H) = R · tr ([ρ, d] H),

™

tr (ρH) = (5.86)

because tr ([ρ, H] H) = 0, since tr (ρHH) = tr (HρH).8 The second part

of (5.85) can be written out as follows:

™ ™ ™

tr (ρH) = R · tr (ρ∇R H) + P · tr (ρ∇P H). (5.87)

The tricky term is ∇R H:

ˆ ˆ ˆ

∇R Hnm = ∇R φn |H|φm + φn |∇R H|φm + φn |H|∇R φm .

The middle term equals ’Fnm , according to (5.81); in (5.87) it produces a

™™

term ’R · P , which exactly cancels the last term in (5.87). The ¬rst term

ˆ

can be rewritten with the use of the Hermitian property of H:

dr ∇R φ— Hφm = dr (H∇R φn )— φm

ˆ ˆ ˆ

∇R φn |H|φm = n

= (Hd)— = [(Hd)† ]nm ,

mn

while the third term equals (Hd)nm . Realizing that

(Hd)† = d† H† = ’dH,

the ¬rst and third term together are equal to [H, d]nm . Thus (5.87) reduces

to

™ ™ ™

tr (ρH) = R · tr (ρ [H, d]) = ’R · tr ([ρ, d] H),

which exactly cancels the term (5.86) left over from the ¬rst part of (5.85).

8 The trace of a matrix product is invariant for cyclic permutation of the matrices in the product.

5.3 Embedding in a classical environment 133

5.3.2 Forces in the adiabatic limit

The considerations in the previous section allow an extrapolation to the

pure Born“Oppenheimer approximation (or adiabatic limit) and allow an

analysis of the proper forces that act on the classical degrees of freedom in

the adiabatic limit, which have been the subject of many discussions. The

leading principle, again, is the conservation of total energy.

In the Born“Oppenheimer approximation it is assumed that for any point

in classical phase space (R, P ) the time-independent Schr¨dinger equation

o

ˆ = Eψ has been solved exactly for the ground state. The classical coor-

Hψ

dinates R are only parameters in this solution. The system is assumed to be

and remain in the pure ground state, i.e., in terms of a density matrix with

the exact solutions of the Schr¨dinger equation as basis functions, numbered

o

™

n = 0, 1, . . ., ρ00 = 1 and ρ = 0. No transitions to excited states are allowed.

It is this assumption that forms the crucial approximation. It also implies

that the system always remains in its exact ground state, i.e., that it follows

adiabatic dynamics.

The assumption ρ = 0 implies not only that the Hamiltonian is diagonal,

™

including whatever small perturbations there are, but also that the second

™

term in (5.35): R[ρ, d], is completely neglected.

The ground-state energy E0 (R) functions as the potential energy for the

Hamiltonian dynamics of the classical degrees of freedom. Hence the forces

must be the negative gradients of the ground-state energy in order to con-

serve the total energy:

—ˆ

F = ’∇ dr ψ0 Hψ0 . (5.88)

Since all three elements in the integral depend on R, this force is not equal

to the Hellmann-Feynman force:

— ˆ

F HF = ’ dr ψ0 ∇Hψ0 . (5.89)

The di¬erence is the Pulay force due to the dependence of the wave function

on the nuclear coordinates:

—ˆ —ˆ

F Pulay = ’ dr ∇ψ0 Hψ0 ’ dr ψ0 H∇ψ0 = [d, H]00 , (5.90)

where the last equality follows from the same reasoning as was followed in

the proof of (5.84) on page 132. It is not surprising that the Pulay force

reappears, because it only canceled in (5.84) against the now neglected term

in the density matrix evolution.

The Pulay force seems to be a nasty complication, but it isn™t. When

134 Dynamics of mixed quantum/classical systems

the Hamiltonian is really diagonal, the term [d, H]00 is equal to zero and

the Pulay force vanishes. So, for a pure adiabatic process, the Hellmann-

Feynman forces su¬ce.

5.3.3 Surface hopping dynamics

The method of surface hopping (SH), originating from ideas of Pechukas

(1969a, 1969b), was introduced by Tully and Preston (1971) and speci¬ed

with the minimum number of switches by Tully (1990). The method was

designed to incorporate the in¬‚uence of excited electronic states into the

atomic motion. The basic notion is that there is no single best atomic

trajectory subject to the in¬‚uence of electronic transitions (which would

lead to mean-¬eld behavior), but that trajectories must split into branches.

The splitting is accomplished by making sudden switches between electronic

states, based on the diagonal elements of the quantum density matrix. The

density matrix is evolved as usual by (5.35). The sudden switches are de-

termined in each time step ”t by a random process based on the transition

probabilities from and to the state before the time step. More precisely, if

the present state is n, then consider the rate of change of the population

ρnn from (5.35):

i

(ρnm Hnm ) ’ R · (ρnm dnm ’ ρ— dnm )

— ™

ρnn =

™ (5.91)