Nonadiabatic Hybrid Quantum and Molecular Mechanic Simulations of

J. Phys. Chem. A 2010, 114, 745–754

745

Nonadiabatic Hybrid Quantum and Molecular Mechanic Simulations of Azobenzene Photoswitching in Bulk Liquid Environment Marcus Bo¨ckmann, Nikos L. Doltsinis,*,† and Dominik Marx Lehrstuhl fu¨r Theoretische Chemie, Ruhr-UniVersita¨t Bochum, 44780 Bochum, Germany ReceiVed: October 22, 2009

A nonadiabatic hybrid quantum and molecular mechanical (na-QM/MM) molecular dynamics scheme has been implemented recently combining the nonadiabatic Car-Parrinello molecular dynamics method by Doltsinis and Marx [Phys. ReV. Lett. 2002, 88, 166402] with the QM/MM coupling approach by Laio et al. [J. Chem. Phys. 2002, 116, 6941]. Here an extensive validation of the underlying, density functional theory based, electronic structure methods by comparison to CASPT2 ab initio data is presented for the case of azobenzene. The “on the fly” na-QM/MM method is then applied to study ZfE and EfZ photoisomerization of azobenzene in a bulk liquid environment. The isomerization mechanism is found to be a pedal motion of the central CNdNC group in both cases. While the ZfE reaction is barely affected by the environment, EfZ photoisomerization is slowed down considerably in the liquid compared to the gas phase. This effect is due to the fact that reorientation of the phenyl rings is significantly hindered in the liquid by steric nearest neighbor interactions. Nonradiative decay is found to be substantially faster for Z-AB (subpicosecond regime) than for E-AB (picosecond regime). The main molecular motions responsible for nonadiabatic coupling have been identified as the oscillations in the NN and CN bond lengths, the CNN bond angles, and the CNNC dihedral angle. I. Introduction Azobenzene, H5C6sNdNsC6H5 (AB), is one of the most widely used photoswitches in physics, (bio)chemistry, and materials science1-19 owing to its property that it can be reversibly switched between the Z (also known as “cis”) and E (“trans”) isomers by exposing it to light of a suitable wavelength. The isomerization from Z to E is accompanied by a substantial expansion of length of the AB molecule (by ∼2.4 Å), a fact that has been exploited in optomechanical applications. For instance, AB has been integrated into synthetic peptides and foldamers in order to photocontrol their conformational dynamics.6,19,20 The same principle has underlied attempts to design molecular nanomachines and optomechanical cycles that convert light into mechanical action.15,16 In materials science, AB functional units are incorporated in polymers and liquid crystals to manipulate reversibly their shape or phase order by irradiation,11,21 properties which may, for instance, be exploited in photoaddressable image storage media.4 Furthermore, photochromic molecules such as AB can also be used as optoelectronic building blocks in molecular electronics due to their photoswitchable, isomer-dependent conductivity.22,23 Despite the widespread use of AB-based photoswitches in practical applications, the details of the underlying photoisomerization process had been poorly understood for decades. Only a recent surge of ab initio quantum chemical studies finally seemed to establish a consensus on the basic photoswitching mechanism of isolated AB molecules in the gas phase.24-31 Broadly speaking, it was believed that photoexcitation to the S1 (nfπ*) excited state results in out-of-plane torsional isomerization, while in-plane inversion governs isomerization in the * To whom correspondence should be addressed. E-mail: nikos.doltsinis@ kcl.ac.uk. † Present address: Department of Physics, King’s College London, Strand, London WC2R 2LS, United Kingdom.

S2 (πfπ*) state. In a recent first principles molecular dynamics study of the ZfE photoisomerization in the liquid phase it was proposed that the switching is governed by a small-amplitude pedal motion of the central CNNC moiety.32 This mechanism is largely unaffected by spatial restrictions of the phenyl rings, as encountered, for example, in the liquid phase. Similar scenarios have been put forward for rotation-restricted AB derivatives33,34 and related systems.35,36 While understanding ZfE photoisomerization of AB in the liquid and in the gas phase was certainly an important step, most applications such as the aforementioned ones occur in condensed phases consisting predominantly of E-AB isomers. This presents an additional challenge to first principles computational modeling due to the significantly longer isomerization times in the EfZ direction compared to ZfE. Computational studies of photoswitching in complex condensed matter are far beyond the capabilities of current ab inito methods due to their computational demand. A viable option is the use of hybrid quantum classical (“QM/MM”) two-scale techniques37-40 in which the photoactive chromophore is treated using a quantum mechanical (QM) electronic structure method whereas the photoinactive environment is approximated by classical force fields, i.e., molecular mechanics (MM). Here we present a powerful and versatile nonadiabatic QM/MM (na-QM/ MM) simulation method that involves not only the excited state as such but includes a proper nonadiabatic coupling to the ground state in a fully dynamical framework. In this approach, the electronic structure of the QM subsystem is treated within density functional theory41,42 and nonadiabatic effects are included in a mixed quantum-classical manner through Tully’s surface hopping algorithm.43 The approach can also be viewed as a QM/MM extension of the nonadiabatic ab inito molecular dynamics (AIMD) method41,42 by Doltsinis and Marx41,42,44,45 using the efficient CPMD/Gromos QM/MM interface46 (see the monograph40 for detailed and unifying discussions of the

10.1021/jp910103b  2010 American Chemical Society Published on Web 11/23/2009

746

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

Bo¨ckmann et al.

underlying techniques). We would like to point out that the procedure applied here crucially differs from that described by Groenhof et al.47 in that the latter does not involve the computation of nonadiabatic coupling elements. This na-QM/MM method is applied here to photoisomerization in bulk liquid AB based on our preliminary work,32 which is a prototype of photoswitchable condensed matter systems.8,13,48 The present study now gives important qualitative but detailed insights into the effects of molecular environments and confinement on the photoisomerization mechanism as well as rough trends that affect the quantum yield, which provides valuable clues to understand photoaddressable polymeric or liquid crystalline materials used in actual applications.2,4,6,8,9,11,13,15,16,19-21 In particular, transcending our short technical report32 of Z-AB, we focus here on the technologically relevant trans photoswitch, E-AB, in a bulk environment and find significant differences of the EfZ switching process compared to the previously reported ZfE scenario using refined order parameters. Finally, we provide comprehensive benchmarking of the quality of the restricted open-shell Kohn-Sham approach for AB photoisomerization in comparison to reference methods such as CASPT2. II. Computational Methods A. Na-QM/MM Molecular Dynamics. Inclusion of nonadiabatic effects into otherwise classical nuclear dynamics is achieved using a mixed quantum classical approach43 in which the nuclei follow a single classical trajectory, R(t), whereas the electrons in the QM subsystem are described by ΦQM/MM satisfying the time-dependent Schro¨dinger equation (TDSE). The QM/MM coupling is established via a Hamiltonian46 HQM/MM which is a function of all the nuclear coordinates, i.e., both of the QM and the MM subsystems. Likewise, the total wave function, ΦQM/MM, depends on the entire set of nuclear coordinates and is expanded

ΦQM/MM(r, R, t) )

∑ ai(t)φiQM/MM(r, R)

(1)

i

in terms of known electronic state functions, φQM/MM (r, R). The i time-dependent expansion coefficients, ai(t), are determined by inserting this ansatz into the TDSE, resulting in a system of coupled differential equations

i a˙i(t) ) - ai(t)EiQM/MM p

∑ aj(t)CijQM/MM

where EQM/MM is the energy of electronic state i, and i

〈

CijQM/MM ) φiQM/MM

||

(2)

j

〉 〈

| |

〉

∂ QM/MM ˙ φiQM/MM ∂ φjQM/MM φ )R ∂t j ∂R (3)

are the nonadiabatic couplings between states i and j. The present two-state implementation couples nonadiabatically the closed-shell Kohn-Sham ground state, φQM/MM , to the 0 reorthonormalized restricted open-shell Kohn-Sham (ROKS) representation49,50 of the S1 first singlet excited state, φ1QM/MM, following the successful single-scale na-QM technique.41,51-56 As a two-determinant representation,40,42,57 the ROKS S1 state provides an improved reference to compute nonadiabatic couplings,58 and yields reliable S1 nonradiative lifetimes and decay mechanisms when nonadiabatically coupled to the KS ground state41,51-56 the level of accuracy obtained by this

efficient approach for the system of interest to this study, i.e., AB photoisomerization, will demonstrated in detail in section IIIA. Along the classical trajectory, the probability Πij for a transition between S1 and S0 is determined by Tully’s fewest switches surface hopping algorithm43 as

Πij ) max(0, Pij)

(4)

with the hopping parameter

Pij ) -∆t

|| ||

d 2 a dt i 2

(5)

ai

where i, j ) 0, 1 and ∆t is the time step. In addition to this QM component, a force field parametrization suitable for condensed phase simulations48 must be used to define the MM part of HQM/MM. Importantly, the consistent electron density in the excited state, which deviates from that of the ground state, has to be used to include the QM T MM electrostatic coupling within HQM/MM. B. Technical Details. The na-QM/MM simulation method has been implemented in the ab initio molecular dynamics package CPMD40,57,59 extended by the CPMD/Gromos QM/MM interface46 in conjunction with a consistent AB force field for the MM part.48 The simulations were carried out in the Born-Oppenheimer propagation mode employing the PBE exchange and correlation functional and a plane wave basis set truncated at 70 Ry together with dual-space pseudopotentials for the QM part. QM atoms carry the same nonbonded parameters as the MM atoms and the cutoff for the explicit interaction of the remaining MM atoms with the QM charge density, i.e., the nearest-neighbor cutoff, between the charge group centers was set to 20.0 au (10.6 Å); the nearest-neighbor list was updated every 10 time steps. Charge groups at a distance greater than 20.0 au interact with the QM subsystem through a multipole expansion of the QM charge density. Liquid AB comprises 343 AB molecules in the all-Z or all-E configuration in a cubic 45.2 Å box subject to periodic boundary conditions. Therein a single AB molecule was treated quantummechanically in a 18 Å cubic box using cluster boundary conditions.57 Initial conditions for na-QM/MM simulations were sampled from a standard QM/MM ground state run at 400 K48 and extended, after vertical S0fS1 photoexcitation at t ) 0, in the nonadiabatic propagation mode. In addition, reference gas phase nonadiabatic QM runs of a single AB molecule have been performed analogously for the ZfE and EfZ photoisomerization using the same setup as for the QM part of the na-QM/ MM calculations together with a time step of 2 au. Ten initial conditions have been sampled for each type of setup, i.e., EfZ and ZfE photoswitching in the liquid and in the gas phase for reference, resulting into 40 independent nonadiabatic trajectories. The trajectory length was at least 500 fs but for the slower EfZ condensed phase process most of the trajectories are around 2 ps or longer and one has been extended up to 3.2 ps. In total, more than 40 ps of na-QM/MM simulations have been carried out “on the fly” including full electronic coupling according to Tully’s fewest switches surface hopping criterion on which the presented analysis is based. Clearly, this amount of statistics is far from being sufficient if the aim is to extract quantitative results such as excited-state lifetimes or quantum yields. Still, a wealth of qualitative insights into the photoinduced isomerization mechanisms and, in par-

Azobenzene Photoswitching in Bulk Liquid Environment

J. Phys. Chem. A, Vol. 114, No. 2, 2010 747

TABLE 1: Comparison of Ground State Structures of Isolated AB: Bond Lengths in Ångstroms and Angles in Degrees (r ) ∠CNN, O ) ∠CNNC, θ ) ∠CCNN) E-AB method 25

CAS(14,12)/6-31G* CAS(10,10)/6-31G63 CAS(10,8)/6-31G*65 CAS(6,6)/6-31G70 CAS(6,6)/4-31G27 CAS(6,5)24 MP2/cc-pVTZ62 MP2/6-31G+*71 DFT/BP86/TZVP62 DFT/LDA28 DFT/B3LYP/6-31G63 DFT/B3LYP/6-31G+*71 HF/3-21G72 AM173 AM1+V74 this work exp. (X-ray)68,75 exp. (electr. diffr)76 a

Z-AB

Erel

rNN

rCN

R

φ

θ

E

Erel

rNN

rCN

R

φ

0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0

1.243 1.260 1.243 1.259 1.253 1.24 1.268 1.279 1.267 1.278 1.278 1.258 1.240 1.231 1.239 1.271 1.259a 1.25

1.422 1.424 1.423 1.483 1.428 1.42 1.417 1.423 1.420 1.408 1.424 1.420 1.427 1.437

115.1 116.8 115.0 114.2 116.8 115.0 113.7 113.4 114.8

180.0 180.0 180.0 180.0 180.0 180.0 180.0 177.6 180.0 177.2 180.0 180.0 180.0 179.3 180.0 180.0 180.0 180.0

0.0 0.1 0.0 0.0 0.0 0.0 0.0 18.5 0.0 –8.8 0.2 0.0 0.0 15.3 0.0 0.0 0.0 0.0

0.71

1.242

1.435

122.9

4.2

62.0

0.87

1.242 1.240

1.437 1.472

115.0 121.9

3.2 0.0

62.8 60.8

1.24 1.261 1.274 1.255 1.270

1.43 1.432 1.440 1.437 1.421

123.0 120.8 120.4 124.1

0.0 7.3 6.4 11.4 13.2

56.0

1.249 1.235 1.204 1.221 1.262 1.253

1.437 1.448 1.442

124.1 124.8 129.2 124.3 121.1 121.9

9.6 0.0 2.3 4.1 9.9 8.0

46.2 51.0 45.8

1.426 1.428 1.43

115.6 115.3 116.3 119.6 117.5 115.0 114.1 114.1

0.79 0.59 0.69 0.90

0.61

1.449 1.449

52.2 52.9

58.0

A lower value of 1.247 Å given in ref 75 is related to disorder effects.

ticular, into the significant differences between EfZ vs ZfE photoswitching in the condensed phase can be gained. In order to illustrate our findings representative trajectories have been selected which are found to be typical based on the available ensemble of trajectories. The single point state-averaged CASSCF(14,12) calculations have been performed with the 6-31G* basis set employing the Gaussian03 and MOLCAS60 program packages. In line with Cembran et al.,25 the active space comprises lone pairs, π and π* of the NdN double bound, and the π system of the two phenyl rings omitting the lowest two bonding and the highest two antibonding orbitals. CASPT2 calculations on top of the CASSCF wave functions were done in MOLCAS with the lowest 14 orbitals kept frozen. III. Results and Discussion A. Potential Energy Surfaces and Stationary Points. Let us begin by examining the S0 and S1 potential energy surfaces (PESs) that underlie the nonadiabatic surface hopping simulations discussed below. The ground-state-optimized structures of E-AB and Z-AB obtained from our DFT calculations are listed in Table 1 together with data from literature. We do not wish to explicitly discuss every single entry in this table and in the subsequent two tables of this sort; however this offers a comprehensive overview of previous work and demonstrates in addition the spread and accuracy of results obtained for AB using different state-of-the-art ab initio methods such as CASSCF, CASPT2, CC2, or MRCI. Expectedly, the groundstate data are all in good general agreement, differences in the degree of nonplanarity being the only noteworthy issue. Taking the experimental data and highest accuracy ab initio calculations (i.e., CASSCF(14,12)/6-31G*25) as a reference; however, we conclude that E-AB is a planar, C2h symmetrical structure, while Z-AB is nonplanar with a small central CNNC dihedral angle of less than 10° and a large CCNN dihedral around 60° indicating a tilt in the phenyl rings. The Z-AB isomer is predicted to be 0.61 eV higher in energy than E-AB by our calculations, lying in between the CASSCF value of 0.71 eV and the CASPT2 value of 0.51 eV.25 Table 2 collects the vertical and adiabatic S1 excitation energies from this work and the literature. The gas-phase experimental values61 serve as a reference to calibrate the various

theoretical methods. Again, a comprehensive literature overview is given without discussing every work individually. To correct for the well-known49 underestimation of the vertical excitation gap by ROKS we have added a constant shift of 0.7 eV to all our ROKS energies, which we will maintain in all dynamics simulations discussed below. We note that the reported CASPT2 vertical excitation energies25 are consistently lower than the (corrected) ROKS and experimental numbers by about 0.3-0.4 eV, while the CASSCF(14,12) values25 are too high by the same amount. The vertical de-excitation, i.e., fluorescence, and adiabatic excitation energies can easily deviate from these rules as the underlying structures exhibit large differences as we will see in Table 3. It should also be noted that the popular TDDFT method underestimates the vertical excitation energies by 0.6-0.7 eV when used in conjunction with a GGA functional,62 and still by ca. 0.5 eV with the hybrid functional B3LYP,63 see Table 2. In particular, the DFT/ROKS data for vertical excitation from both E-AB and Z-AB are found to compare best to the available CC2 data and to gas phase experiment. The excited-state-optimized structures in Table 3 exhibit a widespread difference between the different methods, illustrating the challenge at hand. Although all calculations consistently predict a contraction of the CN bonds by about 0.05 Å, there is controversy about the degree of torsion in the CNNC group with dihedrals ranging from 45 to 180°. In particular those CASSCF calculations with larger active spaces and basis sets all predict planar S1 global minima of the E-AB type.25,64,65 This is in contrast to the ROKS calculations whose S1 global minimum structure is a rotamer with a twisted CNNC group, explaining why comparison of adiabatic excitation and fluorescence energies in Table 2 can be misleading. From the large differences in excited-state global minimum structures one is tempted to conclude that the underlying S1 potential energy surfaces (PESs) are vastly different. However, as we shall see in the following, this is not the case, but rather the large deviations arise from the pronounced flatness of the S1 potential in the vicinity of the global minimum. To assess the quality of the ROKS S1 PES, we have calculated ROKS, CASSCF, and CASPT2 energy profiles along the CNNC dihedral angle, which obviously is an important internal coordinate for photoisomerization in the S1 state. On the lefthand side of Figure 1 we present DFT, CASSCF, and CASPT2

748


Bo¨ckmann et al.

TABLE 2: Vertical and Adiabatic (De)Excitation Energies in eV of E-AB, Z-AB, and the Optimized S1 Rotamer R; Adiabatic Excitation Energies Are Relative to the E-AB Ground State method

ref

εvert(E)

εvert(Z)

this work MP2/cc-pVTZ//TDDFT MP2/cc-pVTZ//CC2 MP2/cc-pVTZ//CCSD CAS(14,12) CAS(14,12)//PT2 AM1*+Vcorr CAS(6,6) CAS(6,6)//PT2 CAS(6,6)//CIPSI-EN CAS(6,6)//CIPSI-MBP CAS(6,5) CAS(6,5)//MRSDCI CAS(6,5)//CIPSI CAS(6,6) CAS(6,6)//cc-pVDZ CAS(10,10) CAS(10,10)//PT2 B3LYP//TDDFT DFT//CDFT exp. (gas phase) exp. (ethanol) exp. (n-hexane)

calibrating blue-shift of 0.7 eV TDDFT, BP86/aug-TZVP; CC2, aug-cc-pVTZ; CCSD, [3s2p1d/s]62 TDDFT, BP86/aug-TZVP; CC2, aug-cc-pVTZ; CCSD, [3s2p1d/s]62 TDDFT, BP86/aug-TZVP; CC2, aug-cc-pVTZ; CCSD, [3s2p1d/s]62 6-31G*25,31 6-31G*25,31 74 6-31G(*)70 6-31G(*)70 6-31G(*)70 6-31G(*)70 small basis set (larger basis set)24 small basis set (larger basis set)24 small basis set (larger basis set)24 4-31G27 4-31G27 6-31G63 6-31G63 6-31G63 double ζ + polarization28 61 66 77

2.87 2.19 2.84 2.95 3.18 2.50 (2.53) 2.94 4.24 2.67 3.02 2.82 2.85 (3.48) 3.11 2.81 2.59 2.75 3.11 2.34 2.33 2.2 2.82 2.88 2.76

2.96 2.34 3.00 3.17 3.38 2.71 (2.72) 3.32

εfluor(R)

εadiab(R)

1.20

2.30

1.92 1.25 (1.27)

2.61 1.95 (2.01)

3.64 (3.16) 2.95 (2.76) 3.65 (4.53) 3.95 2.94

2.1 2.92 2.88

TABLE 3: Comparison of S1 Excited State Structures of the Isolated AB (Bond Lengths in Ångstroms and Angles in Degrees) rCN

∠CNN

CASSCF(14,12)/6-31G*25a CASSCF(14,12)/6-31G*64a CASSCF(10,8)/6-31G*65 CASSCF(6,6)/4–31G27b ROKS/BLYP34 CASSCF(6,6)/6-31G70 AM1*+Vcorr74

optimized 1.253 1.253 1.253 1.263 1.281 1.251 1.231

1.366 1.368 1.366 1.398 1.370 1.394

CASSCF(6,5)/small24 ROKS/PBE

1.26 1.280

1.38 1.368

128.8 128.5 128.7 127.9 126.0 126.3 137.0 126.8 126.0 124.7

method literature

this work

rNN

∠CNNC

∠CCNN

180.0 180.0 180.0 128.5 129.4 125.0 95.0

0.0 0.0 0.0 –5.7 –6.4 5.6

135.0 116.0

–3.0 –1.7

118.0 100–110 100

–9.0

estimated from MEP scan literature

a

CDFT/caPZ28 B3LYP//TD-DFT64 B3LYP//CASPT2/6-31G*64

C2h. b Optimization not converged.

S0 and ROKS, CASSCF, and CASPT2 S1 potential energy curves for the DFT ground-state minimum energy path (MEP) along the CNNC dihedral angle, θ. All ground-state curves have a maximum at 90°, the DFT curve being very close to the CASPT2 curve, while CASSCF is seen to overestimate by about 0.2 eV compared to CASPT2 over a wide range of θ. The agreement between the ROKS and CASPT2 excited-state curves throughout the entire range of the isomerization coordinate θ is remarkable, whereas the CASSCF S1 curve is very similar in shape but shifted upward by about 0.5 eV. Next, the right-hand side of Figure 1 shows ground- and excited-state energy profiles along the ROKS optimized S1 excited-state MEP along the CNNC dihedral. It is seen that the ROKS curve has a shallow S1 minimum around θ ≈ 120° which can be reached by a barrierless path from both the Z and E Franck-Condon (FC) points at θ ) 10° and 180°, respectively. As for the S0 MEP, there is again striking agreement between the ROKS and CASPT2 excited-state curves, while the CASSCF energies are higher than the CASPT2 data by 0.3-0.5 eV. With regards to the differences between the EfZ and ZfE photoisomerization dynamics which we shall discuss below, it is important to realize that there is a considerably larger potential

Figure 1. S0 and S1 energy profiles along the CNNC dihedral angle using structures optimized using DFT (PBE) for the S0 ground state (left) and ROKS (PBE) for the S1 first excited state (right). The ROKS energies (2) are compared to the state-averaged CAS(14,12) (4) and CASPT2 (]) data. S0 ground-state energies are shown for DFT (b), state-averaged CAS(14,12) (O), and CASPT2 (]). All energies are relative to the respective S0 energies of E-AB optimized with DFT in the S0 state; 0.7 eV have been added to the ROKS energies as explained in section IIIA.



Figure 3. (a) Visualization of plane normal vectors nR (cyan), nR′ (orange), nN (blue), and nN′ (red), together with reference coordinate system (x, purple; y, pink; z, green) and numbering. (b) Visualization of the projection onto the xz and yz planes for arbitrary vector (see text).

Figure 2. Analysis of the S0/S1 energy gap for a sample trajectory of the ZfE (left) and EfZ (right) photoisomerization in the gas phase (top) and the liquid (bottom). For selected structures obtained from the ROKS/DFT (filled black symbols) approach, state average CASSCF(14,12)/6-31G* energies have been calculated for S0 (open triangles) and S1 (open circles). The ROKS/DFT values for the S1 state are blue-shifted by 0.7 eV as explained in section IIIA, the CASSCF curves are red-shifted by 0.3 eV according to Figure 1.

energy difference between the FC point and the S1 global minimum upon vertical excitation of the Z isomer as compared to E-AB. At this stage it is concluded based on the direct comparison to CASPT2 reference data in Figure 1 that the aforementioned constant blue-shift applied to the ROKS data indeed corrects consistently the gap not only at 90° but along the entire MEPs along the CNNC dihedral angle, which is a major component of the reaction coordinate of AB photoisomerization, both in the S0 and in the S1 in the full range between 0° and 180°. Despite the convincing performance of the DFT/ROKS approach along the S0 and S1 MEPs shown in Figure 1 we have sought additional reassurance that the DFT gound state and the ROKS excited state are indeed reliable representations in all regions of configuration space visited during photoisomerization dynamics. For this reason we have calculated CASSCF S0 and S1 energies at regular intervals along typical examples of nonadiabatic EfZ and ZfE trajectories in the gas phase and also in the liquid. The comparison between these CASSCF and DFT/ROKS data shown in Figure 2 demonstrates that the overall shapes of the PESs obtained with the different methods are very similar. The fact that the CASSCF energies generally lie above the DFT/ROKS values is in line with the data presented in Figure 1. This it is not a shortcoming of the DFT/ROKS approach but rather must be seen as a deficiency of the CASSCF

method in view of the higher accuracy CASPT2 data which are virtually on top of the DFT/ROKS results according to the MEP scans in both ground and excited state, S0 and S1, depicted in Figure 1. In conclusion, the presented DFT/ROKS approach provides a useful representation of the potential energy landscape in the two relevant states including the (de)excitation gap. This is supported by both MEP scans along the isomerization angle and configurations sampled along dynamical trajectories undergoing photoisomerization. In particular, the quantitative agreement with CASPT2 data along the isomerization angle both in the ground and first excited state for θ from 0° to 180° suggests that the nuclear gradients as well as the electronic gap provided by DFT/ROKS for the subsequent nonadiabatic “on the fly” dynamics are effectively at about this level of accuracy along the main component of the reaction coordinate. B. photoisomerization Mechanism: Conformational Analysis. 1. Definition of Internal Coordinates. To analyze in detail the isomerization mechanism and possible differences between the gas and the liquid phase, we describe internal motion in terms of the plane normal vectors nR and nR′, of the two aromatic rings at their geometric centers, R and R′, together with the normal vectors nN and nN′ of the C(1)NN′ and NN′C(1′) coordination planes at N and N′ (see Figure 3). We measure torsion of the C(1)NN′ and NN′C(1′) coordination planes by the change in angle of the projection of the normal vector nN onto the yz plane (see Figure 3), ψR(t) ) ∠(nRyz(t), nRyz(0)), while ψRβ(t) ) ∠(nR(t), nβ(t)), R, β ∈ {N, N′, R, R′}, gives relative changes, e.g., ψNN′ is the C(1)NN′C(1′) dihedral angle and ψRN captures rotation of the phenyl rings. For this purpose we define an intrinsic (right-handed) coordinate system whose x axis is parallel to the NdN′ bond, the z axis is parallel to the arithmetic mean of nN and nN′, and the origin is at the geometric midpoint of the two nitrogen atoms. 2. ZfE photoisomerization. In this section, we analyze in detail the mechanism of Z-AB to E-AB photoisomerization. Condensed-phase effects will be discussed by comparing results for the liquid and for the gas phase. In Figure 4 we present the time evolution of ψNN′ (i.e., the CNNC dihedral) for typical ZfE trajectories in the liquid and in the gas phase together with the corresponding time evolution of ψN and ψN’. It is seen that for the liquid as for the gas phase ψNN′ changes rapidly, in about 30 fs, to a value of ∼90° upon photoexcitation at t ) 0; after the S1fS0 transition to the ground state, ψNN′ reaches a value of ∼180°, thus indicating a successful ZfE isomerization in both cases. Inspection of the order parameters ψN and ψN′ (cf.

750


Figure 4. Time evolution of ψNN′ (solid line), ψN (dashed line), and ψN′ (b) for typical ZfE trajectories in the liquid (top) and in the gas phase (bottom). Vertical lines indicate S1fS0 (dashed line) and S0fS1 (solid line) hops. Gray bars indicate the ψNN′ E reference range. The insets show the rmsds of N, N′ (solid) and R, R′ (dashed) relative to t ) 0.

Figure 5. Time evolution of ψRN (solid line), ψR (O), ψR′N′ (dashed line), and ψR′ (b) for typical ZfE trajectories in the liquid (top) and in the gas phase (bottom). Vertical lines indicate S1fS0 (dashed line) and S0fS1 (solid line) hops. Gray bars indicate the ψRN E reference range.

Figure 4) reveals that the total change in ψNN′ is due to equal contributions from the two coordination planes at N and N′ in opposite directions. Note that the analysis presented in Figure 4 shows no significant differences between the liquid and the gas phase. The liquid environment obviously does not impose any major constraints on the dynamics of the CNNC moiety. This is due to the fact that photoisomerization is dominated by a pedal motion of the CNNC group and not by large amplitude rotation of the phenyl rings. This is illustrated by the insets of Figure 4, which show the root-mean-square deviations (rmsd) of the N atoms and the phenyl ring centers R. We can see that during the first 30 fs it is mainly the translocation of the N atoms that is responsible for the change in ψNN′ by ∼90°, while the phenyl rings remain largely fixed in space. The average time is takes a molecule to reach ψNN′ ) 90° is 42 fs in the gas phase and 46 fs in the liquid, further underpinning the environment independence of this motion. However, the drastic impact of the condensed phase becomes apparent by the data shown in Figure 5. For the gas phase, a smooth reorientation of the phenyl rings from a ψR (ψR′) angle of about 10° to about 60° is observed shortly after the S1fS0 hop (t > 180 fs). Compared to the reorientation of the nitrogen coordination planes (cf. Figure 4), the relaxation of the phenyl rings is much slower. This implies that the phenyl rings follow

Bo¨ckmann et al.

Figure 6. Time evolution of ψNN′ (solid line), ψN (dashed line), and ψN′ (b) for typical EfZ trajectories in the liquid (top) and in the gas phase (bottom). Gray bars indicate the ψNN′ Z reference range. See caption of Figure 4.

up the reoriented CNNC moiety in order to reach their zero equilibrium angles ψRN and ψR′N′, indicating planarity, by rotation about the C-N bonds. Similarly, the large, rapid changes in ψRN and ψR′N′ immediately after excitation and deexcitation are not caused by a torsion of the rings (ψR and ψR′ are rather stable during these periods) but rather by the fast motion of the nitrogen coordination plane (cf. inset of Figure 4), which tend to orient parallel to the neighboring phenyl ring (i.e., ψRN ≈ 0°). For the liquid (Figure 5, top panel), in stark contrast, the reorientation that occurs after the S1fS0 transition is hindered to a large extent (open and filled circle lines); the rotation angles ψRN and ψR′N′ tend to stay far away from their equilibrium reference values (solid and dashed lines) over the entire time period shown (t e 500 fs). It is noted in passing that the ultrafast decay times extracted from our trajectories are in accordance with a time constant of ∼200 fs that is reported for the ZfE photoisomerization of AB in ethanol.66,67 3. EfZ photoisomerization. Analogous to the above discussion of the ZfE photoisomerization mechanism, we now analyze the EfZ photoisomerization in the gas phase and the liquid. For the gas phase, the time evolution of ψNN′ is diplayed in the bottom panel of Figure 6. It shows a more or less smooth decrease from 180° at t ) 0 to ∼90° at t ≈ 300 fs, where a S1fS0 transition occurs. After the hop ψNN′ rapidly falls to a value close to 0°, thus indicating a successful EfZ isomerization. As observed for the ZfE isomerization, the total change of ∼180° in ψNN′ is eventually accomplished by equal contributions (∼90°) from ψN and ψN′ (dashed and circled lines), in opposite directions. Interestingly, however, initially both nitrogen coordination planes rotate in the same direction, thus producing no net change of ψNN′ (e.g., at t ) 50 fs). This is a crucial difference from the ZfE photoisomerization mechanism. In the liquid (see upper panel of Figure 6), we observe a much slower decrease in ψNN′. For the trajectory shown a value of ∼120° is reached after 450 fs, still far from the 90° value where the S1fS0 hop took place in the gas phase. In fact none of the trajectories in the liquid reached 90°, and the average time to reach 120° is 672 fs, compared to 259 fs in the gas phase. In contrast to the gas phase, there is no sustained sign change in either ψN or ψN′, explaining the fact that changes in ψNN′ are smaller in the liquid. As for the ZfE isomerization, the change in ψNN′ is produced by the translocation of the N atoms, again indicating a pedal-motion-like mechanism (see insets of Figure 6). Interestingly, a similar translocation at fixed ψNN′ has been found in X-ray diffraction experiments of azobenzene crystals.68



Figure 9. Time evolution of the energy gap between ground and excited state, ∆E (black line), and the absolute CNNC dihedral angle relative to 90° (red line) for a typical ZfE trajectory.

Figure 7. Time evolution of ψRN (solid line), ψR (O), ψR′N′ (dashed line), and ψR′ (b) for typical EfZ trajectories in the liquid (top) and in the gas phase (bottom). Gray bars indicate the ψRN Z ground-state reference range. See caption of Figure 5.

Figure 8. Time at which a given CNNC angle is reached for the first time after photoexcitation of E-AB in the liquid (left) and in the gas phase (right). Values for individual trajectories are indicated by circles, ensemble averages are connected by a solid line, and the standard deviation is marked as error bars.

Inspection of the reorientation of the phenyl rings (see Figure 7) reveals for the gas phase (bottom panel) a slow but steady increase of ψR and ψR′ starting at t ≈ 200 fs and leading up to ∼20° at the S1fS0 hop. Considering, on the other hand, the rapid time evolution of ψRN and ψR′N′, this implies that this is chiefly due to the reorientation of its neighboring nitrogen coordination planes which each phenyl ring follows up. After the S1fS0 hop, ψRN and ψR′N′ are seen to smoothly approach their Z equilibrium value of 50° accompanied by an overall steady, but more rapidly fluctuating, increase of ψR and ψR′, before the system is pushed back to the excited state by a S0fS1 hop. For the liquid (Figure 7, top panel), we observe a qualitatively different behavior of ψR and ψR′. Showing a pronounced increase for t < 150 fs in opposite directions, ψR′ continuously increases (as in the gas phase) up to 30° at t ) 500 fs, whereas ψR is seen to fluctuate with a maximum value of ca. 25° at t ) 250 fs. Again, the rapid fluctuations observed for ψRN and ψR′N′ are chiefly due to the changes seen for ψN and ψN′ (see Figure 6), and in contrast to the ZfE direction it is a rare occurrence that ψRN and ψR′N′ are close to 0° simultaneously. The impact of the condensed phase for the EfZ isomerization can also be seen from Figure 8, where we have plotted the time needed to reach a specific CNNC value in the S1 excited state for all EfZ trajectories performed. While it takes roughly 230 fs on average to reach a CNNC angle of 140° both in the liquid and in the gas phase, it apparently takes increasingly more time on average to reach smaller CNNC dihedral angles in the liquid. The longer decay times as compared to our ZfE trajectories are in accordance with a larger time constant of ≈ 3 ps that has

been reported for the EfZ photoisomerization of AB in ethanol67 recently. Being the characteristic feature also of the EfZ photoisomerization, the pedal motion of the two N atoms accounts for the experimental finding of fast isomerization dynamics in rotationrestricted E-AB33 and clearly rules out an inversion type mechanism as deduced from resonance Raman intensity analysis.67 Note, that there is no large-amplitude rotation of the phenyl rings involved as has been suggested for the interpretation of fluorescence anisotropy data.69 C. Nonradiative Decay. In the context of Tully’s fewest switches surface hopping approach, nonradiative transitions occur with the hopping probability Πij (eq 4), which depends on fluctuations in both the S1 and S0 wave functions. As these wave functions are functions of the molecular configuration, their fluctuations are directly related to vibrational motion. It should therefore be possible to identify certain structural changes that have a strong influence on the time evolution of Πij. In the previous sections we have seen that during photoisomerization the most important geometrical changes are related to a translocation of the nitrogen atoms of the CNNC moiety and progress can be monitored by means of the CNNC dihedral angle. To this end we would like to demonstrate that the CNNC dihedral is indeed the main variable determining the S0-S1 energy gap, ∆E, as can be seen from Figure 9. Thus in the following discussion the two quantities should be seen as interchangeable; we will carry out the analysis in terms of the CNNC angle. Further candidates for important molecular motions for radiationless decay can be deduced from the main structural changes upon photoexcitation. This suggests to consider internal coordinates that are directly related to the nitrogen atoms, i.e., the CN and NN bond lengths, the CNN bond angle, as well as the related order parameters ψN, ψN′, ψRN, and ψR′N′. 1. Photoexcited Z-AB. A comparison of the time evolution of the hopping parameter Pij together with that of the timederivatives of various structural variables is presented in Figure 10 for gas-phase Z-AB. We observe that the time-derivative of the NN bond length, V(NN), shares many extrema and zero points with the envelope of Pij (top panel). Over long periods these are seen to coincide with the corresponding curves for one or both of the CN bond lengths (second panel). The third panel illustrates that some of the lower frequency modulations of Pij are well described by temporal changes in the CNN bond angles. A particularly good example is seen between 600 and 700 fs. The correlation between V(CNNC) and Pij (bottom panel) is less pronounced. However, the combination of a small energy gap (i.e., CNNC angle ≈ 90°) and a high CNNC velocity is seen to produce high amplitudes in Pij, e.g., around the first S1fS0 hop at 180 fs (see Figure 9 for the corresponding ∆E curve). It is important to realize that there is no single vibrational

752


Figure 10. Time evolution of transition parameter Pij (gray shaded background) and the time derivative of the NN and CN bond lengths, the CNN bond angle, and the CNNC dihedral (red and green curves) for a characteristic ZfE trajectory in the gas phase. Vertical lines indicate S1fS0 (dashed line) and S0fS1 (solid line) hops.

Figure 11. Time evolution of transition parameter Pij (gray shaded background), and the time derivative of the NN and CN bond lengths, the CNN bond angle, and the CNNC dihedral (red and green curves) for a characteristic ZfE trajectory in the liquid phase. Vertical lines indicate S1fS0 (dashed line) and S0fS1 (solid line) hops.

mode that can model the evolution of Pij. Rather a combination of many structural parameters is required, and we believe that the ones selected by us make the most dominant contributions. An indication of this is provided by the fact that whenever the extrema of all or many of the selected quantities coincide the corresponding amplitude of Pij is particularly large, as for instance at 180 or 640 fs. An analogous analysis for the liquid phase is shown in Figure 11. We observe a similar match between the time derivatives of the selected structural variables and the hopping parameter Pij as in the gas phase. Obvious examples are around 200, 300, and 650 fs where the extrema of all four molecular quantities coincide, producing a large amplitude in Pij. There are intervals in which the fluctuations in Pij appear slightly less regular then in the gas phase. This may be due to the liquid environment, but our limited statistics does not permit a definite conclusion on this matter. Figures 10 and 11 further show that there are generally a number of nonadiabatic transitions in both directions. Back transitions from the ground to the excited state are most likely to occur shortly after a S1fS0 hop, since the nonadiabatic transition probability dramatically increases due to the large amount of kinetic energy pumped into the coupling modes. Over

Bo¨ckmann et al.

Figure 12. Time evolution of transition parameter Pij, bond lengths and bond angles for characteristic ZfE trajectories in the liquid (left) and gas phase (right). E reference values are shown as horizontal lines (solid, average; dashed, ( standard deviation).

Figure 13. Time evolution of transition parameter Pij (gray shaded background) and the time derivative of the NN and CN bond lengths, the CNN bond angle, and the CNNC dihedral (red and green curves) for a characteristic EfZ trajectory in the gas phase. P01 has been scaled down by a factor of 10. Vertical lines indicate S1fS0 (dashed line) and S0fS1 (solid line) hops.

a longer period this energy dissipates into other, noncoupling modes and/or the environment (in the liquid). Figure 12 illustrates the large amplitude fluctuations of the NN and CN bond lengths and the CNN bond angle following a hop to the ground state. These are seen to be far from the equilibrium range determined from separate ground-state CP-MD simulations. 2. Photoexcited E-AB. We have applied the same procedure as on the previous section to analyze the nonadiabatic coupling modes in the case of photoexcited E-AB. For a characteristic gas phase trajectory, Figure 13 shows the comparison of the surface hopping parameter Pij with the selected structural quantities. As for Z-AB, the many extrema of the NN and CN bond lengths, the CNN bond angles, and the CNNC dihedral angle are seen to match the peaks in the envelope of Pij. The picture is very similar for the liquid phase, see Figure 14. Again, the precise shape of the Pij signal will be determined by a superposition of many coupling modes; we expect the ones selected here to make the dominant contributions. Beyond the primitive structural variables discussed so far in the context of nonadiabatic transitions, we have also constructed additional reaction coordinates based on the order parameters ψN, ψN′, ψRN, and ψR′N′. A comparison of the coupling modes ˙ RN - ψ ˙ R′N′| and ∆ψ ˙ N ) |ψ ˙N-ψ ˙ N′| with P10 is shown ˙ RN ) |ψ ∆ψ in Figure 15. The two modes are seen to be a good approximation to the envelope of P10 over wide ranges both in the gas


Figure 14. Time evolution of transition parameter Pij (gray shaded background) and the time derivative of the NN and CN bond lengths, the CNN bond angle, and the CNNC dihedral (red and green curves) for a characteristic EfZ trajectory in the liquid phase.

Figure 15. Comparison of the hopping criterion P10 (gray lines, left ˙ RN and ∆ψ ˙ N (see text) for EfZ axis) to the coupling modes ∆ψ trajectories in the liquid (top) and the gas phase (bottom).

˙ RN serves as an envelope of phase and in the liquid. While ∆ψ ˙N P10 for both EfZ and ZfE trajectories, the parameter ∆ψ appears to be more suitable for the EfZ than the ZfE case. The latter is better described by the coupling mode ∆V(CNN) ) |V(C(1)NN′) - V(NN′C(1′))|. IV. Conclusions We have investigated the photoisomerization of azobenzene in the bulk liquid phase at finite temperatures using a recently developed computationally efficient nonadiabatic surface hopping QM/MM simulation method. To be able to discern any effects of the condensed phase environment on the photoisomerization we have also carried out reference nonadiabatic simulations in the gas phase at the same temperature. An

J. Phys. Chem. A, Vol. 114, No. 2, 2010 753 extensive assessment of the quality of the underlying DFT/ ROKS potential energy surfaces of the ground/first-excited states was made by carrying out CASSCF and CASPT2 ab initio calculations along relevant minimum energy paths and photoisomerization trajectories. We found remarkable agreement between DFT/ROKS and CASPT2 for both ground and excited state thus convincingly validating our computationally efficient approach. In particular, our nonadiabatic “on the fly” simulations show that the key mechanistic feature of the azobenzene photoisomerization in the first excited state is a pedal-motion-like translocation of the innermost nitrogen atoms accounting for a net change in the central dihedral angle. We have found no hint at an inversion at the nitrogen atoms. ZfE photoisomerization was seen to be only mildly affected by the liquid environment. This is traced back to the fact that isomerization occurs by a small amplitude pedal motion of the central CNNC group and not by large amplitude rotation of the phenyl rings. Orientational relaxation of the rings after isomerization was found to be slowed down in the condensed phase. Much larger environmental effects were observed for EfZ photoisomerization, where the CNNC dihedral was seen to change much more slowly in the liquid compared to the gas phase. As a consequence no successful EfZ isomerization event occurred in the liquid during the simulation period. The difference between the dynamics of photoexcited E-AB and Z-AB mainly comes from significantly higher energy of the Z-AB FC region relative to E-AB. This is also the reason for the shorter S1 lifetime of Z-AB determined from an exponential fit to the S1 population decrease as a function of time. As a further difference, we point at the initial orientation of the phenyl rings being nearly coplanar for E-AB and roughly perpendicular for Z-AB. Here, we might expect a chemically modified E-AB with tilted phenyl rings to show faster dynamics. We have further studied the mechanism of nonradiative decay by identifying molecular vibrational modes which are strongly related to the time-dependent modulation of the surface hopping transition probabilities as a function of time. The most important structural parameters are the NN and CN bond lengths, the CNN bond angles, and the CNNC dihedral angle. For each trajectory numerous nonadiabatic transitions generally took place, as the coupling modes are far from their equilibrium immediately after a S1fS0 transition, thus producing high probabilities for a back transition. Having been validated in this detailed condensed phase study of photoswitching of both Z-azobenzene and E-azobenzene, we plan to apply our nonadiabatic QM/MM simulation method to model photoaddressable bulk materials such as AB-containing liquid crystals that are known to undergo light-induced ordered-disordered phase transitions. Acknowledgment. We are grateful to the Volkswagen Stiftung for supporting our project “Adaptive Multiscale Simulation: Connecting the Quantum to the Mesoscopic Level” within the framework of the program ”New Conceptual Approaches to Modeling and Simulation of Complex Systems Computer Simulation of Molecular and Cellular Biosystems as well as Complex Soft Matter”. M.C.B. would like to thank J. Ribas-Arino for helpful discussions of the MOLCAS calculations. The simulations were performed using resources from NIC Ju¨lich, Bovilab@RUB, and Rechnerverbund-NRW. References and Notes (1) Kumar, G. S.; Neckers, D. C. Chem. ReV. 1989, 89, 1915–1925. (2) Liu, Z. F.; Hashimoto, K.; Fujishima, A. Nature 1990, 347, 658.

754


(3) Sekkat, Z.; Dumont, M. Appl. Phys. B: Laser Opt. 1992, 54, 486. (4) Ikeda, T.; Tsutsumi, O. Science 1995, 268, 1873. (5) Hagen, R.; Bieringer, T. AdV. Mater. 2001, 13, 1805. (6) Spo¨rlein, S.; Carstens, H.; Satzger, H.; Renner, C.; Behrendt, R.; Moroder, L.; Tavan, P.; Zinth, W.; Wachtveitl, J. Proc. Nat. Acad. Sci. 2002, 99, 7998. (7) Wachtveitl, J.; Spo¨rlein, S.; Satzger, H.; Fonrobert, B.; Renner, C.; Behrendt, R.; Oesterhelt, D.; Moroder, L.; Zinth, W. Biophys. J. 2004, 86, 2350. (8) Sekkat, Z.; Knoll, W. PhotoreactiVe organic thin films; Academic Press: San Diego, 2002. (9) Natansohn, A.; Rochon, P. Chem. ReV. 2002, 102, 4139–4175. (10) Shibaev, V.; Bobrovsky, A.; Boiko, N. Prog. Polym. Sci. 2003, 28, 729–836. (11) Yu, Y.; Nakano, M.; Ikeda, T. Nature 2003, 425, 145. (12) Banerjee, I.; Yu, L.; Matsui, H. J. Am. Chem. Soc. 2003, 125, 9542. (13) Feringa, B. L. Molecular Switches; Wiley-VCH: Weinheim, 2001. (14) Stolow, A. Annu. ReV. Phys. Chem. 2003, 54, 89. (15) Browne, W. R.; Feringa, B. L. Nat. Nanotechnol. 2006, 1, 25. (16) Hugel, T.; Holland, N. B.; Cattani, A.; Moroder, L.; Seitz, M.; Gaub, H. E. Science 2002, 296, 1103. (17) Yesodha, S. K.; Pillai, C. K. S.; Tsutsumi, N. Prog. Polym. Sci. 2004, 29, 45. (18) Yager, K. G.; Barrett, C. J. J. Photochem. Photobiol. A 2006, 182, 250. (19) Renner, C.; Moroder, L. ChemBioChem 2006, 7, 869. (20) Khan, A.; Kaiser, C.; Hecht, S. Angew. Chem., Int. Ed. 2006, 45, 1878. (21) Finkelmann, H.; Nishikawa, E.; Pereira, G. G.; Warner, M. Phys. ReV. Lett. 2001, 87, 015501. (22) Li, J.; Speyer, G.; Sankey, O. F. Phys. ReV. Lett. 2004, 93, 248302. (23) Valle, M. D.; Gutierrez, R.; Tejedor, C.; Cuniberti, G. Nat. Nanotechnol. 2007, 2, 176. (24) Ishikawa, T.; Noro, T.; Shoda, T. J. Chem. Phys. 2001, 115, 7503. (25) Cembran, A.; Bernardi, F.; Garavelli, M.; Gagliardi, L.; Orlandi, G. J. Am. Chem. Soc. 2004, 126, 3234. (26) Toniolo, A.; Ciminelli, C.; Persico, M.; Martinez, T. J. Chem. Phys. 2005, 123, 234308. (27) Diau, W.-G. J. Phys. Chem. A 2004, 108, 950. (28) Tiago, M. L.; Ismail-Beigi, S.; Louie, S. G. J. Chem. Phys. 2005, 122, 094311. (29) Crecca, C. R.; Roitberg, A. E. J. Phys. Chem. A 2006, 110, 8188– 8203. (30) Granucci, G.; Persico, M. Theor. Chem. Acc. 2007, 117, 1131. (31) Conti, I.; Garavelli, M.; Orlandi, G. J. Am. Chem. Soc. 2008, 130, 5216. (32) Bo¨ckmann, M.; Doltsinis, N. L.; Marx, D. Phys. ReV. E 2008, 78, 036101. (33) Pancur, T.; Renth, F.; Temps, F.; Harbaum, B.; Kru¨ger, A.; Herges, R.; Na¨ther, C. Phys. Chem. Chem. Phys. 2005, 7, 1985–1989. (34) Nonnenberg, C.; Gaub, H.; Frank, I. ChemPhysChem 2006, 7, 1455. (35) Liu, R. S. H.; Hammond, G. S. Proc. Nat. Acad. Sci. 2000, 97, 11153–11158. (36) Andresen, M.; Wahl, M. C.; Stiel, A. C.; Grater, F.; Schafer, L. V.; Trowitzsch, S.; Weber, G.; Eggeling, C.; Grubmuller, H.; Hell, S. W.; Jakobs, S. Proc. Nat. Acad. Sci. 2005, 102, 13070–13074. (37) Sherwood, P. Hybrid Quantum Mechanics/Molecular Mechanics Approaches. In Modern Methods and Algorithms of Quantum Chemistry; Grotendorst, J., Ed.; NIC: Ju¨lich, 2000. (38) Senn, H. M.; Thiel, W. Top. Curr. Chem. 2007, 268, 173–290. (39) Senn, H. M.; Thiel, W. Angew. Chem., Int. Ed. 2009, 48, 1198– 1229. (40) Marx, D.; Hutter, J. Ab Initio Molecular Dynamics: Basic Theory and AdVanced Methods; Cambridge University Press: Cambridge, 2009. (41) Doltsinis, N. L.; Marx, D. Phys. ReV. Lett. 2002, 88, 166402. (42) Doltsinis, N. L.; Marx, D. J. Theor. Comp. Chem. 2002, 1, 319– 349. (43) Tully, J. C. J. Chem. Phys. 1990, 93, 1061. (44) Doltsinis, N. L. Nonadiabatic dynamics: mean-field and surface hopping. In Quantum Simulations of Complex Many-Body Systems: From Theory to Algorithms; Grotendorst, J., Marx, D., Muramatsu, A., Eds.; NIC: FZ Ju¨lich, 2002; www.fz-juelich.de/nic-series/volume10/doltsinis.pdf.

Bo¨ckmann et al. (45) Doltsinis, N. L. Molecular dynamics beyond the Born-Oppenheimer approximation: mixed quantum-classical approaches. In Computational Nanoscience: Do it Yourself!; Grotendorst, J., Bl¨ugel, S., Marx, D., Eds.; NIC: FZ Ju¨lich, 2006; www.fz-juelich.de/nic-series/volume31/doltsinis1.pdf. (46) Laio, A.; VandeVondele, J.; Rothlisberger, U. J. Chem. Phys. 2002, 116, 6941. (47) Groenhof, G.; Scha¨fer, L. V.; Boggio-Pasqua, M.; Goette, M.; Grubmu¨ller, H.; Robb, M. A. J. Am. Chem. Soc. 2007, 129, 6812. (48) Bo¨ckmann, M.; Peter, C.; Delle Site, L.; Doltsinis, N. L.; Kremer, K.; Marx, D. J. Chem. Theory Comput. 2007, 3, 1789. (49) Frank, I.; Hutter, J.; Marx, D.; Parrinello, M. J. Chem. Phys. 1998, 108, 4060. (50) Grimm, S.; Nonnenberg, C.; Frank, I. J. Chem. Phys. 2003, 119, 11574. Grimm, S.; Nonnenberg, C.; Frank, I. J. Chem. Phys. 2003, 119, 11585. (51) Doltsinis, N. L. Mol. Phys. 2004, 102, 499. (52) Langer, H.; Doltsinis, N. L. Phys. Chem. Chem. Phys. 2004, 6, 2742. (53) Langer, H.; Doltsinis, N. L.; Marx, D. ChemPhysChem 2005, 6, 1734. (54) Markwick, P. R. L.; Doltsinis, N. L. J. Chem. Phys. 2007, 126, 175102. (55) Doltsinis, N. L.; Markwick, P. R. L.; Nieber, H.; Langer, H. In Radiation Induced Molecular Phenomena in Nucleic Acid; Shukla, M. K., Leszczynski, J., Eds.; Springer: Netherlands, 2008. (56) Nieber, H.; Doltsinis, N. L. Chem. Phys. 2008, 347, 405. (57) Marx, D.; Hutter, J. Ab Initio Molecular Dynamics: Theory and Implementation. In Modern Methods and Algorithms of Quantum Chemistry; Grotendorst, J., Ed.; NIC: Ju¨lich, 2000; www.theochem.rub.de/go/cprev. html, (accessed June 2007). (58) Billeter, S. R.; Egli, D. J. Chem. Phys. 2006, 125, 224103. (59) Hutter, J. Car-Parrinello Molecular Dynamics: An Ab Initio Electronic Structure and Molecular Dynamics Program, see www.cpmd.org, (accessed June 2007). ˜ .; Roos, B. O.; Ryde, (60) Karlstro¨m, G.; Lindh, R.; Malmqvis, P.-A U.; Veryazov, V.; Widmark, P.-O.; Cossi, M.; Schimmelpfennig, B.; Neogrady, P.; Seijo, L. Comput. Mater. Sci. 2003, 28, 222. (61) Andersson, J.-Å.; Petterson, R.; Tegne´r, L. J. Photochem. 1982, 20, 17–32. (62) Fliegl, H.; Ko¨hn, A.; Ha¨ttig, C.; Ahlrichs, R. J. Am. Chem. Soc. 2003, 125, 9821. (63) Schultz, T.; Quenneville, J.; Levine, B.; Toniolo, A.; Martı´nez, T. J.; Lochbrunner, S.; Schmitt, M.; Shaffer, J. P.; Zgierski, M. Z.; Stolow, A. J. Am. Chem. Soc. 2003, 125, 8098. (64) Gagliardi, L.; Orlandi, G.; Bernardi, F.; Cembran, A.; Garavelli, M. Theor. Chem. Acc. 2004, 111, 363. (65) Wang, L.; Xu, W.; Yi, C.; Wang, X. J. Mol. Graph. Model. 2009, 27, 792–796. (66) Na¨gele, T.; Hoche, R.; Zinth, W.; Wachtveitl, J. Chem. Phys. Lett. 1997, 272, 489. (67) Stuart, C. M.; Frontiera, R. R.; Mathies, R. A. J. Phys. Chem. A 2007, 111, 12072–12080. (68) Harada, J.; Ogawa, K.; Tomoda, S. Acta Cryst. B 1997, 53, 662– 672. (69) Chang, C.-W.; Lu, Y.-C.; Wang, T.-T.; Diau, E. W.-G. J. Am. Chem. Soc. 2004, 126, 10109–10118. (70) Cattaneo, P.; Persico, M. Phys. Chem. Chem. Phys. 1999, 1, 4739. (71) Kurita, N.; Ikegami, T.; Ishikawa, Y. Chem. Phys. Lett. 2002, 360, 349–354. (72) Cimiraglia, R.; Hofmann, H.-J. Chem. Phys. Lett. 1994, 217, 430– 435. (73) Kurita, N.; Tanaka, S.; Itoh, S. J. Phys. Chem. A 2000, 104, 8114– 8120. (74) Ciminelli, C.; Granucci, G.; Persico, M. Chem.-Eur. J. 2004, 10, 2327. (75) Bouwstra, J. A.; Schouten, A.; Kroon, J. Acta Crystallogr., Sect. C 1983, 39, 1121. (76) Traetterberg, M.; Hilmo, I.; Hagen, K. J. Mol. Struct. 1977, 39, 231. (77) Fujino, T.; Arzhantsev, S. Y.; Tahara, T. J. Phys. Chem. A 2001, 105, 8123.

JP910103B

Nonadiabatic Hybrid Quantum and Molecular Mechanic Simulations of

Recommend Documents