0Quantum dynamics study of fulvene double bond photoisomerization: The role of intramolecular vibrational energy redistribution and excitation energy
Llms Blancafort, Fabien Gatti, and Hans-Dieter Meyer
Citation: The Journal of Chemical Physics 135, 134303 (2011); doi: 10.1063/1.3643767 View online: http://dx.doi.org/10.1063/1.3643767
View Table of Contents: http://scitation.aip.org/content/aip/journal/jcp/135/13?ver=pdfcov Published by the AIP Publishing
Articles you may be interested in
Azole energetic materials: Initial mechanisms for the energy release from electronical excited nitropyrazoles J. Chem. Phys. 140, 034320 (2014); 10.1063/1.4861670
Assessing computationally efficient isomerization dynamics: SCF density-functional theory study of azobenzene molecular switching
J. Chem. Phys. 135, 224303 (2011); 10.1063/1.3664305
Theoretical investigation of highly excited vibrational states in DFCO: Calculation of the out-of-plane bending states and simulation of the intramolecular vibrational energy redistribution J. Chem. Phys. 126, 024302 (2007); 10.1063/1.2402920
Vibrational spectroscopy and intramolecular dynamics of 1-butyne J. Chem. Phys. 121, 5860 (2004); 10.1063/1.1786923
On the photoisomerization of 5-hydroxytropolone: An abinitio and nuclear wave function study J. Chem. Phys. 107, 6275 (1997); 10.1063/1.474344
Quantum dynamics study of fulvene double bond photoisomerization: The role of intramolecular vibrational energy redistribution and excitation energy
Lluís Blancafort,1,a) Fabien Gatti,2 and Hans-Dieter Meyer3
1 Institut de Química Computacional, Department de Química, Universitat de Girona, Campus de Montilivi, 17071 Girona, Spain
2CTMM, Institut Charles Gerhardt Montpellier (UMR 5253), CC 1501, Université Montpellier 2, 34095 Montpellier Cedex 05, France
3Theoretische Chemie, Ruprecht-Karls-Universität, Im Neuenheimer Feld 229, 69120 Heidelberg, Germany (Received 22 July 2011; accepted 8 September 2011; published online 3 October 2011)
The double bond photoisomerization of fulvene has been studied with quantum dynamics calculations using the multi-configuration time-dependent Hartree method. Fulvene is a test case to develop optical control strategies based on the knowledge of the excited state decay mechanism. The decay takes place on a time scale of several hundred femtoseconds, and the potential energy surface is centered around a conical intersection seam between the ground and excited state. The competition between unreactive decay and photoisomerization depends on the region of the seam accessed during the decay. The dynamics are carried out on a four-dimensional model surface, parametrized from complete active space self-consistent field calculations, that captures the main features of the seam (energy and locus of the seam and associated branching space vectors). Wave packet propagations initiated by single laser pulses of 5-25 fs duration and 1.85-4 eV excitation energy show the principal characteristics of the first 150 fs of the photodynamics. Initially, the excitation energy is transferred to a bond stretching mode that leads the wave packet to the seam, inducing the regeneration of the re-actant. The photoisomerization starts after the vibrational energy has flowed from the bond stretching to the torsional mode. In our propagations, intramolecular energy redistribution (IVR) is accelerated for higher excess energies along the bond stretch mode. Thus, the competition between unreactive decay and isomerization depends on the rate of IVR between the bond stretch and torsion coordinates, which in turn depends on the excitation energy. These results set the ground for the development of future optical control strategies. © 2011 American Institute of Physics. [doi:10.1063/1.3643767]
I. INTRODUCTION
One of the main challenges for chemists today is to use laser pulses to control reactions.1 In some cases, mode selectivity has provided a way to guide these processes successfully, so that desired products are obtained preferentially. This field also offers new possibilities for the understanding of fundamental phenomena involving the conversion of light into mechanical motion such as the elementary steps of vision, photosynthesis, protein dynamics, and electron and proton transport in DNA.1 However, the control of chemical reactions by laser pulses cannot be done without knowing and controlling the time scales and the properties of the intramolecular energy redistribution in molecules. Much theoretical effort has thus been directed toward the investigation of the underlying quantum dynamics.
In former papers (see Chapter 21 in Ref. 2 for a review), some of us have demonstrated that the multi-configuration time dependent Hartree (MCTDH) algorithm3,4 is an efficient tool to investigate the energy redistribution in molecules in their electronic ground state after excitation by a laser pulse. We have, for instance, studied the cis-trans isomerization5-7 of HONO or the inversion8 of NHD2 in full dimensionality.
a)Electronic mail: lluis.blancafort@udg.edu.
In these works, the processes were guided through quantum mechanical effects, namely, tunneling, but the molecules were considered in their electronic ground state only (Born-Oppenheimer approximation). In contrast to this, photochemistry is often governed by strong non-Born-Oppenheimer couplings between the electronic structure and the nuclear dynamics, at regions of conical intersection. In general, an intersection enables fast decay to the ground state, and it acts as a bifurcation between different reaction paths that can lead to the initial reactant and a photoproduct, or to different photoproducts.9-13 It is therefore amenable to optical control, at least in principle. A further possibility of control is given by the extended nature of the conical intersection seam.14 Several potential energy surface studies on different molecules have shown that the seam of intersection can be composed of different segments, each one associated to reaction channels leading to different species.15-19 In this case, the outcome of the excitation process depends on the seam segment encountered in the decay to the lower state. Therefore, the present work can be seen as a first step to extend our studies of the intramolecular energy redistribution with the MCTDH approach to a reaction where the key step is the passage through a conical intersection seam. In such a challenging scenario, our mechanistic knowledge of the reaction shall provide the guiding principles for future control strategies.
0021-9606/2011/135(13)/134303/11/$30.00
135, 134303-1
© 2011 American Institute of Physics
The idea of optical control at a seam has been tested before in direct dynamics studies, both at the quantum and mixed quantum classical dynamics level. In these studies the wave packet propagations or mixed quantum classical trajectories are guided to different seam segments varying the initial position or momentum.15,16,18,20 While this scheme provides a proof of principle of this type of control, our aim here is to implement it in a quantum dynamics formalism that allows to include explicitly the laser pulse required for a realistic simulation. In particular, in this paper we derive a suitable four-dimensional model to simulate the photoisomer-ization of a polyatomic molecule, fulvene, and carry out propagations with different laser pulses. These propagations prove the suitability of the model surface for the dynamics, and provide a first approach to understand the factors on which the photoisomerization depends.
Our idea of control is illustrated in Figure 1 for the case of fulvene, an unsaturated, non-fluorescent hydrocarbon with a delocalized n system. The photophysics of fulvene is characterized by a broad, weak, and diffuse S1 absorption band, corresponding to a valence excited state of B2 symmetry.21-23 The shape of the band is thought to arise from large geometry changes that occur immediately after the excitation, and the absence of fluorescence suggests the presence of a conical intersection seam between the ground and lowest excited state that enables the fast decay of the excited species. These experimental features have been explained by ab initio calculations that show that there is an energetically accessible S1/S0 seam of conical intersection.10,24 The potential energy surface along the two relevant coordinates, a totally symmetric bond alternating mode and the CH2 torsional mode, is sketched in Figure 1. Along the torsional coordinate there are two symmetric minima corresponding to two chemically equivalent species (coded by the different color of the methylene hy-
FIG. 1. Potential energy surface sketch for the S1 and So states of fulvene along the torsional and symmetric bond alternation coordinates, displaying the paths for photostability (1) and photoisomerization (2).
drogen atoms in Figure 1). The two isomers are connected by a double bond isomerization. They are labeled FC, which stands for Franck-Condon geometry, and FC'. Every point on the seam plotted in Figure 1 is a conical intersection if another coordinate (the interstate coupling coordinate, see below) is taken into account. When a wave packet is promoted to the electronically excited state from one of the ground state minima, FC, there are two possibilities: (1) decay in the region of the seam close to FC, reached along the bond alternating coordinate, which will regenerate the reactant, or (2) double bond twisting before the decay, which will allow for the decay at the seam segment close to FC' and induce the double bond isomerization. Based on this picture, fulvene is a good system to test optical control strategies which aim to govern the reactivity by controlling the segment of the seam where decay of the wave packet takes place.
The photodynamics of fulvene has been the subject of several previous theoretical dynamics studies. The decay along the seam was studied in an early example with mixed quantum classical dynamics (classical propagation of the nuclei combined with a trajectory surface hopping algorithm) which showed that different trajectories decay to the ground state at different points on the seam.10 Recently there have been several quantum dynamics studies on two- and three-dimensional surfaces with the goal of separating nuclear spin isomers.25-27 In one of these studies, it was found that the occurrence of the isomerization depended on the shape of the excitation laser pulse.26 The control of the pho-toisomerization of fulvene has also been addressed recently in a direct full dimensional variational multi-configurational Gaussian (vMCG) dynamics study.20 In this case, the isomer-ization depended on the initial momentum and position of the wave packet. In the present paper we carry out a quantum dynamics study on a four-dimensional model surface, using the MCTDH approach.2-4,28 We use the idea of control at the conical intersection seam to guide our choice of coordinates, centering on the coordinates corresponding to seam segments of low energy. As we will show in Sec. II, a set of four coordinates provides a very good model for the potential energy surface along these segments. This happens at the expense of taking some approximations, but it is largely compensated by the good description of the surface and the relative simplicity of the potential. On this model we have carried out propagations with single laser pulses as a first attempt of a realistic simulation of the control. Initially, the vibrational energy is located on a bond stretching mode, and the photoisomeriza-tion depends on the intramolecular vibrational redistribution (IVR) to the torsional mode, which in turn depends on the energy of the excitation pulse.
II. ELECTRONIC STRUCTURE CALCULATIONS AND PARAMETRIZATION OF THE ELECTRONIC HAMILTONIAN
Our choice of coordinates aims to capture the lowest energy region of the seam and its most important features, including the branching space vectors, and the minimum energy path from the FC structure to the seam. The excited state relaxation and intersection space of fulvene have been studied
1.2/1.1 v(Si)m
-0.3 /-0.2
^^plan 0.0/ 0.0
.349 (FC)
1479 (FC) 1 1.498 (SiU
1 399 X1.578 (CIplan)
1.372 (О^Ж. W
FC -2.9/-2.1 (a)
.355 (FC) 1471 (S^n _ .531 (CIplan)
.483 (
1.360 (S^' 1.320 (CIplan)
FIG. 2. (a) Potential energy profile for the initial excited state relaxation coordinate along planar structures (energies in eV with respect to S1 energy at FC structure; ab initio data in blue and model data in red). (b) Bond lengths in A for optimized structures along the relaxation coordinate.
in detail previously with electronic structure calculations at the complete active space self-consistent field (CASSCF) level of theory.10,24 We base our choice on these results, which are summarized briefly. The initial part of the decay (see Figure 2) retains the planarity of the molecule and involves a bond inversion coordinate that leads to a minimum on the excited state surface, (S1)Min, and further to a conical intersection with the ground state, CIplan (see optimized ab initio energies in blue in Figure 2(a) and bond lengths in Figure 2(b)). CIplan belongs to an extended intersection space that contains two low energy segments along the CH2 torsion and pyramidalization coordinates. Overall, the seam contains four relevant conical intersection stationary points (see Table I for a summary of the energies). These structures are shown in Figure 3, together with the displacements corresponding to the branching space vectors. These vectors are the gradient difference (gd) and interstate coupling (ic) vectors that lift the degeneracy at the intersection at first order (see Eq. (1)):
gd = Vq (E2 - Ei), ic = mVQHele |Ф2).
(1a) (ib)
TABLE I. Energies of the optimized critical points relative to CIptan at the ab initio level (CASSCF and CASPT2) compared to the four-dimensional model Hamiltonian W.
Structure ecasscf Erel (eV)a E CASPT2 Erel (eV)b pmodel Erel (eV)
FC (So) - 2.88 - 2.77 - 2.11
FC (Sl)c 1.25 - 0.97 1.12
(S1)Min - 0.28 - 0.17 - 0.19
CIplan 0.00 0.00(0.3l)d 0.00
CIperp - 0.33 - 0.l6(0.06)d - 0.33
CI63 - 0.07 - 0.01 (0.28)d - 0.06
CIpyr - 0.42 - 0.34 (0.20)d - 0.40
aEnergies of CASSCF(6,6)/cc-pvdz optimized points using GAUSSIAN 03 (Ref. 29) (see Ref. 24 for details).
bSingle point CASPT2(6,6)/cc-pvdz energies on CASSCF optimized structures. Energies calculated with Molcas 7.2 (Ref. 30) state averaging over the two lowest states with equal weights, a real level shift (Ref. 31) of0.1 a.u. and an ionization potential electron affinity shift (Ref. 32) of 0.25 a.u. cExperimental vertical excitation: 3.44 eV (Ref. 21). dAverage of CASPT2 S1 and S0 energies; S1-S0 energy gap in brackets.
(a) CIpi,„ (gd) = Qx, (b) CIpi,„ (ic) = Qi2 (c) С1рит (gd) (d) CIp„p (ic)
(j) в
FIG. 3. (a)-(h) Relevant conical intersection structures for fulvene, showing the mass-weighted displacements corresponding to the gd and ic vectors. (i) and (j) Torsion and pyramidalization coordinates.
In Eq. (1b), is the adiabatic electronic wave function for state i. The four conical intersections differ in the orientation of the CH2 group with respect to the ring: at CIperp, the CH2 group lies perpendicular to the plane (torsion angle of 90° around the double bond), CI63 has an intermediate structure between CIplan and CIperp (torsion angle of 63°), and CIpyr has the CH2 group pyramidalized with respect to the plane. From the point of view of energy, the energy of the seam with respect to CIplan decreases along the torsion and pyramidalization coordinates. CI63 coincides to a very good approximation with the global energy minimum of the seam, labeled CImin in Ref. 24, which is less than 1 meV lower in energy than CI63. As we discuss in this section, this space of intersection can be described by four coordinates (Figure 3): a totally symmetric and a non-totally symmetric bond stretch coordinate (QX1 and Qx2, respectively), which correspond to the branching space vectors at CIplan (Figures 3(a) and 3(b)), and the torsion angle p and the pyramidalization angle 0 (Figures 3(i) and 3(j)). QX1 and QX2 are treated as mass-weighted Cartesian displacements, while p and 0 are curvilinear coordinates corresponding to rotation of the rigid CH2 group around the axes ep and e0 in Figures 3(i) and 3(j), respectively (ep is the axis through the C1-C6 bond, while e0 is the axis through C6 in the plane of the ring, perpendicular to ep).
From the point of view of the electronic structure, the seam involves three electronic configurations, as shown by the analysis of the CASSCF wave function summarized in Figure 4 (orbital correlation diagram and configurations of the three lowest states along p). At planar structures of C2v symmetry, the ground state has A1 and the excited state B2 symmetry. Along the torsional coordinate the symmetry is lowered to A and B, respectively. At p = 90°, the structure has C2v symmetry, and the two lowest states have A2 and B1 symmetry. The 1A2 state at perpendicular geometries correlates diabatically with S2 at the planar structures, while the ground state at planar structures correlates with S2 at the structure with p = 90°. In principle such a situation could be described
as a three state problem where the two states of A symmetry are highly coupled along the torsional coordinate. However, the two lowest states are separated from S2 by at least 1.9 eV at the CASSCF level, and a two-state treatment is satisfactory for the present dynamics.
For the dynamics, the two-state potential energy matrix W was parametrized with respect to CASSCF(6,6)/cc-pvdz calculations over two states, state averaging with equal weights. This level of theory only gives non-dynamic electronic correlation, but complete active space second order perturbation (CASPT2) calculations on the critical points indicate that the effect of dynamic correlation is small (see Table I), with a maximum deviation of 0.3 eV between the CASSCF and CASPT2 relative energies. The CASPT2 energy gaps at the CASSCF conical intersection structures are
also smaller than 0.3 eV, and in all cases we have located the CASPT2 intersections in the vicinity of the CASSCF optimized structures. Therefore, the CASSCF level of theory is satisfactory for our purpose.
The potential energy matrix W along the four-dimensional vector of coordinates Q = {QX1, Qx2, y, 0} is given by a modified version of the diabatic vibronic coupling (2 x 2) Hamiltonian,33 where the diabatic basis is used to avoid singularities at the points of conical intersection. The electronic Hamiltonian W has the following form:
Vl,2 = W =
' Wi Wi2 N Wi2 W2
= W 0I +
Wii Wi2
Wi2 W22
Wii = -
f) O OOO
W0 = (A.0 + sin2 ((p))QXi + -(a>Xi + a>Xup sin2 (p))Q2 + sin2 (p)
i 4 (i i i 2 \ 2 + 2Çp sin4 (p) +1 2me + 2Vxi,e Qxi + 2mpß sin2 (p) I e2
( i i i 2 \ 4 i 2
+ 12 qe + 2 ÇxioQxi + 2 Çp,e sin2 (pn e4 + 2 ^ q^ ,
2 i 2 O i 2
[8k + 8Kxup sin2 (p)]Qxi + 2[5kxi + Syxi,p sin2 (p)]QXi + 2sYp sin2 (p)
i 4 ( i i i 2 \ 2
+ 2Xp sin4(p) + I 28Ye + 28Vxi,eQxi + 28Yp,e sin2 (p) I e2
( i i i 2 \ 4 i
+ ( 2Xe + öXxi,e Qxi + ^Xp,e sin2 (p) J e4 + -8yx2 Qi2
W22 Wi
W12 = «B sin2 (p) + kAb62 sin4 (p)) 9 + k£b QX2. (2e)
In Eq. (2a), V1>2 are the adiabatic state energies. For the diagonal elements of W, the W0 expansion (Eq. (2b)) is the average energy between the diabatic states, and W11 is half the energy difference. In our choice of the
3b, 10b 7b,
la2 12a la2
2b, 9b 3b, V=0° tp=45° (p=90° (a)
2b,2,1a22,3b,0
2b,2,1a22,3b,0
FIG. 4. Electronic wave function analysis along the torsional coordinate. (a) Correlating orbitals involved in the excitations. (b) Electronic states and dominant configurations.
terms contained in W0 and W11 we have aimed at obtaining a reasonably accurate fit with the smallest possible number of terms. Thus, the W0 expansion contains a linear and a quadratic term along QX1 (A.0 and &>X1 parameters), a quadratic term along QX2 (parameter &>X2), and a quadratic and a quartic term along 0 (parameters w0 and Z0). The torsional coordinate y is described with the sin2 and sin4 functions (wy and Zy parameters, respectively) because the potential is periodic in n along y. W0 also includes six bimodal terms between the different coordinates, which are the most relevant ones to obtain a good fitting of the surface. The W11 expansion has an analogous structure to W0. Turning to the coupling element W12, we include a single linear coupling term along QX2 (kA2b) because the length and direction of the ic vector is approximately constant along the surface (see below for details). Two additional terms are included in W12 to describe the bimodal coupling along 0 and y, since its value changes significantly in the space spanned by the two angles.
Before the details of the parametrization are discussed, a brief consideration on symmetry is helpful to understand how the elements of the diabatic potential energy matrix are approximated. At planar configurations of C2v symmetry, the four coordinates QX1, QX2, y, and 0 have a1, b2, a2, and b1 symmetry, respectively. The 1A1 and 1B2 states provide a natural choice for the diabatic states, and the diagonal
elements of W coincide with the adiabatic energies for cuts of the surface where Qx2 = 0 and p0 = 0. Therefore, all parameters can be fitted directly to the ab initio data, except those that involve terms that depend on Qx2 or couple p and 0. For the latter parameters we use the regularized diabatic states approach,34 where the couplings are expressed through offdiagonal potential terms of the diabatic Hamiltonian (W12). These terms are obtained using the unitary transformation that converts the adiabatic ab initio energies into the elements of W, as we explain below.
The origin of W lies at CIplan, of C2v symmetry, and the totally symmetric stretching coordinate QX1 corresponds to the gd vector at the intersection (Eq. (1a)). In our model, QXl also connects CIplan with the FC structure, since it is the only coordinate of a1 symmetry. The energy profile along this coordinate (surface cut for (Qx2,p,0) = (0,0,0)) in the model is compared with the ab initio energies of optimized points in Figure 2(a) (model energies in red and ab initio data in blue; see also the data in Table I). There is good agreement between the CASSCF data and the model for the energy of (S1)Min and the excited state energy at the FC structure, but there is a deviation of almost 1 eV for the vertical excitation energy (4.1 eV at the ab initio level and 3.2 eV in the model). This discrepancy is not due to the fit, but to the fact that we use the gd vector at CIplan as the totally symmetric coordinate to obtain a good description of the seam. A different coordinate might be chosen to improve the agreement between the relative energies, but this would cause a worse description of the seam. Therefore, inclusion of more totally symmetric modes would be the best way to improve the vertical excitation energy of the model. However, the early stages of the dynamics depend mainly on the S1 surface, which is well reproduced by the present model.
The Qx1 and p coordinates are highly coupled because of the high coupling between the S0 and S2 electronic states along the torsional mode (S2 is not included in the Hamilto-nian). To obtain the parameters for terms of QX1 and p we have followed a two-step procedure: in the first step, the parameters for the terms that depend only on p (wp, ZP,&YP,XP) are fitted with respect to ab initio calculations for geometries in the surface cut along p for (QX1,QX2,0) = (0,0,0). The parameters that couple Qxx and p (&xup, Zxup, Y ,p, Xxu p) are then fitted to ab initio calculations subject to the constrain that the energy of CIperp relative to CIplan and the second derivatives of the S0 and S1 energies at CIperp in the model match the ab initio values. The resulting surface cut in the (Qx1,p) plane for (Qx2,0) = (0,0) is shown in Figure 5. The root mean square error (RMSE) for this surface cut is 0.16 eV and 0.18 eV for V1 and V2, respectively (parametrization based on 102 ab initio points). The agreement between the model and the ab initio data is further assessed in Figure 6. In Figure 6(a) we compare the energy of the model seam along the torsion coordinate with the ab initio data, and in Figure 6(b) we compare the locus of the model seam in the (Qx1 ,p) plane (surface cut for (Qx2,0 ) = (0,0)) with the ab initio data. In both cases there is excellent agreement between the ab initio data (filled diamonds) and the model (continuous line). The agreement is confirmed by comparing the relative energies of CIplan, CIperp, and CI63, optimized at the CASSCF level, with the
E [eV]
3.8n 2,8 1,8-^ 0,8; -0,2-; -1,2 -2,2-J
FIG. 5. Plot of the model surfaces in the (Qx1,p) plane (cuts for (Qx2,0) = (0,0)). Dark surface: W11; light surface: W22.
energies of the stationary points of the model (see Table I). For the terms in 0, the parameters of W0 are obtained fitting the model to the average of the adiabatic energies. Following the regularized diabatic states approach, the parameters of W11 and W12 are fitted to the ab initio data using the equality:
AV2 = 4W22 + (W22 - Wn)2
Overall, the parameters corresponding to 0 have been fitted to 82 ab initio points, with RMSE of 0.16 eV and 0.15 eV for V1 and V2, respectively. The parametrization provides a good description of the seam segment along the pyramidaliza-tion mode, and the relative energy of CIpyr is well reproduced by the model (see Table I).
-Model
♦ ab initio
^[rad]
-10 --20 --30 --40 -50
^[rad] (b)
FIG. 6. Seam energy and locus of the seam in the (Qx1 ,p) plane (a and b, respectively; cuts for (Qx2,0) = (0,0)), comparing the model (continuous line) with the ab initio data (full diamonds).
-Model
♦ ab initio
TABLE II. Summary of propagation results.
Pulse Excitation Electric field Time for
duration t energy hrn strength Eq Maximal V2 Final V2 Final Wii product from Qxp
(fs) (eV) (a.u.) population P^ax population Pj population pp rj tivr (f
5 2 1.5 x 10-1 0.17 0.02 0.03 72
5 2.5 1.0 x 10-1 0.37 0.05 0.12 40
5 3 7.5 x 10-2 0.39 0.07 0.15 38
5 3.5 7.5 x 10-2 0.34 0.07 0.14 36
5 4 1.0 x 10-1 0.30 0.07 0.12 36
25 1.85 1.2 x 10-1 0.12 0.01 0.01 132
25 1.9 1.2 x 10-1 0.20 0.02 0.02 124
25 2 1.2 x 10-1 0.38 0.04 0.04 60
25 2.5 5.0 x 10-2 0.27 0.06 0.12 54
25 3 3.5 x 10-2 0.31 0.10 0.20 46
25 3.5 3.5 x 10-2 0.25 0.09 0.17 44
25 4 7.5 x 10-2 0.24 0.10 0.16 42
Finally, for QX2 the parameter KqB is the length of ic at CIplan (see Eq. (1b)). Similarly, the parameters mX2 and SyX2 are obtained from the curvatures of the adiabatic states
at CIplan.
In addition to the energy and locus of the conical intersection seam segments, the model also describes correctly the branching space along the seam. First, analysis of the conical intersection stationary points shows that the interstate coupling vector ic is similar in length and direction for all critical points on the seam (see Figures 3(b), 3(d), 3(f), and 3(h)), i.e., Qx2 provides a good approximation to the interstate coupling vector ic along the seam segments included in the model. The case of the gradient difference vector gd is more complicated, since this vector changes its direction along the seam.14 In general, the seam segments lie along combinations of two or more coordinates, and the gd vector is spanned by the same space. In the case of fulvene, for the segments along the rotation and pyramidalization coordinates the seam coordinate is a combination of QXl and y, and QXl and 0, respectively, and gd along each segment should be spanned by the same coordinate pair. Our model agrees with this analysis. Thus, according to Eq. (1), gd in the (QXi,y) plane (surface cut for (QX2,0) = (0,0)) is a linear combination of the vectors that span the plane, QX1 and y:
d (V2 - V1)
Qx2 =0,0=0
= [&K + 8KX1,y sin2 (y) + (SyX1 + SYxp,y sin2 (y))QxJQX1 + [ (2SKX1 ,y QX1 +SYxuy QX1 +Syp) sin (y) + 2/y sin3 (y) ]y■
A similar expression can be derived for gd in the (QX1,0) plane (surface cut for (QX2,y) = (0,0)):
d (V2 - V1)
Qx2 =0,y=0
1 2 1 4
Sk + SYX! Qxp + 2 SYx!,00 2 + 2 Xx!,004
These expressions are consistent with the ab initio results. For example, gd at CI63, which lies in the (QX1 ,y) plane, has a component along the CH2 torsion angle y (Figure 3(e)). Similarly, gd at CIpyr has a component along the CH2 pyra-midalization angle 0 (Figure 3(g)).
To summarize, the parameters of W have been fit to a total of 233 points, with RMSE values of 0.15 eV and 0.16 eV for y1 and V2, respectively. More importantly, the set of four coordinates {QX1, QX2, y, 0} provides a good description of the energies of the critical points, including the relaxation path on S1 along QX1, and the seam segments along y and 0. In future work, the four-dimensional model could be improved by including more terms in the expansion of W (Eq. (2)). However, it seems more important to include more coordinates to improve the agreement between the ab initio data and the model, in particular the vertical excitation energy (see Figure 2).
III. MCTDH PROPAGATION DETAILS
The propagations are carried out integrating the time dependent Schrodinger equation of the form:
d /K (t)
= (W + T + Hint)
K (t))
K (t))
where W is the potential energy matrix described in Sec. II. T is approximated by a diagonal kinetic energy operator (KEO) with terms:
t = E fQi, + fy + f0,
i = 1,2
Tq_ = -
2№9QX
(i = 1, 2),
Tk =---T a = y, 0),
? 2ik dk2
+ [(SY0 + SYx1,0QX1)0 + 2(X0 + Xxi,0 Qxi )0 310. (5)
where m is the reduced mass of QXi and is the moment of inertia corresponding to coordinate k. The KEO given above is obviously approximate. First, 0 involves a displacement of
the center of mass (see Figure 3(j)) which is neglected in this approximation, but the effect is small. We have also neglected the coupling between the normal coordinates and the two angles. The derivation of this coupling would lead to a very involved expression of the KEO, in a form which would not be apt for an optimal use of MCTDH, i.e., a direct product form in terms of the coordinates. Moreover, we assume that the omission of the coupling leads to an error that is smaller than the error due to the fact that we take into account only four degrees of freedom in the potential. In fact, using a mixture of normal and curvilinear coordinates along with the above approximate expression of the KEO is the price to pay for a correct description of the reaction path and the branching space with only four degrees of freedom and a simple expression of the KEO. A further improvement of our approach would be to reformulate the description of the system in terms of the so-called "polyspherical coordinates".35 A rigorous derivation of the corresponding KEO could then be performed, and the final expression would be perfectly adapted to a direct implementation in MCTDH. However, more than four degrees of freedom would be required to describe the seam and branching space correctly.
Hint describes the interaction between the wave packet and an external electric field (laser pulse) and has the form:
M42-E (t)
-^12E (t) 0
where ^12 is the transition dipole moment between So and S1. The value of /¿12 is estimated as 0.2222 a.u. at the CASSCF level (FC geometry) and is assumed to be independent of the coordinates, following the Condon approximation.
The dynamics were carried out with the MCTDH approach.3,4,28 For Qx1 and Qx2 we used a primitive basis of 48 and 25 harmonic oscillator functions, respectively, for p a fast fourier transform collocation of 128 functions, and for 0 a primitive basis of 25 sine functions. For the single particle basis, we have used a combined mode for Qx1 and p (50 single particle functions (SPFs) per diabatic state), as explained in Section 4.5 of Ref. 4, and 25 and 20 SPF per state for QXl and 0, respectively (i.e., a total of 25 000 configurations per state). The propagations, which were run over 150 fs, are well converged with respect to the number of SPF, since in all runs the maxima over time of the lowest natural weight for all modes did not exceed 0.001. Regarding the propagation length, test runs show that phase space saturation starts to appear at times longer than 150 fs. Inclusion of more degrees of freedom in our model seems to be necessary to avoid this problem and propagate for longer times.
The propagations were initiated exciting the wave packet from the vibrational ground state to the electronic excited state with a single laser pulse of the form:
E (t) = E0 cos (at) sin2 (—^ ,
where E0 is the intensity of the pulse, ha is the excitation energy, and t is the pulse duration (the pulse was set to zero outside the [0,t] interval). We ran two sets of calculations with a pulse duration of 5 fs and 25 fs, varying the excitation energy from 1.85 eV (0-0 transition) to 4.0 eV. For every
propagation E0 was tuned such that the maximum population of V2 (adiabatic excited state) lied approximately between 0.2 and 0.4. The resulting intensities lie between 0.03 and 0.15 a.u., which corresponds to 6 x 1013-2 x 1015 W cm-2, approximately. These intensities may not be practical because they may induce Coulomb explosion, but in the present work it is preferable to use these values to get significant population transfer to the excited state and facilitate the population analysis. In future work, it should be possible to apply less intense fields and extend the length of the pulse. Longer propagations will become reliable when more degrees of freedom are included in the model, as mentioned above.
Fulvene is symmetric with respect to the torsion coordinate, which implies that the nuclear wave function should be symmetric with respect to p .25-27 Thus, the ground state V1 has a double well potential along p with a high barrier, and the vibrational ground state has two doubly degenerate, symmetry adapted solutions (the energy difference is smaller than 10-6 a.u.). To monitor the double bond isomerization it is more convenient to propagate a symmetry broken, localized wave packet centered around p = 0, obtained from a linear combination of the symmetry adapted ones. In this case, we calculate the evolution of the expectation value of (p) in two intervals, pe[-n/2,n/2] (initial configuration) and pe[n/2,3n/2] (isomerized configuration). The expectation values were calculated for the adiabatic states following the procedure described in Ref. 36. In this procedure, the expectation values of the adiabatic state Vi can be calculated from the diabatic wave functions using Eq. (10):
pi* = (*d P'Wi*d
In Eq. (10), Cpstep is a step operator for the desired interval, and P(,) is the projection matrix in diabatic representation onto the ith adiabatic state. The matrix elements of the projection matrix are first generated on the product grid and then transformed to a product representation with the so-calledpotfit algorithm.4,37,38 The resulting adiabatic expectation values are interpreted as populations of the reactant configuration, p].ct and p2rct (pe[-n/2,n/2]) and product populations, pxpr and p2pr (pe[n/2,3n/2]).
IV. DYNAMICS RESULTS AND DISCUSSION A. Vertical excitation propagation
As a first approach to the photodynamics, and for comparison with the previous quantum dynamics studies,20,25-27 we carried out one propagation where the initial ground state wave packet is transferred to the excited state with an instantaneous vertical excitation. The results are plotted in Figures 7(a) and 7(b). These figures show the key parameters for the isomerization, the adiabatic populations of the reactant, and product configurations obtained from (p). Initially the wave packet is on V2 in the reactant configuration (p?ct = 1). The propagation has two phases. During the first 30 fs, approximately, there is a substantial population transfer to V1 (decrease of p2ct and increase of p1ct), without appreciable torsion around the double bond (ppr and ppr stay zero). The efficient deactivation can be readily understood from the
0.75 -
-— Pp
— pp
50 100
t [fs]
u* n 0.2 -
50 100
t [fs]
0.75 -
i>r p 2
—- pp
— p2
t [fs] (c)
FIG. 7. (a) Reactant and product populations for the two adiabatic states after instantaneous vertical excitation. p.ct (dotted line); p^ct (short hashed line); pp,r (long hashed line); ppr (continuous line). (b) Detailed plot of p and ppr . (c) Same as (a) for the diabatic states.
energy profile along Qx1 (see Figure 2). Thus, the wave packet acquires more than 1 eV initial vibrational energy along this coordinate, and the seam of intersection is energetically accessible. Initially the decay occurs in the vicinity of CIplan, since no population builds up for ppr. The second phase starts after around 30 fs, when the double bond starts to twist (rise of ppr and p2pr). For the initial value p = 0 the V2 surface has a negative curvature along p, and the inset of the torsion gives a measure of the time that is required for IVR from Qx1 to p. The value of 30 fs determined in the present case agrees approximately with the value given in Ref. 27 for propagation on a two-dimensional surface, where the torsion of the methylene groups takes ~50 fs. The part of the wave packet that spreads along p on V2 decays to V1 almost instantaneously, as seen by the rise of pp1 r, because the seam has a peaked topology around p = n/2, i.e., all the population around that
region is funneled almost immediately from the upper to the lower state. After 150 fs, 81% of the total population has been transferred to the ground state, where the reactant and product populations are 0.45 and 0.36, respectively. Thus, there is no selectivity in favor of unreactive or photochemical decay. Moreover, at that time the reactant and product populations on V2 have converged to the same value (~0.10), i.e., the wave packet has spread over p on V2 because of the flat energy profile along that coordinate. This implies that any control scheme must achieve the selectivity in the early stages of the dynamics, before the wave packet spreads on V2 and the selectivity is lost.
Our discussion is based on the adiabatic populations (see Eq. (10)). For comparison, in Figure 7(c) we show the reac-tant and product populations for the diabatic states. The differences are remarkable. For instance, the total population transferred to the ground state after 150 fs is 63% in the diabatic basis, compared to 81% in the adiabatic one. The adiabatic populations require an additional computational step and most publications therefore discuss only diabatic ones. However, given the substantial differences, we prefer to concentrate on the more physical adiabatic populations and will discuss only those in the following.
Overall, the results of the instantaneous excitation propagation are different from those described with previous approaches, where no double bond twist is observed after instantaneous vertical excitation of the wave packet. In the case of Ref. 26, where a three-dimensional model was used, the difference is probably be due to the different choice of the totally symmetric coordinate. In the present model, this mode (Qx1) contains a C1-C6 bond stretch component that is necessary for a good description of the surface, as seen from the changes in that distance along the relaxation coordinate (see Figure 2). Regarding Ref. 20, the different results obtained with the reduced dimensionality and the direct quantum dynamics approaches must be due to the different ways in which the IVR from the stretching to the torsional mode is treated, since both approaches use a similar electronic structure method for the potential. Regarding IVR, our reduced dimensionality propagations have two advantages: first, the dynamics is well converged (see Sec. III), and second, our curvilinear coordinates p and 0 are well adapted to the present problem. One may still question the fact that our model includes only four degrees of freedom. Although these are certainly the most relevant coordinates, inclusion of more degrees of freedom may slow down the IVR to the torsional mode, and this may account, in part, for the difference with respect to the direct dynamics approach. This issue will be addressed in future studies. On the other hand, the direct dynamics vMCG approach of Ref. 20 accounts for the full dimensionality of the problem but suffers from more approximate dynamics.
B. Pulse induced propagations
To approach experimental conditions, we have carried out propagations initiated with a single laser pulse. We have used pulses of 5 fs or 25 fs duration, and excitation energies
t [fs]
— - Pp
50 100 150
t [fs] (b)
FIG. 8. (a) Reactant and product adiabatic populations for the two states after excitation with a pulse of 25 fs and 3.0 eV (see Eq. (9)). p].ct (dotted line); plct (short hashed line); p.r (long hashed line); ppr (continuous line). (b) Detailed plot of ppr and ppr.
0.75 -
ft o ft
0.25 -
t [fs]
— Pp
t [fs] (b)
FIG. 9. (a) Reactant and product adiabatic populations for the two states after excitation with a pulse of 25 fs and 2.0 eV (see Eq. (9)). p].ct (dotted line); p2rct (short hashed line); ppr (long hashed line); ppr (continuous line). (b)
Detailed plot of ppr and ppr.
between 1.85 eV and 4.0 eV. Table II shows a summary of the propagation results. The competition between unreactive decay and isomerization to the product is monitored by determining the time required for 1% of product population to appear, which corresponds to the minimal time required for IVR from QX1 to y, tvr. This value is comparable between the different runs because the amount of population transferred to the excited state during the pulse is similar for all runs. This is monitored with the help of P^ax, the maximum value of the adiabatic excited state population during the propagation, which ranges from 0.23 to 0.43 (see Table II).
First we describe a representative propagation initiated with a pulse of 25 fs and 3.0 eV energy (Figure 8). During the pulse, the population on v2 starts to build up and pp2r reaches a maximum value of 0.31 (P2ax). After that, the propagation is similar to the vertical excitation case. Thus, a substantial amount of population is rapidly transferred back to v1, and after 150 fs the excited state population is reduced to approximately one third of P^ax (Pj = 0.10). Turning to the selectivity, the product configuration is only populated after ~40 fs, which corresponds to tvr. At that time, a large amount of population has already been transferred back from the excited to the ground state. As soon as the product population appears on V2, it also builds up on V., similar to the vertical excitation case. At the end of the propagation, the population of product on the ground state is 0.20. In turn, p2rct and ppr are approximately equal, i.e., the wave packet has spread along y on V2. Excitation with energies between 2.5 eV and 4.0 eV produces similar results to the ones described for the 3.0 eV
pulse. The maximal amount of population transfer to Vr is observed at 3.0 eV excitation energy, which is also the closest energy value to the vertical excitation energy from the FC structure in the model (3.2 eV).
Turning to the selectivity, the results for the runs with a 25 fs excitation pulse show that tivr becomes shorter as the excitation energy increases. This trend implies that IVR from QXi to y is accelerated when the amount of excess energy along Qx1 increases. More importantly, the increase in tivr becomes more pronounced as the excitation energy approaches the 0-0 limit (hrn < 2.0 eV). As a consequence, the isomerization is almost completely suppressed at low excitation energies because it cannot compete with the decay near CIplan, which is very efficient at all wave lengths. Such behavior is exemplified in Figure 9 for a run initiated by a 25 fs pulse of 2.0 eV excitation energy. During the propagation the wave packet is confined almost exclusively to the region around the initial value of y, and the total product population hardly reaches 0.04. The dependence of tivr on the excitation energy is also observed for the propagations initiated by short 5 fs pulses, although tivr is somewhat shorter compared to the 25 fs pulse for the same excitation energy. In contrast to this, the efficiency of the population transfer from the excited to the ground state does not depend so critically on the excitation energy (compare the values of Pj in Table II and Figures 8(a) and 9(a)). This is probably due to the fact that the barrier to reach CIPlan from (S1)Min is only 0.2 eV on our model surface (0.3 eV on the reference ab initio surface, see Table I). With such a small barrier, even the
propagations started close to the 0-0 transition (1.85 eV) show a very fast population decay (shorter than 10 fs). This differs from the results described using the vMCG approach on a full dimensional surface. In that case, propagation of a wave packet set initially in the region of (S1)Min induces a damped decay. The population transfer to the ground state is slowed down to about 70 fs, leaving time for the isomerization.
V. CONCLUSIONS
We have presented a four-dimensional model for the pho-todynamics of fulvene to model the competition between un-reactive return to the ground state and double bond isomer-ization. The photoreactivity depends on a seam of conical intersection along a torsional and a pyramidalization coordinate, and our model reproduces the energies and locus of the seam along these coordinates and the branching space vectors associated to the seam. The model involves several approximations concerning the KEO, the number of modes, and the level of theory for the electronic structure calculations, but wave packet propagations on this model have allowed us to determine the principal characteristics of the early stage of the photodynamics. Initially, the excitation energy is transferred to a bond stretching mode that leads the wave packet straight to a region of conical intersection around planar or quasi-planar nuclear configurations, and this leads to regeneration of the reactant. The photoisomerization only takes place after some energy flows from the bond stretching to the torsional mode, and therefore the competition between unreactive decay and isomerization depends on the rate of IVR between these two coordinates. In our propagations, IVR is accelerated for higher excess energies along the bond stretch mode. Thus, when the excitation energy is sufficiently high (2.5 eV or more in the present case, which is ~0.6 eV above the 0-0 transition energy), IVR takes less than 50 fs and the photoiso-merization competes with the unreactive decay. In contrast to this, for excitation energies near the 0-0 threshold (2.0 eV or less) IVR is slowed down by up to a factor of three, and the isomerization is virtually suppressed. Given the importance of IVR in the outcome of the decay, in future work we will inspect if the relationship between IVR and excitation energy that we have observed here is preserved when more degrees of freedom are included in the model.
From a general point of view, the photodynamics of ful-vene can be described as a two-step, two-mode process. Initially the excitation energy flows into the first mode, which leads to the unreactive decay region of the seam. If IVR is sufficiently fast, decay along the second mode can take place, which leads to the reactive region of the seam. This feature is quite general for photochemical processes. For example, in the polyenes used to model the photoisomerization of retinal, the chromophore of rhodopsin proteins,39 the first mode is a bond inversion mode of the polyene chain, and the second mode the cis-trans isomerization coordinate (double bond rotation).40 In benzene, the first mode is a bond stretching coordinate of the ring, and the second one a ring puckering mode (strictly speaking, a combination of two out-of-plane normal modes) that leads to formation of a prefulvene product
(channel 3 of benzene).18 In both cases, the photochemical event requires IVR from the first to the second mode. Based on our results, it will be interesting to investigate these processes focusing on how IVR affects the product formation rate and how this depends on the excitation energy.
Turning to our aim of controlling the photoreactivity, we can compare the present results with our previous MCTDH studies on intramolecular energy redistribution in molecules in their electronic ground state (Chapter 21 in Ref. 2). The basis of this approach was to create specific wave packets as linear combinations of converged vibrational eigenstates. In addition to the propagations initiated by single pulses we have tried a similar approach to control the photodynamics of fulvene, but due to the enormous density of states it has not been possible to converge any eigenstates in the vicinity of the conical intersection. Therefore, control of chemical processes involving conical intersections seems to be a much more difficult task. In such a context, the present work highlights the role of IVR in the photodynamics, and it suggests that controlling IVR will be necessary to control the photoreactivity.
ACKNOWLEDGMENTS
L.B. thanks the Generalitat de Catalunya (Spain) for project 2008-BE2-00316 and the Ministerio de Ciencia e Innovación (Spain) for project CTQ2008-06696BQU. Financial support from the German and French science foundations through the common project DFG-Me623/17 and ANR-09-BLAN-0417 is gratefully acknowledged by H.D.M. and F.G. L.B. thanks Oriol Vendrell for helpful discussions.
1 A. H. Zewail, in Femtochemistry, edited by F. C. De Schryver, S. De Feyter, and G. Schweitzer (Wiley VCH, Weinheim, Germany, 2001), pp. 1-86. 2H.-D. Meyer, F. Gatti, and G. A. Worth, Multidimensional Quantum Dynamics: MCTDH Theory and Applications (Wiley VCH, Weinheim, Germany, 2009).
3H. D. Meyer, U. Manthe, and L. S. Cederbaum, Chem. Phys. Lett. 165, 73
(1990).
4M. H. Beck, A. Jackle, G. A. Worth, and H. D. Meyer, Phys. Rep. 324, 1 (2000).
5F. Richter, M. Hochlaf, P. Rosmus, F. Gatti, and H. D. Meyer, J. Chem. Phys. 120, 1306 (2004).
6F. Richter, P. Rosmus, F. Gatti, and H. D. Meyer, J. Chem. Phys. 120, 6072 (2004).
7 F. Richter, F. Gatti, C. Leonard, F. Le Quere, and H. D. Meyer, J. Chem. Phys. 127, 164315 (2007).
8R. Marquardt, M. Sanrey, F. Gatti, and F. Le Quere, J. Chem. Phys. 133, 174302 (2010).
9G. J. Atchity, S. S. Xantheas, and K. Ruedenberg, J. Chem. Phys. 95, 1862
(1991).
10M. J. Bearpark, F. Bernardi, M. Olivucci, M. A. Robb, and B. R. Smith, J.
Am. Chem. Soc. 118, 5254 (1996). 11L. Blancafort, F. Ogliaro, M. Olivucci, M. A. Robb, M. J. Bearpark, and A. Sinicropi, in Computational Methods in Photochemistry, edited by A. G. Kutateladze (Taylor & Francis, Boca Raton, FL, 2005), vol. 13, pp. 31-110.
12M. Klessinger and J. Michl, Excited States and Photochemistry of Organic
Molecules (VCH Publishers, Inc., New York, 1995). 13D. R. Yarkony, J. Phys. Chem. A 105, 6277 (2001). 14L. Blancafort, B. Lasorne, M. J. Bearpark, G. A. Worth, and M. A. Robb, in The Jahn-Teller Effect: Fundamentals and Implications for Physics and Chemistry, edited by H. Köppel, D. R. Yarkony, and H. Barentzen (Springer, Heidelberg/Berlin, 2009), pp. 169-200. 15M. Araujo, B. Lasorne, A. L. Magalhaes, M. Bearpark, and M. A. Robb, J. Phys. Chem. A 114, 12016 (2010).
16D. Asturiol, B. Lasorne, G. A. Worth, M. A. Robb, and L. Blancafort, Phys.
Chem. Chem. Phys. 12, 4949 (2010). 17M. Kobylecka, A. Migani, D. Asturiol, J. Rak, and L. Blancafort, J. Phys.
Chem. A 113, 5489 (2009). 18B. Lasorne, M. J. Bearpark, M. A. Robb, and G. A. Worth, J. Phys. Chem.
A 112, 13017 (2008). 19A. Migani, L. Blancafort, A. D. DeBellis, and M. A. Robb, J. Am. Chem.
Soc. 130, 6932 (2008). 20D. Mendive-Tapia, B. Lasorne, G. A. Worth, M. J. Bearpark, and
M. A. Robb, Phys. Chem. Chem. Phys. 12, 15725 (2010). 21K. R. Asmis, M. Allan, O. Schafer, and M. Fulscher, J. Phys. Chem. A 101, 2089 (1997).
22P. J. Domaille, J. E. Kent, and M. F. Odwyer, Chem. Phys. 6, 66 (1974). 23 J. E. Kent, P. J. Harman, and M. F. Odwyer, J. Phys. Chem. 85, 2726 (1981). 24F. Sicilia, M. J. Bearpark, L. Blancafort, and M. A. Robb, Theor. Chem.
Acc. 118, 241 (2007). 25S. Alfalah, S. Belz, O. Deeb, M. Leibscher, J. Manz, and S. Zilberg, J.
Chem. Phys. 130, 124318 (2009). 26S. Belz, T. Grohmann, and M. Leibscher, J. Chem. Phys. 131, 034305 (2009).
27O. Deeb, M. Leibscher, J. Manz, W. von Muellern, and T. Seideman,
ChemPhysChem 8, 322 (2007). 28G. A. Worth, M. H. Beck, A. Jackle, and H. D. Meyer, The MCTDH Package, Version 8.2, University of Heidelberg, Heidelberg, Germany, 2000;
H.-D. Meyer, The MCTDH Package, Version 8.3, 2002 and Version 8.4, 2007; see http://mctdh.uni-hd.de. 29M. J. Frisch, G. W. Trucks, H. B. Schlegel et al., GAUSSIAN 03, Revision
B.02, Gaussian, Inc., Pittsburgh, PA, 2003. 30G. Karlstrom, R. Lindh, P. A. Malmqvist, B. O. Roos, U. Ryde, V. Verya-zov, P. O. Widmark, M. Cossi, B. Schimmelpfennig, P. Neogrady, and L. Seijo, Comput. Mater. Sci. 28, 222 (2003). 31B. O. Roos, and K. Andersson, Chem. Phys. Lett. 245, 215 (1995). 32G. Ghigo, B. O. Roos, and P. Ä. Malmqvist, Chem. Phys. Lett. 396, 142
(2004).
33H. Köppel, W. Domcke, and L. S. Cederbaum, Adv. Chem. Phys. 57, 59 (1984).
34H. Köppel, J. Gronki, and S. Mahapatra, J. Chem. Phys. 115, 2377 (2001).
35F. Gatti and C. Iung, Phys. Rep. 484, 1 (2009).
36M. R. Brill, F. Gatti, D. Lauvergnat, and H. D. Meyer, Chem. Phys. 338, 186(2007).
37A. Jäckle and H. D. Meyer, J. Chem. Phys. 104, 7974 (1996). 38A. Jäckle and H. D. Meyer, J. Chem. Phys. 109, 3772 (1998). 39 Y. Shichida and T. Yoshizawa, in CRC Handbook of Organic Photochemistry and Photobiology, 2nd ed., edited by W. Horspool and F. Lenci (CRC, Boca Raton (FL), 2004), vol. 2. 40M. Olivucci, A. Lami, and F. Santoro, Angew. Chem. Int. Ed. 44, 5118
(2005).
