Nonadiabatic Molecular Dynamics Study of the cis–trans

ARTICLE pubs.acs.org/JPCA

Nonadiabatic Molecular Dynamics Study of the cis
trans Photoisomerization of Azobenzene Excited to the S1 State Marek Pederzoli and Ji
rí Pittner* J. Heyrovsky Institute of Physical Chemistry, Academy of Sciences of the Czech Republic, v.v.i., Dolej
skova 3, 18223 Prague 8, Czech Republic

Mario Barbatti Max-Planck-Institut f€ur Kohlenforschung, Kaiser-Wilhelm-Platz 1, D-45470 M€ulheim an der Ruhr, Germany

Hans Lischka* Institute for Theoretical Chemistry, University of Vienna, Waehringerstrasse 17, A 1090 Vienna, Austria

bS Supporting Information ABSTRACT: Ab initio nonadiabatic dynamics simulations of cis-to-trans isomerization of azobenzene upon S1 (n
π*) excitation are carried out employing the fewest-switches surface hopping method. Azobenzene photoisomerization occurs purely as a rotational motion of the central CNNC moiety. Two nonequivalent rotational pathways corresponding to clockwise or counterclockwise rotation are available. The course of the rotational motion is strongly dependent on the initial conditions. The internal conversion occurs via an S0/S1 crossing seam located near the midpoint of both of these rotational pathways. Based on statistical analysis, it is shown that the occurrence of one or other pathway can be completely controlled by selecting adequate initial conditions.

I. INTRODUCTION Azobenzene and its derivatives belong to most widely studied molecules.1
4 Its clean and reversible photoisomerization between the cis and trans form is unquestionably one of the most fundamental photochemical reactions. Photochromic compounds based on azobenzene and its derivatives can be used to trigger reversible changes of various properties upon irradiation, which is the basis for many valuable applications including optical switches,5
7 high density data storage,8,9 switchable superhydrophobic/superhydrophilic surfaces,10 and other materials with photomodulable properties.11
13 It is also widely used in molecular engineering, where it serves as an important compound of artificial molecular motors powered by light14,15 and other proposed light-driven nanoscaled devices.16,17 Despite considerable work in the field, the precise mechanism for the photoisomerization of azobenzene is still not completely clear. Originally, two pathways were suggested. An in-plane inversion of the NdN double bond and a rotation around the CNNC dihedral angle (see Figure 1). Based on the assumptions that rotational motion is not allowed in sterically hindered azobenzene derivatives and that after S2 excitation the molecule can relax to the S1 state, Rau deduced that the inversion path should be assigned to the S1(n,π*) excitation, while the rotational pathway should be assigned to the S2(π,π*) excitation.2 His assumption about sterically hindered azobenzenes has been, however, found incorrect, as even in stilbenophanes, where the inversion pathway r 2011 American Chemical Society

is not available, the photoisomerization was observed.18 The possibility of rotational mechanism in azobenzophanes was shown later by a semiempirical simulation.19 Rau’s picture of inversion mechanism seemed to be confirmed by the first theoretical calculation performed by Monti et al.,20 using a small STO-3G basis set. They found a substantial barrier on the S1 surface following the rotational pathway, while there was no barrier found for the inversion. These conclusions have been used for interpretation of experimental results (e.g., refs 21 and 22) for a long time. More accurate calculations by Cattaneo and Persico23 showed no barrier for the rotational motion and suggested that the rotational mechanism might also be possible. This suggestion was supported by the work of Ishikawa et al.,24 which reports a conical intersection near the midpoint of the rotational pathway and indicates that the rotational pathway should be the main channel for the isomerization process after S1 excitation. These results have also been confirmed by static calculations on complete active space self-consistent field (CASSCF) and complete active

Special Issue: Pavel Hobza Festschrift Received: February 9, 2011 Revised: May 4, 2011 Published: June 20, 2011 11136

dx.doi.org/10.1021/jp2013094 | J. Phys. Chem. A 2011, 115, 11136–11143

The Journal of Physical Chemistry A

ARTICLE

Figure 1. Conventional inversion versus rotation mechanisms.

space with second-order perturbation theory (CASPT2) levels,25 as well as using the density functional theory (DFT).26 Simulations based on direct semiempirical calculation of the electronic wave functions and potential energy surfaces (direct trajectories with surface hopping (DTSH) method)27 using both surface hopping method28 and full multiple spawning (FMS) method29 have been performed, concluding that the mechanism of the isomerization is purely rotational by both approaches. Similar results have been obtained employing the Car
Parrinello dynamics30 and molecular dynamics on the DFT level.31,32 Only very recently the first ab initio dynamics in the cis
trans direction has been performed by Ootani et al.33 employing the STO-3G basis set. In this work, for the first time the actual direction of the rotation was studied. The two nonequivalent rotational pathways are denoted as clockwise and counterclockwise based on the direction of the rotation of the central CNNC moiety in the initial phase of the isomerization process (see Figure 2). In the simulation,33 three different pathways have been observed: clockwise, counterclockwise, and mixed “counterclockwise S1 þ clockwise S0” rotation. Even though most of the recent studies agree that the favored mechanism of the cis to trans isomerization of azobenzene after excitation to the S1(n,π*) state is rotational, many important questions remain to be reliably answered. Most of the dynamical results so far come from semiempirical calculations or ab initio simulations using a small basis set. In this work we aim at verifying these previous results, expanding them and providing a clear picture of the dynamics of this important photoisomerization reaction while paying special attention to the aspects important for future applications, mainly (i) various isomerization pathways and direction of the rotational motion, (ii) importance of the inversion coordinates, and (iii) dependence on the initial conditions. For these purposes, we employ nonadiabatic dynamics based on high-level ab initio potential energy surfaces and detailed statistical analysis of the obtained data.

II. COMPUTATIONAL DETAILS Complete active space self-consistent field (CASSCF) and multireference configuration interaction with single excitations (MR-CIS) were used to describe the electronic structure of azobenzene. For the CASSCF method, an active space of 10 electrons in 8 orbitals has been selected (CAS(10,8)) as a compromise between computational cost and accuracy of the description of excited states of interest. This selection includes two n, three π, and three π* orbitals in the Franck
Condon region. The state averaging procedure has been employed, averaging over three states (SA-3) with equal weights for each state. Subspaces of the CAS(10,8) were used as reference spaces

Figure 2. Clockwise versus counterclockwise rotation mechanisms.

for subsequent supporting MR-CIS(6,5) and MR-CIS(4,3) calculations. The 24 lowest orbitals were kept frozen in the MR-CIS calculations. The basis sets 6-31G* for the nitrogen and hydrogen atoms and the 6-31G for the carbons were used. Our choice of methods and basis sets has been tested by a series of single point calculations and geometry optimizations, including the optimization of the conical intersection (see Results and Discussion). The energies, energy gradients, and nonadiabatic coupling vectors were calculated analytically34
38 using the COLUMBUS package.39
41 Ab initio nonadiabatic surface hopping dynamics42 has been carried out. In this mixed quantum-classical approach, the motion of the nuclei was represented by classical trajectories computed by numerical integration of Newton’s equations using the velocityVerlet algorithm with a 0.5 fs time step.43 All required quantities were interpolated between two time steps, and the time-dependent Schr€odinger equation was integrated using the Butcher algorithm with a 0.025 fs time step.44 The nonadiabatic effects were taken into account by means of Tully’s fewest switches method.45 Decoherence corrections were included using the scheme described in ref 46. Energies, gradients, and nonadiabatic coupling terms were computed on the fly at each time step at the SA-3-CASSCF(10,8) level of theory in the aforementioned basis set. In total, 70 trajectories have been computed for a maximal simulation time of 500 fs. In 50 trajectories, the nonadiabatic couplings were computed with the finite difference methods described in refs 47
49. In the remaining 20 trajectories (started from a common subset of initial conditions), nonadiabatic couplings were computed using the analytical procedures for comparison. The initial geometries and momenta were selected from a Wigner (quantum harmonic oscillator) distribution in the 11137

dx.doi.org/10.1021/jp2013094 |J. Phys. Chem. A 2011, 115, 11136–11143

The Journal of Physical Chemistry A

ARTICLE

Table 1. Geometric Parametersa geometry cis S0

N
N

N
C

N
N
C

C
N
N
C

MR-CIS(4,3)

125.3

144.1

122

CASSCF(10,8)

124.6

143.7

123

3

exp.(X-ray)55

125.3 ( 0.4

144.9 ( 0.8

121.9 ( 0.3

8 ( 0.3

MR-CIS(4,3)

126.0

142.9

114

180

CASSCF(10,8)

126.0

142.9

114

180

exp.(GE)54

126.0

142.7

113.6

180

exp.(X-ray)55

124.7 ( 0.4

142.8 ( 0.8

N/A

180

trans S1

MR-CIS(4,3) CASSCF(10,8)

127.2 125.5

137.8 137.3

126 128

180 180

CI S0/S1

CASSCF(10,8)

126.3

135.8/140.4

118/136

93

trans S0

a

method

2

Bond lengths are given in pm, angles in deg.

Table 2. Vertical Excitation Energies method

trans

cis

ref.

CAS(10,8)

3.24

3.36

present work

MR-CISD CIPSI

3.11 2.81

3.95 2.94

24 23

in solution

2.8

2.87

57

in gas phase

2.82

2.92

58

Experiment

ground state.50 The trajectories were started in the S1 state of these geometries. The dynamics simulations, including the initial conditions generation were carried out using the NEWTON-X program51,52 employing the COLUMBUS program for the ab initio calculations. After the statistical analysis of the results obtained from the dynamics, an additional set of 34 shorter trajectories has been calculated to confirm our hypotheses about the initial conditions dependence of the isomerization direction. The statistical analysis of the initial conditions and hopping geometries has been performed using the tools and programs included in the NEWTON-X program package and the program Statistica 7.53

III. RESULTS AND DISCUSSION A. Static Calculations: Investigation of Potential Energy Surfaces. Geometry optimizations at CASSCF and MR-CIS

levels have been performed to locate the optimal geometries of relevant structures. The ground state minima of both isomers, the S1 minimum of the trans form and the S0/S1 conical intersection have been determined. We could not locate a S1 cis local minimum; thus, there is no barrier between the vertically excited geometry and the conical intersection. The resulting geometries are shown in Figure 2. The Cartesian coordinates and more detailed illustrations of the optimized structures are provided in the Supporting Information. The ground state minimum of the trans isomer is planar and has C2h symmetry. The minimum of the cis isomers lies about 0.7 eV higher in energy, it is nonplanar and has a C2 symmetry. Geometry parameters and symmetry of the molecules are in good agreement with experimental data.54,55 The most important geometric parameters are presented in Table 1. The conical intersection is located in the midway of the rotational pathway at the CNNC angle close to 90°, in agreement with previous investigations.24,25,27,56

Figure 3. LIIC curve for S1 and S0 states at SA3-CAS-SCF(10,8) level.

Table 2 presents vertical excitation energies, obtained at the same level as used for the dynamics, in comparison with experimental data and previous results obtained earlier by computationally more demanding methods, multireference configuration interaction with single and double excitations (MR-CISD),24 and multireference perturbation CI method (CIPSI).23 Several other active spaces and basis sets, including basis sets with diffuse functions, have been tested prior to the dynamics (see Supporting Information for more data). It can be seen that CASSCF overestimates the excitation energies by about 0.4
0.5 eV in comparison to CIPSI. An approximate isomerization pathway between the cis ground state minimum, the conical intersection and the trans ground state minimum has been computed by means of linear interpolation of internal coordinates (LIIC). The potential energy curve (Figure 3) confirms that there is no barrier in the rotational pathway starting from the cis isomer. Also, the curve is much steeper on the cis side, which should lead to shorter times for the cis
trans isomerization and considerably longer for the trans
cis, which was indeed confirmed experimentally22 and is also observed in the potential curves for azomethane.59 The character of both S0 and S1 states has been examined. The main contributions for these states are the ground state configuration and the n1(π*)1 excited configuration, respectively. We found out that a double excited configuration S(n,n;π*,π*) plays an important role in the isomerization process. As the CNNC angle increases, the double excited S(n,n;π*,π*) configuration becomes more important in both S0 and S1 states. Near the conical intersection, both S0 and S1 are a mixture of n2(π*)0, n1(π*)1, and n0(π*)2 configurations. B. Photodynamics of the Isomerization. The nonadiabatic dynamics of azobenzene has been carried out, starting from the 11138

dx.doi.org/10.1021/jp2013094 |J. Phys. Chem. A 2011, 115, 11136–11143

The Journal of Physical Chemistry A

ARTICLE

S1 state of the cis isomer. A total of 70 trajectories have been computed. The resulting average adiabatic populations are shown in Figure 4. The trajectories using analytical and overlap-based nonadiabatic coupling terms gave very similar results. A detailed comparison of these methods has been published earlier.48 The isomerization lifetime was obtained by fitting the S1 population with the function f(t) = e
(t
td)/τ, where td is the onset of the nonadiabatic events and τ is the time-constant for the exponential decay. The lifetime is obtained as τ þ td and yields 67 fs, which is in good agreement with most recent experimental values (see Table 3). The quantum yield, estimated as a ratio of the number of trajectories where the isomerization to the trans form occurred and the total number of trajectories is somewhat higher, ϕ = 0.65 ( 0.11, than the experimental values obtained in the condensed phase given in Table 3. Both direct observation of trajectories and the analysis of most important modes of motion obtained using the principal component analysis (PCA)60,61 have shown that the isomerization occurs as a pedal motion of the central CNNC moiety, while the (much heavier) benzene rings remain more or less parallel to each other during the actual isomerization. This type of motion, which was previously described for the isomerization of the related stilbene molecule,62 also allows isomerization in some azobenzene derivatives where the motion of the rings is restricted.19

Upper right panel of Figure 5 shows the evolution of potential energies of the S1 and S0 states in a typical trajectory that involved isomerization to the trans form (reactive trajectory); the upper left panel shows the same time evolution for a trajectory that remained in the cis form (nonreactive trajectory). The current state at each time is indicated by the thick line. One can see that in both cases the trajectories quickly approach the crossing seam. After passing the conical intersection, the molecules remain in the ground state for the rest of the simulation. The lower panels of Figure 5 show the evolution of the most important geometry parameters in the aforementioned trajectories. The first two of these parameters are the CNN and NNC angles. No inversion motion is observed. The isomerization does not depend strongly on these parameters. The next part of Figure 5 shows the time dependence of the central CNNC dihedral angle and CCNN and NNCC dihedrals, further denoted as ring rotations. As expected, the CNNC angle shows the largest difference between the reactive and nonreactive trajectories and can be used to keep track of the progress of the isomerization. The CNNC dihedral angle is also closely related to the length of the molecule, which is an important property for certain applications.63 The change in the length of the molecule during the isomerization is shown in Figure 6. The length, computed as the maximum distance between any two atoms, is smallest near the cis minimum. During the isomerization process it gradually increases to its maximum, which corresponds to the planar trans form. Because of relatively large momenta of the benzene rings and rather flat potential energy surface near the trans minimum, the molecule passes through the trans geometry and continues further in its rotational motion. Only after about 400 fs the molecule begins to return to its planar form. The motion of the benzene rings shows another interesting feature. The rings motions stay synchronized at the early stages of all trajectories. For nonreactive trajectories this synchronization is only slightly distorted. On the other hand, in reactive trajectories, the synchronization disappears completely after passing the CI (an illustration of this behavior is given in the Supporting Information). In the later stages, the rings start to move in an opposite fashion, as can be seen also in the third panel in Figure 5. The lowest part of Figure 5 shows the development of o the NN bond length, which has also been found very important. In the trajectories shown in Figure 5 it can be seen that the NN distance at the hopping time is larger in the reactive trajectory. This result holds for most of the trajectories. The subsequent statistical analysis (see also the next section) confirmed that the reactive trajectories have significantly (p = 0.013, the so-called p-value is a standard statistical measure for probability

Figure 4. Average adiabatic populations in the S1 and S0 states as a function of time.

Table 3. Quantum Yields and Lifetimes (in fs) of cis to trans Photoisomerization after S1 Excitation method

quantum yield

lifetime (fs)

ref.

ab initio SA3-CASSCF(10,8), 6-31G*/6-31G

0.65

67

ab initio SA5-CASSCF(6,4), STO-3G

0.45

60

33

semiempirical, SH

0.61

≈50

23

semiempirical, FMS

0.68

≈50

29

≈100 100

21 64

Molecular Dynamics present work

Experiment femtosecond transient absorption (in ethanol) transient absorption/steady state emission (in DMSO) measurement in low polarity, low viscosity solvent

0.40
0.56 11139

65
70 dx.doi.org/10.1021/jp2013094 |J. Phys. Chem. A 2011, 115, 11136–11143

The Journal of Physical Chemistry A

ARTICLE

Figure 7. CNNC dihedral angle as a function of time (hopping events are indicated with circles, and the horizontal lines indicate the CNNC dihedral angle of the conical intersection).

Figure 5. Potential energies and the most important geometrical parameters in a typical nonreactive (left) and reactive (right) trajectory.

Figure 6. Evolution of the length of the molecule during isomerization.

of obtaining a result at least as extreme as the one that was actually observed, assuming that the null hypothesis is true) larger NN bond length than nonreactive trajectories. This important result indicates that it might be possible to directly influence the quantum yield of the cis
trans photoisomerization through the excitation of the NN vibration during the initial stage of the isomerization process. We will investigate this possibility in our future work.

The isomerization process in individual trajectories is most easily followed by investigating the evolution of the CNNC dihedral angle in time, as shown in Figure 7. Two important points can be readily seen from this figure. First, none of the geometries remained planar in the process of isomerization, which excludes the purely inversion mechanism. Second, the isomerization proceeds mainly through the counterclockwise direction channel; only a small fraction of trajectories proceeds through the clockwise channel. In contrast to the work of Ootani et al.,33 no mixed counterclockwise þ clockwise motion has been observed. The sample trajectories for all four pathways (reactive and nonreactive in both clockwise and counterclockwise directions) are given in the Supporting Information. To see whether inversion plays any role in the isomerization process, we analyzed the values of the CCN and NNC angles during the dynamics. The differences between time evolution of these parameters in the reactive and nonreactive trajectories are very small, and the consequent analysis showed that they are not statistically significant (for more details, see Supporting Information and also the next section). Based on these results, we conclude that the isomerization process is basically independent of the values of angles characteristic for inversion and confirm the previous conclusions that rotational mechanism dominates after the n-π* excitation.25,28
30 Aiming at explaining the preference for the counterclockwise path, new potential curves have been computed. To account for the reduced motion of the benzene rings, average hopping geometries (obtained by averaging the coordinates of hopping geometries either in clockwise or in counterclockwise direction) have been taken as the final points for the LIIC procedure. The resulting LIIC curves for clockwise and counterclockwise rotations have been merged into curve shown in Figure 8. The clockwise direction shows a notable barrier, which accounts for the much lower rates of this isomerization channel. This barrier is mainly caused by obstructing benzene rings, as can also be seen from Figure 2. C. Analysis of the Influence of Initial Conditions. The statistical analysis has been performed using the one-factor and two-factor analysis of variance (ANOVA) method.71 The trajectories have been divided into several groups, based on whether the isomerization occurred and also on the direction of the rotation (clockwise, counterclockwise). Five trajectories that finished 11140

dx.doi.org/10.1021/jp2013094 |J. Phys. Chem. A 2011, 115, 11136–11143

The Journal of Physical Chemistry A

ARTICLE

Figure 9. Analysis of variance; illustration of confidence intervals.

Figure 8. LIIC curve for clockwise and counterclockwise rotations at SA3-CAS-SCF(10,8) level.

too early to clearly distinguish whether they would end up in cis or trans isomer have been excluded from the total of 50 trajectories starting with different initial conditions. For each trajectory the following set of variables has been tested: the most important internal coordinates, two velocity projections, the initial kinetic energy and the vertical excitation at the initial geometry. The two velocity projections were defined as: (i) projection of the initial velocities into the direction toward the average hopping geometry for all atoms and (ii) the same projection for nitrogen atoms only. All these variables should have close to normal distribution among the initial conditions, which is necessary for the ANOVA method to be applicable. The validity of the statistical model has been also confirmed with nonparametric tests. The tests have been done in three stages. First, we tested whether there is some connection between the variables mentioned above and the fact whether the isomerization occurred or not. As expected, no significant connection between the initial conditions and the occurrence of the isomerization has been found. Next, the influence of all the variables on the isomerization direction has been tested. Statistically significant influence has been observed in several internal coordinates, most notably the ring rotations, and CNNC dihedral angle. The effect of these variables is shown in the left part of Figure 9. Several other internal coordinates (mostly ring deformations) also proved to be relevant. The most likely explanation is that the motion against the particular deformation in some cases helps to push the molecule to either direction. The total projection of velocities of all atoms in the molecule was not found to be statistically significant (p = 0.22), but for the projection of the velocities of nitrogen atoms only, we obtained a significant (p = 0.03) correlation. The visual comparison for values of the nitrogen-velocity projections are shown on the right-hand side of Figure 9. This result is in agreement with previous findings19,30 that only the motion of the central NdN region is important for the first stage of the isomerization. In the last stage, we have considered both factors (occurrence of isomerization and direction of the rotation) at once. The results for ring rotations and CNNC dihedral angle did not show any new information. The results for projection of velocity of nitrogens, however, are different and suggest that this parameter might have an influence on the quantum yield. In other words, if we can separate the initial conditions that favor clockwise and counterclockwise motion, then, we probably can also, to some extent, influence the quantum yield of the isomerization by controlling the velocity of the nitrogen atoms.

To confirm our results, we generated a new set of 13 trajectories, with initial conditions generated in a way that should promote the clockwise mechanism. Six of these geometries are based directly on the individual initial geometries of trajectories that proceed via the clockwise channel and another six geometries were based on their average; the last geometry used was the average itself. New geometries have been generated from reference ones (six individual geometries and the average) by adding random numbers having a normal distribution with variance of either onehalf or the same as value observed among the individual trajectories. Using these 13 geometries, two new sets of trajectories have been created. For the first set the initial velocities have been set to zero. For the second set the velocities obeying the Wigner distribution were used. These velocities have been taken from one of the earlier trajectories that proceeded via the clockwise channel and have a positive value of projections to the clockwise isomerization direction for nitrogen atoms. The trajectories have been computed and stopped 10 fs after the hopping occurred. The maximum simulation time has been set to 100 fs. Figure 10a shows the time evolution of the CNNC dihedral angle for the first set of new trajectories. Because in this set the initial velocity was set to zero, apparently even the change in the geometry alone is sufficient to greatly modify the ratio of counterclockwise and clockwise pathways. By controlling both initial geometry and initial velocity, see Figure 10b, it is possible to obtain exclusively clockwise rotations. The results are summarized in Table 4. A similar procedure could have been used to maximize counterclockwise rotations. D. Conclusions. An ab initio nonadiabatic molecular dynamics simulation of azobenzene isomerization has been carried out at SA3CASSCF(10,8) level. The results of static calculations and geometry optimizations of local minima and crossing seam indicate that the selected methods and basis sets are capable of describing the structures and mechanisms correctly at a computational cost, which allowed us to run enough molecular dynamics trajectories. A total of 20 trajectories were calculated employing the analytical procedure for calculating the nonadiabatic couplings. These trajectories have been successfully reproduced with the overlap method, and 30 more (for the total of 50 trajectories using the overlap method) trajectories have been carried out. Two additional sets of 13 trajectories each have been computed to test the effect of the initial conditions bias on the rotation directions. Our results show that the photoisomerization of azobenzene after excitation to the S1 state follows the rotational mechanism with the quantum yield of 0.65 and the lifetime of 67 fs. The inversion coordinates do not significantly influence the isomerization process. Two distinct rotational pathways, denoted as clockwise and counterclockwise, have been observed. The counterclockwise pathway is dominant in most cases (80
95%), while the clockwise pathway occurs only rarely. The reason behind this fact was found 11141

dx.doi.org/10.1021/jp2013094 |J. Phys. Chem. A 2011, 115, 11136–11143

The Journal of Physical Chemistry A

ARTICLE

Table 4. Percentage of Trajectories counter(S1) þ clockwise normal initial

13

counterclockwise

clockwise(S0)

87

0

conditions geometry biased

≈50

≈50

0

geometry and

100

0

0

14.5

71

14.5

’ AUTHOR INFORMATION Corresponding Authors

velocity biased Ootani et al.33

supporting figures for MD results (comparison of average adiabatic populations for the approximate and the analytical couplings, illustration of correlation of the motion of benzene rings, example trajectories for all pathways, and figures supporting the analysis of the inversion angles). This material is available free of charge via the Internet at http://pubs.acs.org.

*E-mail: [email protected] (J.P.), [email protected] (H.L.).

’ ACKNOWLEDGMENT This work has been supported by the Granting Agency of the Academy of Sciences of Czech Republic (Project No. IAA400400810) and Austrian Science Fund within the framework of the Special Research Program and F41 Vienna Computational Materials Laboratory (ViCoM). The calculations were performed in part at the Vienna Scientific Cluster (project nos. 70019 and 70151) of the Vienna University Computer Centre. ’ REFERENCES

Figure 10. CNNC dihedral angle as a function of time for the additional trajectories: (a) Trajectories with zero initial velocity; (b) Trajectories with initial velocity with positive projection into the clockwise direction.

out to be a slight barrier on the potential energy surface in the clockwise direction, mostly caused by sterical hindrance of benzene rings. Statistical analysis showed a nontrivial dependence of the isomerization direction on the initial conditions. The most important factors are the dihedral angles corresponding to rotations of the benzene rings and projection of the initial velocities to the direction of the isomerization, especially those of the nitrogen atoms. Based on these results, the two additional sets of 13 trajectories each have been constructed, which proved that the direction of the rotational motion can be fully controlled by selecting appropriate initial conditions.

’ ASSOCIATED CONTENT

bS

Supporting Information. Cartesian coordinates of the optimized structures, the table of vertical excitation energies, and

(1) Tamai, N.; Miyasaka, H. Chem. Rev. 2000, 100, 1875. (2) Rau, H. Photochromism: Molecules and Systems. Studies in Organic Chemistry; Elsevier: Amsterdam, 1990; Vol. 40, p 165. (3) Granucci, G.; Persico, M. Theor. Chem. Acc. 2007, 117, 1131. (4) Bockmann, M.; Doltsinis, N. L.; Marx, D. Angew. Chem. 2010, 122, 3454. (5) Ikeda, T.; Tsutsumi, O. Science 1995, 268, 1873. (6) Bach, H.; Anderle, K.; Fuhrmann, T.; Wendorff, J. H. J. Phys. Chem. 1996, 100, 4135. (7) Liu, Z.; Zhao, C.; Tang, M.; Cai, S. J. Phys. Chem. 1996, 100, 17337. (8) Liu, F. Z.; Hashimoto, K.; Fujishima, A. Nature 1990, 347, 658. (9) Rasmussen, P. H.; Ramanujam, P. S.; Hvilsted, S.; Berg, R. H. J. Am. Chem. Soc. 1999, 121, 4738. (10) Lim, H. S.; Han, J. T.; Kwak, D.; Jin, M.; Cho, K. J. Am. Chem. Soc. 2006, 128, 14458. (11) Hagen, R.; Bieringer, T. Adv. Mater. 2002, 13, 1805. (12) Natansohn, A.; Rochon, P. Chem. Rev. 2002, 102, 4139. (13) Shibaev, V.; Bobrovsky, A.; Boiko, N. Prog. Polym. Sci. 2003, 28, 729. (14) Credi, A. Aust. J. Chem. 2006, 59, 157. (15) Hugel, T.; Holland, N. B.; Cattani, A.; Moroder, L.; Seitz, M.; Gaub, H. E. Science 2002, 296, 1103. (16) Harada, A. Acc. Chem. Res. 2001, 34, 456. (17) Ballardini, R.; Balzani, V.; Credi, A.; Gandolfi, M. T.; Venturi, M. Acc. Chem. Res. 2001, 34, 445. (18) Tanner, D.; Wennerstrom, O. Tetrahedron Lett. 1981, 22, 2313. (19) Ciminelli, C.; Granucci, G.; Persico, M. J. Chem. Phys. 2005, 123, 174317. (20) Monti, S.; Orlandi, G.; Palmieri, P. Comp. Phys. 1982, 104, 1616. (21) N€agele, T.; Hoche, R.; Zinth, W.; Wachteitl, J. Chem. Phys. Lett. 1997, 272, 148. (22) Lednev, I. K.; Ye, T.-Q.; Matousek, P.; Towrie, M.; Foggi, P.; Neuwahl, F. V. R.; Umapathy, S.; Hester, R. E.; Moore, J. N. Chem. Phys. Lett. 1998, 290, 68. (23) Cattaneo, P.; Persico, M. Phys. Chem. Chem. Phys. 1999, 1, 4739. (24) Ishikawa, T.; Noro, T.; Shoda, T. J. Chem. Phys. 2001, 115, 7503. 11142

dx.doi.org/10.1021/jp2013094 |J. Phys. Chem. A 2011, 115, 11136–11143

The Journal of Physical Chemistry A (25) Conti, I.; Garavelli, M.; Orlandi, G. J. Am. Chem. Soc. 2008, 130, 5216. (26) Gagliardi, L.; Orlandi, G.; Bernardi, F.; Cembran, A.; Garavelli, M. Theor. Chem. Acc. 2004, 111, 363. (27) Granucci, G.; Persico, M.; Toniolo, A. J. Chem. Phys. 2001, 114, 10608. (28) Ciminelli, C.; Granucci, G.; Persico, M. Chem.—Eur. J. 2004, 10, 2341. (29) Toniolo, A.; Ciminelli, C.; Persico, M.; Martinez, T. J. J. Chem. Phys. 2005, 123, 234308. (30) B€ockmann, M.; Doltsinis, N. L.; Marx, D. Phys. Rev. E 2008, 78, 36101. (31) Shao, J.; Lei, Y.; Wen, Z.; Dou, Y.; Wang, Z. J. Chem. Phys. 2008, 129, 164111. (32) Sauer, P.; Allen, R. E. Chem. Phys. Lett. 2008, 450, 192. (33) Ootani, Y.; Satoh, K.; Nakayama, A.; Noro, T.; Taketsugu, T. J. Chem. Phys. 2009, 131, 194306. (34) Lischka, H.; Dallos, M.; Shepard, R. Mol. Phys. 2002, 100, 1647. (35) Shepard, R. Modern Electronic Structure Theory; World Scientific: Singapore, 1995; Vol. 1, p 345. (36) Shepard, R.; Lischka, H.; Szalay, P. G.; Kovar, T.; Ernzerhof, M. J. Chem. Phys. 1992, 96, 2085. (37) Lischka, H.; Dallos, M.; Szalay, P. G.; Yarkony, D. R.; Shepard, R. J. Chem. Phys. 2004, 120, 7322. (38) Dallos, M.; Lischka, H.; Shepard, R.; Yarkony, D. R.; Szalay, P. G. J. Chem. Phys. 2004, 120, 7330. (39) Lischka, H.; Shepard, R.; Shavitt, I.; Pitzer, R. M.; Dallos, M.; T. M€uller, Szalay, P. G.; Brown, F. B.; Ahlrichs, R.; H. J. B€ohm, Chang, A; Comeau, D. C.; Gdanitz, R.; Dachsel, H.; Ehrhardt, C.; Ernzerhof, M.; P. H€ochtl, Irle, S.; Kedziora, G.; Kovar, T.; Parasuk, V.; Pepper, M. J. M.; Scharf, P.; Schiffer, H.; Schindler, M.; M. Sch€uler, Seth, M.; Stahlberg, E. A.; Zhao, J.-G.; Yabushita, S.; Zhang, Z.; Barbatti, M.; Matsika, S.; Schuurmann, M.; Yarkony, D. R.; Brozell, S. R.; Beck, E. V.; Blaudeau, J.P.; Ruckenbauer, M.; Sellner, B.; Plasser, F.; Szymczak, J. J.; COLUMBUS, an ab-initio electronic structure program, release 5.9.2; 2008; www. univie.ac.at/columbus. (40) Lischka, H.; Shepard, R.; Pitzer, R. M.; Shavitt, I.; Dallos, M.; M€uller, T.; Szalay, P. G.; Seth, M.; Kedziora, G. S.; Yabushita, S.; Zhang, Z. Phys. Chem. Chem. Phys. 2001, 3, 664. (41) Lischka, H.; Shepard, R.; Brown, F. B.; Shavitt, I. Int. J. Quantum Chem. 1981, S 15, 91. (42) Tully, J. C. Faraday Discuss. 1998, 110, 407. (43) Swope, W. C.; Andersen, H. C.; Berens, P. H.; Wilson, K. R. J. Chem. Phys. 1982, 76, 637. (44) Butcher, J. J. Assoc. Comput. Mach. 1965, 12, 124. (45) Tully, J. C. J. Chem. Phys. 1990, 93, 1061. (46) Granucci, G.; Persico, M. J. Chem. Phys. 2007, 126, 134114. (47) Hammes-Schiffer, S.; Tully, J. C. J. Chem. Phys. 1994, 101, 4657. (48) Pittner, J.; Lischka, H.; Barbatti, M. Chem. Phys. 2009, 356, 147. (49) Barbatti, M.; Pittner, J.; Pederzoli, M.; Werner, U.; Mitric, R.; Bona
cic-Koutecky, V.; Lischka, H. Chem. Phys. 2010, 375, 26. (50) Barbatti, M.; Aquino, A. J. A.; Lischka, H. Phys. Chem. Chem. Phys. 2010, 12, 4959. (51) Barbatti, M.; Granucci, G.; Ruckenbauer, M.; Plasser, F.; Pittner, J.; Persico, M.; , Lischka, H. NEWTON-X, a package for Newtonian dynamics close to the crossing seam; 2008
2011; http:// www.univie.ac.at/newtonx. (52) Barbatti, M.; Granucci, G.; Persico, M.; Ruckenbauer, M.; Vazdar, M.; Eckert-Maksic, M.; Lischka, H. J. Photochem. Photobiol., A 2007, 190, 228. (53) StatSoft, Inc. STATISTICA, data analysis software system, version 7; 2004; www.statsoft.com. (54) Tsuji, T.; Takashima, H.; Takeuchi, H.; Egawa, T.; Konaka, S. J. Phys. Chem. A 2001, 105, 9347. (55) Bouwstra, J. A.; Schouten, A.; Kroon, J. Acta Crystallogr., Sect. C 1983, 39, 1121. (56) Diau, E. W.-G. J. Phys. Chem. A 2004, 108, 950.

ARTICLE

(57) Birnbaum, P. P.; Linford, J. H.; Style, D. W. G. Trans. Faraday Soc. 1953, 49, 735. (58) Beveridge, D. L.; Jaffe, H. H. J. Am. Chem. Soc. 1966, 88, 1948. (59) Sellner, B.; Ruckenbauer, M.; Stambolic, I.; Barbattti, M.; Aquino, A. J. A.; Lischka, H. J. Phys. Chem. A 2010, 114, 8778. (60) Amadei, A.; Linssen, A. B. M.; Berendsen, H. J. C. Protein: Struct., Funct., Genet. 1993, 17, 412. (61) Plasser, F.; Barbatti, M.; Aquino, A. J. A.; Lischka, H. J. Phys. Chem. A 2009, 113, 8490. (62) Toniolo, A.; Levine, B.; Thompson, A.; Quenneville, J.; M. Ben-Nun, Owens, J.; Olsen, S.; L.Manohar, Martinez, T. J. Computational Methods in Organic Photochemistry; Marcel-Dekker: New York, 2005; p 167. (63) Fujiwara, M.; Akiyama, M.; Hata, M.; Shiokawa, K.; Nomura, R. ACS Nano 2008, 2, 1671. (64) Satzger, H.; Sporlein, S.; Root, C.; Wachtveitl, J.; Zinth, W.; Gilch, P. Chem. Phys. Lett. 2003, 372, 216. (65) Rau, H.; L€uddecke, E. J. Am. Chem. Soc. 1982, 104, 1616. (66) Zimmerman, G.; Chow, L.-Y.; Paik, U.-J. J. Am. Chem. Soc. 1958, 80, 3528. (67) Gegiou, D.; Muszkat, K. A.; Fischer, E. J. Am. Chem. Soc. 1968, 90, 12. (68) Ronayette, J.; Arnaud, R.; Lebourgeois, P.; Lemaire, J. J. Am. Chem. Soc. 1974, 52, 1848. (69) Bortolus, P.; Monti, S. J. Phys. Chem. 1979, 89, 648. (70) Siampiringue, N.; Guyot, G.; Monti, S.; Bortolus, P. J. Photochem. 1987, 37, 185. (71) Fisher, R. A. Trans. R. Soc. Ed. 1918, 52, 399.

11143

dx.doi.org/10.1021/jp2013094 |J. Phys. Chem. A 2011, 115, 11136–11143