m=n

2

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

™

= (5.92)

nm

m=n

which can be written as

ρnn =

™ bnm . (5.93)

m=n

Now de¬ne a switching probability gnm within a time step ”t from the

current state n to other states m as

”tbmn

gnm = . (5.94)

ρnn

If gnm < 0, it is set to zero (i.e., only switches to states with higher prob-

abilities are allowed; this in fact corresponds to the condition of a minimal

number of switches). The cumulative probabilities hm = m gnk are de-

k=1

termined. A uniform random number between 0 ¤ ζ < 1 is generated and

a switch to state m occurs when hm < ζ < hm+1 .

In the vast majority of cases no switch will occur and the system remains

5.3 Embedding in a classical environment 135

in the same state n. The classical forces are now calculated as Hellmann-

Feynman forces from the nth state, not from the complete density matrix. If

a switch occurs, one must take ad hoc measures to conserve the total energy:

scale the velocity of the classical degrees of freedom (in the direction of the

nonadiabatic coupling vector). If the kinetic energy of the particle does not

su¬ce to make up for the increase in energy level after the switch, the switch

is not made.

There are many applications of surface hopping. Mavri (2000) has made

an interesting comparison between mean-¬eld DME, SH and exact quantum

calculation of a particle colliding with a quantum oscillator and found in

general that SH is closer to exact behavior than mean-¬eld DME, but also

concluded that both have their advantages and disadvantages, depending

on the case under study. There is nothing in the Schr¨dinger equation that

o

compels a treatment one way or another.

5.3.4 Other methods

The situation regarding the consistent treatment of the back reaction from

quantum-dynamical subsystems is not satisfactory. Whereas mean-¬eld

DME fails to select a trajectory based on quantum probabilities, SH contains

too many unsatisfactory ad hoc assumptions and fails to keep track of the

coherence between various quantum states. Other approaches have similar

shortcomings.9 A possible, but not practical, solution is to describe the sys-

tem as an ensemble of surface-hopping classical trajectory segments, keeping

simultaneously track of the trajectories belonging to each of the states that

mix into the original state by DME (Ben-Nun and Martinez, 1998; Kapral

and Ciccotti, 1999; Nielsen et al., 2000). The e¬ect of quantum decoherence

was addressed by Prezhdo and Rossky (1997a, 1997b).

A promising approach is the use of Bohmian particle dynamics (quantum

hydrodynamics, see section 3.4 on page 64). The quantum particle is sam-

pled from the initial distribution ψ 2 (r, 0) and moves as a classical particle

in an e¬ective potential that includes the quantum potential Q (see (3.93)

on page 68). The latter is determined by the wave function, which can be

computed either by DME in the usual way or by evaluating the gradient

of the density of trajectories. Thus the quantum particle follows a single

well-de¬ned trajectory, di¬erent for each initial sample; the branching oc-

curs automatically as a distribution of trajectories. The back reaction now is

9 There is a vast literature on non-adiabatic semiclassical dynamics, which will not be reviewed

here. The reader may wish to consult Meyer and Miller(1979), Webster et al. (1991), Laria et

al. (1992), Bala et al. (1994), Billing (1994), Hammes-Schi¬er (1996), Sun and Miller (1997),

M¨ller and Stock (1998).

u

136 Dynamics of mixed quantum/classical systems

dependent on the Bohmian particle position. The method has been applied

by Lepreore and Wyatt (1999),10 Gindensperger et al. (2000, 2002, 2004)

and Prezhdo and Brooksby (2001). The latter two sets of authors do not

seem to agree on whether the quantum potential should be included in the

classical back reaction. There are still problems with energy conservation11

(which is only expected to be obeyed when averaged over a complete en-

semble), and it is not quite clear how the Bohmian particle approach should

be implemented when the quantum subsystem concerns some generalized

coordinates rather than particles. Although at present the method cannot

be considered quite mature, it is likely to be the best overall solution for the

trajectory-based simulation of mixed quantum/classical dynamics.

Mixed quantum-classical problems may often be approximated in a practi-

cal way, when the details of the crossing event are irrelevant for the questions

asked. Consider, for example, a chromophore in a protein that uses light

to change from one conformation to another. Such systems are common

in vision (rhodopsin), in energy transformation (bacteriorhodopsin) and in

biological signaling processes. After excitation of the chromophore, the sys-

tem evolves on the excited-state potential surface (generally going downhill

along a dihedral angle from a trans state towards a cis state), until it reaches

the conical intersection between excited state and ground state.12 It then

crosses over from the excited state to either a trans or a cis ground state,

proceeding down-hill. The uphill continuation of the excited state is unlikely,

as its probability in the damped, di¬usion-like multidimensional motion is

very small; if it happens it will lead to re-entry into the conical intersection

and can be disregarded as an irrelevant process. The fate of the system

upon leaving the conical intersection is of much more interest than the de-

tails during the crossing event. Groenhof et al. (2004), in a simulation study

of the events after light absorption in the bacterial signaling protein PYP

(photoactive yellow protein), used a simple approach with a single surface

hop from excited to ground state after the excited state was found to cross

the ground state. The potential surface of ground and excited states were

determined by CASSCF (complete active space SCF, see page 105). After

each time step the con¬guration-interaction vector is determined and it is

seen whether the system crossing has occurred. If it has, the classical forces

are switched from the excited to the ground state and the system continues

10 See also Wyatt (2005).

11 See comment by Salcedo (2003).

12 A conical intersection is the multidimensional analog of the two-state crossing, as depicted in

Fig. 5.2. The main coordinate in retinal-like chromophores is a dihedral angle or a combination

of dihedral angles, but there are other motions, called skeletal deformations, that aid in reaching

the intersection.

Exercises 137

on one of the descending branches. Thus the crossing details are disre-

garded, and hopping between states before or after they cross are neglected.

Still, the proper system evolution is obtained with computed quantum yield

(fraction of successful evolutions to the cis state) close to the experimentally

observed one.

Exercises

5.1 Derive the adiabatic wave functions and energies for two states with

diabatic energy di¬erence ”E0 and o¬-diagonal real coupling ener-

gies C (see (5.1)).

5.2 How do the results di¬er from those of the previous exercise when

the o¬-diagonal coupling is purely imaginary?

Show that(5.45) and (5.46) imply that the quantity 4ρ12 ρ— + z 2 is

5.3 12

a constant of the motion.

5.4 Prove that (5.77) follows from (5.76).

6

Molecular dynamics

6.1 Introduction

In this chapter we consider the motion of nuclei in the classical limit. The

laws of classical mechanics apply, and the nuclei move in a conservative

potential, determined by the electrons in the Born“Oppenheimer approxi-

mation. The electrons are assumed to be in the ground state, and the energy

is the ground-state solution of the time-independent Schr¨dinger equation,

o

with all nuclear positions as parameters. This is similar to the assumptions

of Car“Parrinello dynamics (see Section 6.3.1), but the derivation of the

potential on the ¬‚y by quantum-mechanical methods is far too compute-

intensive to be useful in general. In order to be able to treat large systems