Photochemical Ring-Opening of Cyclohexadiene: Quantum

Mar 8, 2010 - show that after irradiation the molecule arrives in the second excited state, the ..... second glance, the deviations of the experimenta...
0 downloads 0 Views 2MB Size
Download PDF
4036

J. Phys. Chem. A 2010, 114, 4036–4044

Photochemical Ring-Opening of Cyclohexadiene: Quantum Wavepacket Dynamics on a Global Ab Initio Potential Energy Surface J. B. Scho¨nborn, J. Sielk, and B. Hartke* Institut fu¨r Physikalische Chemie, Christian-Albrechts-UniVersita¨t, Olshausenstrasse 40, 24098 Kiel, Germany ReceiVed: September 29, 2009; ReVised Manuscript ReceiVed: February 6, 2010

We have assembled a global CASSCF potential energy surface for the excited 2A state of the cyclohexadiene -hexatriene system, in two degrees of freedom, with full relaxation in all other degrees of freedom. Quantum wavepacket dynamics on this surface yields simple interpretations of recent experimental data on the ultrafast photochemical ring-opening of cyclohexadiene as well as predictions on preferred product configurations. The feasibility of this system as a model for fulgide molecular switches is discussed. 1. Introduction The photochemical dynamics of cyclohexadiene is regarded to be a key element for understanding the photochemistry of more complicated compounds in the field of photoswitchable molecules such as furylfulgides. The photodynamics of cyclohexadiene have therefore been at the center of interest of various experimental and theoretical groups in the last few years, and research in this area is still ongoing. During various experimental studies using femtosecond mass spectrometry, Fuss et al.1 and Kosma et al.2 were able to gain major insights into the time constants of the electrocyclical ringopening of cyclohexadiene to hexatriene. Experimental data show that after irradiation the molecule arrives in the second excited state, the 1B state. It leaves this state driven by steep gradients within approximately 50 fs, reaching the dark 2A state in which the population has a lifetime of 80 fs. After leaving the 2A state, the molecule finishes the reaction in the ground state. To explain the experimental data, Celani et al.3,4 and Garavelli et al.5 performed detailed ab initio quantum-chemical investigations into the stationary points on the potential energy hypersurface (PES) of the electronic states involved, which resulted in minimum energy pathways for the reaction. On the basis of these calculations, Hofmann et al.6,7 were able to derive two reactive coordinates and set up a twodimensional (2D) PES consisting of two states, mainly the first excited and the ground state that contain the calculated reaction paths and rigid displacements orthogonal to these. By quantum wavepacket propagation they were able to show the importance of the incorporation of the conical intersection (CI) seam in the simulation. In further studies, the effects of barriers in the reaction pathways were investigated, and shaped laser pulses were recommended for optimal control schemes. Most of these studies were performed using the modified PES developed for cyclohexadiene but using kinetic energy Hamiltonians for related molecules and not for cyclohexadiene itself.8,9 Up to now, the quantum dynamics simulations on cyclohexadiene are highly influenced by the assumption that the reaction follows some kind of reaction path between precalculated stationary points around which the potential energy surface was set up. This * To whom correspondence should be addressed. E-mail: hartke@ phc.uni-kiel.de.

makes it hard to discriminate between true dynamical propensities of the system and influences of these a priori assumptions. In the present work we have used an alternative approach for describing the dynamical processes occurring during the photochemically triggered electrocyclic ring-opening of cyclohexadiene. Instead of limiting the description to the region close to the reaction path, we have constructed a PES that globally covers the full configuration space spanned by the two coordinates most relevant for the reaction, as far as it is accessible in typical experimental setups. All other degrees of freedom are taken into account via full geometry relaxations. Although the full relaxation carries the implicit assumption that all inactive degrees of freedom relax on a much faster time scale than the active ones, it is worth following this model, since it is the opposite limit of the reaction path approach used before, in which all coordinates except the reaction path are assumed to relax infinitesimally slowly. In contrast to other models used before, our approach does not contain any assumptions about locations important for the reaction, especially regions of high nonadiabatic couplings. CI seams arise naturally as points of lowest energy for a given combination of active coordinate values and are covered completely in our model. On the resulting PES, we have performed quantum-mechanical wavepacket propagations for the motions of the nuclei. From these, we extract data that can be directly compared to the experimental findings by Kosma et al.2 This comparison establishes direct links between the experimental measurements and the wavepacket movements, providing very simple explanations of the former and also predictions for the outcome of possible future experiments. This indicates that ultrafast chemical dynamics experiments should be accompanied not only by static analyses of the PESs involved but also by dynamcis simulations. 2. Theoretical Background A full-dimensional treatment of cyclohexadiene at any ab initio level of theory that reliably describes multiple coupled excited electronic states is not possible at present. Therefore, one has to settle for a lower-dimensional treatment, selecting both the active coordinates and the approximate model for the inactive coordinates in advance. We have chosen the same two active degrees of freedom (DOFs) as in previous work in the literature,3,4,6,7 with the following justifications: The CI consensus geometry in the literature as well as the most representative geometries in our CI region (see below) both feature a strong

10.1021/jp909362c  2010 American Chemical Society Published on Web 03/08/2010

Photochemical Ring-Opening of Cyclohexadiene

Figure 1. The coordinates q ) (r, φ) explicitly spanning our potential energy surface (for further explanations see text). One of the relaxed coordinates, ξ, is also shown; its importance is discussed in Section 3.2.

distortion of the 6 carbon atoms away from planarity, toward a screw-like arrangement. The single coordinate that embodies most of this distortion is the dihedral angle φ (involving carbon atoms No. 1, 4, 3, 6, and Figure 1) that symmetrically opens the CH2-CH2 σ-bond. Likewise, the interatomic distance r (between carbon atoms No. 1 and 4) is obviously involved in the same distortion since it serves to change the distance between the two “loose ends” of the screw-like arrangement just mentioned and hence it modulates the interaction between these two ends. At the same time, it offers the possibility to distort the molecule asymmetrically, serving as a counterpart to the symmetrical φ coordinate. Regarding the treatment of the other DOFs, it should be noted that keeping them fixed is not a good model for the whole course of the reaction. For example, the two carbon atoms next to the single bond that is about to be opened during the reaction have sp3 character in cyclohexadiene but then switch to sp2 hybridization in hexatriene. Therefore, at these C-atoms bond angles and dihedral angles involving CH bonds cannot be kept fixed without incurring unphysical energy penalties. Similar considerations apply to several other internal coordinates. In reaction path approaches of the system-bath style, such problems are avoided by relaxation of all passive DOFs. Typically, however, the active coordinate is constructed as lowest-energy path linking certain points on the PES deemed to be important (stationary points, CI minima, Franck-Condon points, etc.). Treating all other DOFs as harmonic oscillator baths during the dynamics is equivalent to forcing the system onto this preconceived reaction path, which may or may not be relevant for the true dynamics. Our present approach has similarities to and differences from the above cases. Similar to reaction surface approaches, we relax all other DOFs but we avoid further biases of the dynamics by performing a complete scan through a grid of the two active DOFs. Clearly, this corresponds to the assumption that motion in the passive DOFs is faster than in the active ones. This is likely true in some cases (e.g., high-frequency vibrations of CH bond lengths and bond angles) but can be questioned in others. The ultimate justification could only come from a higher-dimensional treatment, which currently is prohibitive. A weaker justification would be the a posteriori finding of agreement with experimental data. As we will show below, this is the case here; therefore, we have indications that our model may actually describe the reaction correctly. If some of the passive DOFs are approximately as slow as the active DOFs, it would be better to relax them only partially. However, in the traditional two-stage Born-Oppenheimer separated treatment, it is impossible to guess in advance just how far to relax which DOFs before doing the actual dynamcis. At the stage of solving the electronic structure problem, the only well-defined limits are the completely rigid model, the reaction path model, and the completely relaxed model. We take the latter approach for the reasons described above. Thus, full relaxation has the advantages of being well-defined and of avoiding unphysical energy penalties. In the presence

J. Phys. Chem. A, Vol. 114, No. 12, 2010 4037 of several electronic states, however, a disadvantage is that the relaxation can be performed consistently only for one of these electronic states; for all others we will not obtain meaningful PESs. For the present work, we focus on the 2A state where the reaction is in its decisive middle phase. Hence, we have to introduce assumptions about what happens directly before this middle stage (on the 2B state) and after it (in the ground state). Nevertheless, it should be emphasized that we have done stateaveraged calculations including these three states (1A, 2A, and 1B). Therefore, we can make reliable statements about energy differences between these three states. In particular, as depicted in Section 3.1, we do see large parts of the CI seam between 2A and 1A on our PES; hence, no additional assumptions are needed to include this CI into our description. Since we cannot directly include the actual passage of the molecular system from one PES to another in our dynamics (as described above, only one PES is present explicitly), we had to model the effects of the relevant CIs. The CI between the 1B and 2A electronic states determines how the system arrives at the 2A state; hence, the basic data of this CI are used to fix the initial conditions of our wavepacket on the 2A PES. This strategy of starting a dynamics simulation at a CI has recently been used by Sellner et al.10 for classical-mechanical trajectories. Similarly, the CI between the 2A and 1A electronic states determines the departure of the system from the 2A state; hence, we substitute its influence by the action of an imaginary absorbing potential (IAP), which we place in the coordinate space regions where we find 2A/1A degeneracy. This idea is better justified than it may seem at first sight: Our own experience with wavepacket dynamics at PES crossings as well as standard formulas (like the Landau-Zener approximation) tell us that for finite velocities (nonzero and noninfinite) we should expect neither complete absorption of the wavepacket (as a model for complete crossover to the electronic ground state) nor zero absorption (complete retention in the excited state). Therefore, what is needed here is an imperfect absorbing potential, which is of course much easier to set up than a (nearly) perfectly absorbing one. Also, the Landau-Zener formula and the adiabatic theorem tell us that wavepacket components with small velocities are less likely to cross over into the ground state. This feature is nicely mirrored by an imperfect absorbing potential, which is well-known to absorb the fast components better than the slow ones.11 Technical data for these settings (initial conditions and final absorption) are given in Section 4.1. Further justification for using an absorbing potential as a model for the effects of a CI is provided in the Appendix (Section 8.1). 3. Potential Energy Surface 3.1. Computational Details. The relaxed 2A surface was set up by a geometry optimization of the first exited state of a state-averaged CASSCF calculation (including the three states 1A, 2A, 1B) in Ahlrich’s TZVP basis12 using the Molpro program suite.13 The active space consisted of six electrons in six orbitals containing two pairs of π and π* orbitals in cyclohexadiene and -hexatriene, and the σ and σ* orbitals describing the bond cleavage in cyclohexadiene, which adopted π character in hexatriene. A comparison of a fully relaxed 2A/ 1A CI geometry with previous results is presented in Section 8.3. To further check the accuracy of the results, additional single-point calculations in a diffuse basis and larger active space were performed with the Molcas program,14 namely the 6-311 g(2df,p)15,16 basis was augmented by an additional diffuse p-function on each carbon atom, and the active space was

4038

J. Phys. Chem. A, Vol. 114, No. 12, 2010

Figure 2. The raw CASSCF data used to generate the 2A PES, depicted as final energy after relaxation in all other coordinates, at each (r,φ)-coordinate value pair.

Scho¨nborn et al.

Figure 4. Smoothed version of the PES depicted in Figure 2. The shaded region correponds to the CI seam in the (r,φ)-coordinates used here.

the initial excitation from the ground state to the 1B state (“FC”) and the location of the 1B/2A CI (“CI”). 4. Wavepacket Dynamics 4.1. Computational Details. The initial wavepacket for the quantum dynamics calculation was set up to optimally represent the physical conditions met in the experiment. The wavepacket enters the 2A potential surface through the 1B/2A CI. The initial wavepacket is a complex 2D Gaussian function (in au) Figure 3. Smoothed version of the PES shown in Figure 2, together with three representative molecular conformations, marking characteristic points on this PES.

2

Ψ(q, t0) )

j)1

extended to a CAS(8,12). These extensions did not lead to any significant differences for the first excited state. To take dynamic correlation into account, CASPT2 calculations on the results of both active spaces were done, also using Molcas. Although these CASPT2 data show larger deviations from the CASSCF results, they are still qualitatively very similar. These differences could therefore be further minimized by a geometry optimization at the CASPT2 level. Due to the lack of analytical gradients for this method in the used programs this was not feasible. 3.2. Results. The PES of the first excited state resulting from the calculations described in Sections 2 and 3.1 is shown in Figure 2. It exhibits several discontinuities, the most obvious one being the triangular-shaped step at intermediate values of φ and large values of r. The reason for this step lies in the relaxed dihedral angle ξ (cf. Figure 1), that is, at (r, φ)combinations corresponding to the step edge, an infinitesimal step in (r, φ) suffices to “flip” the molecule from one possible relaxed ξ value to the other one. Although the energy of the molecule is sufficient to enable this rotation, it is clearly unphysical in this instantaneous form. Furthermore, on the expected time scale of below 100 fs, this large-amplitude motion is presumably not important for the desired reaction. To arrive at a meaningful PES in this region, we have therefore eliminated this step, by effectively keeping ξ frozen. The remaining isolated discontinuities can be traced back to minor convergence problems in the CASSCF calculations and hence have been smoothed out, employing a bicubic spline approximation to our regularly gridded data. This smoothed PES is depicted in Figure 3. Further important features of our PES are indicated in Figure 4. In the shaded area, the ground state and the first excited state are nearly degenerate and feature large nonadiabatic coupling elements. This classifies that region of the surface to be part of a conical intersection seam. This seam is energetically lowest at this part of our 2D coordinate subspace. Additionally, Figure 4 shows the location of the Franck-Condon point for

(

∏ exp -

(qj - q˜j)2 σj2

)

+ ipj(qj - q˜j)

(1)

with its center at the above-mentioned CI, which has the coordinates q˜1 ) 2.84 Å and q˜2 ) 24.2° in our coordinate system q ) (r,φ). The width of the Gaussian was set to σ1 ) 0.15 Å and σ2 ) 12.0°, respectively. The energy of the CI is 184 kJ/ mol above the 2A state minimum. Additionally, the FC point for the initial excitation from the ground state to the 1B state is 37 kJ/mol above the 1B/2A CI. Thus, if all of this excess energy was channelled exclusively into our two active coordinates, the amplitude of the initial momentum should be chosen such that the total energy of our initial wavepacket is 221 kJ/mol. This, however, is both very unlikely and quite inconsistent with our model assumptions, which include relaxation of the other degrees of freedom on a time scale commensurate with the reaction. Hence, the total energy should be smaller, but without a full-dimensional treatment on all three electronic states it is very difficult to arrive at a trustworthy number. Therefore, we have adjusted the amplitude of the initial momentum to p1 ) -56 au and p2 ) 60 au, respectively, so that total energies of about 200 kJ/mol are accomplished. We have verified that our main results are insensitive to the exact choice of these numbers, for total energies in the broad range between 100 and 260 kJ/ mol (see below). For the direction of the initial momentum, it seems natural to take the direction from the FC point (optimized to r ) 2.87 Å and φ ) 15.6° in our coordinate system) to the CI point (Figure 4). This direction, however, has almost no component in the r-coordinate, it is almost exactly parallel to the edge of the 2A/1A CI. Given the featureless shape of the PES in the φ-coordinate, this would lead to repeated, only weakly hindered rotations in this φ-coordinate, not touching the CI seam. Such a long retention in the 2A state, however, is incompatible with all experiments. Therefore, we have to assume an initial momentum with a significantly larger r-coordinate component. In other words, we deduce from experiment that the motion of the wavepacket when entering the 2A surface through the 1B/

Photochemical Ring-Opening of Cyclohexadiene

J. Phys. Chem. A, Vol. 114, No. 12, 2010 4039

2A CI is directed largely toward the 2A/1A CI. One possible explanation for this could be that the initial gradients and/or the movement of the wavepacket on the 1B state (between the FC point and the 1B/2A CI) provide this additional momentum component. Another possible explanation could be that the notion of a single CI point is not really correct for molecules of more than three atoms. As shown for the 2A/1A CI in Figure 4, the 1B/2A CI should also be thought of as a rather extended region in configuration space, and correspondingly it is easily conceivable that larger parts of this CI seam (not just the lowestenergy point) are active in funneling the wavepacket from the 1B state to the 2A state. In fact, as we show in Section 4.2, this is exactly the case also for the 2A/1A CI. Therefore, we have chosen an initial momentum direction pointing more directly into the 2A/1A CI seam. This leads to wavepacket dynamics signatures in perfect agreement with experiment. As shown below, the essential characteristics of this dynamics are roboust against changes in the starting conditions of the wavepacket. In the Appendix (Section 8.2), we demonstrate that even an initial momentum pointing away from the CI seam (i.e., in the direction opposite to the one assumed here) is not changing our main results and conclusions. This initial wavepacket was propagated with a symplectic integrator of fourth order,17 using a converged time step of 0.01 fs and a total propagation time of 350 fs. The kinetic part of the Hamiltonian (kinetic energy operator, KEO) was calculated numerically using the Tnum package,18 as in earlier work (cf. refs 19-21 for details). Reduced masses for the KEO were calculated using the rigid model at the 1B/2A CI-structure. Spatial first and second derivatives of the wave function were calculated using the traditional fast Fourier transform (FFT) technique22,23 on a rectangular grid with 256 points in the r-coordinate and 384 points in the φ-coordinate, spanning a range from 1.9 to 5.0 Å in r and a range from 0 to 360° in φ, respectively. All dynamics calculations were performed with our own wavepacket propagation program suite MrPropa.24,25 As already mentioned in Section 2, the 2A/1A CI is modeled by an imaginary absorbing potential. We used the form given by Manolopoulos and Zhang:11,26

Vcap(r) ) iVIexp(-RIκ)

(2)

Figure 5. Negative imaginary absorbing potential to simulate the CI seam, in its (r, φ)-dependent form.

Figure 6. Snapshots of the wavepacket at t ) 0 fs (top), t ) 40 fs (middle) and t ) 180 fs (enlarged by factor 10). See text for further explanation.

where

rmax - r κ) r - r0

(3)

In eq 3, rmax denotes the end of the grid. The amplitude is RI ) 3.071, while the corresponding coefficient is set to VI ) -0.05ETmax with ETmax ) 350 kJ/mol. The edge of the seam is not a straight line in the (r,φ)-coordinates; therefore a φ-dependency was introduced in the absorbing potential by making the “starting point” of the absorbing potential r0 a function of φ according to the CI-seam obtained from our quantum chemistry calculations. The resulting absorbing potential is shown in Figure 5. The effect of this adjustment, however, turned out to be minor. We have not included a negative real part11 into this absorbing potential since its acceleration effect would be unphysical in this case. Of course, without having calculated all nonadiabatic coupling matrix elements, we cannot set the amplitude of the imaginary absorbing potential with any confidence; instead we treat it as an adjustable parameter in our comparison to the experimental data.

To avoid numerical instabilities, complex absorbing potentials are located on all other borders of the PES, extending the φ-grid by 20° at both ends. In contrast to the CAP simulating the CI, these absorbing potentials also employ a real part to accelerate slow parts of the wavepacket, because they are more difficult to absorb. This leads to a similar analytical form

Vcap(Q) )

∑ VR,jexp(-RR,jκj) + iVI,jexp(-RI,jκj) j

(4) with equivalent definition of κ as in eq 3. The already mentioned parameters VI, RI, and ETmax are unchanged, while RR ) 0.739 and VR ) -1.27ETmax. 4.2. Results. The initial wavepacket (Figure 6, top) is rather compact. In our setting, this is dictated by the initial conditions explained in Sections 2 and 4.1, but this also corresponds to the experimental observations: The system is very quickly (within 20 fs)2 accelerated out of the FC region on the 1B state, corresponding to the action of rather steep gradients there, and

4040

J. Phys. Chem. A, Vol. 114, No. 12, 2010

reaches the 1B/2A CI shortly thereafter (after less than additional 35 fs). Therefore, the wavepacket has had only little time to spread before it reaches the 2A state. Due to the large acceleration in the 1B state, the wavepacket initially moves very rapidly on the 2A state. This rapid movement is directed straightly into the 2A/1A CI seam, which starts not far away from the 1B/2A CI. In this initial phase of the wavepacket dynamics, we can corroborate the diagnosis by Fuss et al.2 who called this movement “ballistic” (i.e., largely unimpeded and uninfluenced by any features on the 2A PES). As we will see shortly, however, in later stages the character of the movement changes and the 2A PES shape does become important. In the 2A/1A CI seam, a relatively large part of the wavepacket is absorbed (crosses to the ground state); this first part is larger than the second part (see below) mainly because of the initially high wavepacket velocity. As shown in Figure 4, the 2A state PES rises up in the CI seam (parallel to the energetically degenerate ground state). Therefore, the wavepacket goes through a turning point and its remaining part exits again from the CI seam, toward larger values of the r-coordinate (Figure 6, center). From then on, further movement of the wavepacket is dictated by the shape of the 2A PES, with two deciding features: (1) In the φ-direction the PES is nearly flat, with only two comparatively shallow minima, whereas (2) in the r-direction the PES has the shape of a comparatively steep and rather anharmonic oscillator. Feature (2) leads to a marked oscillation of the molecule, sending it into and out of the CI seam while keeping the wavepacket fairly compact, while simultaneously feature (1) leads to a rapid spreading of the packet in the φ-direction (roughly parallel to the CI seam edge), so strongly that it covers the whole φ-space within our short overall propagation time. Therefore, at a later stage, the wavepacket re-enters the CI seam again (Figure 6, bottom). In contrast to the first time, it is now not compact anymore but very spread-out, and due to this internal vibrational energy redistribution (IVR) it also has a lower velocity in the r-direction. Therefore, in this second crossing event, only a smaller part of the remaining packet crosses over to the ground state. Additionally, we predict a very different conformational distribution of the reaction products: At the first crossover to the ground state, the still compact packet enters the CI with a narrow molecular geometry distribution that obviously resembles the reactant (cyclohexadiene) and even more the similar (helical) cZc-form of the hexatriene product. Therefore, we put forward the hypothesis that the product will be formed preferentially in this helical form. Of course, this hypothesis hinges upon the assumption that the ensuing ground state dynamics (which is not included in our modeling) only has little influence on the product conformation. At the second crossover event, however, the very broad molecular geometry distribution will also lead directly to the more “open” cEc-form of the hexatriene product.1 Obviously, a small part of the wavepacket will remain in the 2A state after this second crossover event with the ground state. It continues its dynamics as described above (spreading in φ and oscillation in r; not shown in Figure 6), leading to a third encounter with the CI seam. In order to compare to experiment we record the population, that is, the norm of the wavepacket, on the 2A surface. In the absence of more definite information, this is our model of the experimental detection process.2 For the mentioned starting conditions this results in the data shown in Figure 8, which are discussed in the next section. To substantiate the picture and

Scho¨nborn et al. time constants read from this graph, the propagation was repeated several times with different starting conditions of the wavepacket: The initial wavepacket width σ2 was varied by (8.0°, its momentum q2 by (20. σ1 and q1 were adjusted respectively in order to keep the total energy at around 200 kJ/mol. Its center was varied between 2.5 and 3.0 Å in q˜1, while q˜2 was varied between 16 and 32°. Additionally, the total energy was varied between 100 and 260 kJ/mol. Apart from this, the strength of the CAP was varied in order to simulate different sizes of nonadiabatic coupling elements. In all these tests, the main characteristics of the resulting time-dependent 2A population did not change, in particular the times of the two absorption steps, where the wavepacket leaves the 2A surface through the CI to the electronic ground state, remained almost unaffected. The relative heights of the steps, however, do change somewhat, mainly when the absorption strength is changed (which in our model stands for the 2A-1A coupling). Further tests are reported in the Appendix, Section 8.2. 5. Comparison to Experiment There are several observations in the experimental femtosecond mass spectrometry work by Kosma et al.2 which we can check without much effort. For example, these authors state that product formation starts already in the 1B state, goes on in 2A, but is not complete until some time after the molecule has reverted back to 1A. This is of course implicitly present in all of our wavepacket propagations, but it is also illustrated by several data features shown here, for example, in Figure 4, where the molecular structures representative for the 2A/1A CI obviously are rather far from both the reactant and the product configurations. Kosma et al.2 also mention that the frequencies they observe for what they call the “L3” stage of the photodynamical process (which is what we model here) must have a small restoring force and/or a large effective mass. We can corroborate a small gradient on our 2A PES in φ-direction (c.f., Figures 3 and 4) and also a fairly large effective mass of 107 amu (in the 1B/2A CI), which is not surprising given that increments in this coordinate have to move around large parts of the molecule. The FC-active coordinates identified by Kosma et al.2 (largely based upon the ab initio quantum-chemical calculations by Celani et al.3,4) are (1) the C-C stretching of the π-system, (2) a symmetric (conrotatory) CdC twist (together with some C-C torsion), and (3) the “erection” and stretching of the CH2-CH2 bond. To a large extent, these movements are captured by the coordinates r and φ used here. Since we have adopted these coordinates from the work of Hofmann et al.6,7 and have not attempted higher-dimensional treatments, this should be viewed as a priori input rather than as output of our work. What we simulate here essentially is the dynamics in region L3 of Kosma et al., that is, the red curve (in the time domain) of Figure 2b in ref 2, reproduced here as Figure 7. This graph should be almost directly comparable to the population in the 2A state in our wavepacket dynamics simulation (Figure 8). In this comparison, one should note that “time zero” is not identical. In the experiment, this is the pump pulse, that is, the first arrival of excited molecules in the 1B state. In our simulation, time zero is the advent of the wavepacket on the 2A state, which is (21 + 35) fs ) 56 fs later, according to the analysis by Kosma et al. This offset of 60 fs nicely matches the first maximum of the experimental signal in Figure 7. With this different zero of time taken into account, there is a very good agreement in the overall shape of the experimental and theoretical curves: Both consist of three downward steps,

Photochemical Ring-Opening of Cyclohexadiene

Figure 7. Experimental signature of the dynamics in the 2A state from femtosecond mass spectrometry (the black dots show the experimental ion signal that is assigned to the population of the first excited state) (Figure 2(b) of ref 2 reproduced by permission of the PCCP Owner Societies).

Figure 8. Population in the 2A state in our simulation, to be compared with the experimental data in Figure 7. See text for detailed discussion.

at similar times and of similar height. As mentioned at the end of Section 4.2, the relative heights of these steps are modified if the 2A-1A coupling strength (modeled by an absorbing potential) is changed. However, the times of all steps turn out to be insensitive to the initial conditions. This is easy to understand: As described in Section 4.2 (c.f., Figure 6), each downward step corresponds to the wavepacket moving into the CI region, while at the beginning of a flat part of the graph in Figure 8 the wavepacket leaves the CI region again for another excursion in the r-oscillation. Therefore, the timing of these steps is an intrinsic feature, and also an experimentally visible signature, of the shape of the 2A PES. In contrast, the height of the steps is a measure for the crossover probability from the 2A state to the ground state. As such, it is influenced by several factors, including the wavepacket velocity (distribution). This velocity depends on the preceding history of the wavepacket journey, in particular on the initial gradient in the 1B state and on the additional acceleration effects in the 1B/2A CI, both of which are not explicitly included in our simulation. Further discussion of these aspects is deferred to the Appendix in Section 8.2. We would like to emphasize that, in contrast to earlier simulations of the photochemical ring-opening of cyclohexadiene in the literature, the decisive features (“steps”) in Figure 8 occur on the correct time scale, as observed in experiment (in fact, they even match almost quantitatively). This provides some support for the claim that our modeling indeed has captured some decisive physical characteristics of this photoreaction. Additionally, we would also like to point out that our

J. Phys. Chem. A, Vol. 114, No. 12, 2010 4041 direct dynamical interpretation of Figure 7 differs in some respects from that given in ref 2. There is no need to invoke several frequencies building up the experimental signal; all of its main features are easily understood by the wavepacket movement depicted in Figure 6, which in turn is dictated in a simple way by the shape of the 2A PES. Of course, there is also no need to link these deduced frequencies to any normal modes of the molecule, obtained at stationary points that may or may not play a major role in the dynamics. On the other hand, the simple picture of a single crossover event from the 2A state to the ground state needs to be substituted by a periodic succession of several such events, directly linked to the observed oscillation pattern. Of course, there are also obvious differences between the experimental curve shown in Figure 7 and our theoretical one depicted in Figure 8. Most importantly, the regions between the steps are perfectly flat in the theoretical simulation while they appear more like oscillations in the experimental trace. The most likely reason for this difference is that the theoretical figure depicts nothing but the population in the 2A electronic state, without any influence of the probing scheme. In contrast, the experimental curve clearly is the outcome of the experimental detection mechanism. In other words, our simulated curve is tantamount to the assumption that the experimental probe mechanism is equally effective everywhere on the 2A PES but ceases to function as soon as the system has crossed over into the electronic ground state. While the latter may be true, given the experimental energetics, the former appears to be unlikely at first sight. In fact, similar oscillatory traces have often been interpreted with a strong emphasis on the wavepacket moving in and out of a probe window, and rather independent of the actual de-excitation process, both recently, as for example, in ref 27 and already in the beginning years of femtochemistry.28,29 Therefore, we have also attempted to model the experimental curve by assuming different probe efficiencies at different locations of the 2A PES (essentially by calculating crosscorrelation functions with target functions of suitable shapes). This, however, turned out to be far less convincing than the simple total 2A population shown in Figure 8. After all, at second glance, the deviations of the experimental steps from perfect flatness are rather small. This indicates that while a movement of the wavepacket into and out of regions of higher probe effectiveness may have some residual influence, the main effect visible in the experiment is a stepwise de-excitation from the 2A state to the ground state. This is also in line with the rather indiscriminate nature of the experimental probe process.1 Therefore, we are quite confident to have captured essential aspects of the photodynamics of the cyclohexadiene system in the 2A state, in spite of strong model assumptions (relaxation of all other degrees of freedom, on the time scale of the reactive events; neglect of details in the probe process; and absence of the initial dynamcis from the model: the excitation process, the 1B state, and the 1B/2A CI). 6. Cyclohexadiene As a Model System for Fulgides In previous studies, cyclohexadiene has frequently been viewed as a model system for larger molecules, in particular for fulgide molecular switches. From our experience, we would like to put in a word of caution here. In several stages of the electrocyclic ring-opening studied here, several C-C bonds acquire rather low-energy rotational barriers, and the available total kinetic energies after photochemical excitation certainly are sufficient to overcome these barriers. One case in point is the dihedral angle ξ mentioned in

4042

J. Phys. Chem. A, Vol. 114, No. 12, 2010

Section 3.2, which occasionally “flipped around” already during our relaxation of the inactive degrees of freedom. Another case in point is the low barrier for rotations corresponding to the active φ-coordinate. In compounds such as furylfulgides, some or all of these dihedral rotations are hindered or forbidden by additional ring systems fused onto the central cyclohexadiene moiety. Furthermore, the electronic excitation characteristics are qualitatively different: although CHD is excited to its second excited statesthe first one is dark at the Franck-Condon pointsthis is not true for fulgides, for example, for certain furylfulgides30 that can easily be excited to the first excited state, as can other fulgides. This excitation might lead to significantly different dynamics. Other important electronical differences such as a lowering in the intersection seam and the opening of new relaxation pathways have been reported by Nenov et al.31 Therefore, we do not view cyclohexadiene as a particularly good model system for fulgides. 7. Conclusions With this work, we were able to show that quantum dynamical wavepacket propagation on a fully relaxed potential energy surface for the 2A excited state of cyclohexadiene yields theoretical data in good agreement with the experimental results by Fuss et al.2 on the ultrafast photochemical ring-opening of this compound. Specifically, we are able to directly match the essential features of the raw 2A component of the experimental signal, namely the relative delay times of the steps/oscillations observed there. From our wavepacket dynamics, we can also give a simple explanation of these features: The initial “ballistic” movement on the 2A state drives the wavepacket directly into the CI seam where most of it crosses over into the ground state; the signature of this behavior is the initial rapid fall-off of the signal. This is followed by an anharmonic oscillation in the r-coordinate, driving the packet out of the CI seam but then back into it again several times, causing the second and third fall-off visible in the experimental data. Simultaneously, we also observe rapid wavepacket spreading in the φ-coordinate. On this basis we predict that the “first wave” of products will lead to closed-ring and helical conformations only (both of reactant and product), while at later times a much broader distribution should be observable, including various more “open” forms. In current quantum chemical and dynamical work, we are examining suitably substituted cyclohexadienes that may turn out to be better models for furylfulgide switches, or even better molecular switches by themselves. 8. Appendices 8.1. Imaginary Absorbing Potential as Model for the CI. We have used an imaginary absorbing potential (IAP) to model the (partial) population transfer away from a given potential energy surface, occurring in a region of strong nonadiabatic coupling to another potential energy surface. Although this is a simple idea, we are not aware of any previous instances of it in the literature. Therefore, this idea is further illustrated and justified in the following paragraphs. For this illustration purpose, we have arbitrarily employed the 1D model system of ref 32. Specifically, we have propagated a wavepacket on two diabatic potential energy surfaces

Scho¨nborn et al.

V(1, 1)(X) ) AX2,

V(2, 2)(X) ) B/X,

V(1, 2)(X) ) C (5)

with A ) B ) 0.1 and a constant diabatic coupling C ) 0.012 (explicit coupling). The resulting adiabatic potentials exhibit an avoided crossing as shown in Figure 9. The propagated wavepacket is transformed back to its adiabatic representation in order to record the population on the upper adiabatic potential energy surface. In a second calculation, we have propagated a wavepacket only on the upper adiabatic potential energy surface and have modeled the region of strong nonadiabatic couplings with an IAP, which is also shown in Figure 9 (simulated coupling). The absorbing potential was set up in order to simulate the nonadiabatic couplings in this region. It should be pointed out that none of the results shown are strongly sensitive to any features of this absorbing potential. Hence, it was a quick and simple matter to find an absorbing potential that leads to almost perfect agreement between the results of the full calculation and of the model calculation, as shown below. Similarly, in a situation where the exact result is not known (as in the main body of this article) it can be expected that the results will not depend sensitively on the features of the IAP. In general, when a wavepacket enters a region of strong nonadiabatic couplings, essentially two things happen: (1) A part of the wavepacket is transferred to the other potential energy surface, that is, it disappears from the first surface; (2) these two parts of the original wavepacket interact with each other as long as the nonadiabatic coupling is strong. Clearly, by modeling the coupling events by an IAP on a single surface, we are able to mimick the first effect. For both dynamics simulations (explicit coupling and simulated coupling), Figures 10 and 11 show the population on the upper adiabatic potential and two snapshots of the wavepackets, respectively. Obviously, the populations and the wavepackets before (t ) 0 fs) and after (t ) 15 fs) passing through the region of high nonadiabatic coupling are in

Figure 9. Adiabatic potential energy curves (red) and IAP (blue), which was used in the modeling example.

Figure 10. Comparison of the population of the wavepacket on the upper adiabatic potential.

Photochemical Ring-Opening of Cyclohexadiene

Figure 11. Comparison of two snapshots of the initial wavepacket (t ) 0 fs) and the wavepacket after passing the region of strong nonadiabatic coupling (t ) 15 fs). Both lines (red and blue) are almost perfectly on top of each other, hardly showing visible differences.

(almost) perfect agreement between the two simulations. The only differences are finer details of the interaction between the two parts of the wavepacket inside the region of strong coupling, specifically, the simulated coupling cannot reproduce the half-period of a Stueckelberg oscillation visible around t ) 7 fs. Nonetheless, we claim that we can establish an “effective” model also for the second effect, as we can still capture the end result, that is, not only the final distribution over the two PESs after the coupling between them has become weak (Figure 10) but even the shape of the wavepacket at that stage (Figure 11). Thus, in our perception, modeling a region of strong nonadiabatic coupling between two surfaces by a local IAP on just one surface is an appropriate model for the common situation of a localized region of strong surface coupling, as it is also present in the molecular system studied in the main body of this article. 8.2. Initial Momentum in Different Direction. It may seem that our ad-hoc assumption of an initial momentum pointing toward the 2A/1A CI seam is critical for a successful comparison to experiment. In this appendix, we present additional simulation data to demonstrate that the central features of this comparison remain unaffected, namely the relative delay times of the steps/oscillations in the 2A population (or in the experimental signal, respectively). Here, “relative” is meant to indicate that “time zero” is not welldefined (and, in fact, is affected by the choice of initial momentum), while the delay time differences between sucessive population falloff steps are occurring on the correct, absolute scale. To support these claims, we have performed additional dynamics simulations, with all settings as explained in the main text but with the initial momentum in the r-coordinate changed from -56 au to values between +50 and +60 au. Figure 12 shows how this affects our main result, that is, the populationvs-time graph of Figure 8. Obviously, the choice of initial momentum has a decisive influence on how long the wavepacket needs to first reach the 2A/1A CI. As we already discussed in the main body of the text, however, this time lag is not well-defined in our simulations anyway, since we do not (and by construction cannot) include the dynamics on the 1B state and the crossover from 1B to 2A. Similarly obviously, the choice of initial momentum has a decisive influence on the size of the falloff “steps”, that is, on how much population is transferred from 2A to 1A during each encounter with the 2A/1A CI seam. This is a simple realization of the standard expectation that the amount of population transfer scales directly with the size of the wavepacket velocity. Again,

J. Phys. Chem. A, Vol. 114, No. 12, 2010 4043

Figure 12. Population in the 2A state for test simulations with opposite direction and different sizes of initial momentum, to be compared with the data in Figure 8. See text for detailed discussion.

Figure 13. Position expectation value of the propagating wavepacket in the r-coordinate, for an initial momentum pr of -56 au (as in the main body of the text; red line) and of +55 au (blue line). See text for detailed discussion.

as already discussed, also this item it not well-defined in our model setup. Upon closer inspection, however, it becomes obvious that the time differences between successive falloff steps are hardly affected by the choice of initial momentum. The reason for this becomes obvious if the position expectation value of the wavepacket is followed over time. As Figure 13 shows, different initial momenta (including radically different directions) merely lead to a phase shift (and possibly minor amplitude changes) of what essentially remains the same oscillation, which thus reveals itself as an inherent feature of our 2A PES. Since these delay time differences are the items we are claiming to capture correctly in the main part of this article, we can also claim that initial momentum only has a very small influence here. Of course we have to rule out the “pathological” case that the initial momentum points exactly parallel to the (φ-averaged) CI seam edge, but this obviously is rather unlikely. Also, if (again more realistically) a distribution over an interval of initial momenta is assumed, this pathological case looses its destructive consequences. 8.3. Comparison of Stationary Points Obtained with the Method Used in this Study to Published Structures. Due to the different basis sets, a comparison of characteristic data calculated on the level used in this paper to previously published data is necessary to see how well the new data fits within the known framework. The most important feature of the PES presented in the main part of the paper (section 3) is the part of the conical intersection seam of the 2A/1A CI, since this part of the seam has not been subject to any investigations up to now. The nonsymmetric minimum structure obtained on the CASSCF(6,6)/TZVP level is compared to the one given in

4044

J. Phys. Chem. A, Vol. 114, No. 12, 2010

Figure 14. Comparison C-C interatom-distances of the 2A/1A CI in Å on CASSCF(6,6)/TZVP level to the values given by Garavelli et al.5 (latter in brackets).

the literature,5 where the same active space but a smaller one-electron basis set was used. Detailed analysis of these two structures, as partially illustrated in Figure 14, shows differences well within the limits of the method and basis sets, and no systematic deviations. Acknowledgment. Financial support by the German Science Foundation (DFG) via grant Ha2498/6 and via the Collaborative Research Center SFB 677 “Function by Switching” is gratefully acknowledged. References and Notes (1) Fuss, W.; Schmid, W. E.; Trushin, S. A. J. Chem. Phys. 2000, 112, 8347–8362. (2) Kosma, K.; Trushin, S.; Fuss, W.; Schmid, W. Phys. Chem. Chem. Phys. 2009, 11, 172–181. (3) Celani, P.; Ottani, S.; Olivucci, M.; Bernardi, F.; Robb, M. A. J. Am. Chem. Soc. 1994, 116, 10141–10151. (4) Celani, P.; Bernardi, F.; Robb, M.; Olivucci, M. J. Phys. Chem. 1996, 100, 19364–19366. (5) Garavelli, M.; Celani, P.; Fato, M.; Bearpark, M.; Smith, B.; Olivucci, M.; Robb, M. J. Phys. Chem. A 1997, 101, 2023–2032. (6) Hofmann, A.; de Vivie-Riedle, R. Chem. Phys. Lett. 2001, 346, 299–304. (7) Hofmann, A.; de Vivie-Riedle, R. J. Chem. Phys. 2000, 112, 5054– 5059.

Scho¨nborn et al. (8) Geppert, D.; Seyfarth, L.; de Vivie-Riedle, R. Appl. Phys., B 2004, 79, 987–992. (9) Voll, J.; Kerscher, T.; Geppert, D.; de Vivie-Riedle, R. J. Photochem. Photobiol., A 2007, 190, 352–358. (10) Sellner, B.; Barbatti, M.; Lischka, H. J. Chem. Phys. 2009, 131, 024312. (11) Ge, J.-Y.; Zhang, J. Z. H. J. Chem. Phys. 1998, 108, 1429–1433. (12) Scha¨fer, A.; Huber, C.; Ahlrichs, R. J. Chem. Phys. 1994, 100, 5829–5835. (13) Werner, H.-J. MOLPRO, version 2008.1, a package of ab initio programs; 2008. See http://www.molpro.net. (14) Karlstro¨m, G.; Lindh, R.; Malmqvist, P.-Å; Roos, B. O.; Ryde, U.; Veryazov, V.; Widmark, P.-O.; Cossi, M.; Schimmelpfennig, B.; Neogrady, P.; Seijo, L. Comput. Mater. Sci. 2003, 28, 222–239. (15) McLean, A. D.; Chandler, G. S. J. Chem. Phys. 1980, 72, 5639– 5648. (16) Krishnan, R.; Binkley, J. S.; Seeger, R.; Pople, J. A. J. Chem. Phys. 1980, 72, 650–654. (17) Gray, S. K.; Manolopoulos, D. E. J. Chem. Phys. 1996, 104, 7099– 7112. (18) Lauvergnat, D.; Nauts, A. J. Chem. Phys. 2002, 116, 8560–8570. (19) von Horsten, H. F.; Rauhut, G.; Hartke, B. J. Phys. Chem. A 2006, 110, 13014–13021. (20) von Horsten, H. F.; Hartke, B. Chem. Phys. 2007, 338, 160–167. (21) Sielk, J.; von Horsten, H. F.; Krüger, F.; Schneider, R.; Hartke, B. Phys. Chem. Chem. Phys. 2009, 11, 463–475. (22) Kosloff, R. In Numerical Grid Methods and Their Applications to Schro¨dinger Equation; Cerjan, C., Ed.; Kluwer: Dordrecht, 1993; Vol. 412; Ch The Fourier Method, p 175. (23) Kosloff, D.; Kosloff, R. J. Comput. Phys. 1983, 52, 35. (24) von Horsten, H. F. Ph.D. Thesis, Quantum dynamics of larger molecules using reduced-dimensional models; University of Kiel: 2008. (25) von Horsten, H. F.; Sielk, J.; Hartke, B. MrPropa - A Quantum Dynamical Wavepacket Propagation Program. See http://www.mrpropa.de. (26) Manolopoulos, D. Talk at the “Charles Coulson Summer School in Theoretical Chemistry”; University of Oxford: Oxford, 1996. (27) Hudock, H.; Levine, B.; Thompson, A.; Satzger, H.; Townsend, D.; Gador, N.; Ullrich, S.; Stolow, A.; Martı´nez, T. J. Phys. Chem. A 2007, 111, 8500. (28) Gerdy, J.; Dantus, M.; Bowman, R.; Zewail, A. Chem. Phys. Lett. 1990, 171, 1. (29) Hartke, B. Chem. Phys. Lett. 1990, 175, 322–326. (30) Siewertsen, R.; Renth, F.; Temps, F.; So¨nnichsen, F. Phys. Chem. Chem. Phys. 2009, 11, 5952–5961. (31) Nenov, A.; Ko¨lle, P.; Robb, M.; de Vivie-Riedle, R. J. Org. Chem. 2010, 75, 123–129. (32) Horenko, I.; Salzmann, C.; Schmidt, B.; Schu¨tte, C. J. Chem. Phys. 2002, 117, 11075–11088.

JP909362C