Extended Hückel

Nov 3, 2016 - The self-consistent coupling between quantum and classical degrees of freedom is achieved by nonadiabatic Hellmann–Feynman–Pulay ...
1 downloads 0 Views 5MB Size
Subscriber access provided by University of Otago Library

Article

A Nonadiabatic Excited State Molecular Mechanics/Extended Hückel Ehrenfest Method Robson da Silva Oliboni, Graziele Bortolini, Alberto Torres, and Luis G. C. Rego J. Phys. Chem. C, Just Accepted Manuscript • DOI: 10.1021/acs.jpcc.6b09606 • Publication Date (Web): 03 Nov 2016 Downloaded from http://pubs.acs.org on November 8, 2016

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

The Journal of Physical Chemistry C is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 51

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

A Nonadiabatic Excited State Molecular Mechanics/Extended Hückel Ehrenfest Method Robson da Silva Oliboni, Graziele Bortolini, Alberto Torres, and Luis G. C. Rego∗ Department of Physics, Universidade Federal de Santa Catarina, SC, CEP 88040-900, Brazil E-mail: [email protected]



To whom correspondence should be addressed

1

ACS Paragon Plus Environment

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

November 2, 2016 Abstract In this work we study intramolecular vibronic relaxation effects and photoisomerization dynamics in atomistic molecular systems by means of a mixed self-consistent quantum-classical Ehrenfest formalism for nonadiabatic molecular dynamics (NA-MD) in the excited state. The quantum mechanical part of the method is based on the extended Hückel formalism whereas the nuclei are treated with the molecular mechanics method. The self-consistent coupling between quantum and classical degrees of freedom is achieved by nonadiabatic Hellmann-Feynman-Pulay forces that conserve the total (quantum-classical) energy. Moreover, this work demonstrates the capabilities of a method that combines two simple though efficient computational schemes to describe complex excited state effects such as photoisomerization dynamics of stilbene and azobenzene molecules. In particular, for stilbene we point out the key importance of hydrogen out-of-plane (HOOP) modes to trigger the isomerization of the phenyl rings. For azobenzene, we identify the inversion-assisted rotation and the associated pedal-cycle motion of the CNNC dihedral as the driving isomerization mechanism. The method can also be applied for NA-MD studies of molecules adsorbed on extended solid surfaces and charge transfer processes.

2

ACS Paragon Plus Environment

Page 2 of 51

Page 3 of 51

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Introduction The understanding of the mechanisms that drive structural changes in the photoexcited state of molecular systems is essential for the development and design of new functional materials, with applications on energy conversion, optoelectronic devices like switches and sensors, biophysics, etc. 1–5 That is so because, on the atomic scale, structural dynamics strongly influences function and electronic behavior. There is, thereby, a compelling need for new theoretical methodologies for describing excited state electron-nuclear dynamics, specially for large scale molecular systems. Such computational approaches also provide theoretical support to the experimental counterpart, which consist of the ever increasing number of ultrafast time-resolved techniques 6–9 available for the study of charge and energy transfer processes that occur in supramolecular structures, 10 clusters, 11,12 quantum dots, 13–15 and various heterostructures. Conventional adiabatic molecular dynamics methods, based on the Born-Oppenheimer approximation, 16 cannot fully describe certain quantum mechanical effects, like vibrational relaxation, photoisomerization, phonon assisted tunneling and nonadiabatic effects in general. 17 For capturing such phenomena nonadiabatic molecular dynamics formalisms (NAMD) are necessary, because they allow for electronic quantum transitions between different potential energy surfaces (PES). Particularly, in the field of photochemistry the nonadiabatic approaches are paramount, since the structural motions occur on an excited potential energy surface and electronic transitions to nearby potential energy surfaces can happen during reaction time. 18,19 The natural way to deal with nonadiabatic transitions is a full quantum mechanical approach, 5,20 however, the high computational cost hinders the use of such formalisms for real interest systems. Alternatives for small systems include approaches such as the ab-initio multiple spawning (AIMS) 21–24 and full multiple spawning (FMS) techniques 21,25,26 based on Heller’s semiclassical approach, 27 and the multiconfigurational time-dependent Hartree (MCTDH) method. 28–30 Most common NA-MD methods separate the electron and nuclear 3

ACS Paragon Plus Environment

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

dynamics propagation by treating the nuclei via classical mechanics in mixed quantum– classical approaches. 18,19 In these methods the key issue is the self-consistency between the classical and quantal degrees of freedom. The nuclear dynamics can induce nonadiabatic transitions between electronic levels, and electronic transitions can alter the potential energy surface. Most used mixed quantum–classic NA-MD methods are the surface hopping 17,31,32 and Ehrenfest 18,33 methods. The two methods have much in common but differ in how the nonadiabatic transitions are taken into account in order to define the classical path. In surface hopping, the nuclei propagates on a single adiabatic PES but hops between surfaces occur with a certain rate. The most used rate criteria is given by the fewest switching surface hopping (FSSH) algorithm, 17 but recently several other surface hopping algorithms have been proposed and implemented. 34–39 In the Ehrenfest mean-field approach the potential energy function that drives the classical subsystem (nuclear dynamics) is obtained as an average over adiabatic potential energy surfaces, weighted by the coefficients of the quantum mechanical (electronic dynamics) wavepacket. 40 The general advantages and limitations of the two methods are well documented. 18,19 Both provide accurate quantum transition probabilities and conserve total (quantum plus classical) energy. The Ehrenfest method, however, does not obey microscopic reversibility 18,41 and cannot describe small polaronic effects, 42 among other limitations. 19 Surface hopping methods, on the other hand, are not invariant to the choice of quantum representation (adiabatic or diabatic) but they satisfy microscopic reversibility depending on the hopping algorithm and classical-quantum correlation is included in surface hopping methodologies. 18,19 Muller and Stock have compared self-consistent meanfield (Ehrenfest) with surface-hopping (SH) descriptions. 43 Despite the advantages of SH, though, they concluded that both methods are capable of describing nonradiative relaxation in molecules. In principle, the method used to calculate the instantaneous potential energy surfaces determines the computational cost of the NA-MD method. Several ab-initio electronic struc-

4

ACS Paragon Plus Environment

Page 4 of 51

Page 5 of 51

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

ture methods are used for surface hopping and Ehrenfest applications, 23,44–47 but the large number of degrees of freedom hinders the use of first principle methods, particularly if there is need to describe excited state dynamics of condensed systems. Therefore, demand for computationally feasible NA-MD methods for systems of practical interest has promoted the use of semiempirical methods for mixed quantum–classical dynamics. 48,49 This work presents a novel quantum-classical Ehrenfest approach for nonadiabatic excitedstate molecular dynamics. The quantum mechanical part of the method is based on the extended Hückel formalism whereas the nuclei are treated via molecular mechanics (MM). Moreover, we propose a formalism to extend molecular mechanics calculations to the realm of excited states. The self-consistent coupling between quantum and classical degrees of freedom is calculated on-the-fly by nonadiabatic Hellmann- Feynman-Pulay forces that conserve the total (quantum-classical) energy. The focus of this study is the excited-state dynamics relaxation effects and the underlying mechanisms for responsible isomerization. We demonstrate the capabilities of the method by studying the excited-state intramolecular vibrational relaxation (IVR) in benzene molecule as well as the photoisomerization dynamics of stilbene and azobenzene molecules (shown in Figure 2). The method is implemented in the package DynEMol. 50 The method can also be applied for NA-MD studies of molecules adsorbed on extended solid surfaces and charge transfer processes. The paper is organized as follows. Section II presents theoretical and computational details of the method. Section III presents results and discussions of nonadiabatic Ehrenfest simulations of vibronic relaxation effects and photoisomerization dynamics for prototypical systems. Section IV presents conclusions and prospects for future applications of the method.

5

ACS Paragon Plus Environment

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 6 of 51

The Mixed Molecular Mechanics/Quantum Mechanics Nonadiabatic Ehrenfest Method The method combines semiempirical – within the formalism of the extended Hückel theory – quantum dynamics calculations for the electronic degrees of freedom with molecular mechanics calculations for the nuclei. The procedure for performing the quantum dynamics calculations of electronic excitations has been used in studies for several molecular systems – gas phase, 51 molecules in solution, 52,53 dye sensitized semiconductor interfaces 54–56 – and it is described in detail elsewhere. 57,58 Herein we just point out that the quantum mechanP ical wavepackets that describe such excitations, |Ψi = φ Cφ (t)|φ(t)i, are written in terms of the adiabatic molecular orbitals {φ} and they evolve according to the time-dependent

Schrödinger equation (TDSE)

i~

∂ ˆ el (Rt ) |Ψ(r; R, t)i , |Ψ(r; R, t)i = H ∂t

(1)

where r designates the electronic coordinates, Rt ≡ R(t) are the time-dependent nuclear ˆ el (Rt ) is the time-dependent extended Hückel hamiltonian operator for an coordinates and H instantaneous molecular configuration. The TDSE is thus converted to an equation for the adiabatic coefficients Cφ (t), X i d d Cφ + Cγ hφ| γi + Cφ Eφ = 0 , dt dt ~ γ

(2)

that is solved sequentially for each time slice of dynamics. As for the nuclei, all of them in the system are treated with the molecular mechanics method, differently from the usual partitioning of the system between quantum and classical regions. Apart from that our implementation of MM follows standard procedures. We use a reliable force field to drive the nuclear dynamics; for the sake of completeness we describe

6

ACS Paragon Plus Environment

Page 7 of 51

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

below the OPLS force field (FF) 59,60 used in our calculations

MM VGS ({R}) =

X

X

Kb (R − Ro )2 +

bonds

+

X i,j6=i

4εij

"

angles !12

σij Rij

Kθ (θ − θo )2 +

X X

torsions



σij Rij

!6 #

+

X i,j6=i

Cn (cosφ)n

n

qj qi , 4πǫo Rij

(3)

where R is the atomic position and Rij = |Ri − Rj | is the distance between atoms i and j, which have partial charges qi and qj . Likewise, θ and φ are the angular and dihedral torsional variables. The parameters Ro and θo are the equilibrium bond length and angle; Kb , Kθ and Cn are the intramolecular FF parameters, and εij and σij are the Lennard-Jones parameters. The charges in Eq. (3) are kept fixed, since fluctuating electrostatic terms give rise to poor overall (quantum-classical) energy conservation. The intramolecular forces are calculated explicitly. However, it is well known that semiempirical hamiltonian operators, despite being capable of providing very good description of electronic energy levels, are not as good for describing interatomic potentials. On the other hand, MM force fields are built mainly to describe the interatomic classical potentials in a very efficient manner. Moreover, the standard force fields are parametrized for the electronic ground state (GS) and, most importantly, do not take into account the electronic back-reaction as the electron-hole (e-h) excitations evolve according to the TDSE (Eq. 1). To treat the NA dynamics we, therefore, extend the groundstate molecular mechanics to the excited-state and take into account the electronic backreaction. That is done by using nonadiabatic Hellmann-Feynman-Pulay forces in addition to those derived from the ground-state force field potential of Eq. (3). In the following, we describe the proposed self-consistent coupling between the electronic and nuclear degrees of freedom. We start by assuming that the electron-hole excitation consists of a Hartree product comprised of a wavepacket that describes the electronic excitation, Ψel (t), and another that describes the hole, Ψhl (t). A particular HOMO-LUMO excitation is depicted in Figure 1.

7

ACS Paragon Plus Environment

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 ACS Paragon Plus Environment

Page 8 of 51

Page 9 of 51

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

the unoccupied manifold, a perturbation potential is produced that is a functional of the electron and hole wavefunctions. This second term is responsible for the coupling between quantum and classical degrees of freedom in the excited state. It is calculated on-the-fly from the extended Hückel hamiltonian 61 as     X VEH Ψel (R, t), Ψhl (R, t) = T r ρEH (R, t)H(Rt ) = Eφ ρEH φφ ,

(7)

φ

el hl where ρEH ψφ = ρψφ − ρψφ is the electron-hole density matrix, written in terms of the adiabatic

states ψ and φ, with

|Ψel ihΨel | =

X

Bφ Bψ∗ |φihψ| = ρel ψφ |φihψ|

(8)

Aφ A∗ψ |φihψ| = ρhl ψφ |φihψ| .

(9)

φψ

hl

hl

|Ψ ihΨ | =

X φψ

The generalized Hellmann-Feynman-Pulay (HFP) force produced on atom N by an e-h excitation is given by   FN = −∇N VEH Ψel (R, t), Ψhl (R, t) ,

(10)

where ∇N ≡ ∇RN . This force is responsible for the electronic back-reaction on the classical degrees of freedom and, moreover, it gives rise to adiabatic and nonadiabatic effects. It is possible to derive an expression in closed form for the HFP force (the details of the derivation are provided as Supporting Information). For that matter, let us consider a wavepacket that satisfies Eq. (2), written as a time-dependent linear combination of adiabatic orbitals. In ˆ this case, E = hΨ|H({R})|Ψi is the energy of the wavepacket – not the variational energy eigenstate – and we have X X dE dEφ = |Cφ |2 + Eφ dRN dRN φ φ



dCφ∗ dCφ Cφ + Cφ∗ dRN dRN



,

(11)

where the first term on the r.h.s. is the adiabatic component of the force whereas the second 9

ACS Paragon Plus Environment

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 10 of 51

term yields the nonadiabatic components. To derive the adiabatic component we make explicit use of the extended Hückel hamiltonian, given by Hαβ = ξαβ Sαβ , with ξαβ ≡ Kαβ (hα + hβ )/2, whose eigenstates are obtained from HQφ = Eφ SQφ , where Sαβ = hα|βi are the overlap matrix elements. In the foregoing definitions and for the expressions thereafter we use the notation fα ≡ fa (r−RA ) to designate a Slater-type atomic orbital fa centered at position RA of atom A, thereby {α} ≡ {a, A}, and likewise for fβ ≡ fb (r − RB ). Following the standard procedures (described in Supporting Information) we obtain X φ

|Cφ |2

X X X ˆ dEφ ∂hφ|H|φi ∂Sαβ φ = = Qβ . |Cφ |2 |Cφ |2 (ξαβ − Eφ ) Qφα dRN ∂R ∂R N N φ φ αβ

(12)

The nonadiabatic force terms are obtained from the TDSE (Eq. 2) and can be written as X φ





dCφ∗ dCφ Cφ + Cφ∗ dRN dRN



φ6=ψ

=

X

ˆ Cφ∗ Cψ hφ|∇N H|ψi .

(13)

φ,ψ

The above equation is associated with the nonadiabatic coupling term dN ψφ = hφ|∇N ψi via equation dN ψφ

ˆ hφ|∇N H|ψi = = Eφ − Eψ

P

αβ

ˆ Qφα Qψβ hα|∇N H|βi , Eφ − Eψ

(14)

ˆ ˆ ˆ where hα|∇N H|βi = ∇N Hαβ − h∇N α|H|βi − hα|H|∇ βi. To proceed we write the extended N Hückel hamiltonian in operator representation as ˆ = H

X αβ

  |αi S−1 HS−1 αβ hβ|

(15)

and make use of the extended Hückel generalized eigenvalue equation. Here we present the

10

ACS Paragon Plus Environment

Page 11 of 51

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

ensuing expression for the nonadiabatic component of the force on atom N X

Eφ ∇N |Cφ |2 = 2

φ

X

ψ>φ

 Re Cφ∗ Cψ ×

XX

n∈N ∀α

Qφn Qψα (ξnα − Eψ ) + Qφα Qψn (ξαn − Eφ ) hfα |∇N fn i . (16)

For the sake of clarity, we use the notation fη ≡ fn (r − RN ) that stands for atomic orbital n centered at atom N and likewise for fα ≡ fa (r − RA ). Finally, we can write in closed form the complete force produced on atom N by the electron-hole excitation, including adiabatic and nonadiabatic terms, within the extended Hückel theory formalism

= −2 FEH N

XX

hfα |∇N fn (r − RN )i

n∈N ∀α

+

X

ψ>φ

(

X

φ φ ρEH φφ (ξnα − Eφ ) Qα Qn

φ

 φ ψ  Qn Qα (ξnα − Eϕ ) + Qφα Qψn (ξαn − Eφ ) Re ρEH φψ

)

,

(17)

el hl φ recalling that ρEH φψ = ρφψ − ρφψ , and that Q and Eφ are the adiabatic eigenstate vector and

its energy, and ξαβ = Kαβ (hα + hβ )/2 is the extended Hückel coupling element. It is important to notice that the total energy of the photoexcited system is conserved during the dynamics simulation

  MM E = EM M + EQM = K + VGS + T r ρEH H ,

(18)

MM where EM M = K + VGS is the overall classical energy (kinetic plus potential) of the ions   in the molecular mechanics formalism and EQM = T r ρEH H is the quantum energy of the

photoexcited electron-hole pair.

11

ACS Paragon Plus Environment

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 ACS Paragon Plus Environment

Page 12 of 51

Page 13 of 51

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 ACS Paragon Plus Environment

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 ACS Paragon Plus Environment

Page 14 of 51

Page 15 of 51

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

intersection, which is more realistic than static calculations carried along constrained reaction paths. The present NA-MD simulations show that, as the geometry of the molecule becomes non-planar, the LUMO and HOMO states collapse into a conical intersection and a population transfer occurs between S1 and S0 states. The internal conversion is also evinced     by the energy of the wavepackets, Eel = T r ρel H and Ehl = −T r ρhl H , as shown in

Figure 4-a). To compare the results obtained with the mixed quantum-classical Ehrenfest method with ab-initio formalisms we have calculated the energy of the S1 and S0 states

with the complete active space self-consistent field method (CASSCF(6,6)/6-31G∗ ) for the geometries obtained from the NA-MD simulations. The results are presented in Figure 4-b). The agreement between the semiempirical and the ab-initio method is excellent for benzene, including the position of the conical intersection as well as the gap of the molecule; for the sake of comparison, the CASSCF(6,6) total energy of the S0 state (-6280 eV) is set to zero at t=0. To make the comparison between CASSCF and EHT more direct, Figure S1 (in the Supporting Information) shows the LUMO and HOMO energies obtained from EHT vis-à-vis CASSCF energies.

Figure 5: Energy conservation during IVR in benzene. The overall QM-MM energy of the quantum and classical degrees of freedom; details provided in the text. The energy is conserved within numerical precision. The time stept of the mixed quantum-classical simulation is dt = 0.005 fs; results improve as the time step of simulation decreases. We now analyze QM-MM energy balance during the nonadiabatic internal conversion process (Figure 5). As previously determined by the initial conditions, the simulation started from rest with the S0 → S1 vertical excitation at the geometry optimized in the ground15

ACS Paragon Plus Environment

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

  MM state. Therefore, E(t = 0) = EQM = T r ρEH H , since VGS = 0. It is important to notice   that the existence of the e-h pair implies that T r ρEH H > 0. During the excited-state

MD, part of the e-h pair energy is converted into vibrational energy in the molecule, thus   MM E = EM M + EQM = K + VGS + T r ρEH H . It is noticeable in Figures 3 and 4-a) that the system dwells for a short period in the conical intersection, but it soon decays to the S0

state and remains there. If the vibrational relaxation is effective the excess energy is quickly distributed throughout the vibrational modes that can not drive the system back through the conical intersection. In the ground state (S0 ) the energy is computed as E = EM M =   MM K + VGS , because T r ρEH H < 0, and thus the force on a given atom N is simply given M MM by FM = −∇N VGS . As evinced by Figure 5, the method conserves QM-MM total energy N

within numerical precision. The results for energy conservation improve as the time step of the simulation decreases; the time step used in the present simulation was dt = 0.005 fs. It is also evident from the graph that it is more challenging to conserve the energy during the excited state dynamics, however, the overall energy drift is smaller that 10 meV whereas the energy exchange from QM to MM degrees of freedom is of the order of a few eV. The small energy drift indicates that the ansatz of Eq. (6) is indeed a good approximation for the excited state PES.

Mechanism for Photoisomerization of Stilbene Photoisomerization of stilbene (Figure 2-a) has been the subject much debate. 65–67 The difficulties to describe it in terms of standard Rice-Ramsperger-Kassel-Marcus (RRKM) theory 68 rendered it the character of a textbook case of a non-RRKM process. 66 Common reasons for discrepancy between theory and experiment in cases like this can be ascribed to uncertainties on vibrational modes at both reactant and transition states, which give rise to the reaction path, as well as inaccuracy of barrier heights. Though rates of isomerization are obtained with statistical mechanics procedures, like in the RRKM theory, such calculations must be based on sound models that include the correct dynamical mechanisms. Understanding of 16

ACS Paragon Plus Environment

Page 16 of 51

Page 17 of 51

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

the underlying dynamics is important to build reliable statistical models to describe the effect under thermodynamical conditions. Therefore, we apply the nonadiabatic Ehrenfest method to investigate the molecular mechanisms that trigger isomerization in stilbene after photoexcitation. In the following, we evince that the torsion of the ethylenic bond (C = C) in isomerization consists actually of a concerted two step process: first, the torsion of the H − C = C − H dihedral, followed by the rotation of the phenyl rings around the ethylenic bond. For the NA-MD simulations on stilbene, the preparation of the system differed from the case of benzene. Before carrying out the excited state NA dynamics, we performed a ground state classical MD run on stilbene, putting the molecule in contact with a NoséHoover thermostat at 300 K for a period of 100 ps (details provided in the Supporting Information). During the entire GS dynamics the molecule remained in the same (initial) conformation, that is, the trans geometry. Moreover, the ethylenic bond at the bridge did not rotate significantly (torsion angle < 10°). After the initial thermalization, we took the ensuing molecular geometry to start the mixed quantum-classical NA-MD simulation. We thus produced in it an electron-hole excitation by placing an electron in the LUMO and a hole in the HOMO. Therefore, we start from the S0 → S1 Franck-Condon excitation at the geometry of the thermalized ground-state; the atomic positions and velocities at the end of the classical MD were used as initial conditions for the NA-MD in the excited state. Figure 6 shows the isomerization dynamics of stilbene in terms of the H − C = C − H (red) and C − C = C − C (black) dihedrals as a function of time, after photoexcitation. The striking feature of the dynamics is the abrupt trans → cis rotation of the H − C = C − H dihedral, thus setting up the reaction path for isomerization of the phenyl rings. The effect happens twice, first ∼ 150 fs and later, around 700 ps, when it rotates back from cis to trans. Molecular geometries at instantaneous configurations show the role of the atomic mechanism. For the back rotation (cis → trans), the isomerization starts likewise by the rotation of the H − C = C − H dihedral, being followed by the phenyl rings. The effect is

17

ACS Paragon Plus Environment

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 ACS Paragon Plus Environment

Page 18 of 51

Page 19 of 51

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

cm−1 ) and the ethylenic HOOP wag (956 cm−1 ) vibrational modes that can be associated with the bridge. For a similar ethylenic structure, the retinal molecule in rhodopsin, femtosecond time-resolved spectroscopy 10 has detected the following vibrational modes that can also be associated with the bridge: the ethylenic C = C stretch at

1550 cm−1 , C − C

single-bond stretching and C − H rocking modes in the region from 1100 to 1300 cm−1 , and hydrogen-out-of-plane (HOOP) modes between 800 and 1000 cm−1 . Signatures of these modes appear at the velocity density of states (VDOS), shown in Figure 7, which is obtained from the short-time Fourier transform of the velocity auto-correlation functions performed for the the atoms that comprise the ethylenic bridge along the NA-MD trajectory (for details see Supporting Information). At the beginning the vibrational energy is concentrated around 1300 cm−1 , due to excitation of the C − C stretch mode by the e-h impulsive force as well as the C − H rocking modes evidenced in Figure 6. That energy is transferred afterwards to the other vibrational modes of the bridge. It is interesting to notice that the vibrational energy of the bridge decays considerably after 500 ps and it is negligible at 1 ps.

Photoisomerization of Azobenzene Azobenzene is a prototypical photochromic molecule. Its reversible photoisomerization properties rendered it various applications in photoswitches and nanoscale devices, for both biological and industrial areas. 1,3 Azobenzene derivatives are the most used compounds among photoswitches owing to their good light absorption, photoisomerization efficiency and, specially, the large relative conformational change as compared to other molecules. Besides the technological interest, knowledge of the underlying molecular mechanisms driving its photochromism has prompted several studies, 1,4,8,9,70–75 apart from questions of quantum yield and influence of environmental conditions. Despite the wide applicability of azobenzene and its reliability as a photochromic compound, the excited state dynamics of azobenzene remains a subject of debate, as indicated by recent studies. 1,4,73 Conclusions gathered from various investigations lead to four general 19

ACS Paragon Plus Environment

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

isomerization mechanisms: 1,76 (1) rotation, (2) inversion, (3) concerted inversion and (4) inversion-assisted rotation. In short, the rotation mechanism (1) involves the change of the N = N bond order, followed by torsion of the C − N − N − C dihedral with the N − N − C angles fixed at ∼ 120°. The inversion mechanism (2) has the C − N = N − C dihedral angle fixed at 0° whereas one of the N = N − C increases to 180°. In the concerted inversion (3) a transition state is formed when both N = N − C bond angles increase to ∼ 180°. Finally, the inversion-assisted rotation (4) is a dynamic mechanism that arises when a complete rotation of the C − N − N − C dihedral angle is assisted by large changes in N = N − C bond angles (100° . ∡N N C . 140°). Unfortunately, the nomenclature found in literature to designate the isomerization mechanisms is not consistent, particularly for describing the combination of rotation and inversion motions. Therefore, we adopt the above definition. In the following we apply the Nonadiabatic Ehrenfest method to study the photoisomerization dynamics of azobenzene, concentrating our interest at the isomerization mechanism and the timescales of the process, which shall be compared with the available experimental and ab-initio results for the molecule in gas phase. Before carrying out the excited-state simulations, we performed a classical MD in the ground state, with the system in contact with a Nosé-Hoover thermostat for a period of 100 ps. During this period the molecule remained in the original conformation, either trans or cis, and did not evince any hint of isomerization. The same set of parameters were used for the stilbene and azobenzene simulations. More details about the system preparation are provided in the Supporting Information. After thermalization, the system is prepared in the S2 state, assuming a π → π ∗ FranckCondon excitation at the geometry of the thermalized ground-state. Atomic positions and velocities at the end of the classical MD are used to initiate the excited state NA-MD discussed next. Figure 8 shows the ensuing S2 → S0 IVR process, as the electronic wavepacket (blue) and the hole wavepacket (red) follow a reaction path to the conical intersection. For t . 600 fs the electronic wavepacket evolves almost exclusively on the LUMO (π ∗ ) whereas the hole is described by a mixture of the nonbonding (n) and bonding (π) orbitals. The

20

ACS Paragon Plus Environment

Page 20 of 51

Page 21 of 51

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Figure 8: Photoisomerization dynamics of the azobenzene molecule (Figure 2-b) after the S2 state (π → π ∗ ) Franck-Condon excitation. Curves describe time evolution of the wavepacket coefficients |Cφ (t)|2 ≡ |C(Eφ (t), t)|2 for electron (blue) and hole (red). The vertical axis, Eφ (t), is the energy of the time evolving adiabatic molecular orbital φ(t). The color code for |C(Eφ (t), t)|2 is so defined: blue (positive) describes the electron wavepacket and red (designated by negative values) the hole wavepacket. behavior of the hole is consistent with experimental reports 8,9,70,73 of a fast S2 → S1 internal conversion that maintains the planar geometry, thus producing a hot S1 species that can overcome occasional energy barriers for isomerization. Reported lifetimes for the ππ ∗ state vary from ∼ 110 to 170 fs, 8,9 in addition to ∼ 500 fs for the S1 state. 9 The graph clearly shows the time scales of the process: intramolecular vibrational energy distribution and relaxation put the system at the isomerization path in ∼ 600 fs; once the edge of the conical intersection is reached the system quickly goes through it within ∼ 300 fs, and isomerization is complete within 1 ps. The average energies of the electron (blue) and hole (red), shown     in Figure 9, given by Eel = T r ρel H and Ehl = −T r ρhl H , corroborate the wavepacket

dynamic description. The total QM-MM energy drift for this simulation is 13 meV (shown in Supporting Information). We now analyze the structural dynamics of azobenzene during isomerization, for the previous NA-MD simulation. For that matter, Figure 10 shows two graphs: above, the temporal behavior of the C − N − N − C dihedral and, below, the C − N − N (green) and N − N − C 21

ACS Paragon Plus Environment

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

 el  Figure 9: The average  hl  energies of electron (blue) and hole (red), given by Eel = T r ρ H∗ and Ehl = −T r ρ H , as a function of the time, for NA-MD simulation upon S2 (π → π ) excitation in azobenzene. (orange) angles as a function of time. The NA-MD starts with the azobenzene in the trans geometry and, after internal conversion, it isomerizes to cis conformation at ∼ 0.9 ps. At this point the system decays to S0 , abruptly transferring the energy of the e-h pair to the vibrational modes of the reaction coordinate. This tremendous vibrational energy is not, however, transferred to other vibrational modes as efficiently, thus the azobenzene continues to isomerize along the path cis → trans → cis until the simulation ends, as shown in Figure 10. As observed in other simulations, the vibrational energy contained in the reaction coordinate is eventually transferred to other vibrational modes and the molecule stabilizes in a specific conformation. The lack of anharmonic terms in the MM-force field, Eq. (3), is a factor that undermines intramolecular vibrational energy distribution. Looking at the structural data, it evinces that the C − N − N − C dihedral rotates easily around its axis. In regard to the C − N − N (green) and N − N − C (orange) angles, they never approach 180°, as supposed to characterize the inversion or the concerted inversion mechanisms. Figure 10 shows that these angles oscillate vigorously between 100° and 130° though, particularly ∼ 0.9 ps, which characterizes the inversion-assisted mechanism. Moreover, by analyzing the simulation frames during the isomerization process we observe that the C − N − N − C dihedral

22

ACS Paragon Plus Environment

Page 22 of 51

Page 23 of 51

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 ACS Paragon Plus Environment

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 ACS Paragon Plus Environment

Page 24 of 51

Page 25 of 51

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

ticular, the pedal-like mechanism is also observed in previous studies of azobenzene 4,74,75 and other molecular systems. 77 Fewer recent articles, 8,9 to our knowledge, support the plain inversion mechanism. Therefore, we conclude that the isomerization mechanism of azobenzene is a multidimensional one, with rotation dominating the process but assisted by inversion.

Conclusions We have developed a new mixed quantum-classical Ehrenfest self-consistent method for nonadiabatic molecular dynamics in the excited state and used it to investigate the molecular mechanisms of photoisomerization. The method combines two simple though efficient computational schemes: the quantum mechanical part is based on the extended Hückel formalism whereas the nuclei are treated with the molecular mechanics method. The quantum-classical self-consistent coupling is based on the approach

  MM hΨ(r; R, t)|V (r, R)|Ψ(r; R, t)ir ≈ VGS + VEH Ψel (t), Ψhl (t) , which assumes that the excited state PES can be approximated as the ground state MM MM force field potential, VGS , plus a correction that is due to the electron-hole excitation,  el  VEH Ψ (t), Ψhl (t) . Moreover, the resulting force on atom N, which is given by the closed-

form expression in Eq. (17), includes adiabatic and nonadiabatic effects. The viability of the

approximation is proved by the conservation of the overall quantum-classical energy within numerical precision, as evinced by energy drifts of the order of a few meV in comparison to an energy exchange of ∼ 4 eV between QM and MM degrees of freedom during nonradiative relaxation. The method conserves energy since the semiempirical EH hamiltonian is written as a function of the overlap matrix, thereby the nuclear dynamics affects directly the electronic energy and its quantum dynamics. Since this is a mean-field description some effects are neglected, such as quantum-classical correlation effects that are described by surfacehopping. It is an open question, though, whether the present approach can be adapted to 25

ACS Paragon Plus Environment

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

surface-hopping. We have described the mechanisms of internal conversion in benzene, evincing the very good agreement of our semiempirical calculations with ab-initio CASSCF results. Furthermore, the method has been used to gain insights into the atomistic mechanisms that control the photoisomerization dynamics for stilbene and azobenzene molecules. For stilbene, our simulations demonstrate that the isomerization mechanism consists of a two step process, which initiates with an abrupt torsion of the H − C = C − H dihedral, triggered by hydrogen out-of-plane (HOOP) vibrations, followed by the rotation of the phenyl rings around the ethylenic bond. The short-time vibrational density of states (VDOS) that is calculated from the excited state NA-MD trajectory shows strong signatures of C − C stretching, C − H rocking and HOOP modes. Calculations performed on the azobenzene molecule identified the inversion-assisted rotation as the driving isomerization mechanism for both the ππ ∗ and nπ ∗ excitations in isolated molecules. Moreover, the C − N − N − C dihedral exhibits a pedal-cycle motion (see Figure 11) during the isomerization. For all the molecules investigated, the nonradiative decay through the conical intersection happened in less than 200 fs whereas the excess vibrational energy stored in the modes of the reaction path decayed with a longer time scale. To conclude, we point out that this work has demonstrated the capabilities of an approach that could be used to extend molecular mechanics methods to the excited state realm with reduced computational costs. However, in the present form the method is designed for simulation periods around ten picoseconds, which is the time that coherent quantum effects can be relevant. For longer simulations, conventional MD simulations might be more appropriate.

Supporting Information Available This document contains: description of the classical (MM) parameterization; description of quantum (extended Hückel, EH) parameterization; description the system preparation for NA-MD simulations; brief introduction to the extended Hückel (EH) method and derivation 26

ACS Paragon Plus Environment

Page 26 of 51

Page 27 of 51

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

of Hellmann-Feynman-Pulay force for the EH method; description of the procedures for calculating the velocity density of states (VDOS). Three movie files displaying the nonadiabatic dynamics of benzene, stilbene and azobenzene. This material is available free of charge via the Internet at http://pubs.acs.org/.

Acknowledgement The authors are indebted to the Brazilian National Counsel of Technological and Scientific Development (CNPq) and the Coordenação de Aperfeiçoamento de Pessoal de Nível Superior (CAPES) for funding the project.

References (1) Bandara, H. M. D.; Burdette, S. C. Photoisomerization in different classes of azobenzene. Chem. Soc. Rev. 2012, 41, 1809–1825. (2) Galperin, M.; Nitzan, A. Molecular optoelectronics: the interaction of molecular conduction junctions with light. Phys. Chem. Chem. Phys. 2012, 14, 9421–9438. (3) Garcia-Iriepa, C.; Marazzi, M.; Frutos, L. M.; Sampedro, D. E/Z Photochemical switches: syntheses, properties and applications. RSC Adv. 2013, 3, 6241–6266. (4) Bockmann, M.; Marx, D.; Peter, C.; Site, L. D.; Kremer, K.; Doltsinis, N. L. Multiscale modelling of mesoscopic phenomena triggered by quantum events: light-driven azomaterials and beyond. Phys. Chem. Chem. Phys. 2011, 13, 7604–7621. (5) Curchod, B. F. E.; Rothlisberger, U.; Tavernelli, I. Trajectory-Based Nonadiabatic Dynamics with Time-Dependent Density Functional Theory. ChemPhysChem 2013, 14, 1314–1340. (6) Zewail, A. H. Femtochemistry: Atomic-scale dynamics of the chemical bond. J. Phys. Chem. A 2000, 104, 5660–5694. 27

ACS Paragon Plus Environment

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(7) Radloff, W.; Stert, V.; Freudenberg, T.; Hertel, I.; Jouvet, C.; Dedonder-Lardeux, C.; Solgadi, D. Internal conversion in highly excited benzene and benzene dimer: femtosecond time-resolved photoelectron spectroscopy. Chem. Phys. Lett. 1997, 281, 20 – 26. (8) Schultz, T.; Quenneville, J.; Levine, B.; Toniolo, A.; Martínez, T. J.; Lochbrunner, S.; Schmitt, M.; Shaffer, J. P.; Zgierski, M. Z.; Stolow, A. Mechanism and dynamics of azobenzene photoisomerization. J. Am. Chem. Soc. 2003, 125, 8098–8099. (9) Fujino, T.; Arzhantsev, S. Y.; Tahara, T. Femtosecond time-resolved fluorescence study of photoisomerization of trans-azobenzene. J. Phys. Chem. A 2001, 105, 8123–8129. (10) Kukura, P.; McCamant, D. W.; Mathies, R. A. Femtosecond Stimulated Raman Spectroscopy. Annu. Rev. Phys. Chem. 2007, 58, 461–488. (11) Gessner, O.; Gühr, M. Monitoring Ultrafast Chemical Dynamics by Time-Domain Xray Photo- and Auger-Electron Spectroscopy. Acc. Chem. Res. 2016, 49, 138–145. (12) Chen, J.; Wu, K.; Rudshteyn, B.; Jia, Y.; Ding, W.; Xie, Z.-X.; Batista, V. S.; Lian, T. Ultrafast Photoinduced Interfacial Proton Coupled Electron Transfer from CdSe Quantum Dots to 4,4’-Bipyridine. J. Am. Chem. Soc. 2016, 138, 884–892. (13) Long, R.; English, N. J.; Prezhdo, O. V. Minimizing Electron–Hole Recombination on TiO2 Sensitized with PbSe Quantum Dots: Time-Domain Ab Initio Analysis. J. Phys. Chem. Lett. 2014, 5, 2941–2946. (14) Liu, J.; Neukirch, A. J.; Prezhdo, O. V. Non-radiative electron–hole recombination in silicon clusters: ab initio non-adiabatic molecular dynamics. J. Phys. Chem. C 2014, 118, 20702–20709. (15) Zhu, H.; Yang, Y.; Wu, K.; Lian, T. Charge Transfer Dynamics from Photoexcited Semiconductor Quantum Dots. Annu. Rev. Phys. Chem. 2016, 67, 259–281. 28

ACS Paragon Plus Environment

Page 28 of 51

Page 29 of 51

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

(16) Born, M.; Oppenheimer, R. Zur quantentheorie der molekeln. Annalen der Physik 1927, 389, 457–484. (17) Tully, J. C. Molecular dynamics with electronic transitions. J. Chem. Phys. 1990, 93, 1061–1071. (18) Tully, J. C. Mixed quantum–classical dynamics. Faraday Discuss. 1998, 110, 407–419. (19) Tully, J. C. In Classical and Quantum Dynamics in Condensed Phase Simulations; Berne, B. J., Ciccotti, G., Coker, D. F., Eds.; 1998; Vol. 1; pp 489–514. (20) Guallar, V.; Batista, V. S.; Miller, W. H. Semiclassical molecular dynamics simulations of excited state double-proton transfer in 7-azaindole dimers. J. Chem. Phys. 1999, 110, 9922–9936. (21) Ben-Nun, M.; Quenneville, J.; Martínez, T. J. Ab initio multiple spawning: Photochemistry from first principles quantum molecular dynamics. J. Phys. Chem. A 2000, 104, 5161–5175. (22) Tao, H.; Levine, B. G.; Martínez, T. J. Ab initio multiple spawning dynamics using multi-state second-order perturbation theory. J. Phys. Chem. A 2009, 113, 13656– 13662. (23) Levine, B. G.; Coe, J. D.; Virshup, A. M.; Martinez, T. J. Implementation of ab initio multiple spawning in the Molpro quantum chemistry package. Chem. Phys. 2008, 347, 3–16. (24) Mori, T.; Glover, W. J.; Schuurman, M. S.; Martinez, T. J. Role of Rydberg states in the photochemical dynamics of ethylene. J. Phys. Chem. A 2012, 116, 2808–2818. (25) Martinez, T. J.; Ben-Nun, M.; Levine, R. Multi-electronic-state molecular dynamics: A wave function approach with applications. J. Phys. Chem. 1996, 100, 7884–7895.

29

ACS Paragon Plus Environment

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(26) Martinez, T.; Ben-Nun, M.; Levine, R. Molecular collision dynamics on several electronic states. J. Phys. Chem. A 1997, 101, 6389–6402. (27) Heller, E. J. Time-dependent approach to semiclassical dynamics. J. Chem. Phys. 1975, 62, 1544–1555. (28) Meyer, H.-D.; Manthe, U.; Cederbaum, L. S. The multi-configurational time-dependent Hartree approach. Chem. Phys. Lett. 1990, 165, 73–78. (29) Beck, M. H.; Jäckle, A.; Worth, G.; Meyer, H.-D. The multiconfiguration timedependent Hartree (MCTDH) method: a highly efficient algorithm for propagating wavepackets. Phys. Rep. 2000, 324, 1–105. (30) van Harrevelt, R.; Manthe, U. Multidimensional time-dependent discrete variable representations in multiconfiguration Hartree calculations. J. Chem. Phys. 2005, 123, 064106. (31) Tully, J. C.; Preston, R. K. Trajectory surface hopping approach to nonadiabatic molecular collisions: The reaction of H+ with D2 . J. Chem. Phys. 1971, 55, 562–572. (32) Barbatti, M. Nonadiabatic dynamics with trajectory surface hopping method. Wiley Interdiscip. Rev.: Comput. Mol. Sci. 2011, 1, 620–633. (33) Ehrenfest, P. Bemerkung über die angenäherte Gültigkeit der klassischen Mechanik innerhalb der Quantenmechanik. Zeitschrift für Physik A Hadrons and Nuclei 1927, 45, 455–457. (34) Fabiano, E.; Groenhof, G.; Thiel, W. Approximate switching algorithms for trajectory surface hopping. Chem. Phys. 2008, 351, 111–116. (35) Lasser, C.; Swart, T. Single switch surface hopping for a model of pyrazine. J. Chem. Phys. 2008, 129, 034302.

30

ACS Paragon Plus Environment

Page 30 of 51

Page 31 of 51

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

(36) Wang, L.; Akimov, A.; Prezhdo, O. V. Recent Progress in Surface Hopping: 2011-2015. J. Phys. Chem. Lett. 2016, 7, 2100–2112. (37) Wang, L.; Trivedi, D.; Prezhdo, O. V. Global Flux Surface Hopping Approach for Mixed Quantum-Classical Dynamics. J. Chem. Theory Comput. 2014, 10, 3598–3605. (38) Wang, L.; Sifain, A. E.; Prezhdo, O. V. Fewest Switches Surface Hopping in Liouville Space. J. Phys. Chem. Lett. 2015, 6, 3827–3833. (39) Wang, L.; Sifain, A. E.; Prezhdo, O. V. Communication: Global flux surface hopping in Liouville space. J. Chem. Phys. 2015, 143 . (40) Bornemann, F. A.; Nettesheim, P.; Schütte, C. Quantum-classical molecular dynamics as an approximation to full quantum dynamics. J. Chem. Phys. 1996, 105, 1074–1083. (41) Parandekar, P. V.; Tully, J. C. Mixed quantum-classical equilibrium. The Journal of Chemical Physics 2005, 122 . (42) Wang, L.; Beljonne, D. Charge transport in organic semiconductors: Assessment of the mean field theory in the hopping regime. The Journal of Chemical Physics 2013, 139 . (43) MuÌĹller, U.; Stock, G. Surface-hopping modeling of photoinduced relaxation dynamics on coupled potential-energy surfaces. J. Chem. Phys. 1997, 107, 6230–6245. (44) Barbatti, M.; Granucci, G.; Persico, M.; Ruckenbauer, M.; Vazdar, M.; EckertMaksić, M.; Lischka, H. The on-the-fly surface-hopping program system Newton-X: Application to ab initio simulation of the nonadiabatic photodynamics of benchmark systems. J. Photochem. Photobiol., A 2007, 190, 228–240. (45) Doltsinis, N. L.; Marx, D. Nonadiabatic Car-Parrinello molecular dynamics. Phys. Rev. Lett. 2002, 88, 166402.

31

ACS Paragon Plus Environment

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(46) Castro, A.; Appel, H.; Oliveira, M.; Rozzi, C. A.; Andrade, X.; Lorenzen, F.; Marques, M. A.; Gross, E.; Rubio, A. Octopus: a tool for the application of time-dependent density functional theory. Phys. Status Solidi B 2006, 243, 2465–2488. (47) Du, L.; Lan, Z. An on-the-fly surface-hopping program jade for nonadiabatic molecular dynamics of polyatomic systems: implementation and applications. J. Chem. Theory Comput. 2015, 11, 1360–1374. (48) Akimov, A. V.; Prezhdo, O. V. The PYXAID program for non-adiabatic molecular dynamics in condensed matter systems. J. Chem. Theory Comput. 2013, 9, 4959–4972. (49) Akimov, A. V. Libra: An open-Source “methodology discovery” library for quantum and classical dynamics simulations. J. Comp. Chem. 2016, 37, 1626–1649. (50) Rego, L. G. C.; Oliboni, R. S.; Hoff, D. A.; Torres, A. DynEMol: Tools for studying Dynamics of Electrons in Molecules. 2014; http://sourceforge.net/projects/ charge-transfer/. (51) Rego, L. G. C.; Hames, B. C.; Mazon, K. T.; Joswig, J.-O. Intramolecular Polarization Induces Electron–Hole Charge Separation in Light-Harvesting Molecular Triads. J. Phys. Chem. C 2014, 118, 126–134. (52) Hoff, D. A.; Silva, R.; Rego, L. G. C. Subpicosecond dynamics of metal-to-ligand chargetransfer excited states in solvated [Ru(bpy)3 ]2+ complexes. J. Phys. Chem. C 2011, 115, 15617–15626. (53) Hoff, D. A.; da Silva, R.; Rego, L. G. C. Coupled electron–hole quantum dynamics on D–π–A dye-sensitized TiO2 semiconductors. J. Phys. Chem. C 2012, 116, 21169– 21178. (54) Rego, L. G. C.; da Silva, R.; Freire, J. A.; Snoeberger, R. C.; Batista, V. S. Visible

32

ACS Paragon Plus Environment

Page 32 of 51

Page 33 of 51

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Light Sensitization of TiO2 Surfaces with Alq3 Complexes. J. Phys. Chem. C 2010, 114, 1317–1325. (55) Bairu, S. G.; Mghanga, E.; Hasan, J.; Kola, S.; Rao, V. J.; Bhanuprakash, K.; Giribabu, L.; Wiederrecht, G. P.; da Silva, R.; Rego, L. G. C.; Ramakrishna, G. Ultrafast Interfacial Charge-Transfer Dynamics in a Donor−π−Acceptor Chromophore Sensitized TiO2 2 Nanocomposite. J. Phys. Chem. C 2013, 117, 4824–4835. (56) Monti, A.; Negre, C. F. A.; Batista, V. S.; Rego, L. G. C.; de Groot, H. J. M.; Buda, F. Crucial Role of Nuclear Dynamics for Electron Injection in a Dye-Semiconductor Complex. J. Phys. Chem. Lett. 2015, 6, 2393–2398. (57) Torres, A.; Oliboni, R. S.; Rego, L. G. C. Vibronic and Coherent Effects on Interfacial Electron Transfer Dynamics. J. Phys. Chem. Lett. 2015, 6, 4927–4935. (58) Hoff, D. A.; Rego, L. G. C. In Specialist Periodical Reports – Chemical Modelling; Springboard, M., Joswig, J. O., Eds.; 2013; pp 102–126. (59) Jorgensen, W. L.; Maxwell, D. S.; Tirado-Rives, J. Development and testing of the OPLS all-atom force field on conformational energetics and properties of organic liquids. J. Am. Chem. Soc. 1996, 118, 11225–11236. (60) Kaminski, G. A.; Friesner, R. A.; Tirado-Rives, J.; Jorgensen, W. L. Evaluation and reparametrization of the OPLS-AA force field for proteins via comparison with accurate quantum chemical calculations on peptides. J. Phys. Chem. B 2001, 105, 6474–6487. (61) Hoffmann, R. An extended Hückel theory. I. Hydrocarbons. J. Chem. Phys. 1963, 39, 1397–1412. (62) Parmenter, C. S. Adv. Chem. Phys.; John Wiley & Sons, Inc., 2007; Vol. 22; pp 365–421. (63) Krogh-Jespersen, K.; Rava, R. P.; Goodman, L. Potential energy surface of the 1 B2u

33

ACS Paragon Plus Environment

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

excited state in benzene. Mode forms and spectroscopic consequences. J. Phys. Chem. 1984, 88, 5503–5513. (64) Thompson, A. L.; Martinez, T. J. Time-resolved photoelectron spectroscopy from first principles: Excited state dynamics of benzene. Faraday Discuss. 2011, 150, 293–311. (65) Myers, A. B. Advances in Multiphoton Processses and Spectroscopy; Advances in MultiPhoton Processes and Spectroscopy Series v. 10; World Scientific, 1998; Chapter Resonance Raman intensities: the roundabout way to Franck-Condon analysis. (66) Leitner, D. M.; Levine, B.; Quenneville, J.; Martinez, T. J.; Wolynes, P. G. Quantum Energy Flow and trans-Stilbene Photoisomerization: an Example of a Non-RRKM Reaction. J. Phys. Chem. A 2003, 107, 10706–10716. (67) Quenneville, J.; Martinez, T. J. Ab Initio Study of Cis-Trans Photoisomerization in Stilbene and Ethylene. J. Phys. Chem. A 2003, 107, 829–837. (68) Han, K.; Chi, T. Reaction Rate Constant Computations: Theories and Applications; RSC Theor. Comput. Chem. Ser.; Royal Society of Chemistry, 2013. (69) Ci, X.; Myers, A. B. A resonance Raman study of solvent effects on the excited state potential surface of trans-stilbene. Chem. Phys. Lett. 1989, 158, 263 – 270. (70) Crecca, C. R.; Roitberg, A. E. Theoretical Study of the Isomerization Mechanism of Azobenzene and Disubstituted Azobenzene Derivatives. J. Phys. Chem. A 2006, 110, 8188–8203. (71) Ciminelli, C.; Granucci, G.; Persico, M. The photoisomerization mechanism of azobenzene: A semiclassical simulation of nonadiabatic dynamics. Chem.–Eur. J. 2004, 10, 2327–2341.

34

ACS Paragon Plus Environment

Page 34 of 51

Page 35 of 51

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

(72) Shao, J.; Lei, Y.; Wen, Z.; Dou, Y.; Wang, Z. Nonadiabatic Simulation Study of Photoisomerization of Azobenzene: Detailed mechanism and load-resisting capacity. J. Chem. Phys. 2008, 129, 164111. (73) Tan, E. M.; Amirjalayer, S.; Smolarek, S.; Vdovin, A.; Zerbetto, F.; Buma, W. J. Fast photodynamics of azobenzene probed by scanning excited-state potential energy surfaces using slow spectroscopy. Nat. Commun. 2015, 6, 5860. (74) Neukirch, A. J.; Shamberger, L. C.; Abad, E.; Haycock, B. J.; Wang, H.; Ortega, J.; Prezhdo, O. V.; Lewis, J. P. Nonadiabatic Ensemble Simulations of cis-Stilbene and cis-Azobenzene Photoisomerization. J. Chem. Theory Comput. 2014, 10, 14–23. (75) Tiberio, G.; Muccioli, L.; Berardi, R.; Zannoni, C. How Does the Trans-Cis Photoisomerization of Azobenzene Take Place in Organic Solvents? ChemPhysChem 2010, 11, 1018–1028. (76) Diau, E. W.-G. A New Trans-to-Cis Photoisomerization Mechanism of Azobenzene on the S1(n,π ∗ ) Surface. J. Phys. Chem. A 2004, 108, 950–956. (77) Chung, W. C.; Nanbu, S.; Ishida, T. Nonadiabatic ab Initio Dynamics of a Model Protonated Schiff Base of 9-cis Retinal. J. Phys. Chem. A 2010, 114, 8190–8201.

35

ACS Paragon Plus Environment

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 ACS Paragon Plus Environment

Page 36 of 51

Page 37 of 51

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

F1 for paper

ACS Paragon Plus Environment

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

F2 for paper

ACS Paragon Plus Environment

Page 38 of 51

Page 39 of 51

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

F3 for paper

ACS Paragon Plus Environment

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

F4 for paper

ACS Paragon Plus Environment

Page 40 of 51

Page 41 of 51

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

F5 for paper 294x159mm (97 x 97 DPI)

ACS Paragon Plus Environment

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

F6 for paper

ACS Paragon Plus Environment

Page 42 of 51

Page 43 of 51

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

F7 for paper 458x264mm (72 x 72 DPI)

ACS Paragon Plus Environment

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

F8 for paper 493x306mm (72 x 72 DPI)

ACS Paragon Plus Environment

Page 44 of 51

Page 45 of 51

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

F9 for paper 281x187mm (99 x 99 DPI)

ACS Paragon Plus Environment

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

F10 for paper

ACS Paragon Plus Environment

Page 46 of 51

Page 47 of 51

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

F11 for paper

ACS Paragon Plus Environment

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

F12 for paper

ACS Paragon Plus Environment

Page 48 of 51

Page 49 of 51

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

TOC figure 82x44mm (300 x 300 DPI)

ACS Paragon Plus Environment

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure 1 for SI document 271x176mm (72 x 72 DPI)

ACS Paragon Plus Environment

Page 50 of 51

Page 51 of 51

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

figure S2 for Supportinf Information 84x99mm (300 x 300 DPI)

ACS Paragon Plus Environment