0The Journal of Chemical Physics
Automatic generation of active coordinates for quantum dynamics calculations: Application to the dynamics of benzene photochemistry
Benjamin Lasorne, Fabrizio Sicilia, Michael J. Bearpark, Michael A. Robb, Graham A. Worth, and Lluis Blancafort
Citation: The Journal of Chemical Physics 128, 124307 (2008); doi: 10.1063/1.2839607 View online: http://dx.doi.org/10.1063/12839607
View Table of Contents: http://scitation.aip.org/content/aip/journal/jcp/128/12?ver=pdfcov Published by the AIP Publishing
Articles you may be interested in
Decay dynamics of a,p-carboxylic methyl esters (CH3OCOCH:CHR) in the lower-lying excited states—Resonance Raman and complete active space self-consistent field calculation study J. Chem. Phys. 141, 134312 (2014); 10.1063/1.4896999
Global sampling of the photochemical reaction paths of bromoform by ultrafast deep-UV through near-IR transient absorption and ab initio multiconfigurational calculations J. Chem. Phys. 138, 124501 (2013); 10.1063/1.4789268
Experimental and theoretical studies of the decomposition of new imidazole based energetic materials: Model systems
J. Chem. Phys. 137, 114303 (2012); 10.1063/1.4752654
Quantum dynamics study of fulvene double bond photoisomerization: The role of intramolecular vibrational
energy redistribution and excitation energy
J. Chem. Phys. 135, 134303 (2011); 10.1063/1.3643767
Real-time propagation time-dependent density functional theory study on the ring-opening transformation of the
photoexcited crystalline benzene
J. Chem. Phys. 124, 124507 (2006); 10.1063/1.2181139
Launching in 2016!
he future of applied photonics research is here
Photonics
Automatic generation of active coordinates for quantum dynamics calculations: Application to the dynamics of benzene photochemistry
Benjamin Lasorne,1 Fabrizio Sicilia,1 Michael J. Bearpark,1 Michael A. Robb,1,a) Graham A. Worth,2 and Lluls Blancafort3
1Department of Chemistry, Imperial College London, South Kensington, London SW7 2AZ, United Kingdom 2School of Chemistry, University of Birmingham, Edgbaston, Birmingham B15 2TT, United Kingdom 3Institut de Química Computational and Departament de Química, Universitat de Girona, E-17071 Girona, Spain
(Received 28 September 2007; accepted 10 January 2008; published online 26 March 2008)
A new practical method to generate a subspace of active coordinates for quantum dynamics calculations is presented. These reduced coordinates are obtained as the normal modes of an analytical quadratic representation of the energy difference between excited and ground states within the complete active space self-consistent field method. At the Franck-Condon point, the largest negative eigenvalues of this Hessian correspond to the photoactive modes: those that reduce the energy difference and lead to the conical intersection; eigenvalues close to 0 correspond to bath modes, while modes with large positive eigenvalues are photoinactive vibrations, which increase the energy difference. The efficacy of quantum dynamics run in the subspace of the photoactive modes is illustrated with the photochemistry of benzene, where theoretical simulations are designed to assist optimal control experiments. © 2008 American Institute of Physics. [DOI: 10.1063/1.2839607]
I. INTRODUCTION
The first step in elucidating a photochemical mechanism computationally is the topographical characterization of the potential energy surfaces (PESs). This involves locating critical points on ground and excited PESs as well as any seams of conical intersections between states.1 Mechanisms of thermal reactions are traditionally described with the concept of a static reaction coordinate along a minimum energy path. In contrast, nonadiabatic mechanisms involve ultrafast processes and depend on changes in the kinetic energy in the energized system. A time-dependent picture is, thus, essential and quantum dynamics is preferable for the description of nonadiabatic effects and coherence preservation during the time scale of the photochemical event.
The cost of quantum dynamics simulations grows quickly with the number of nuclear degrees of freedom (3N-6 if N is the number of atoms and the molecule is not linear). Thus, quantum dynamics simulations are often performed within a subspace of active coordinates (see, e.g., Refs. 3-7). In this paper, we describe a method which enables the a priori selection of these important coordinates for a photochemical reaction using an analytic simultaneous representation of both ground and excited states correct to second order in the energy (i.e., using analytic gradients and second derivatives of the energy difference8-10). The efficacy of the method will be illustrated using the channel 3 photochemistry of benzene,11,12 where the energy difference was analyzed at the S0 minimum.
In any nonadiabatic photochemical reaction, there are two special coordinates:1,2 x1, along the gradient difference,
a)Electronic mail: mike.robb@imperial. 0021-9606/2008/128( 12)/124307/10/$23.00
and x2, along the interstate-coupling vector. The S1 and S0 surfaces for benzene are schematically plotted in this space in Fig. 1. If the energy-difference Hessian is evaluated at the ground state PES minimum and projected onto the space orthogonal to the plane (x1,x2), one has the information about the curvature of the energy difference in this subspace. Diagonalization of this projected Hessian gives a set of vi-brational modes which, as we shall show, can be classified as (a) photoactive modes that decrease the energy gap and can lead to the conical intersection (and which we will use to yield a reduced set of coordinates in quantum dynamics), (b) photoinactive modes that increase the energy gap and lead away from the conical intersection, and (c) bath modes where the two surfaces are parallel (and which can be neglected in a quantum dynamics approach).
For the ^ S0 photochemistry of benzene, the so-called channel 3 process represents the well-known decay route along which fluorescence is quenched above a vibrational excess of 3000 cm-1 (see, e.g., Refs. 11 and 12, and references therein). Our purpose is to understand better the principles that control the nonradiative decay by changing the way the S1 /S0 intersection seam is accessed. For dynamics, we have used the direct dynamics variational multiconfigu-ration Gaussian wavepacket (DD-vMCG) method13-16 to simulate radiationless decay. This is a quantum dynamics algorithm based on the propagation of coupled Gaussian functions. Each function follows a quantum trajectory along
which the PESs are calculated on the fly by a quantum chemistry program. This method has been used within the reduced subspace of vibrational modes selected using the energy-difference Hessian.
The plan of this paper is as follows. The conceptual aspects of the theoretical development are discussed in Sec. II.
© 2008 American Institute of Physics
128, 124307-1
TABLE I. The 30 normal modes of S0 benzene labeled following the Wilson scheme of frequency numbering.
FIG. 1. (Color online) S1 ^ So photochemical reaction path from benzene to prefulvene plotted along the gradient-difference coordinate (xi) and the interstate-coupling-vector coordinate (X2).
The mathematical details are then provided in a subsequent section (Sec. III). The reader should be able to skip such details and proceed directly to the results (Sec. IV) on a first reading. Concluding remarks are gathered in Sec. V.
II. CONCEPTUAL DEVELOPMENT
A photochemical reaction path involves (at least) two electronic states, as shown in Fig. 2. The molecule is initially in the ground state around the equilibrium geometry (GSMin in Fig. 2). After the absorption of light, the system is excited to the Franck-Condon point (FC in Fig. 2) and starts off following the steepest descent path (driving force). The energy gap decreases further until the system reaches a region of conical intersection1 (CoIn in Fig. 2) and decays to form the products in the ground state (GSProd in Fig. 2). In this context, the electronic energy gap is the quantity that dictates the photoreactivity. For the photoproduct to be formed, the system must follow a pathway that makes the energy difference diminish. An efficient pathway will, thus, be characterized by a rapid decrease of the energy gap.
In the case of S1 ^ S0 benzene photochemistry, the photochemical event is illustrated in Fig. 1, where the S0 and S1 surfaces are plotted in the space spanned by the coordinates along the gradient difference (x1) and interstate-coupling vector (x2). Such directions are traditionally computed at the
Mode Symmetry Dominant motion Frequency' (cm-1)
1 aig t CC stretching (breathing) 1042
2 a1g s CH stretching 3390
3 a2g / HCC bending 1500
4 b2g S CCCC out-of-plane ring torsion (chair) 724
5 b2g y HCCC out-of-plane wagging 1035
6 e2g a CCC bending (rectangular/rhomboidal) 654
7 e2g s CH stretching 3360
8 e2g t CC stretching 1739
9 e2g 3 HCC bending 1275
10 e1g y HCCC out-of-plane wagging 874
11 a2u y HCCC out-of-plane wagging 711
12 b1u a CCC bending (triangular) 1099
13 b1u s CH stretching 3350
14 b2u 3 HCC bending 1339
15 b2u t CC stretching + a CCC bending (Kékulé) 1185
16 e2u S CCCC out-of-plane ring torsion 431
(boat/twist)
17 e2u y HCCC out-of-plane wagging 991
18 e1u t CC stretching + a CCC bending 1114
(trapezoidal)
19 e1u 3 HCC bending 1630
20 e1u s CH stretching 3379
FIG. 2. (Color online) Schematic picture of a photochemical event. The energies of the electronic ground state and excited state are plotted against a photochemical reaction coordinate. The arrows describe the mechanism leading from the Franck-Condon point (FC) to the ground state product (GSProd) passing through a conical intersection (CoIn).
aThe frequencies were calculated at the S0 equilibrium geometry of benzene at the CASSCF(6,6)/6-31G* level. No scaling was applied for reasons of consistency with subsequent direct dynamics calculations at the same ab initio level. The relative error with respect to experimental data does not exceed 10%.
conical intersection point, where they define the branching
plane. Although the energy gap is not zero, the same vectors can be defined at the Franck-Condon point. Here, x1 and x2 span a "pseudobranching plane," which may be different from the branching plane at the conical intersection point. In the following, we describe x1 and x2 at both points by comparing them to a reference set of vibrational modes computed at the S0 equilibrium geometry of benzene (see Table I, where the Wilson scheme of frequency numbering18 is used).
For benzene, the interstate-coupling vector remains substantially unchanged at both the conical intersection and Franck-Condon points, i.e., mode 15 (Fig. 3). However, the gradient difference varies considerably. At the conical intersection, this vector is a combination of modes 1, 4, and 16x (Fig. 3), whereas at the Franck-Condon point, it corresponds to mode 1 only (Fig. 3). We show here that this variation can be predicted by analyzing the representation of the energy difference through first and second orders at the Franck-Condon geometry.
At the Franck-Condon point, the first-order representation of the energy difference corresponds to the gradient difference (mode 1 in Fig. 3). The second-order contributions are computed by diagonalizing the energy-difference Hessian projected onto the space orthogonal to the pseudobranching plane. The resulting eigenvectors can be classified according to the magnitude and sign of the corresponding eigenvalues. As shown in Fig. 4, three types of modes must be distinguished.
FIG. 3. (Color online) The nine dominant motions for the photochemistry of benzene (the labels refer to the most similar normal modes of S0 benzene listed in Table I).
Three one-dimensional cuts of two PESs are plotted in the space orthogonal to the pseudobranching plane in Fig. 4 for a generic case. It should be emphasized that the vectors defining the pseudobranching plane, as well as those perpendicular to it, are defined at the Franck-Condon point. Thus, the gradient difference simply is the S1 gradient. In Fig. 4, the minima on both the ground and excited state PESs are at the same value of the coordinates, since the S0 gradient is zero and the S1 gradient is orthogonal to these directions at the Franck-Condon point.
FIG. 5. (Color online) Schematic representation of two trajectories starting from the benzene S1 Frank-Condon point (FC). Trajectory (a) oscillates within the S1 minimum valley, whereas trajectory (b) reaches the ground state product passing through a conical intersection point (CoIn).
In Fig. 4, we show that three different types of motions can be distinguished according to the eigenvalues obtained from the energy-difference (excited state minus ground state) Hessian diagonalization. The first class of modes makes the energy difference decrease, i.e., negative eigenvalues, and we call them photoactive modes [Fig. 4(a)]. The modes along which the energy difference increases, i.e., positive eigenvalues, are called photoinactive modes [Fig. 4(b)]. Finally, those eigenvectors where the energy difference does not significantly change, i.e., almost zero eigenvalues, are called bath modes [Fig. 4(c)]. For photochemistry to take place, the system must reduce the energy difference to access the conical intersection. The most important coordinates to describe a photochemical event are, thus, the gradient difference, the interstate-coupling vector, and the additional photoactive modes. Quantitative results of this analysis applied to benzene will be presented and discussed in Sec. IV.
To explore the photochemical behavior of benzene, dynamics simulations were performed with modes 1 and 15, i.e., pseudobranching-space modes, and seven photochemical modes, viz., modes 4, 6 (degenerate pair), 16 (degenerate pair), and 18 (degenerate pair), depicted in Fig. 3. The time evolution of the system is represented in Fig. 5 in a schematic way. Two coordinates are used: the gradient difference
at the Franck-Condon point (Qi) and a generic photoactive mode (Qpa). At the Franck-Condon point (D6h benzene), the gradient difference is identical to the Si gradient. The driving force (opposite of the gradient) leads to the Sl minimum (stretched D6h benzene). The first motion to be activated is, thus, a breathing motion (totally symmetric stretching) that can keep oscillating until symmetry is broken (or until the system eventually fluoresces). Trajectory (a) illustrates such a case with no initial momentum. If some initial momentum is added along the photoactive mode Qpa, as illustrated by trajectory (b), the system can escape this well and reach a point on the seam of conical intersection between the two surfaces. A nonradiative decay will happen if the system crosses the seam and keeps evolving on S0. As further discussed in Sec. IV, stimulating the photochemical modes identified with our approach proved to be a very efficient way for the molecule to reach geometries where the energy gap is small enough to induce strong nonadiabatic transitions from Si to So.
In summary, efficient photochemical reaction paths can be predicted by analyzing the local properties of the energy difference between Si and S0 at the Franck-Condon point. This approach gives the modes that must be stimulated in order to enhance nonradiative decay. In addition, it provides an objective criterion to select active coordinates and run quantum dynamics in reduced dimensionality, as we show for benzene in Sec. IV.
III. THEORETICAL DEVELOPMENT
The analysis of the photochemical activity of nuclear coordinates is now presented in more detail. An integral part
of our approach'
8-10,19
to the general quadratic representation
of conical intersections2,20-28 is an analytic representation of the energy difference between the states. The approach in this paper is based on a second-order analytic representation of the local two-state potential energy matrix to calculate a quadratic expansion of the energy difference around the ground state equilibrium geometry (i.e., Franck-Condon point in the excited state), where the energy difference is not zero. We now recall the features and concepts of this analysis. They are further applied to the case of benzene where D6h symmetry is explicitly taken into account.
The theoretical concepts used in the second-order analysis of the energy difference are illustrated below in the ideal case of a two-state problem. Our purpose is to clarify the meaning of the terms that appear in the second-order expansion of the energy difference as a function of the nuclear coordinates. Practical details about how the corresponding quantum chemistry calculations are carried out can be found elsewhere (see Refs. 8-10 and references therein). The pair of adiabatic states |S0; Q) and |Si; Q) (parametrized by the nuclear geometry Q) are taken as known at a reference geometry Q = Qo, chosen as the S0 equilibrium geometry. In the representation formed by the pair of adiabatic states, the matrix of the clamped-nucleus Hamiltonian operator H/(Q) is diagonal at Q = Q0,
H(Qo) :
S ; QoH(Qo) So ; Qo> S ; Qo|H(Qo) |S ; Qo> ,<Si;Qo|HH(Qo)|So;Qo> <Si;Qo|H(Qo)|Si;Qo> J
Vo(Qo) o o Vi(Qo)J
However, it is no longer diagonal when using the same states |S0; Q0) and |Si; Q0) at a displaced geometry Q = Q0+ SQ, where H (Q) = H(Qo) + SH,
H(Qo + SQ)--
<So;Qo|HH(Qo) + SH|So;Qo> <So;Qo|#(Qo) + S?Si;Qo> <Si;Qo|H(Qo) + SH|So;Qo> <Si;Qo|#(Qo) + SH|Si;Qo> J
Vo(Qo) o Vi(Qo)J
SHoo SHoi LSHio SHn\
The mass-weighted geometrical displacement SQ generates a finite variation of the off-diagonal elements SHi0=SH0i (real-valued electronic wavefunctions) because the adiabatic states frozen at Q = Q0 define trivial diabatic states when Q varies in the Hamiltonian operator. Such states are often re-
ferred to as "crude adiabatic."
In the two-state approximation assumed here, the two electronic wavefunctions, thus, form a complete basis set. Any discrepancy can formally be attributed to neglecting couplings with more energetic states. In this representation, the positive difference between the adiabatic potential energies varies with SQ according to
AV(Qo + SQ) = V[A V(Qo) + S(Hn - Hoo)]2 + 4SHoi2,
where AV(Q) = Vi(Q)- Vo(Q). We now introduce two functions of the coordinates Q
fi(Q) = Hii(Q) - Hoo(Q)
= <Si; Qo|HH(Q)|Si; Qo> - <So; Qo|#(Q)|So; Qo>,
f2(Q) = Hoi (Q) = <So ; Qo|H (Q) |Si ; Qo>. (4)
Note that a conical intersection would correspond to fi(Qo) = o andf2(Qo) = o. Here, fi(Qo) >o andf2(Qo) = o. As a con-
sequence, the second-order variation of the adiabatic energy difference satisfies
coordinate along the gradient difference and Qx the rectilinear coordinate along the interstate-coupling vector, so that
SAV = AV(Qq + SQ) - AV(Qq) « S/1 + 2
(SS2)2
where AV0 = AV(Q0). The Hamiltonian operator H(Q) varies with Q in Eq. (4) but the states |S0; Q0) and |S1; Q0) do not. Functions /1(Q) and /2(Q) are, thus, diabatic quantities playing the role of parameters for the adiabatic energy difference. They depend implicitly on the reference geometry Q0, where the diabatic and adiabatic representations coincide. They can be used to define local curvilinear coordinates adapted to A V. In a quadratic approximation in terms of rectilinear coordinates, second-order contributions from S/1 and squared firstorder contributions from S/2 will define parabolic coordinates.
Since the energy gap is finite, the second term in Eq. (5) characterizes a second-order Jahn-Teller effect, also called pseudo-Jahn-Teller effect (see, e.g., Refs. 30 and 31). Note that a ot/ 4 rotation of the two diabatic states would swap the roles of /1 and /2, changing a second-order Jahn-Teller effect into an avoided crossing between two new diabatic states of the same energy (/2(Q0) = 0) and with a nonzero interstate coupling (/ 1(Q0) > 0). Such rotated states coincide with the adiabatic states at infinity rather than at the origin (Nikitin transformation, as discussed in Ref. 32). In benzene, they are related to the pair of resonant valence-bond Kekule structures.
First-order terms are now discussed. Since |S0; Q0) and |S1 ; Q0) are assumed to form a complete and "exact" basis set, the Hellman-Feynman theorem is valid at Q = Q0. As a consequence, the local gradient of AV (adiabatic representation) is equal to that of /1 (two-state crude-adiabatic, i.e., diabatic, representation),
[djAVl = <S1; Q0II4HI0IS1; Q0) - (S0; Q0|[jH]0|S0; Q0) = [j1]0, (6)
Sf 1 = [dX1(Hu- H00)]0SQX1 +
Sf2 = KH01LSQX, + - •
These coordinates, as well as second-order terms, are discussed below for benzene where D6h symmetry is explicitly taken into account.
In the following, we compare the coordinates to the 30 normal modes of benzene calculated at the S0 equilibrium geometry (D6h) Q = Q0. They are listed in Table I. There are 20 distinct frequencies corresponding to the nondegenerate modes plus ten pairs of twofold degenerate modes. A displacement parallel to the totally symmetric S1 -S0 gradient difference at Q = Q0 (coordinate QXl) involves only the a1g modes (mainly) 1 and 2 (see gradient difference similar to mode 1 in Fig. 3). S0 is A1g and S1 is B2u so a displacement along the interstate-coupling vector (coordinate Qx^ combines b2u modes 14 and (mainly) 15 (see interstate-coupling vector of approximately mode 15 in Fig. 3).
The geometrical basis set is now completed by introducing a set of 28 rectilinear coordinates Qj (j =1,28) orthogonal to Qx and Qx . These generalize the concept of first-order intersection space.33-36 Q1 is specifically chosen as the a1g mode such that Q1 and Qx span the same plane as modes 1 and 2. Similarly, Q2 is chosen as the b2u mode such that Q2 and Qx span the same plane as modes 14 and 15. For j = 3-28, the coordinates Qj are not specified yet but transform as the specific irreducible representations of the remaining normal modes.
The quadratic approximation, in other words the local harmonic approximation, of AV reads, thus,
where [3j]0 stands for the local partial derivative 3/3Qj|q=qq. Note that in the actual implementation, the corrected state-averaged gradients are used. In addition, the derivative-coupling vector in the adiabatic representation is parallel to the interstate-coupling vector in the diabatic representation (gradient of the interstate coupling /2),
[<Sq; Q\âj\S1-, Q>]
<Sp; Qo\[dyH]c|S1; QQ) AV(Qq)
AVq •
We now introduce a pair of mass-weighted rectilinear coordinates QQx and QQx that correspond to the generalization of the concept of a first-order branching space, i e , branching plane,33-36 to a nondegenerate case The pseudobranching-plane directions coincide in both adiabatic and diabatic pictures, as shown in Eqs . (6) and (7) . We call Qx the rectilinear
AV(Q) « AVq + [âx1/1]oqx1 + -[^f^qx
{.¿x^f 1]oQx1Q 1 + 2 2 [djdkf1]oQjQk
-j,k=1
[¿x^if 1]oQx2Q 2 +1 [¿f 1]oQ
[¿x^f2]Q
+ 2—2-Qx ,
AV0 x2
where [3j]0 stands now for the local partial derivative =Q0.
with Eq. (5):
3/3<2j|q=qo. This can be summarized as follows, together
f i(Q) ~ АУо + AQ + 2BSArx1x1ß21
+ BS/ISArxiiQx1Qi + 2 2 *b7jkQjQk
■i,k=i
+ BS/ISA у 2Q Q2 + -BSAy Q2 ,
f2(Q)" K"S°SiQx2, where (see also Refs. 8-i0) A^= [Vi]o,
Ayxixi = Й/i]0, АУх2Х2= Й/ i]0,
BS/ISAyxii = IA^lAL.
А Ух.
: lA^i]0,
ISAyjk = [jf]o (jk e IS).
Superscripts BS and IS refer to branching plane X and x2) and intersection space, respectively. The symbol Д was preferred here for quantities related to the energy difference Д V, while S refers to variations induced by a small geometrical displacement SQ. Note that the quadratic expansions of f1 and Д V differ only by a supplementary term due to f2, which alters the curvature along Qx in the Hessian of Д V. As mentioned above, this term is responsible for the second-order Jahn-Teller effect in benzene along Qx , which involves mostly the Kekule mode 15. It is always positive and leads to
the exaltation of the Sj Kekule frequency. '
In Eqs. (9)—(11), the cross terms are nonzero only between modes of the same irreducible representation (Г <E> Г DA1g). The (3N-8) X (3N-8) matrix (вДуд) is the mass-weighted Hessian of f1 projected out of the branching
plane in the spirit of the reaction-path Hamiltonian method. We now define explicitly the intersection-space coordinates Qj (j =1 to 28) as mass-weighted displacements along the eigenvectors of the projected Hessian38 (КДу,к) with eigen-
values \j. The eigenvectors are the normal modes and the eigenvalues are the normal curvatures (force constants) of the energy difference ДV within the intersection space (3N -8 dimensions). Note that coordinates Q1 and Q2, which had already been defined above, correspond to eigenvectors of (lSДyjk) for symmetry reasons. However, they involve couplings with the branching-space coordinates in the full (non-projected) Hessian of f1. Neglecting these couplings leads to a separate form for the adiabatic energy difference describing a paraboloid with a slope along Qx ,
FIG. 6. (Color online) Eigenvalues of the energy-difference Hessian computed at the Franck-Condon point of benzene in the 28-dimensional space orthogonal to the pseudobranching plane. The labels refer to the most similar normal modes of S0 benzene (see Table I). The dominant local motions are indicated in boxes.
AV(Q) « AV0 + AkQxi + 2BSAyxixiQ2i
4(^S0Si)2
Ayxx +
The magnitude and sign of the eigenvalues, ISXj, are used to classify the 28 eigenvectors of (КДуд) in terms of photochemical activity (see discussion in Sec. 2). Large negative eigenvalues correspond to photoactive modes. Evolution of the system along these modes from the Franck-Condon point is expected to lead more easily to the seam of conical intersection [see Fig. 4(a)]. In a quantum dynamics picture, the wavepacket starting around the Franck-Condon region will tend to spread along such directions. In contrast, large positive eigenvalues correspond to photoinactive modes. The wavepacket will be energized along such directions and will tend to contract. [see Fig. 4(b)]. Bath modes, with a near-zero eigenvalue [see Fig. 4(c)], will not play any significant role in the dynamics. The wavepacket will stay similar to the ground vibrational state along such directions, which can be neglected in a first approach.
IV. RESULTS AND DISCUSSION
The numerical result of this analysis is illustrated in Fig. 6. Calculations were performed with a complete active space self-consistent field (CASSCF) of six electrons spread over six tt molecular orbitals at the 6-31G* level. The new coordinates Qj, obtained from the diagonalization of the energy-difference Hessian (excited state minus ground state), were projected onto the original coordinates Qj, i.e., the normal modes of S0 benzene. Both sets are actually quite similar. This confirms that Duschinsky rotations are not large for them.39-41 On Fig. 6, the Qj coordinates are labeled with the corresponding main components in terms of Qj coordinates (listed in Table I). Note that Qx , Qi, Qx , and Q2 correspond mostly to Qi, Q2, Q15, and Qi4, respectively. A common
TABLE II. Columns 2-5: S0 frequencies of the 28 nontotally symmetric modes of benzene (specified in column 1 with the Wilson scheme of frequency numbering), for four D6h geometries defined by CH = 1.075 A and values of CC given in the second line. Columns 6-9: difference between Si and So frequencies under the same conditions.
"(So) (cm"1)' "(S,) -v(SQ) (cm-1)'
CC (Â) 1.350 1.396 1.434 1.500 1.350 1.396 1.434 1.500
"3 1a2g 1529 1500 1474 1425 6 6 5 5
"4 1b2g 690 724 737 738 -400 -296 -249 -195
"5 2b2g 1078 1035 995 913 -308 -289 -276 -259
"6 1-2e2g 659 654 643 612 -88 -72 -60 -41
"7 7-8e2g 3351 3360 3367 3376 -4 -5 -5 -5
"8 5-6e2g 1922 1739 1616 1468 61 60 56 42
"9 3-4e2g 1331 1275 1218 1086 16 22 32 55
"10 1-2e1g 901 874 849 799 -263 -249 -241 -233
"11 1a2u 726 711 698 671 -147 -149 -151 -156
"12 1b1u 1109 1099 1086 1054 -17 -18 -18 -18
"13 2b1u 3339 3350 3358 3369 -3 -3 -4 -5
"14 2b2u 1564 1339 1284 1239 658 675 577 391
"15 1b2u 1298 1185 962 447 28 110 309 785
"16 1-2e2u 339 431 475 516 -595 -255 -178 -126
"17 3-4e2u 1028 991 957 886 -306 -288 -276 -260
"18 1-2e1u 1195 1114 1043 917 -47 -64 -77 -96
"19 3-4e1u 1716 1630 1573 1493 -38 -27 -20 -11
"20 5-6e1u 3373 3379 3383 3389 -6 -6 -6 -6
aThe frequencies were calculated at the CASSCF(6,6)/6-31G* level and no scaling was applied. "Negative frequencies" are actually imaginary and correspond to negative curvatures (force constants).
feature of the Qj coordinates is that they tend to decouple the H motions (s CH stretching, ( HCC bending, and y CCCH wagging) from the C6-ring motions.
As a result, twelve modes are obviously of no interest (bath modes): the six s CH stretching modes and the six ( HCC bending modes. As a first approximation, they involve independent motion of the six H nuclei with respect to the C6 ring. In contrast, the photoactive modes (large negative value of IS\j) describe deformations of the C6 ring: three out-of-plane 8 CCCC ring-puckering modes (torsions), six out-of-plane y CCCH wagging modes, and nine in-plane modes mixing t CC stretching and a CCC bending. There is only one degenerate pair of photoinactive modes, similar to the pair 8 (e2g t CC stretching).
The predictions from the quadratic analysis of the energy difference were confirmed by two different tests. First, the evolution of the nontotally symmetric-mode frequencies on both states was analyzed along a totally symmetric deformation (same level of calculation as in the energy-difference analysis). Second, quantum dynamics simulations using the DD-vMCG approach14-16 were run within a reduced subspace of active coordinates. This method uses an expansion of the wavepacket on a time-dependent basis set of Gaussian functions. A local harmonic approximation of the PESs is calculated on the fly along the trajectory followed by the center of each Gaussian function. A diabatic picture is used to represent the pair of coupled electronic states. The dynamics code is implemented in a development version of the
Heidelberg MCTDH package and is currently interfaced with a development version of the GAUSSIAN program.43 The same theoretical level as in the static analysis was used.
We compare here the normal frequencies on S0 and S,.
At a given geometry, a negative/positive frequency difference indicates that motion along the mode decreases/ increases the energy gap ("negative frequencies" are to be understood as imaginary frequencies corresponding to negative curvatures). This should, thus, correspond to a negative/ positive value of IS\j, i.e., a photoactive/photoinactive mode when the Duschinsky rotation is small. To analyze the coupling between the gradient-difference mode (coordinate Qxi) and the nontotally symmetric modes, the S0 and S1 frequencies were calculated at four D6h geometries around the S0 and S1 equilibrium geometries: CC= 1.350, 1.396 (S0 minimum), 1.434 (S1 minimum) or 1.500, and CH= 1.075 A (S0 minimum). The results are listed in Table II. A single value was chosen for CH in order to compare frequency calculations in both states along the same cut line, namely, at the S0 minimum, which differs only slightly from that of the S1 minimum (1.073 A). Such frequency calculations away from the stationary point are meaningful only for the 28 nontotally symmetric modes. Modes presenting a strong dependence of their frequencies along the a1g deformation are expected to play a more significant role as they couple with the gradient difference to reduce the energy difference. This is a semiquantitative way of getting third-order effects.
Three out-of-plane modes dominate: the b2g mode 4 and the e2u pair 16 (S CCCC motions). The six other modes of this type also show a frequency lowering but it is less dependent on the value of CC. Some of the remaining in-plane skeletal deformations of the C6 ring also show a frequency weakening. The four most significant are the e2g pair 6 and the e1u pair 18. In addition, this confirmed that the e2g pair 8 is a photoinactive mode, giving here the largest positive fre-
TABLE III. Frequencies of the two b2u modes of benzene in the S0 and S electronic states and the differences v(Si) — v(Sq) for six D6h geometries defined by CH = 1.075 A and values of CC given in the first column.
^(So)(cm-1)a v(Si) (cm—1)a // — v(So) (cm-1)
CC (A)
Diabatic lebel V14 j15 j14 j15 /14 j15
1.350 1298/) 1564/) 1326/y 2222/) 28 658
1.396 1339/) 1185/) 1295/y 2014/) —44 829
1.434 1284/) 962/) 1271/) 1861/) —13 900
1.500 1239/) 447/) 1232/y 1631/) —8 1184
1.550 1212/) —460/) 1203/y 1479/) —10 1939
1.600 1188/) —722/) 1171 (/¡5) 1348/) —17 2069
aThe frequencies were calculated at the CASSCF(6,6)/6-31G* level and no scaling was applied. "Negative
frequencies" are actually imaginary and correspond to negative curvatures (force constants).
quency difference. This frequency analysis is consistent with the quadratic analysis of the energy difference. Among the 14 photoactive modes previously identified, seven modes (see Fig. 3) are shown here to be of greatest importance: 4, 6 (pair), 16 (pair), and 18 (pair). They were, thus, included in quantum dynamics calculations.
The only modes for which a large Duschinsky rotation has been shown are the b2u modes 14 and 15.39-41 As a consequence, they are not easy to treat separately. They compose the interstate-coupling vector at the Franck-Condon point. Note that they still dominate at the minimum of the conical intersection seam. It makes things clearer to correlate frequencies and geometries in a vibrational diabatic-like way. This is illustrated in Table III, where the frequencies are correlated with respect to the nature of the vibrational modes rather than the order of their values. Data for CC= 1.550 and 1.600 A have been added here. In this picture, the frequency that varies most with CC (third, fifth, and seventh columns in Table III) is associated with frequency v'5 in S0 but v'14 in S1 when CC ^ 1.396 A, and with v''4 in S0 and vj4 in S1 when CC ^ 1.350 A. This exalted frequency corresponds to the Kekule mode (t CC stretching + a CCC bending), while the other (second, fourth, and sixth columns) describes a ( HCC bending. The Kekule mode correlates with v'5 (lower frequency of the two b2u modes in the S0 state) at the reference geometry and, thus, the nomenclature mode 15 (v^) is retained for this mode even though it is labeled as v[4 in the S1 state. This emphasizes the corresponding second-order JahnTeller effect30,31 by showing how the subsequent exaltation effect of the S1 Kekule frequency is enhanced when the cycle expands. For rather stretched geometries, the twin Kekule structures appear even more stable (v15 < 0) leading, thus, to a valley-ridge inflection point. However, the second-order Jahn-Teller effect is known to be artificially overestimated at the CASSCF level.44,45
Finally, the relevance of the subset of active coordinates was tested with quantum dynamics simulations. To this end, DD-vMCG quantum dynamics simulations were used with a nine-dimensional (9D) model including the nine most important modes displayed on Fig. 3, namely, 1, 4, 6 (pair), 15, 16 (pair), and 18 (pair).
Simulations were started with a Franck-Condon Gaussian wavepacket placed on S1 at t=0 and approximated by a harmonic product of one-dimensional Gaussian functions with parameters based on a normal frequency analysis at Q0.
We focused on enhancing nonadiabatic transitions by stimulating photoactive modes. In order to control specifically the amount of excess energy deposited in the system, an additional mean momentum p=hk was given to the initial wave-packet, where the k vector has components with higher or lower magnitude along the normal coordinates Q. This corresponds to a shifted momentum distribution and is achieved by multiplying the real-valued multidimensional Gaussian function by the spatial phase factor exp(;k ■ Q).
Here, the effect of our choice of photoactive coordinates is illustrated in Fig. 7 by one case made up of two Gaussian functions in the wavepacket, one on each electronic state. The initial momentum is shown by an arrow at the Franck-Condon point. Components along modes 4, 6x, 16x, and 18x were chosen to target the S1 / S0 conical intersection minimum projected onto the 9D subspace (initial momentum: k4=6.7, k6x=-4.4, k16x =15.1, and k18x=2.2; components of k not explicitly mentioned are initially set to zero). Note that y components do not need to be considered for symmetry reasons. The evolution of the wavepacket is represented by the trajectories followed by the center of the two Gaussian functions over 27 fs. These are shown in three different planes of projection. Panel (a) of Fig. 7 corresponds to the geometrical plane spanned by coordinates Q4 and Q16x. Panel (b) corresponds to the (Q1, Q16x) plane and panel (c) to the (Q1, Q4) plane. The three panels together give a three-dimensional picture (Q1, Q4, Q16x) of the trajectories. The projections of the main stationary points are also plotted to exhibit where the trajectories go. The wavepacket starts off from the Franck-Condon point on S1 and splits into two components when it crosses the seam of intersection (dashed line in Fig. 7). The first component (red line in Fig. 7, see online version) represents the diabatic crossing from S1 to S0 (nonra-diative decay). The second component (blue line in Fig. 7, see online version) represents the part of the wavepacket that ends up trapped in S1. The trajectory of the latter is represented after 6.5 fs, once the corresponding Gaussian function is on S1 with a population larger than 1%.
In this example, mode 1 (gradient difference) is not initially activated (no initial momentum component). However, it gets stimulated because the system tends to follow the driving force (opposite of the gradient difference) from the Franck-Condon point toward the S1 minimum. This is shown by the trajectory going initially to the left (Q1 < 0, i.e., out-
FIG. 7. (Color online) Quantum trajectory of photoexcited benzene decaying from Si to S0 through a seam of conical intersection (dashed line) represented in three planes of projection: (Q4, Qi6x) subspace (a), (Qi, Qi&) subspace (b), and (Qi, Q4) subspace (c). The projections of the Si minimum (min*), the So minimum or Franck-Condon starting point (FC), the Si transition state (TS*), and the So prefulvene are represented with open circles.
ward breathing) on panels (b) and (c) of Fig. 7. However, it is not trapped in an oscillatory motion in the well of the Si minimum [see trajectory (a) in Fig. 5] because the photoactive modes (for example, 4 and i6x, see Fig. 7) were initially stimulated [see trajectory (b) in Fig. 5]. The trajectory starts tangent to the initial momentum (thick arrow in Fig. 7) and evolves toward positive values of Q4 and Qi6x. This concerted chair and boat deformation leads to half-chair geometries, similar to the prefulvenoid geometry of the conical intersection minimum. This can be seen on panel (a) of Fig. 7. A nonadiabatic transition is observed early on (~6-7 fs). The splitting of the wavepacket into two components characterizes the crossing of the conical intersection seam. This nonradiative decay occurs at a similar geometry to the conical intersection minimum but "dilated": about the same values of Q4 and Qi6x [see panel (a) of Fig. 7] and a negative value of Qi [see panels (b) and (c) of Fig. 7]. The component of the wavepacket transferred to S0 tends to go toward the prefulvene intermediate. The other is trapped in Si. The total population transferred to S0 is i2% and stays constant after
i0 fs. The same case was also tried with nine Gaussian functions on each state. The mean trajectory is similar but the population transfer is amplified. This example shows that stimulating photoactive modes is necessary and sufficient to reach efficiently the seam of conical intersection and, thus, induce nonradiative decay.
V. CONCLUSIONS
The optimal control of photochemical reactions is based on shaped laser pulses designed to generate photoproducts selectively. They are based on closed-loop techniques to achieve maximum efficiency through a "black-box" mechanism.46-50 Optimal control experiments are not easy to implement. One often needs to construct special experimental apparatus, and it is often not clear how to achieve the target in practice because there is no direct mapping between the control parameters and the molecular properties. We have, thus, reached a situation where it is desirable to perform theoretical computations before experiments. Simulations based on optimal control theory5i,52 monitor the wave-packet evolution induced by a time-dependent Hamiltonian. However, the effect of the resultant, often complicated pulse shapes is not easy to decipher. In contrast, we are interested here in an "open-loop" approach, where theoretical rationalization can precondition the laser pulse. This strategy is best illustrated by the first optimal control experiments ever performed on a cis-trans isomerization on cyanine dyes.53 Recently, we published results on a model cyanine dye54 that suggested that the behavior of an extended conical intersection seam could be used to control the product ratio by stimulating the skeletal deformations orthogonal to the seam. The experimentalists analyzed the frequency composition of the laser pulse ultimately used and showed conclusively that it was exactly such high-frequency modes that were active.
In the context of theory-assisted optimal control, it is thus essential to establish systematic methods to select such active coordinates. In this work, we proposed a new approach based on the local properties of the energy difference rather than the sole energy of the excited state. A second-order expansion of the energy difference as a function of nuclear coordinates around the Franck-Condon geometry is used to distinguish vibrational modes in terms of their importance for the photoreactivity. This analysis was applied to investigate the channel 3 process of benzene.
The nine most significant modes were automatically identified at the Si Franck-Condon point of benzene. Two of them are the gradient of the energy difference, i.e., gradient difference, which is a totally symmetric breathing mode, and the gradient of the Si-S0 coupling, i.e., interstate-coupling vector, which is a Kekule mode. These two directions define the pseudobranching plane. The remaining modes, called photoactive modes, lower the energy difference at the second order. They consist mainly in out-of-plane C6-ring torsions and in-plane skeletal deformations that lead to prefulvenoid geometries.
Quantum dynamics simulations were performed within the subspace defined by this set of active coordinates. If no mode is excited in the initial Franck-Condon wavepacket,
benzene is trapped in an ineffective breathing oscillation around the Si minimum. In contrast, stimulating photoactive modes proved to be an efficient way of driving the system to the seam of conical intersection, thus inducing nonradiative decay from Si to S0. Our analysis can be, therefore, used in the context of optimal control of photochemical reactivity for predicting which vibrations have to be stimulated by a laser pulse optimized to enhance radiationless decay.
ACKNOWLEDGMENTS
Special thanks are due to Dr. I. Burghardt for a fruitful collaboration in theoretical and algorithmic aspects of the DD-vMCG work. B.L. gratefully acknowledges financial support by EPSRC (Grant No. GR/T203ll/0l). F.S. is pleased to acknowledge the support of Gaussian, Inc. L.B. acknowledges support from Project No. CTQ2005-04563 of the Spanish/Ministerio de Educación y Ciencia/(MEC).
1M. A. Robb, M. Garavelli, M. Olivucci, and F. Bernardi, Rev. Comput. Chem. 15, 87 (2000). 2 Conical Intersections: Electronic Structure, Dynamics & Spectroscopy, edited by W. Domcke, D. R. Yarkony, and H. Koppel (World Scientific, Singapore, 2004).
3D. Lauvergnat, E. Baloitcha, G. Dive, and M. Desouter-Lecomte, Chem. Phys. 326, 500 (2006). 4F. Gatti, Y. Justum, M. Menou, A. Nauts, and X. Chapuisat, J. Mol. Spectrosc. 181, 403 (1997). 5F. Gatti, Chem. Phys. Lett. 373, 146 (2003). 6D. Lauvergnat and A. Nauts, J. Chem. Phys. 116, 8560 (2002). 7 E. Gindensperger, I. Burghardt, and L. S. Cederbaum, J. Chem. Phys.
124, 144103 (2006). 8M. J. Paterson, M. J. Bearpark, M. A. Robb, and L. Blancafort, J. Chem. Phys. 121, ii562 (2004). 9M. J. Paterson, M. J. Bearpark, M. A. Robb, L. Blancafort, and G. A. Worth, Phys. Chem. Chem. Phys. 7, 2i00 (2005). i0 F. Sicilia, L. Blancafort, M. J. Bearpark, and M. A. Robb, J. Phys. Chem.
A 111, 2182 (2007). ilI. J. Palmer, I. N. Ragazos, F. Bernardi, M. Olivucci, and M. A. Robb, J.
Am. Chem. Soc. 115, 673 (1993). 12A. L. Sobolewski, C. Woywod, and W. Domcke, J. Chem. Phys. 98, 5627 (i993).
13 G. A. Worth and I. Burghardt, Chem. Phys. Lett. 368, 502 (2003).
14 G. A. Worth, M. A. Robb, and I. Burghardt, Faraday Discuss. 127, 307 (2004).
15B. Lasorne, M. J. Bearpark, M. A. Robb, and G. A. Worth, Chem. Phys. Lett. 432, 604 (2006).
16 B. Lasorne, M. A. Robb, and G. A. Worth, Phys. Chem. Chem. Phys. 9, 32i0 (2007).
17 G. A. Worth and M. A. Robb, Adv. Chem. Phys. 124, 355 (2002). 18E. B. Wilson, Phys. Rev. 45, 706 (1934).
i9F. Sicilia, M. J. Bearpark, L. Blancafort, and M. A. Robb, Theor. Chem. Acc. 118, 24i (2007).
20 C. A. Mead, J. Chem. Phys. 78, 807 (1983).
21T. A. Barckholtz and T. A. Miller, Int. Rev. Phys. Chem. 17, 435 (1998).
22 T. A. Barckholtz and T. A. Miller, J. Phys. Chem. A 103, 2321 (1999).
23 D. R. Yarkony, J. Chem. Phys. 123, 134106 (2005).
24 D. R. Yarkony, J. Chem. Phys. 123, 204101 (2005).
25M. S. Schuurman and D. R. Yarkony, J. Chem. Phys. 127, 094104 (2007).
26H. Muller, H. Köppel, and L. S. Cederbaum, New J. Chem. 17, 7 (1993).
27H. Köppel and B. Schubert, Mol. Phys. 104, 1069 (2006).
281. B. Bersuker, The Jahn-Teller effect (Cambridge University Press, Cambridge, 2006).
29 H. C. Longuet-Higgins, Advances in Spectroscopy, edited by W. B. Thompson (Interscience, New York, 1961), Vol. 2, p. 429.
30 S. Shaik, A. Shurki, D. Danovich, and P. C. Hiberty, J. Am. Chem. Soc. 118, 666 (1996).
31L. Blancafort and M. Sola, J. Phys. Chem. A 110, 11219 (2006).
32 M. Desouter-Lecomte, C. Galloy, J. C. Lorquet, and M. V. Pires, J. Chem. Phys. 71, 3661 (1979).
33E. R. Davidson, J. Am. Chem. Soc. 99, 397 (1977).
34 M. Desouter-Lecomte, D. Dehareng, B. Leyhnihant, M. T. Praet, A. J. Lorquet, and J. C. Lorquet, J. Phys. Chem. 89, 214 (1985).
35 G. J. Atchity, S. S. Xantheas, and K. Ruedenberg, J. Chem. Phys. 95, 1862 (1991).
36 D. R. Yarkony, Rev. Mod. Phys. 68, 985 (1996).
37 W. H. Miller, N. C. Handy, and J. E. Adams, J. Chem. Phys. 72, 99 (1980).
38A. G. Baboul and H. B. Schlegel, J. Chem. Phys. 107, 9413 (1997).
39 F. Metz, M. J. Robey, E. W. Schlag, and F. Dorr, Chem. Phys. Lett. 51, 8 (1977).
40A. Bernhardsson, N. Forsberg, P. A. Malmqvist, B. O. Roos, and L. Serrano-Andres, J. Chem. Phys. 112, 2798 (2000).
41G. Orlandi, P. Palmieri, R. Tarroni, F. Zerbetto, and M. Z. Zgierski, J. Chem. Phys. 100, 2458 (1994).
42 G. A. Worth, M. H. Beck, A. Jackle, and H. D. Meyer, The MCTDH package, development version 9.0, University of Heidelberg, Heidelberg, Germany, 2006 (see http://www.pci.uni-heidelberg.de/tc/usr/mctdh/).
43 M. J. Frisch, G. W. Trucks H. B. Schlegel et al. GAUSSIAN development version, revision F.02, Gaussian, Inc., Wallingford, CT, 2003.
^R. Berger, C. Fischer, and M. Klessinger, J. Phys. Chem. A 102, 7157 (1998).
45 S. Schumm, M. Gerhards, and K. Kleinermanns, J. Phys. Chem. A 104, 10648 (2000).
46 S. A. Rice and M. Zhao, Optical Control of Molecular Dynamics (Wiley, New York, 2000).
47 M. Shapiro and P. Brumer, Principles of the Quantum Control of Molecular Processes (Wiley, New York, 2003).
48 V. May and O. Kühn, Charge and Energy Transfer Dynamics in Molecular Systems, 2nd ed. (Wiley-VCH, Berlin, 2004).
49 J. L. Herek, J. Photochem. Photobiol., A 180, 225 (2006).
50 Coherent Control of Molecules, edited by B. Lasorne and G. A. Worth (CCP6, Daresbury, 2006).
51D. J. Tannor and S. A. Rice, J. Chem. Phys. 83, 5013 (1985).
52R. S. Judson and H. Rabitz, Phys. Rev. Lett. 68, 1500 (1992).
53 B. Dietzek, B. Bruggemann, T. Pascher, and A. Yartsev, Phys. Rev. Lett. 97, 258301 (2006).
54 P. A. Hunt and M. A. Robb, J. Am. Chem. Soc. 127, 5720 (2005).
