Baer M., Billing G.D. (eds.) - The role of degenerate states in chemistry (Adv.Chem.Phys. special issue, Wiley, 2002)

applying direct molecular dynamics to non-adiabatic systems 383

Importantly, this term is a derivative (nonlocal) operator on the nuclear coordinate space.

The matrix G can in fact be expressed in terms of F (see Appendix A),

X

Gij ¼ $ Fij þ Fik Fkj ð57Þ k

And using this equation the nuclear Schro¨dinger equation can be written in matrix form [54,179]

	h2	ð58Þ
	2 ð$ þ FÞ2þV v ¼ ih qqt v

The matrix of vectors F is thus the defining quantity, and is called the nonadiabatic coupling matrix. It gives the strength (and direction) of the coupling between the nuclear functions associated with the adiabatic electronic states.

The elements of these vectors can be evaluated using an off-diagonal form of the Hellmann–Feynmann theorem, which in Cartesian coordinates, xa, is

a	1		1		^
					ad qHel ad
Fij					hci j qxa jcj i			ð59Þ
	¼ pMa Vj			Vi

Nonscaled coordinates are used here to explicitely include the mass to show that the coupling is modulated by two factors. The first is the mass associated with the coordinate, Ma (atomic mass in Cartesian coordinates, reduced mass in normal mode coordinates, etc.), and the larger the mass the smaller the coupling. This is the basis of the justification for the BO approximation: The mass of the electron is so much smaller than the mass of the nuclei that the motion of the electrons is effectively independent of the nuclear motion, and the electrons instantaneously adjust to the nuclear geometry. The second factor, however, depends inversely on the energy gap between the adiabatic surfaces. This will overwhelm the mass factor as the surfaces approach one another, until at a degenerate point the coupling is infinitely large.

As written, Eq. (52) depends on all the (infinite number of) adiabatic electronic states. Fortunately, the inverse dependence of the coupling strength on energy separation means that it is possible to separate the complete set of states into manifolds that effectively do not interact with one another. In particular, Baer has recently shown [54] that Eq. (57), and hence Eq. (58) also holds in the subset of mutually coupled states. This finding has important consequences for the use of diabatic states explored below.

Choosing a basis set for the nuclear functions ffag allows us to write Eq. (52) in a matrix form, similar to Eq. (10) for the single-surface case, now as

384	g. a. worth and m. a. robb

a matrix of matrices

0

H11

B@ H21

...

H22	10 v2 1					ih0 v2 1			60
H12	.	. . CB	v1	C	¼	B	v1	C	ð Þ
...			...				...
		A@		A		@		A

where the indexes of Hij and vi, which relate to the adiabatic electronic states, run over the states included in the manifold.

The supermatrix notation emphasizes the structure of the problem. Each diagonal operator drives a wavepacket, just as in the adiabatic case of Eq. (10), but here the motion of the wavepackets in different adiabatic states is mixed by the off-diagonal non-adiabatic operators. In practice, a single matrix is built for the operator, and a single vector for the wavepacket. The operator matrix elements in the basis set ffag are

ðHiiÞab ¼

fa

^

^

ii

fb

ð61Þ

TN þ Vi

¼

ð Þ

Hij

f

f

62

ab

a

^ij

b

which are arranged in the blocks. All the methods mentioned in Section II.B for wavepacket dynamics can then be used.

B.The Diabatic Picture

The adiabatic picture is the standard one in quantum chemistry for the reason that, not only is it mathematically well defined, but it is also that used in ab initio calculations, which solve the electronic Hamiltonian at a particular nuclear geometry. To see the effects of vibronic coupling on the potential energy surfaces one must move to what is called a diabatic representation [1,65,180, 181].

In a diabatic representation, the electronic wave functions are no longer eigenfunctions of the electronic Hamiltonian. The aim is instead that the functions are so chosen that the (nonlocal) non-adiabatic coupling operator matrix, ^ in Eq. (52), vanishes, and the couplings are represented by (local) potential operators. The nuclear Schro¨dinger equation is then written

^ jw i þ X jw i ¼ q jw i ð Þ

Tn i Wij j ih q i 63

j

t

where W is the new (nondiagonal) potential matrix, and the coupling between states is now achieved by the off-diagonal elements of this matrix. The adiabatic surfaces are the eigenfunctions of W.

applying direct molecular dynamics to non-adiabatic systems 385

The diabatic electronic functions are related to the adiabatic functions by unitary transformations at each point in coordinate space

/ðRÞ ¼ SðRÞwðRÞ

ð64Þ

From the time-dependent Schro¨dinger equation in the matrix form of Eq. (58), it can be shown [54] that if SðRÞ is chosen so that

Syð$ þ FÞS ¼ 0

ð65Þ

the basis set transformation changes the operators in Eq. (58) to

h2

Sy

F

2S

h2

2

T

66

2

ð$ þ

Þ

¼

2

$

¼

ð

Þ

67

SyVS

¼

W

ð

Þ

and the diabatic representation is rigorously equivalent to the adiabatic representation in the subspace of the coupled states. Baer [53,54] has obtained solutions to this equation, and analyzed the validity of the transformation in regions where the non-adiabatic coupling becomes singular. Note that derivative operators in both terms act on SðRÞ, and so this is not a local transformation. Note also that the diabatic basis is only defined up to a constant rotation. As a result, it is possible to select a point at which the diabatic and adiabatic functions are identical. This simplifies various mathematical manipulations.

Assuming that the diabatic space can be truncated to the same size as the adiabatic space, Eqs. (64) and (65) clearly define the relationship between the two representations, and methods have been developed to obtain the transformation matrices directly. These include the line integral method of Baer [53,54] and the block diagonalization method of Pacher et al. [179]. Failure of the truncation assumption, however, leads to possibly important nonremovable derivative couplings remaining in the diabatic basis [55,182].

Difficulties in obtaining the non-adiabatic coupling elements for polyatomic molecules have lead to the development of alternative approaches to provide the diabatic representation, typically using states that are smooth in a molecular property [183]. Although there is no formal justification for this approach, it seems to work well in practice. Properties used include the dipole moment [184], or retention of the configurational character from an MCSCF wave function [185], or maximization of the overlap between wave functions at neighbouring sites [186]. It has also been shown that the CASSCF method provides a good framework for the definition of diabatic states [187]. A simple scheme that removes the leading terms of the non-adiabatic coupling matrices

386	g. a. worth and m. a. robb

using information from the adiabatic surfaces only has also been recently proposed [188]. For a fuller account of ways to construct diabatic states see [1].

Despite the difficulties in obtaining diabatic states they provide an extremely useful picture for many descriptive purposes. They are in many ways the natural choice for dynamics calculations as the kinetic energy operator is diagonal in this basis, and the singularities associated with the non-adiabatic operator where the adiabatic surfaces meet are not present. In principle, as the electronic basis set is only weakly dependent on R, the electronic character of a state is preserved. The range of properties used to define diabatic states shows the sort of properties that are conserved within them. A typical example is in electrontransfer theory, which uses smooth diabatic states to define the donor–acceptor and charge-transfer states. A more important example in photochemistry is that photoexcitation in the Condon approximation should be modelled as taking place vertically to a single diabatic state, as in this picture the transition matrix element is relatively constant with respect to the nuclear geometry [1].

C.Conical Intersections

For a two-state system, the eigenfunctions of the diabatic potential matrix of Eq. (63) in terms of its elements are

¼	1		þ		Þ q	ð Þ
V	2	ðW11		W22	W2 þ W122	68

where W ¼ 12 ðW22 W11Þ. The functions Vþ and V are the adiabatic PES, and they will meet when

W ¼ 0	ð69Þ
W12 ¼ 0	ð70Þ

In Section III.D, we shall investigate when this happens. For the moment, imagine that we are at a point of degeneracy. To find out the topology of the adiabatic PES around this point, the diabatic potential matrix elements can be expressed by a first order Taylor expansion.

Setting the diabatic basis equal to the adiabatic basis at the degenerate point, R0, the expansions can be written in vector notation as

W ¼ x1	Q	ð71Þ
W12 ¼ x2	Q	ð72Þ

where Q is the vector of nuclear displacements away from the intersection. Note that the constants in both expansions are zero due to the adiabatic–diabatic

applying direct molecular dynamics to non-adiabatic systems			387
correspondence at R0. The vectors that form the first-order coefficients are
x1a ¼	q	W	ð73Þ
	qQa
x2a ¼	q	W12	ð74Þ
	qQa

evaluated at R0. Also due to the adiabatic–diabatic correspondence at this point, these can be written as the gradient difference (GD) vector

x1a ¼

q

ðVþ V Þ

ð75Þ

qQa

and the derivative coupling (DC) vector

^

qHel

x2a

¼ ciad

qQa

cjad

ð76Þ

This latter relationship is obtained by evaluating

q

qcj

qc

^

qHel

i

qQa

hcijH^eljcji ¼

ci

H^el

qQa

þ

qQa

H^el

cj

þ hcij

qQa

jcji ð77Þ

at R0. In the diabatic basis, the first two terms on the right are zero. Due once again to the adiabatic–diabatic correspondence at R0 the third term on the righthand side is

^				ad
		ad		qcj
qHel
hcij	qQa	jcji ¼ ðVj ViÞ ci		qQa		ð78Þ

[see Eqs. (130)–(132) in Appendix A]. And so the expansion coefficient vector lies in the same direction as the adiabatic coupling vector.

The major features of the PES around the degenerate point can now be easily analysed if we write the vector Q in the basis of ðx1; x2; . . .Þ where the unspecified n 2 basis vectors are orthogonal to the ðx1; x2Þ plane, which is called the branching space. First, moving in the n 2-dimensional space orthogonal to the branching space the degeneracy is retained. Second, moving in the plane of the branching space, the degeneracy will be lifted. Ignoring the term 12 ðW11 þ W22Þ, which is the same for both surfaces, the adiabatic PES have the form

q
V x12 þ x22	ð79Þ

388	g. a. worth and m. a. robb

Figure 5. Sketch of a conical intersection. The vectors x1 and x2 are the GD and DC respectively, that lift the degeneracy of the two adiabatic surfaces. The plane containing these vectors is known as the branching space.

where x1; x2 are the components of Q along the respective vectors, and so the topology of the surfaces around the degenerate point is that of a double cone. Hence, this is called a conical intersection. A sketch of such a point is given in Figure 5.

Conical intersections can be broadly classified in two topological types: ‘‘peaked’’ and ‘‘sloped’’ [189]. These are sketched in Figure 6. The peaked case is the classical theoretical model from Jahn–Teller and other systems where the minima in the lower surface are either side of the intersection point. As indicated, the dynamics of a system through such an intersection would be expected to move fast from the upper to lower adiabatic surfaces, and not return. In contrast, the sloped form occurs when both states have minima that lie on the same side of the intersection. Here, after crossing from the upper to lower surfaces, recrossing is very likely before relaxation to the ground-state minimum can occur.

A final point to be made concerns the symmetry of the molecular system. The branching space vectors in Eqs. (75) and (76) can be obtained by evaluating the derivatives of matrix elements in the adiabatic basis

q	ad			ad
	ci	Hel	cj			ð80Þ
qQa

with i ¼ j required for x1 and i ¼6 j for x2. These elements are only nonzero if the product of symmetries of the adiabatic functions i, j, and the symmetry of the nuclear coordinate, Qa contains the totally symmetric irrep

i Qa j Ag

ð81Þ

applying direct molecular dynamics to non-adiabatic systems 389

Figure 6. Two-dimensional (top) and 3D (bottom) representations of a peaked (a) and sloped (b) conical intersection topology. There are two directions that lift the degeneracy: the GD and the DC. The top figures have energy plotted against the DC while the bottom figures represent the energy plotted in the space of both the GD and DC vectors. At a peaked intersection, as shown at the bottom of (a), the probability of recrossing the conical intersection should be small whereas in the case of a sloped intersection [bottom of (b)], this possibility should be high. [Reproduced from [84] courtesy of Elsevier Publishers.]

If the symmetries of the two adiabatic functions are different at R0, then only a nuclear coordinate of appropriate symmetry can couple the PES, according to the point group of the nuclear configuration. Thus if Q are, for example, normal coordinates, x1 will only span the space of the totally symmetric nuclear coordinates, while x2 will have nonzero elements only for modes of the correct symmetry.

D. The Vibronic-Coupling Model Hamiltonian

A more general description of the effects of vibronic coupling can be made using the model Hamiltonian developed by Ko¨ppel, Domcke and Cederbaum [65]. The basic idea is the same as that used in Section III.C, that is to assume a quasidiabatic representation, and to develop a Hamiltonian in this picture. It is a useful model, providing a simple yet accurate analytical expression for the coupled PES manifold, and identifying the modes essential for the non-adiabatic effects. As a result it can be used for comparing how well different dynamics methods perform for non-adiabatic systems. It has, for example, been used to perform benchmark full-dimensional (24-mode) quantum dynamics calculations

390	g. a. worth and m. a. robb

on the photoabsorption of pyrazine into the S2/S1 manifold [109]. Here it will be used to demonstrate the effect of a conical intersection on a systems dynamics using the butatriene radical cation as an example.

The Hamiltonian again has the basic form of Eq. (63). The system is described by the nuclear coordinates, Q, which are relative to a suitable nuclear configuration Q0. In contrast to Section III.C, this may be any point in configuration space. As a diabatic representation has been assumed, the kinetic energy operator matrix, T, is diagonal with elements

X
f	h2 q2
Tii ¼	2 qQa2	ð82Þ
a¼1

The potential matrix elements are then obtained by making Taylor expansions

around Q0, using suitable zero-order diabatic potential energy functions,

Vað0ÞðQÞ.

	X
	f	q
0
Wij Við Þdij ¼ hcijHeljcji þ	a¼1	qQa	hcijHeljcjiQa þ	ð83Þ

where the integrals and derivatives are evaluated at the point Q0. The diabatic functions are again taken to be equal to the adiabatic functions at Q0, and so

f
Wii ¼ Við0Þ þ Ei þ kaðiÞQa þ	ð84Þ
a¼1
X
f
Wij ¼ laðijÞQa þ	ð85Þ
a¼1
X

The model is that the ground-state PES is first altered by the electronic excitations (on-diagonal coupling leads to a change in equilibrium geometry and frequency), and these smooth diabatic states are then further altered by vibronic (off-diagonal) coupling.

The eigenvalues of this matrix have the form of Eq. (68), but this time the matrix elements are given by Eqs. (84) and (85). The symmetry arguments used to determine which nuclear modes couple the states, Eq. (81), now play a crucial role in the model. Thus the linear expansion coefficients are only nonzero if the products of symmetries of the electronic states at Q0 and the relevant nuclear mode contain the totally symmetric irrep. As a result, on-diagonal matrix elements are only nonzero for totally symmetric nuclear coordinates and, if the electronic states have different symmetry, the off-diagonal elements will only

applying direct molecular dynamics to non-adiabatic systems 391

be nonzero for a suitable nonsymmetric mode, as given by the product of the electronic-state symmetries.

For states of different symmetry, to first order the terms W and W12 are independent. When they both go to zero, there is a conical intersection. To connect this to Section III.C, take Q0 to be at the conical intersection. The gradient difference vector in Eq. (75) is then a linear combination of the symmetric modes, while the non-adiabatic coupling vector in Eq. (76) is a linear combination of the appropriate nonsymmetric modes. States of the same symmetry may also form a conical intersection. In this case it is, however, not possible to say a priori which modes are responsible for the coupling. All totally symmetric modes may couple onor off-diagonal, and the magnitudes of the coupling determine the topology.

A conical intersection needs at least two nuclear degrees of freedom to form. In a 1D system states of different symmetry will cross as Wij ¼ 0 for i ¼6 j and so when Wii ¼ 0 the surfaces are degenerate. There is, however, no coupling between the states. States of the same symmetry in contrast cannot cross, as both Wij and Wii are nonzero and so the square root in Eq. (68) is always nonzero. This is the basis of the well-known non-crossing rule.

If the states are degenerate rather than of different symmetry, the model Hamiltonian becomes the Jahn–Teller model Hamiltonian. For example, in many point groups E E E and so a doubly degenerate electronic state can interact with a doubly degenerate vibrational mode. In this, the E E Jahn– Teller effect the first-order Hamiltonian is then [65]

Qx	Qy		ð86Þ
H ¼ ðT þ V0Þ1 þ k Qy	Qx

where x; y denote the two components of the degenerate vibrational mode, 1 is the 2 2 unit matrix, and the zero-order Hamiltonian

X		oi		q2	2
T þ V0 ¼ i	x;y 2			qQi2	þ Qi		ð87Þ
¼

is the unperturbed harmonic state (written here in mass-frequency scaled coordinates). This model results in the splitting of the degeneracy to form a symmetrical moat around a central conical intersection.

The Hamiltonian provides a suitable analytic form that can be fitted to the adiabatic surfaces obtained from quantum chemical calculations. As a simple example we take the butatriene molecule. In its neutral ground state it is a planar

molecule with D2h symmetry. The lowest two states of the radical cation,

~2 responsible for the first two bands in the photoelectron spectrum, are X B2g and

~2B2u. The vibronic coupling model Hamiltonian is set up using the ground-state

A

392	g. a. worth and m. a. robb

~2	~2
Figure 7. The PES of the X	B2g and A B2u states of the butatriene radical cation. (a) Diabatic

surfaces. (b) Adiabatic surfaces. The surfaces are obtained as eigenfucations of the vibronic coupling model Hamiltonain that fitted to reproduce quantum chemical calculations. The coordinates are shown in Figure 1c. See Section III. D for further details.

normal modes, of which there are 15, expanding around the neutral equilibrium geometry using the harmonic ground state as a zero-order surface. Taking symmetry into account, to first order these states can only be coupled by a nuclear degree of freedom with Au symmetry, of which there is only one, the