Automatic generation of active coordinates for quantum dynamics calculations : Application to the dynamics of benzene photochemistry

calculations: Application to the dynamics of benzene photochemistry Benjamin Lasorne, Fabrizio Sicilia, Michael J. Bearpark, Michael A. Robb, Graham A. Worth, and Lluìs Blancafort Department of Chemistry, Imperial College London, South Kensington, London SW7 2AZ, United Kingdom School of Chemistry, University of Birmingham, Edgbaston, Birmingham B15 2TT, United Kingdom Institut de Química Computacional and Departament de Química, Universitat de Girona, E-17071 Girona, Spain


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. 2 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 difference [8][9][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 S 0 minimum.
In any nonadiabatic photochemical reaction, there are two special coordinates: 1,2 x 1 , along the gradient difference, and x 2 , along the interstate-coupling vector.The S 1 and S 0 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 ͑x 1 , x 2 ͒, one has the information about the curvature of the energy difference in this subspace.Diagonalization of this projected Hessian gives a set of vibrational 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 S 1 ← S 0 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 S 1 / S 0 intersection seam is accessed.For dynamics, we have used the direct dynamics variational multiconfiguration Gaussian wavepacket ͑DD-vMCG͒ method [13][14][15][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 17 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.
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 ͑GS Min 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 intersection 1 ͑CoIn in Fig. 2͒ and decays to form the products in the ground state ͑GS Prod 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 S 1 ← S 0 benzene photochemistry, the photochemical event is illustrated in Fig. 1, where the S 0 and S 1 surfaces are plotted in the space spanned by the coordinates along the gradient difference ͑x 1 ͒ and interstate-coupling vector ͑x 2 ͒.Such directions are traditionally computed at the conical intersection point, where they define the branching plane. 1,2Although the energy gap is not zero, the same vectors can be defined at the Franck-Condon point.Here, x 1 and x 2 span a "pseudobranching plane," which may be different from the branching plane at the conical intersection point.In the following, we describe x 1 and x 2 at both points by comparing them to a reference set of vibrational modes computed at the S 0 equilibrium geometry of benzene ͑see Table I, where the Wilson scheme of frequency numbering 18 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.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 S 1 gradient.In Fig. 4, the minima on both the ground and excited state PESs are at the same value of the coordinates, since the S 0 gradient is zero and the S 1 gradient is orthogonal to these directions at the Franck-Condon point.
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 FIG. 3. ͑Color online͒ The nine dominant motions for the photochemistry of benzene ͑the labels refer to the most similar normal modes of S 0 benzene listed in Table I͒.at the Franck-Condon point ͑Q 1 ͒ and a generic photoactive mode ͑Q pa ͒.At the Franck-Condon point ͑D 6h benzene͒, the gradient difference is identical to the S 1 gradient.The driving force ͑opposite of the gradient͒ leads to the S 1 minimum ͑stretched D 6h 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 Q pa , 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 S 0 .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 S 1 to S 0 .
In summary, efficient photochemical reaction paths can be predicted by analyzing the local properties of the energy difference between S 1 and S 0 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][9][10]19 to the general quadratic representation of conical intersections 2,20-28 is an analytic representation of the energy difference between the states. Th 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 D 6h 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 ͉S 0 ; Q͘ and ͉S 1 ; Q͘ ͑parametrized by the nuclear geometry Q͒ are taken as known at a reference geometry Q = Q 0 , chosen as the S 0 equilibrium geometry.In the representation formed by the pair of adiabatic states, the matrix of the clamped-nucleus Hamiltonian operator

͑1͒
However, it is no longer diagonal when using the same states ͉S 0 ; Q 0 ͘ and ͉S 1 ; The mass-weighted geometrical displacement ␦Q generates a finite variation of the off-diagonal elements ␦H 10 = ␦H 01 ͑real-valued electronic wavefunctions͒ because the adiabatic states frozen at Q = Q 0 define trivial diabatic states when Q varies in the Hamiltonian operator.Such states are often referred to as "crude adiabatic." 29n 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 ␦Q according to

͑3͒
where ⌬V͑Q͒ = V 1 ͑Q͒ − V 0 ͑Q͒.We now introduce two functions of the coordinates Q Note that a conical intersection would correspond to f 1 ͑Q 0 ͒ = 0 and f 2 ͑Q 0 ͒ = 0. Here, f 1 ͑Q 0 ͒ Ͼ 0 and f 2 ͑Q 0 ͒ = 0.As a con-sequence, the second-order variation of the adiabatic energy difference satisfies where ⌬V 0 = ⌬V͑Q 0 ͒.The Hamiltonian operator H ˆ͑Q͒ varies with Q in Eq. ͑4͒ but the states ͉S 0 ; Q 0 ͘ and ͉S 1 ; Q 0 ͘ do not.Functions f 1 ͑Q͒ and f 2 ͑Q͒ are, thus, diabatic quantities playing the role of parameters for the adiabatic energy difference.They depend implicitly on the reference geometry Q 0 , where the diabatic and adiabatic representations coincide.They can be used to define local curvilinear coordinates adapted to ⌬V.
In a quadratic approximation in terms of rectilinear coordinates, second-order contributions from ␦f 1 and squared firstorder contributions from ␦f 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 / 4 rotation of the two diabatic states would swap the roles of f 1 and f 2 , changing a second-order Jahn-Teller effect into an avoided crossing between two new diabatic states of the same energy ͑f 2 ͑Q 0 ͒ =0͒ and with a nonzero interstate coupling ͑f 1 ͑Q 0 ͒ Ͼ 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 Kékulé structures.
First-order terms are now discussed.Since ͉S 0 ; Q 0 ͘ and ͉S 1 ; Q 0 ͘ are assumed to form a complete and "exact" basis set, the Hellman-Feynman theorem is valid at Q = Q 0 .As a consequence, the local gradient of ⌬V ͑adiabatic representa-tion͒ is equal to that of f 1 ͑two-state crude-adiabatic, i.e., diabatic, representation͒, where ‫ץ͓‬ j ͔ 0 stands for the local partial derivative ‫ץ͉‬ / ‫ץ‬Q j ͉ Q=Q 0 .Note that in the actual implementation, the corrected stateaveraged gradients are used.In addition, the derivativecoupling vector in the adiabatic representation is parallel to the interstate-coupling vector in the diabatic representation ͑gradient of the interstate coupling f 2 ͒, We now introduce a pair of mass-weighted rectilinear coordinates Q ¯x1 and Q ¯x2 that correspond to the generalization of the concept of a first-order branching space, i.e., branching plane, [33][34][35][36] to a nondegenerate case.The pseudobranchingplane directions coincide in both adiabatic and diabatic pictures, as shown in Eqs.͑6͒ and ͑7͒.We call Q ¯x1 the rectilinear coordinate along the gradient difference and Q ¯x2 the rectilinear coordinate along the interstate-coupling vector, so that These coordinates, as well as second-order terms, are discussed below for benzene where D 6h symmetry is explicitly taken into account.
In the following, we compare the coordinates to the 30 normal modes of benzene calculated at the S 0 equilibrium geometry ͑D 6h ͒ Q = Q 0 .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 S 1 − S 0 gradient difference at Q = Q 0 ͑coordinate Q ¯x1 ͒ involves only the a 1g modes ͑mainly͒ 1 and 2 ͑see gradient difference similar to mode 1 in Fig. 3͒.S 0 is A 1g and S 1 is B 2u so a displacement along the interstate-coupling vector ͑coordinate Q ¯x2 ͒ combines b 2u 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 Q ¯j ͑j =1,28͒ orthogonal to Q ¯x1 and Q ¯x2 .4][35][36] Q ¯1 is specifically chosen as the a 1g mode such that Q ¯1 and Q ¯x1 span the same plane as modes 1 and 2. Similarly, Q ¯2 is chosen as the b 2u mode such that Q ¯2 and Q ¯x2 span the same plane as modes 14 and 15.For j = 3 -28, the coordinates Q ¯j 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 ⌬V reads, thus, where ‫ץ͓‬ j ͔ 0 stands now for the local partial derivative ‫ץ͉‬ / ‫ץ‬Q ¯j͉ Q ¯=Q ¯0.This can be summarized as follows, together with Eq. ͑5͒: 28 IS ⌬␥ jk Q ¯jQ ¯k where ͑see also Refs.8-10͒ Superscripts BS and IS refer to branching plane ͑x 1 and x 2 ͒ and intersection space, respectively.The symbol ⌬ was preferred here for quantities related to the energy difference ⌬V, while ␦ refers to variations induced by a small geometrical displacement ␦Q.Note that the quadratic expansions of f 1 and ⌬V differ only by a supplementary term due to f 2 , which alters the curvature along Q ¯x2 in the Hessian of ⌬V.As mentioned above, this term is responsible for the second-order Jahn-Teller effect in benzene along Q ¯x2 , which involves mostly the Kékulé mode 15.It is always positive and leads to the exaltation of the S 1 Kékulé frequency. 30,31n Eqs.͑9͒-͑11͒, the cross terms are nonzero only between modes of the same irreducible representation ͑⌫ ⌫ ʛ A 1g ͒.The ͑3N −8͒ ϫ ͑3N −8͒ matrix ͑ IS ⌬␥ jk ͒ is the mass-weighted Hessian of f 1 projected out of the branching plane in the spirit of the reaction-path Hamiltonian method. 37e now define explicitly the intersection-space coordinates Q ¯j ͑j = 1 to 28͒ as mass-weighted displacements along the eigenvectors of the projected Hessian 38 ͑ IS ⌬␥ jk ͒ with eigenvalues IS 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 Q ¯1 and Q ¯2, which had already been defined above, correspond to eigenvectors of ͑ IS ⌬␥ jk ͒ for symmetry reasons.However, they involve couplings with the branching-space coordinates in the full ͑non-projected͒ Hessian of f 1 .Neglecting these couplings leads to a separate form for the adiabatic energy difference describing a paraboloid with a slope along Q ¯x1 ,

͑12͒
The magnitude and sign of the eigenvalues, IS j , are used to classify the 28 eigenvectors of ͑ IS ⌬␥ jk ͒ 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 nearzero 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 molecular orbitals at the 6-31G * level.The new coordinates Q ¯j, obtained from the diagonalization of the energydifference Hessian ͑excited state minus ground state͒, were projected onto the original coordinates Q j , i.e., the normal modes of S 0 benzene.Both sets are actually quite similar.0][41] On Fig. 6, the Q ¯j coordinates are labeled with the corresponding main components in terms of feature of the Q ¯j coordinates is that they tend to decouple the H motions ͑s CH stretching, ␤ HCC bending, and ␥ CCCH wagging͒ from the C 6 -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 C 6 ring.In contrast, the photoactive modes ͑large negative value of IS j ͒ describe deformations of the C 6 ring: three out-ofplane ␦ CCCC ring-puckering modes ͑torsions͒, six out-ofplane ␥ CCCH wagging modes, and nine in-plane modes mixing t CC stretching and ␣ CCC bending.There is only one degenerate pair of photoinactive modes, similar to the pair 8 ͑e 2g 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 approach [14][15][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 42 and is currently interfaced with a development version of the GAUSSIAN program. 43The same theoretical level as in the static analysis was used.
We compare here the normal frequencies on S 0 and S 1 .
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 Q ¯x1 ͒ and the nontotally symmetric modes, the S 0 and S 1 frequencies were calculated at four D 6h geometries around the S 0 and S 1 equilibrium geometries: CC = 1.350, 1.396 ͑S 0 mini-mum͒, 1.434 ͑S 1 minimum͒ or 1.500, and CH = 1.075Å ͑S 0 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 S 0 minimum, which differs only slightly from that of the S 1 minimum ͑1.073Å͒.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 a 1g 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 b 2g mode 4 and the e 2u pair 16 ͑␦ 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 C 6 ring also show a frequency weakening.The four most significant are the e 2g pair 6 and the e 1u pair 18.In addition, this confirmed that the e 2g pair 8 is a photoinactive mode, giving here the largest positive fre-TABLE II.Columns 2-5: S 0 frequencies of the 28 nontotally symmetric modes of benzene ͑specified in column 1 with the Wilson scheme of frequency numbering͒, for four D 6h geometries defined by CH = 1.075Å and values of CC given in the second line.Columns 6-9: difference between S 1 and S 0 frequencies under the same conditions.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 b 2u modes 14 and 15. [39][40][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 Å 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 15 Љ in S 0 but 14 Ј in S 1 when CC ജ 1.396 Å, and with 14 Љ in S 0 and 14 Ј in S 1 when CC ഛ 1.350 Å.This exalted frequency corresponds to the Kékulé mode ͑t CC stretching+ ␣ CCC bending͒, while the other ͑second, fourth, and sixth columns͒ describes a ␤ HCC bending.The Kékulé mode correlates with 15 Љ ͑lower fre- quency of the two b 2u modes in the S 0 state͒ at the reference geometry and, thus, the nomenclature mode 15 ͑ 15 ͒ is retained for this mode even though it is labeled as 14 Ј in the S 1 state.This emphasizes the corresponding second-order Jahn-Teller effect 30,31 by showing how the subsequent exaltation effect of the S 1 Kékulé frequency is enhanced when the cycle expands.For rather stretched geometries, the twin Kékulé structures appear even more stable ͑ 15 Ͻ 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,45inally, 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 S 1 at t = 0 and approximated by a harmonic product of one-dimensional Gaussian functions with parameters based on a normal frequency analysis at Q 0 .
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 = បk was given to the initial wavepacket, 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͑ik • 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 S 1 / S 0 conical intersection minimum projected onto the 9D subspace ͑initial momentum: k 4 = 6.7, k 6x = −4.4,k 16x = 15.1, and k 18x = 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 Q 4 and Q 16x .Panel ͑b͒ corresponds to the ͑Q 1 , Q 16x ͒ plane and panel ͑c͒ to the ͑Q 1 , Q 4 ͒ plane.The three panels together give a three-dimensional picture ͑Q 1 , Q 4 , Q 16x ͒ 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 S 1 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 ver-sion͒ represents the diabatic crossing from S 1 to S 0 ͑nonradiative decay͒.The second component ͑blue line in Fig. 7, see online version͒ represents the part of the wavepacket that ends up trapped in S 1 .The trajectory of the latter is represented after 6.5 fs, once the corresponding Gaussian function is on S 1 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 S 1 minimum.This is shown by the trajectory going initially to the left ͑Q 1 Ͻ 0, i.e., out-  ward breathing͒ on panels ͑b͒ and ͑c͒ of Fig. 7.However, it is not trapped in an oscillatory motion in the well of the S 1 minimum ͓see trajectory ͑a͒ in Fig. 5͔ because the photoactive modes ͑for example, 4 and 16x, 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 Q 4 and Q 16x .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 Q 4 and Q 16x ͓see panel ͑a͒ of Fig. 7͔ and a negative value of Q 1 ͓see panels ͑b͒ and ͑c͒ of Fig. 7͔.The component of the wavepacket transferred to S 0 tends to go toward the prefulvene intermediate.The other is trapped in S 1 .The total population transferred to S 0 is 12% and stays constant after 10 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.7][48][49][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 theory 51,52 monitor the wavepacket 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. 53Recently, we published results on a model cyanine dye 54 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 secondorder 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 S 1 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 S 1 -S 0 coupling, i.e., interstate-coupling vector, which is a Kékulé 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 C 6 -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 S 1 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 S 1 to S 0 .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.

FIG. 5 .
FIG. 5. ͑Color online͒ Schematic representation of two trajectories starting from the benzene S 1 Frank-Condon point ͑FC͒.Trajectory ͑a͒ oscillates within the S 1 minimum valley, whereas trajectory ͑b͒ reaches the ground state product passing through a conical intersection point ͑CoIn͒.
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 S 0 benzene ͑see TableI͒.The dominant local motions are indicated in boxes.

TABLE I .
The 30 normal modes of S 0 benzene labeled following the Wilson scheme of frequency numbering.

TABLE III .
Frequencies of the two b 2u modes of benzene in the S 0 and S 1 electronic states and the differences ͑S 1 ͒ − ͑S 0 ͒ for six D 6h geometries defined by CH = 1.075Å and values of CC given in the first column.