Direct Non-Adiabatic Dynamics by Mixed Quantum-Classical

Subscriber access provided by Kaohsiung Medical University

Dynamics

Direct Non-Adiabatic Dynamics by Mixed Quantum-Classical Formalism Connected with Ensemble Density Functional Theory Method: Application to trans-Penta-2,4-Dieniminium Cation Michael Filatov, Seung Kyu Min, and Kwang S. Kim J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.8b00217 • Publication Date (Web): 27 Jul 2018 Downloaded from http://pubs.acs.org on July 29, 2018

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 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 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.

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 50 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

Journal of Chemical Theory and Computation

Direct Non-Adiabatic Dynamics by Mixed Quantum-Classical Formalism Connected with Ensemble Density Functional Theory Method: Application to trans-Penta-2,4-Dieniminium Cation Michael Filatov,∗ Seung Kyu Min,∗ and Kwang S. Kim Department of Chemistry, School of Natural Sciences, Ulsan National Institute of Science and Technology (UNIST), Ulsan 44919, Korea E-mail: [email protected]; [email protected]

Abstract In this work, a direct mixed quantum-classical dynamics approach is presented, which combines two new computational methodologies. The nuclear dynamics is solved by the “decoherence-induced surface hopping based on the exact factorization” (DISH-XF) method, which is derived from the exact factorization of the electronicnuclear wavefunction and correctly describes quantum decoherence phenomena. The “state-interaction state-averaged spin-restricted ensemble-referenced Kohn-Sham” (SISA-REKS, or SSR, for brevity) electronic structure method is based on “ensemble den-

1 ACS Paragon Plus Environment

Journal of Chemical Theory and Computation 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

sity functional theory” (eDFT) and provides correct description of real crossings between the ground and excited Born-Oppenheimer electronic states. The new combined approach has been applied to the excited state non-adiabatic dynamics of the trans-penta-2,4-dieniminium cation (PSB3). The predicted S1 lifetime of trans-PSB3, τ = 99 ± 51 fs, and the quantum yield of the cis-conformation, φ = 0.63, agree with the results obtained previously in non-adiabatic molecular dynamics simulations performed with a variety of electronic structure methods and dynamics formalisms. Normal mode analysis of the obtained classical nuclear trajectories suggests that only a few vibrational normal modes contribute to the nuclear wavepacket; where synchronization between several modes plays a dominant role for the outcome of photoisomerization.

1

Introduction

For modeling non-adiabatic molecular dynamics of realistic molecules with multiple nuclear degrees of freedom, it is necessary to calculate accurately and efficiently electronic (or Born-Oppenheimer) ground and excited states including non-adiabatic couplings between them, and to describe coupled electron-nuclear dynamics with the correct account of quantum decoherence phenomena. For the former, it is crucial to describe correct topology of conical intersections, i.e., real crossings between the electronic states of the same space and spin symmetry. 1 As the quantum mechanical methods used for direct dynamics should satisfy this requirement, only a limited repertoire of methods can be used for this purpose. 2 Predominantly these are the methods based on the complete active space self-consistent field (CASSCF) methodology 3 and CASSCF augmented by perturbational description of the dynamic electron correlation, e.g., CASPT2. 4 These method-


Page 2 of 50

Page 3 of 50 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


ologies should be used in the state-averaged (SA) 5 or multi-state (MS) 6 form to provide for the correct description of the double cone topology of conical intersections and this requirement increases further the computational cost. Although SA-CASSCF can treat sufficiently large molecules, especially when ran on the modern hardware, 7 the lack of the dynamic electron correlation (consequently, underestimation of the energies of the chemical bonds) makes predictions obtained with this methodology less reliable. Inclusion of the dynamic correlation, even at the simplest MS-CASPT2 (MSPT2, for brevity) level, immediately increases the computational cost and restricts applicability of the direct non-adiabatic molecular dynamics to small molecules. For the nuclear dynamics, probably the most accurate method would be the full quantum mechanical description of the nuclear motion, such as in the multi-configurational time-dependent Hartree (MCTDH) method 8 or in the variational multi-configurational Gaussian wave packet (vMCG) method. 9 The latter methods require on-the-fly obtaining of the molecular Hessians, which imposes an unbearable burden on the electronic structure methods; hence, they are limited to a relatively small number of nuclear degrees of freedom. The ab initio multiple spawning (AIMS) method 10 employs classical equations of motion to describe evolution of the trajectory basis functions used to expand the nuclear wavepacket. Hence, AIMS simplifies considerably the full quantum description (e.g., the Hessian is no longer needed) and can be applied to bigger molecular systems. 11 Further simplifications of the non-adiabatic quantum molecular dynamics are achieved with the use of mixed quantum-classical approach, 12,13 where the classical nuclear equations of motion (EOM) are combined with the quantum mechanical description of electrons to obtain on-the-fly electronic energies, energy gradients (forces), and non-adiabatic couplings between the electronic states. Among many other proposed approaches, tra-



jectory surface hopping (TSH) method 12,14 and Ehrenfest dynamics (ED) 15 are by far the most popular methods due to the relative ease of obtaining the required information from the electronic structure calculations. However, the TSH and ED methods possess a number of well-known disadvantages which may result in an inaccurate description of the dynamics. In particular, the lack of quantum decoherence in TSH and ED methods often causes problems in describing non-adiabatic phenomena of realistic molecules. Hence, the quest for a computationally economic and quantitatively accurate methodology for non-adiabatic direct dynamics simulations still goes on. In this work, we combine two recent methodologies, the decoherence-induced surface hopping based on exact factorization (DISH-XF) method 16 for describing nuclear dynamics with the state-interaction state-averaged spin-restricted ensemble-referenced Kohn-Sham (SI-SA-REKS, or SSR, for brevity) method 17–21 for obtaining on-the-fly electronic energies, gradients, and non-adiabatic couplings. In the DISH-XF method, the exact factorization 22,23 of a time-dependent molecular wave function Ψ(r, R, t) into a single product of a nuclear wave function χ(R, t) and an electronic wave function ΦR (r, t), i.e., Ψ(r, R, t) = χ(R, t)ΦR (r, t) with electronic (r) and nuclear (R) degrees of freedom, can correctly describe correlation between electrons and nuclei. 24 The inclusion of the electron-nuclear coupling enables one, with the use of multiple nuclear trajectories, to correctly describe quantum decoherence phenomena arising in various non-adiabatic situations. 25 The DISH-XF method exploits the electron-nuclear correlation contribution, derived from the exact factorization, on top of the conventional surface hopping formalism 12 to accurately account for quantum coherence effects. The SSR method, 17–21 used here for obtaining the input information for DISH-XF simulations, employs ensemble density functional theory (eDFT) 26 to accurately take into


Page 4 of 50

Page 5 of 50 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


account the non-dynamic electron correlation, e.g, when the chemical bonds are broken or formed, in the ground and excited electronic states and to describe the excited electronic states within variational, i.e., time-independent formalism. 27 The method includes the dynamic electron correlation through the use of the exchange-correlation DFT functionals thus correctly describing the strength of the chemical bonds in molecules. The SSR method was previously tested for its ability to describe conical intersections between the ground and excited states and it proved to be a very accurate method matching the most sophisticated multi-reference wavefunction theory methods. 18,28,29 In this work, the DISH-XF/SSR direct dynamics approach will be used to simulate the non-adiabatic dynamics of the lowest excited state of the penta-2,4-dieniminium cation, or PSB3, for brevity. PSB3 is often used as a simplified model of the retinal protonated Schiff base chromophore, 2,30–32 the molecule implicated in the process of vision. In the past, the dynamics of the S1 state of PSB3 was investigated 2,31–36 by TSH methods in connection with SA-CASSCF, MSPT2, and MRCIS (multi-reference configurational interaction with single excitations) 37 as well as with AIMS method 10 in connection with MSPT2. 38 There exists a thick layer of literature dedicated to quantum chemical investigation of the ground and excited states potential energy surfaces (PESs) of PSB3. 2,30,39,40 The wealth of data available for this molecule should facilitate comparison with the new direct dynamics approach and let to draw fair conclusions about its perspectives.



Page 6 of 50

2 Computational methods 2.1

DISH-XF formalism

Here only a brief explanation of the DISH-XF method is given; a more detailed description can be found elsewhere. 16 Similar to other mixed quantum-classical approaches, the electronic degrees of freedom are described quantum mechanically, while classical description in terms of trajectories is used for the nuclear degrees of freedom. In DISH-XF approach, using the Born-Oppenheimer (BO) expansion of the electronic wave function (I)

ΦR(I) (r, t) for a given nuclear trajectory R( I ) , i.e., ΦR(I) (r, t) = ∑l Cl (t)φl (r; R( I ) (t)), the (I)

( I )∗

time evolution of an element ρlk (t) = Cl

(I)

(t)Ck (t) of the reduced density matrix ρ( I ) (t)

is governed by o n o i n (I) d (I) (I) (I) (I) (I) (I) (I) ρlk (t) = El (t) − Ek (t) ρlk (t) − ∑ σlj (t)ρ jk (t) − ρlj (t)σjk (t) dt h¯ j n o (I) (I) (I) (I) + ∑ Q jl (t) + Q jk (t) ρlj (t)ρ jk (t)

(1)

j

(I)

where El

and φl (r; R( I ) (t)) are the energy eigenvalue and the wavefunction of the l(I)

th BO state, σjk is a non-adiabatic coupling matrix element between the jth and kth BO (I)

electronic states, and Q jk is an additional electron-nuclear coupling term which describes quantum corrections to nuclear motion, i.e. the effect of spatial distribution of nuclear density beyond the δ-function. In the above equation, the superscript ( I ) indicates a quantity obtained at a nuclear configuration R( I ) . The first and second terms on the right-hand side of Eq. (1) are the same as in the conventional ED or TSH methods, whereas the last term describes the nuclear quantum (I)

effects and is derived from the exact factorization. More specifically Q jk is a coupling 6 ACS Paragon Plus Environment

Page 7 of 50 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


(I)

(I)

between nuclear quantum momenta i¯h∇ν |χ|/|χ| and electronic phases fν,j and fν,k as (I)

Q jk =

∑ ν

i¯h ∇ν |χ| Mν | χ |

(I) (I) · fν,j − fν,k

(2)

R( I ) (t)

where Mν is a mass of the νth nucleus. In DISH-XF method, the nuclear quantum momenta are obtained from multiple auxiliary nuclear trajectories for electronic states other than the running state l. An auxilliarty trajectory Rk (t0 ) = R( I ) (t0 ) is generated when a nonzero ρkk (t0 ) is encountered at a time t0 . The auxilliary trajectory evolves classically with a uniform velocity, which is obtained from the energy conservation law. The elecRt (I) (I) tronic phase term fν,l = − ∇ν Ek (t0 )dt0 is evaluated by time integration of momentum changes at the given BO state l. To calculate the quantum momentum, a fictituous Gaussian nuclear density |χk |2 with a uniform variance σ is associated with the auxilliary trajectory at each BO state. Then, the nuclear quantum momentum is ∇ν |χ|/|χ|(R(I) (t)) = (I)

(I)

(I)

(I) (I)

− 2σ1 2 ( Rν (t) − h Rν (t)i), where h Rν (t)i = ∑k ρkk Rkν (t). The uniform variance σ can be either obtained from the initial distribution of nuclear trajectories or set as a parameter. In this paper, a uniform value σ = 0.1 a.u. is used. A more detailed description of the DISH-XF algorithm can be found elsewhere. 16 When conducting electronic time evolution according to Eq. (1), a stochastic hopping probability from the running state l to another state k (6= l ) at a time interval [t, t + ∆t] is calculated as 2< Pl →k =

h

(I) (I) ρlk (t)σlk (t) (I) ρll (t)

i ∆t .

(I)

(3)

The phase factor of the ρlk (t) matrix element may result in negative hopping probabil7 ACS Paragon Plus Environment


Page 8 of 50

(I)

ity; in such a case, Pl →k is set to zero. 12 In addition, if the kth BO energy, Ek , is greater than the total energy, the hopping transition is forbidden, which is identical to Tully’s fewest-switches algorithm. 12 After a successful hopping event, the nuclear velocities are rescaled to satisfy the total energy conservation. A nuclear trajectory follows Newtonian equations of motion on the running state l potential energy surface, i.e., the force (I)

acting on the νth nucleus, Fν , is given by Fν

(I)

= −∇ν El . To summarize, the DISH-XF

algorithm combines the electronic equations derived from the exact factorization of the electronic-nuclear wavefunction with the conventional surface hopping formalism. The exact factorization enables one to seamlessly incorporate the effect of nuclear quantum momentum, which depends on the shape of nuclear distribution, into the classical equations of motion. In the DISH-XF algorithm, the nuclear quantum momenta are obtained from a set of auxilliary nuclear trajectories, which are propagated together with the running trajectory. The concept of auxilliary trajecotries for decoherence-induced surface hopping dynamics is not new, see e.g., Ref. 41. However, the decoherence correction employed in the DISH-XF algorithm is derived from the exact quantum equations and no further renormalization of the electron density matrix is neccesary. In this work, the nuclear equations of motion are integrated using the velocity-Verlet (I)

algorithm with the time step of 20 a.u. (0.48 fs). The electronic energies El , forces on the (I)

(I)

nuclei −∇ν El , and the non-adiabatic couplings σlk are calculated quantum mechanically (see below) at the end points of the integration intervals. Within these intervals, the electronic equations of motion (1) are integrated by the 4-th order Runge-Kutta method with the time step of 0.002 a.u. (4.8 × 10−5 fs). When integrating the electronic equations (I)

of motion, the electronic energies El

(I)

and the non-adiabatic σlk couplings are linearly

interpolated between the end points of the integration interval of the nuclear equations


Page 9 of 50 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


of motion.

2.2

REKS methodology

The SSR method employed in this work to obtain on-the-fly the ground and excited electronic state energies, forces on the nuclei, and the non-adiabatic couplings has been previously described in the literature 17–20,42 and only a brief account of its most important features will be given here. The SSR methodology employs ground state eDFT 26 to describe the non-dynamic electron correlation stemming from multireference character of the ground state and eDFT for ensembles of ground and excited states 27 to obtain excitation energies from a variational time-independent formalism. The use of ensemble representation for the density and the energy of a strongly correlated molecular system leads to the occurrence of the fractional occupation numbers (FONs) of several frontier Kohn-Sham (KS) orbitals. In this work, the SSR method including two fractionally occupied KS orbitals accommodating two electrons, i.e., SSR(2,2), is employed. 19,20,42 In the SSR(2,2) method, the energies of the ground and excited states are obtained from variational minimization (with respect to the KS orbitals and their FONs) of an ensemble of the REKS(2,2) perfectly spin-paired singlet (PPS) electronic configuration and an open shell singlet (OSS) configuration obtained within the same (2,2) active space 17–20,43 followed by solving a simple 2×2 secular problem 

PPS  E0



SA ∆01

SA ∆01











  a00 a01   E0 0   a00 a01   =   OSS E1 a10 a11 0 E1 a10 a11

(4)

to include possible coupling between the PPS and the OSS electronic configurations. In Eq. (4), E0PPS and E1OSS are the energies obtained in the SA-REKS(2,2) orbital optimiza9 ACS Paragon Plus Environment


Page 10 of 50

SA is calculated using the SA-REKS(2,2) Lation and the interstate coupling parameter ∆01

grangian matrix element εSA ab between the active orbitals φa and φb as

√ √ SA ∆01 = ( n a − nb ) εSA ab

(5)

where n a and nb are the FONs of the active orbitals. 18,20 The SA-REKS(2,2) orbitals are obtained by minimizing an equiensemble of the E0PPS and E1OSS energies under the constraint of orbital orthonormality; simultaneously, the FONs of the active orbitals are optimized. 17–20 The gradients of the SSR(2,2) energies E0 and E1 are obtained using the analytical derivatives formalism presented in Ref. 21; derivation of the formalism is quite lengthy and is not repeated here. The gradients of the (adiabatic) SSR(2,2) states ∇ Ek , with respect to nuclear coordinates, are related to the gradients of the (diabatic) SA-REKS(2,2) states SA . 21 ∇ E0PPS and ∇ E1OSS and of the interstate coupling element ∇∆01

SA ; ∇ Ek = a2kk ∇ E0PPS + a2lk ∇ E1OSS + 2 akk alk ∇∆01

l 6= k

(6)

Using the gradients of the SSR(2,2) states ∇ Ek and of the interstate coupling element SA , the non-adiabatic coupling vector ∇∆01

H01 = hS0 |∇S1 i

=

1 (( a00 a01 − a10 a11 )g01 + ( a00 a11 + a01 a10 )h01 ) E1 − E0


(7)

Page 11 of 50 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


can be expressed as

H01

2 ( a00 a01 − a10 a11 ) 2 ( a00 a01 − a10 a11 )( a00 a10 − a01 a11 ) G01 − h01 (8) 2 2 2 2 a00 − a01 − a10 + a11 a200 − a201 − a210 + a211 + ( a00 a11 + a01 a10 ) h01 , 1 = E1 − E0

where

G01 =

1 (∇ E0 − ∇ E1 ) 2

h01 = ∇∆SA g01 =

1 (∇ E0PPS − ∇ E1OSS ) . 2

(9) (10) (11)

The gradients of the SSR(2,2) states ∇ Ek and the interstate coupling gradient ∇∆SA are R program 44 (v1.92P, release calculated using the beta-testing version of the TeraChem

7f19a3bb8334). The gradients are picked up by an external script, where the non-adiabatic coupling vector is calculated by Eq. (8), and together with the SSR energies Ek parsed to the UNI-xMD program, a standalone code which implements the DISH-XF method. 16

3 Results and discussion 3.1

Characterization of stationary points on S1 ans S0 PESs

Before running the molecular dynamics simulations, the S0 and S1 state stationary points of PSB3 were studied using the SSR(2,2) method in connection with the ωPBEh range separated hybrid XC functional 45 and the 6-31G* basis set, 46 i.e., SSR-ωPBEh/6-31G* method. The geometries of the S0 and S1 minima as well as the S1 /S0 minimal energy conical in11 ACS Paragon Plus Environment


Page 12 of 50

tersections (MECIs) obtained in this work are compared with the results of a recent investigation of the low-lying stationary points on the S0 and S1 PESs of PSB3 by the CASSCF and MSPT2 methods. 38 Table 1: Relative energies (kcal/mol) of stationary point geometries of PSB3 optimized by SSR-ωPBEh/6-31G* (this work), SA3-MS-CASPT2(6,6)/6-31G* (Ref. 38), and SA3CASSCF(6,6)/6-31G (Ref. 38) methods.

a) b)

geometrya

SSR-ωPBEh/6-31G*

SA3-MS-CASPT2(6,6)/6-31G*

SA3-CASSCF(6,6)/6-31G

trans-S0 cis-S0

0.0/95.3b 2.9/96.3b

0.0/98.3b 3.0/97.2b

0.0/113.2b 3.0/113.0b

S1min CenL S1min Cen S1min CN S1min Ter

– 64.0 95.4 74.7

85.1 – – –

95.2 71.3 77.7 104.0

MECICen MECITer MECICenR MECICN

66.7 78.0 98.5 98.0

60.3 82.7 – 87.8

66.1 108.9 99.4 79.1

See text for abbreviations CenL, Cen, CN, Ter, and CenR. Vertical excitation energy of S1 state relative to trans-S0 .

Table 1 lists the relative energies of the stationary points of PSB3 obtained with the SSR-ωPBEh/6-31G* method and compares them with the results of Ref. 38. The geometries of the energy minima were optimized using the DL-FIND module 47 interfaced with R TeraChem and the MECI geometries were obtained by the CIOpt program 48 using the

penalty function formalism with the analytic energy gradients of the intersecting states. The geometries obtained in Ref. 38 were used to initiate the SSR-ωPBEh/6-31G* geometry optimizations. As seen in Figure 1, where the superimposed geometries of the S0 minima obtained in 12 ACS Paragon Plus Environment

Page 13 of 50 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


1.305/1.334

1.305/1.318 1.404/1.406

1.404/1.421

1.372/1.388

1.369/1.393

1.434/1.439 1.431/1.450 1.345/1.364

1.346/1.374

cis-PSB3

trans-PSB3

Figure 1: Geometries of the S0 minima of PSB3 obtained with the SSR-ωPBEh/6-31G* (conventional color scheme; this work) and SA3-MS-CASPT2(6,6)/6-31G* (red; from Ref. 38) methods. The key bondlength are given: SSR – normal font, Ref. 38 – italicized font. this work and in Ref. 38 are shown, there is a good agreement between the geometries of the trans- and cis-conformations optimized with SSR and MSPT2. The root-mean-square deviation (RMSD) of the atomic coordinates between the SSR and MSPT2 geometries is ˚ for trans-PSB3 and 0.014 A ˚ for cis-PSB3. The SSR-ωPBEh/6-31G* method pre0.058 A dicts somewhat greater bond length alternation (BLA, the difference between the average length of formally single bonds and of formally double bonds); BLA of trans-PSB3 and cis˚ for both minima, while it is 0.069 A ˚ for trans-PSB3 and 0.066 A ˚ PSB3 from SSR is 0.078 A for cis-PSB3 as obtained by MSPT2. Hence, the DFT method predicts a somewhat greater disparity between the single and double bonds in the ground state species of PSB3. The cis–trans energy difference and the vertical excitation energies of the S0 conformations obtained by SSR are in good agreement with the MSPT2 values, see Table 1. The vertical excitation energy of cis-PSB3 obtained by SSR, 96.3 kcal/mol, is also in good 13 ACS Paragon Plus Environment


Page 14 of 50

agreement with the energy of 96.8 kcal/mol from the quantum Monte-Carlo (QMC) calculations of Valsson et al. 49 These comparisons demonstrate that the SSR method provides a fair account of the dynamic correlation effects at the Franck-Condon points of PSB3.

1.311/1.350

1.361/1.437 1.403/1.366

1.441/1.446

1.415/1.355

1.324/1.352 1.385/1.359 1.401/1.437

1.379/1.419 1.408/1.401 1.363/1.372

S1min_Cen

1.429/1.416

1.385/1.353

1.348/1.374

S1min_CN

1.405/1.416

S1min_ter

Figure 2: Geometries of the S1 minima of PSB3 obtained with the SSR-ωPBEh/6-31G* (this work) and SA3-CASSCF(6,6)/6-31G (from Ref. 38) methods. See caption to Fig. 1 for legend. Of the four S1 minima of PSB3 identified in Ref. 38 with the use of SA-CASSCF calculations, three (S1min ter , S1min Cen , and S1min CN ) were found using the SSR-ωPBEh/6-31G* optimizations, see Figure 2. These minima correspond to torsion (dihedral angle ca. ±90◦ ) about the terminal C=C bond, the central C=C bond, and the terminal C=N bond of PSB3, respectively. As CASSCF misses dynamic electron correlation, the geometry obtained with this method in Ref. 38 noticeably deviates from the geometry obtained using SSR. In agreement with earlier works, 2,50 CASSCF predicts the S1min Cen minimum at the BLA, ˚ that corresponds to inversion of the single and double bond lengths in PSB3. -0.012 A, The dynamic correlation strongly modifies the shape of the S1 PES and brings BLA to ˚ with the SSR-ωPBEh/6-31G* method. the positive range of values; 2,40,50–52 BLA = 0.040 A ˚ and A similar trend is observed for the S1min CN minimum with BLACASSCF = -0.019 A ˚ BLASSR = 0.053 A. 14 ACS Paragon Plus Environment

Page 15 of 50 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


Interestingly, the MSPT2 method of Ref. 38 did not yield the same S1 minima as the SSR method in this work. The only S1 minimum identified in Ref. 38 by MSPT2 was the S1min CenL minimum which corresponds to simultaneous twisting of the central and terminal C=C bonds. Presumably, the S1 minima were either very close to the respective MECI points (see below), or even coincided with the respective MECIs; in both cases it would be impossible to optimize their geometries. However, if one assumes that the MSPT2 energy of, e.g., MECICen , 60.3 kcal/mol (see Table 1) is sufficiently close to the energy of the missing minimum, the SSR-ωPBEh/631G* S1min Cen relative energy of 64.0 kcal/mol is in reasonable agreement with this value. According to SSR, the S1min CN minimum lies ca. 30 kcal/mol higher than S1min Cen and at the same energy as the FC point. By contrast, CASSCF, which misses dynamic correlation, predicts S1min CN very close to S1min Cen and well below the FC point, see Table 1. Figure 3 displays geometries of the S1 /S0 MECI points, where the geometries from Ref. 38 (shown in red) are superimposed with the geometries obtained in this work (shown in conventional color scheme). With the exception of MECICenR , for which only the SA3CASSCF(6,6)/6-31G geometry was available in Ref. 38, all other reference geometries were obtained by MSPT2. The dynamic electron correlation has considerable implications for the relative stability of the MECIs and, hence, for the S1 state dynamics of PSB3. The two lowest MECI points, see Table 1, predicted by SSR and MSPT2 are MECICen and MECITer , obtained by torsion of the central C=C bond and the terminal C=C bond, respectively. By contrast, CASSCF predicts MECICN obtained by torsion of the C=N bond as the second lowest MECI. This implies that, although the major reaction channel as predicted by SSR and MSPT2 on one hand and by CASSCF on the other hand is the same (torsion about the



1.292/1.305

Page 16 of 50

1.330/1.343 1.372/1.378

1.457/1.439 1.446/1.462

1.409/1.417 1.357/1.356

1.412/1.391

1.375/1.390

1.358/1.396

MECITer

MECICen

1.390/1.451 1.379/(1.412)

1.452/(1.478)

1.451/1.480

1.379/1.384

1.377/(1.401) 1.425/1.454

1.426/(1.397)

1.354/1.371

1.354/(1.389)

MECICN

MECICenR

Figure 3: Geometries of the S1 /S0 MECIs of PSB3 obtained with the SSR-ωPBEh/631G* (this work) and SA3-MS-CASPT2(6,6)/6-31G* (from Ref. 38) methods; the SA3CASSCF(6,6)/6-31G (from Ref. 38) geometry is reported for MECICenR . See caption to Fig. 1 for legend. The SA3-CASSCF bondlengths are given by parenthesized italic font. central C=C bond), dynamic correlation steers the minor reaction channel in the direction of torsion of the terminal C=C bond and not of the C=N bond. 38 16 ACS Paragon Plus Environment

Page 17 of 50 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


There is a significant difference between the MECICN geometries predicted by the MSPT2 and the CASSCF methods. 38 The MSPT2 geometry features torsion about two bonds, the terminal N1 =C2 bond and the neighboring C2 –C3 single bond. The CASSCF MECICN is obtained by torsion about the terminal N1 =C2 bond alone. The SSR method yields both geometries for this CI: When started from the MSPT2 geometry, the SSR optimization converges to MECICN featuring torsion about the two bonds and lying 98.0 kcal/mol above the S0 trans-PSB3 minimum, see Table 1 and Figure 3. If the CASSCF geometry is used to start the SSR optimization, a crossing point with only N1 =C2 bond torsion is obtained (not shown in Figure 3). The latter crossing point, which is likely a local minimum on the CI seam, lies 18.2 kcal/mol above the SSR MECICN point, see Table 1. Compared to the MSPT2 MECICN , which lies ca. 10 kcal/mol below the FC point, the SSR MECICN is located slightly, by ca. 3 kcal/mol, above the respective FC point. There is a reasonable agreement between SSR and MSPT2 for the geometry of the most mechanistically significant MECI point, MECICen , see Figure 3. However, MSPT2 predicts that MECICen occurs at the bottom of the S1 global minimum, while SSR predicts that MECICen lies 2.7 kcal/mol above S1min Cen . The major difference between S1min Cen and ˚ and BLA(MECICen ) MECICen geometries is in the BLA distortion; BLA(S1min Cen ) = 0.040 A ˚ = 0.069 A. The shape of the S1 and S0 PESs of PSB3 near MECICen as predicted by the SSRωPBEh/6-31G* calculations is generally consistent with the recent results obtained by high level multi-reference methods. 2,40 This is illustrated in Fig. 4 where the S1 and S0 relative energies (with respect to the ground state of cis-PSB3) obtained in Refs. 2 and 40 and in this work are shown along the BLA path constructed previously by Gozem et al. by interpolating between the transition states for the homolytic and heterolytic breaking



80

ΔE, kcal/mol

70

60

50 -0.04

-0.02

0.00

0.02

0.04

0.06

0.08

0.10

BLA, A o

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 18 of 50

Figure 4: S1 and S0 relative energies (kcal/mol, with respect to the cis-PSB3 S0 energy) ˚ path from Ref. 51. The geometries along the path and the geometry along the BLA (A) of the cis-PSB3 conformation are taken from Ref. 51. The SSR-ωPBEh/6-31G* energies are given by blue curves; the MRSDCI (green) and LRDMC (red) energies are reproduced from Fig. 5 of Ref. 40. (Copyright 2015, American Chemical Society) of the central C=C bond of PSB3 obtained in the CASSCF calculations. 51 The BLA path intercepts the MECICen obtained at the same CASSCF level and intercepts the crossing seam at the MRSDCI level. 51 The BLA path of Gozem et al. is standardly used to benchmark the ability of various computational methods to correctly describe the conical intersection of PSB3 obtained by torsion of the central C=C bond and the shape of the S1 and S0 PESs in its vicinity. 2,40,50–52 Fig. 4 compares the SSR-ωPBEh/6-31G* energies with the MRSDCI and the LRDMC (lattice regularized diffusion Monte Carlo) energies from Ref. 40; the purpose of this comparison is to evaluate the quality of the PES shape description provided by the SSR method. Along the BLA path, the S1 /S0 crossing point predicted by SSR occurs at the BLA ˚ which is in good agreement with 0.061 A ˚ predicted by the LRDMC value of ca. 0.055 A, 18 ACS Paragon Plus Environment

Page 19 of 50 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


calculations; both methods predict sloped topography of the crossing point. MRSDCI ˚ and the crossing point has almost peaked predicts crossing point at the BLA of ca. 0.028 A topography, see Fig. 4. In this regard, the MRSDCI calculations of Ref. 40 and the MSPT2 calculations of Ref. 38 are consistent with each other; the latter method predicts MECICen ˚ However, at the bottom of the S1 minimum (i.e., peaked topography) at a BLA of 0.027 A. a more complete account of the dynamic electron correlation provided by the LRDMC methodology shifts the intersection point away from the minimum on the S1 curve and towards larger BLA values; this agrees with the picture predicted by SSR. Judging from the relative energy of the MECIs reported in Table 1 the non-adiabatic relaxation of the S1 state of trans-PSB3 most likely proceeds through MECICen and/or MECITer ; thus leading to isomerization about the central C=C bond or terminal C=C bond. To further assess the accessibility of these MECIs from the FC point of trans-PSB3 the minimum energy pathways (MEPs) connecting the FC point with the two MECIs were optimized using the nudged elastic band (NEB) method 53 with fixed endpoints, as implemented in the DL-FIND module; 47 for each MEP 40 points were optimized. The optimized FC→MECICen and FC→MECIter MEPs are shown in Fig. 5, where they are plotted as functions of the respective torsion angle and the BLA parameter. The initial motion along both MEPs is along the BLA coordinate and a noticeable displacement ˚ to ca. 0.010 A, ˚ see along the torsion angle begins only after the BLA drops from 0.078 A Fig. 5. This picture is generally consistent with the previous investigations of the PSB3 isomerization, where the initial motion along the MEP corresponded to BLA displacement. 51 At the initial stage of torsion, the FC→MECIter MEP is noticeably shallower than the FC→MECICen MEP; there is a shallow, ca. 1 kcal/mol minimum near the FC point along



trans-PSB3 MECIter

MECICen

Figure 5: S1 energy along the MEPs connecting the trans-PSB3 FC point with the MECICen (blue line) and MECITer (red) points. The MEPs are drawn as functions of the respective ˚ parameter. torsion angle (deg.) and the BLA (A) the former MEP. It should be noted that the FC point does not correspond to a minimum on the S1 PES of PSB3. Geometry optimization of S1 trans-PSB3 starting from the FC geometry leads to a minimum ca. 4 kcal/mol below the FC point. Hence, if started from this minimum (not from FC), the two MEPs would encounter low barriers on the way to the MECI points; this picture is typically reported for the PSB3 MEPs in the literature. 51 However, starting the MEP from the FC point, rather than the planar S1 minimum, closely resembles an MD trajectory with zero nuclear kinetic energy; hence providing an insight into possible behavior of real (i.e., non-zero kinetic energy) trajectories. As can be surmised from the MEPs in Fig. 5, isomerization about the terminal C=C bond is slightly hindered compared to the central C=C bond; only a minor fraction of all nuclear trajectories can be expected to follow this isomerization path. Hence, it is not only the relative energy of the MECIs, but also the shape of the MEPs along the respective isomerization pathways that suggests splitting between the major (central C=C bond) and minor (terminal C=C bond) excited state reaction channels. 20 ACS Paragon Plus Environment

Page 20 of 50

Page 21 of 50 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


3.2

Non-adiabatic dynamics of PSB3

Non-adiabatic molecular dynamics (NAMD) of the S1 excited state of PSB3 was studied using the DISH-XF/SSR-ωPBEh/6-31G* method described in Section 2. In total, fifty trajectories were propagated. At the start of the trajectories the molecule was brought to R the S1 state at the geometries generated by the TeraChem code by sampling the Wigner

function of a canonical ensemble at T = 300K; 54,55 the Wigner function was calculated at the SSR-ωPBEh/6-31G* S0 trans-PSB3 equilibrium geometry. The vibrational frequencies of trans-PSB3 were obtained by numeric differentiation of the analytic SSR-ωPBEh/631G* gradient, see Section 2.2. The trajectories were propagated using the forces on nuclei and the non-adiabatic coupling vectors obtained on-the-fly by the SSR-ωPBEh/6-31G* method. The time step between the quantum chemical calculations was set to 0.24 fs and the trajectories were propagated up to 1250 steps (300 fs). The propagation employed the NVE (i.e., micro-canonical) ensemble and the total energy (i.e., the electronic plus the nuclear kinetic energy) was conserved with the fluctuations on the order of 0.04 kcal/mol, see the Supporting Information for sample trajectories. Out of fifty trajectories, 9 trajectories underwent the non-adiabatic S1 → S0 transition through torsion about the terminal C5 =C6 bond, and the remaining 41 trajectories underwent the transition by torsion about the central C3 =C4 bond, see Fig. 6 for the atomic numbering and definition of the torsion angle and BLA parameter. Of the latter trajectories, 63% were isomerized to the cis-PSB3 conformation and 37% fell back to the trans conformation. Hence, as suggested by the analysis of the S1 MEPs in Section 3.1, the major reaction channel is isomerization about the central double bond (probability 0.82) and the minor channel (0.18) is isomerization about the terminal C=C bond. The obtained splitting be21 ACS Paragon Plus Environment


Page 22 of 50

C5

CCW

CW 6

θ

5

C4

H 4

C3 H

3

C2

2 1

BLA=(R23+R45)/2 -(R12+R34+R56)/3

Figure 6: Atomic numbering and definition of the BLA parameter, the direction of the torsion, and the central torsion angle used in this work. The central torsion angle θ is defined as a dihedral angle between the normals to the planes spanning the C2 -C3 -H and C5 -C4 -H planes. tween the major and minor reaction channels is consistent with the results of Ref. 38, where the splitting 0.89/0.11 was obtained in the AIMS/MSPT2 simulations. Although in this work the initial conditions for the trajectories were generated by the Wigner function of a canonical ensemble at T = 300K, which provides a broader distribution of the initial geometries and velocities than the zero-temperature function, no other reaction side channels, e.g., the C=N isomerization, were observed. Earlier TSH/CASSCF and TSH/MRCIS NAMD simulations of the PSB3 E → Z photoisomerization 32,33 as well as the AIMS/CASSCF simulations 38 yielded the C=N isomerization as a side channel; ca. 19–22% of the trajectories were diverted to this channel. The difference in the side channel of the reaction, terminal C=C vs. C=N, is caused by the effect of dynamic electron


Page 23 of 50 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


correlation missing in the CASSCF and MRCIS calculations. 38 The quantum yield of isomerization about the central C=C bond obtained here (0.63) is somewhat lower than the quantum yield obtained in the AIMS/MSPT2 simulations (0.79) 38 and is closer to the AIMS/CASSCF result (0.54) 38 as well as the TSH/CASSCF result (0.57). 32,33 Interestingly, the AIMS/MSPT2 quantum yield agrees well with φ = 0.72 obtained by Gozem et al. 2 in the TSH simulations using the single state CASPT2 method on a PSB3 derivative where a pretwist was applied to the central C=C bond.

Figure 7: Histogram showing the number of S1 /S0 hops as a function of time. The blue curve shows bi-Gaussian fitting of the obtained hoptimes. The average S1 /S0 hoptime for the major reaction channel obtained in this work, 99±51 fs, is in reasonable agreement with the S1 lifetimes obtained in the previous NAMD simulations and varying between ca. 100–150 fs. 31–36 The distribution of the obtained S1 /S0 hoptimes is shown in Fig. 7; the obtained distribution can be best fitted by a bi-



Page 24 of 50

Gaussian function " # # ( t − t1 )2 ( t − t2 )2 c2 F (t) = √ Exp − Exp − +√ ; 2σ12 2σ22 2π σ1 2π σ2 c1

"

c1 + c2 = 1

(12)

with the parameters shown in Table 2. Table 2: Parameters of bi-Gaussian distribution, Eq. (12), of the S1 /S0 hoptimes for isomerization of trans-PSB3 about the central C=C bond. The parameters of the overall S1 population evolution and of populations evolving to cis-conformation and to transconformation are shown. c1 all 0.35 trans→cis 0.53 trans→trans 0.15

t1 , fs 68.7 72.2 57.5

σ1 , fs t2 , fs 14.6 106.9 11.4 148.9 4.5 114.5

σ2 , fs 64.7 64.8 34.6

The observed S1 population dynamics displays pronounced splitting into fast and slow components. The forward going (trans→cis) trajectories show a somewhat broader splitting between the fast (t1 = 72.2 fs) and slow (t2 = 148.9 fs) trajectories than the backward going (trans→trans) trajectories (t1 = 57.5 fs, t2 = 114.5 fs), respectively. The obtained distributions of the hop times are in a very good agreement with the electronic (I)

coherences evaluated from the off-diagonal elements ρ01 (t) of the electronic density matrix (1); see the electronic Supporting Information for plots. This underlines the internal consistency of the obtained dynamics. The majority of fast trajectories quickly slide along the S1 PES towards the CI seam where they undergo surface hop to the S0 state. Evolution of the slow trajectories is more complex; some trajectories are slow starters and spend some time near the FC point, others propagate rather quickly towards CI and become trapped within a shallow S1min Cen minimum. Both types of the dynamics were previously reported in the literature; Nikiforov et al. 56 observed both slow start and trapping in an S1 minimum when studying the 24 ACS Paragon Plus Environment

Page 25 of 50

dynamics of light driven molecular motors derived from PSB3, whereas Weingart et al. 57

390

0.30

360

0.25

330

0.20

300

0.15

270

0.10

240

0.05

210

0.00

180

-0.05

150

0

50

100

150

t, fs

200

250

o

θ, deg.

observed trapping in the dynamics of cis-azobenzene.

BLA, A

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


-0.10 300

˚ dotted Figure 8: Evolution of the torsion angle θ (deg., solid lines) and BLA parameter (A, lines) with time (fs) along the averaged trajectories propagating in the forward (red) and backward (blue) directions. The shaded area shows the time interval within which the majority of surface hops occur. Isomerization about the central C=C bond proceeds with equal probability in both directions, CW and CCW, see Fig. 6 for definition. In view of the equivalence between the directions of isomerization, the further analysis is confined to the trajectories propagating in the CW direction of torsion. Fig. 8 shows evolution of the torsion angle θ and the BLA parameter averaged over all the trajectories propagating in the CW direction. At each time step along the trajectories, the instantaneous geometries were averaged over all individual trajectories propagating in the given sense of torsion. In Fig. 8, the forward going and the backward going trajectories are shown separately. As seen from Fig. 8 the initial motion along the trajectories corresponds to the BLA displacement which takes place within first ca. 10 fs. Then, the BLA displacement undergoes an undulatory change and for both forward going and backward going trajectories 25 ACS Paragon Plus Environment


BLA shows very little difference before entering the strong non-adiabatic coupling region marked in Fig. 8 by shaded area. Within the shaded area the BLA distortions of the forward and backward going trajectories slightly depart from each other, however follow similar trends. This suggests that, most likely, the BLA differences are not responsible for the trajectories taking a particular route, forward or backward. The dihedral angles in Fig. 8 vary smoothly from their start to the finish. At first sight it may seem as though the molecule does not reach the CI seam, which occurs at the torsion angles around 270◦ for the CW torsion, especially in the case of the backward going trajectories. However, the averaging is performed over the trajectories for which the hoptimes are spread across a wide interval such that, at any instance of time, there can be only a few trajectories with the torsion angle near 270◦ while the other trajectories have a smaller torsion angle. Before reaching the shaded area, where the majority of surface hops occur, the dihedral angles of the forward and backward trajectories follow similar trends and begin to depart from each other only after ca. 60 fs. At that moment the first surface hops begin to take place; the very first one at 40 fs. Again the small difference between the torsion angles can hardly suggest a clue as to what might have driven some trajectories in the forward direction and others to turn backward. Although here we inspected trajectories averaged over a pool of trajectories propagating in the CW sense, an inspection of the individual trajectories also does not let to single out a factor or two which would be common for all the trajectories and could explain the propensity to complete the isomerization or to abort it and return to the reactant conformation. In the literature, it is often speculated that specific internal molecular degrees of freedom can explain the observed isomerization quantum yield; specifically, the hydrogen out-of-plane (HOOP) vibrational motion around the central C=C double bond


Page 26 of 50

Page 27 of 50 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


of PSB3 (or more complete retinal models) can modulate the direction of propagation of the reaction. 34–36 In the following, we attempt to inspect the vibrational normal modes involved in the displacement of atoms along the NAMD trajectories and to investigate whether certain modes can be responsible for the outcome of the reaction.

3.3

Normal mode analysis of NAMD trajectories

Let R(t) be the mass-weighted coordinates of atoms in a molecule at time t along some trajectory. Let R0 be a geometry at which the normal coordinates Qk , k ∈ [1, Nnmod ] are known. R0 can be, e.g., the geometry (in mass-weighted coordinates) of the ground state minimum used to start the trajectory, or it can be a local minimum on the excited state PES for which the normal modes are known. As the normal coordinates describe all possible internal displacements of atoms in the molecule, the coordinates of atoms along the trajectory can be represented as Nnmod

R ( t ) = R0 +

∑

ak (t) Qk ,

(13)

k

where ak (t) are the time-dependent coefficients. The normal mode amplitudes ak (t) can be found from orthogonality of the normal modes. ak (t) = Q†k · (R(t) − R0 )

(14)

Typically, the trajectories are initiated at the geometries and nuclear momenta obtained by sampling the Wigner function of the ground state minimum. 58 For multiple trajectories, the average values of the normal mode amplitudes h ak (t)i can be straightforwardly found



from 1 h ak (t)i = Ntraj

Page 28 of 50

Ntraj

∑

Q†k · (R j (t) − R0 )

(15)

j

where Ntraj is the number of trajectories. The outlined procedure represents the well-known normal mode analysis (NMA) applied to NAMD trajectories. Currently NMA is gaining popularity in studying large scale motions of macromolecules, 59 however its application to large scale movements of small molecules is deemed sufficiently fruitful for understanding the role of specific vibrational modes for the dynamics. The amplitudes h ak (t)i provide quantitative information on the involvement of specific normal modes into molecular rearrangement and can potentially serve as useful descriptors of the dynamics. For this purpose however the central structure R0 should be chosen in such a way that its normal modes sensibly describe the dynamical rearrangements. In this regard, the normal modes of an excited state minimum, e.g., S1min Cen , occurring near the relevant CI seam, seem better suited for the normal mode analysis than the modes of the ground state equilibrium conformation. The amplitudes of the normal modes h ak (t)i averaged over all the trajectories, Eq. (15), are shown in Fig. 9; the upper panel shows the amplitudes of the normal modes of the ground state trans-PSB3 equilibrium conformation and the lower panel – the normal modes of the S1min Cen minimum. To obtain the latter the trajectories propagating in the CW direction of torsion were taken apart from the ones propagating in the CCW sense, and their amplitudes were calculated using the normal modes of the respective S1 minima. At the start of the trajectories, the averaged amplitudes of the trans-S0 normal modes are clustering around zero value. This is expected behavior, as it indicates that the Wigner function was probed randomly around the central geometry and shows no bias towards 28 ACS Paragon Plus Environment

Page 29 of 50 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


ak2

t, fs

ν,cm-1

ak2

t, fs

ν,cm-1

Figure 9: Squares of the amplitudes h ak (t)i, Eq. (15), in the decomposition (13) of the NAMD trajectories as functions of vibrational frequencies and propagation time. The amplitudes are averaged over all the trajectories. The upper panel shows decomposition in terms of the normal modes of the trans-PSB3 conformation; the lower panel – in terms of the normal modes of the S1min Cen minimum. another conformation. An important observation from Fig. 9 is that only a handful of normal modes during the trajectory propagation time have amplitudes noticeably departing from zero; this implies that only a few normal modes comprise the nuclear wavepacket propagating in time. Turning to the normal modes of the S1min Cen minimum it is seen that only the low frequency skeletal bending modes (in the range below 500 cm−1 ) and the mode #25 (at 1560 cm−1 ) make considerable contributions to the trajectory displacements; contributions of all other modes average to near zero values. The latter mode cor29 ACS Paragon Plus Environment


Page 30 of 50

responds to the HOOP motion and, as was speculated previously in the literature, 31,34–36 this mode should play a decisive role for the outcome of the rearrangement. The HOOP mode (#25, ν25 = 1560 cm−1 ) and the normal mode #3 (simultaneous torsion of single bonds, ν3 = 280 cm−1 ) are the major contributors to the branching plane (BP) vector h, Eq. (10), as shown in Fig. 10. The BP g-vector, Eq. (11), describes predominantly BLA distortion and comprises higher frequency skeletal stretching modes, see middle panel of Fig. 10. Inspection of the projections h · R av (t) and g · R av (t) of the BP vectors on the averaged trajectory R av (t) (see lower panel of Fig. 10) reveals that there is no substantial difference between the projections of the g-vector (usually associated with the tuning mode) of the forward going and backward going reactions; only after most of the hops completed (shown by the shaded area in Fig. 10) the projections of the forward (red dotted curve) and the backward (blue dotted curve) trajectories depart from each other. The h-vector (coupling mode) shows different trend; the projections of the hvector on the forward going (red solid curve) and the backward going (blue solid curve) trajectories start diverging simultaneously with the first hops, at ca. 70 fs. This is however not surprising as the coupling mode (the h-vector) corresponds to torsion about the central C=C bond and its projection on the forward and backward trajectories should show diverging trends. As seen in Fig. 9, the composition of the swarm of trajectories represented by the averaged amplitudes h ak (t)i varies with time. Covariance Cov(k, l ) = h ak (t)ih al (t)i − h ak (t)i · h al (t)i

(16)

between the amplitudes h ak (t)i of various normal modes (overbar denotes averaging over time interval) shows what modes display similar time evolution and what modes display 30 ACS Paragon Plus Environment

Page 31 of 50

h-vector

g-vector

g→ →

ak2

h

normal mode #

projection

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 10: Branching plane vectors of MECICen (upper panel), their decomposition into the normal modes of S1min Cen (middle panel), and projections of the BP vectors onto the averaged trajectory (lower panel). In the lower panel, solid lines show the projections of the h vectors of the forward (red) and backward (blue) reactions; dotted lines show projections of the g vector, respectively. The shaded area in the lower panel shows the time interval during which the majority of surface hops occur. 31 ACS Paragon Plus Environment


#6: C=C #3: single bond scissoring torsion ν6 = 410 cm-1 ν3 = 280 cm-1

#25: hydrogen out-of-plane ν25 = 1560 cm-1

Covariance

Covariance

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 32 of 50

n2

n2 n1

n1

Figure 11: Covariance between the amplitudes h ak (t)i of the normal modes in the decomposition of forward going trajectories (left panel) and backward going trajectories (right panel). The covariance was calculated over the time interval 0 fs ≤ t ≤ 100 fs. opposite time evolution. Fig. 11 shows covariance Cov(k, l ) between the amplitudes of the normal modes in the forward going (left panel) and backward going (right panel) trajectories calculated over the time interval 0 fs ≤ t ≤ 100 fs where the upper limit corresponds to the middle of the shaded area in Figs. 8 and 10, i.e., the time from the beginning of the reaction to the moment when the first half of the hops occurred. Although the covariance in Fig. 11 was calculated up to t = 100 fs, using the full time interval 0 fs ≤ t ≤ 300 fs does not lead to substantial qualitative change of the picture; the purpose of using the short time interval was to analyze the behavior of the amplitudes at the beginning of the reaction and to investigate whether similarities in the initial behavior of the amplitudes can explain the outcome of the reaction.


Page 33 of 50

1.5

1.0

0.5

0.0 -0.5 -1.0 -1.5

normal mode #

normal mode #

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


time, fs

time, fs

Figure 12: Density plot showing the covariance between the normal mode #6 (central C=C scissoring mode) and the rest of the S1min Cen normal modes as a function of time. Left panel shows the covariance for the forward going trajectories, right panel – for the backward going trajectories. It is seen in Fig. 11 that the forward going trajectories afford considerably stronger correlation between the variations of the normal mode amplitudes h ak (t)i; which suggests a greater degree of synchronization between the normal modes contributing to the forward going trajectories. Especially, the anticorrelation (seen as dark spots below the respective dips of the plots in Fig. 11) between the amplitudes of the normal modes #3 and #25 on one hand and of the mode #6 on the other is prominent. The mode #6 describes scissoring motion of the allyl and vinylamine moieties around the central C=C bond (C=C scissoring), the other two modes make the greatest contributions to the h-vector of MECICen ’s branching plane. Time evolution of the anticorrelation of the mode #6 with the modes #3 and #25 is illustrated in Fig. 12 where the covariance Cov(6, k ), k ∈ [1, 36] is shown as a



function of time. From Figs. 10, 11, and 12 the following scenario can be suggested. Initially, the swarm of trajectories begins to slide down the slope on the S1 PES in the direction of torsion about the central C=C bond. This is evidenced by the solid curves in Fig. 10 showing that the average projection of all the trajectories on the BP coupling mode (h-vector, describing the C=C bond torsion) is large at the beginning and is indistinguishable for the forward, i.e., trans → cis (red curves in Fig. 10), propagating trajectories and the backward, i.e., returning to trans (blue curves in Fig. 10), propagating trajectories. After ca. 70 fs, when the first S1 → S0 hops occur, the swarm of trajectories begins to split into two branches, where, as evidenced by Figs. 11 and 12, the trajectories with stronger synchronization between certain normal modes (specifically, #6 and #3 and #25) continue to propagate in the forward direction, whereas less synchronized trajectories start turning back and returning to the trans-conformation. In deciding whether a trajectory will propagate in the forward or the backward direction, the antisynchronization between the C=C scissoring mode, mode #6, with the normal modes comprising the BP h-vector, primarily modes #3 and #25, seems to play the dominant role; those trajectories where the latter modes are stronger synchronized with one another and, simultaneously, are antisynchronized with the C=C scissoring mode have a greater propensity for the forward propagation. As the period of the vibrational motion along the modes #3 (ν3 = 280 cm−1 , τ = 119 fs) and #6 (ν6 = 410 cm−1 , τ = 81 fs) matches reasonably the S1 lifetime (99±51 fs), the cooperation between these low frequency modes determines whether the reaction will propagate forward or will turn backward. This scenario is generally consistent with the analysis of the CASSCF NAMD trajectories of PSB3 by Weingart et al., 31 who proposed that synchronization between the displacements along the C2 –C3 =C4 –C5 dihedral angle and the


Page 34 of 50

Page 35 of 50 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


HOOP dihedral angle results in predominant formation of the reaction product, i.e. the cis-PSB3 conformer.

4 Conclusions In this work, we employed two new computational methodologies, the DISH-XF nonadiabatic molecular dynamics algorithm 16 and the SSR ensemble density functional method, 17–21 to model the non-adiabatic dynamics of the lowest excited singlet state of trans-PSB3 cation in the gas phase. The mixed quantum-classical DISH-XF NAMD algorithm is derived from the exact factorization of the electronic-nuclear wavefunction 22–24 and is capable of correctly describing the decoherence phenomena arising from the electronicnuclear correlations. 25 The SSR method is based on ensemble DFT for the ground 26 and excited 27 electronic states and provides an accurate description of the ground and excited states of molecular systems undergoing bond breaking/bond formation processes. In particular the SSR method is capable of describing the conical intersections between the ground and excited electronic states with an accuracy matching the most sophisticated multi-reference wavefunction theory methodologies. 2,18,20,28,29,50 The geometries of the stationary points on the S1 PES of PSB3 and the S1 /S0 conical intersections were optimized using the SSR-ωPBEh/6-31G* method and compared with the available reference data from the literature. A comparison with the DMC calculations 40 for the S1 /S0 PESs crossing along the standard BLA path 51 revealed that SSR-ωPBEh/631G* reproduces the most accurate multi-reference ab initio results with an outstanding ˚ for the geometry and 0.5 kcal/mol for the energy, see Fig. 4. The SSR accuracy, 0.006 A method has identified most of the local minima and conical intersections obtained in Ref.



38 with the use of the SA-CASSCF and MSPT2 methods and relevant for the non-adiabatic dynamics of the S1 state. The direct dynamics simulations using the DISH-XF/SSR-ωPBEh/6-31G* method produced the following parameters of the non-adiabatic decay of the S1 state of trans-PSB3: the S1 /S0 hoptime τ = 99 ± 51 fs and the quantum yield of the cis-PSB3 conformation φ = 0.63. These parameters are consistent with the previously reported results of the AIMS/MSPT2 and AIMS/CASSCF dynamics simulations. 38 Similar to AIMS/MSPT2, the DISH-XF/SSR-ωPBEh/6-31G* simulations predict that the S1 decay is split into the major, trans → cis isomerization, and a minor, torsion about the terminal C=C bond, reaction channels with the ratio of 0.82/0.18 (AIMS/MSPT2 ratio 0.89/0.11). The obtained quantum yield of the trans → cis isomerization (0.63) is somewhat lower than the AIMS/MSPT2 value (0.79) and is closer to the AIMS/CASSCF value (0.54). However, the AIMS/MSPT2 simulations employed a truncated 6-31G basis set, where the polarization functions were omitted; 38 the DISH-XF/SSR-ωPBEh/6-31G* simulations did not truncate the basis set. Given the accuracy demonstrated by the SSR-ωPBEh/6-31G* method for the description of the most significant mechanistic features of the S1 and S0 PESs, the dynamics obtained in this work seems sufficiently reliable. Although in this work mixed quantum-classical NAMD simulations were carried out, where the nuclear motion was described by the classical formalism, analysis of possible composition of the quantum mechanical nuclear wave packet was attempted using the normal mode analysis of the DISH-XF/SSR-ωPBEh/6-31G* trajectories. The performed NMA revealed that only a few vibrational normal modes make a considerable contribution to the atomic displacements along the trajectories. In particular, the normal modes corresponding to C=C scissoring and C–C single bond torsion and hydrogen out-


Page 36 of 50

Page 37 of 50 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


of-plane motion play a prominent role for the outcome of the non-adiabatic relaxation of the S1 state of trans-PSB3. Synchronized (coherent) displacements along these normal modes drives the molecule towards the cis-conformation during the non-adiabatic S1 decay, as opposed to a dis-synchronous (incoherent) displacement which predominantly returns the molecule to the initial trans-conformation. This result shows that the quantum yield of trans-PSB3 photoisomerization is controlled not only by the HOOP normal mode alone, 34–36 but by its synchronization with other low frequency normal modes. The obtained results also suggest that NMA can be a useful tool in the repertoire of non-adiabatic molecular dynamics simulations of individual molecules.

Acknowledgement This work was supported by Brain Pool Program (2018H1D3A2000493) and National Honor Scientist Program (2010-0020414) through the National Research Foundation of Korea(NRF) funded by the Ministry of Science and ICT. SKM acknowledges financial support by Basic Science Research Program through the NRF funded by the Ministry of Education (2016R1C1B2015103).

Supporting Information Available Cartesian coordinates of the PSB3 stationary points, population analysis of the trajectories, and a few sample trajectories obtained in this work.



References (1) (a) Atchity, G. J.; Xantheas, S. S.; Ruedenberg, K. Potential energy surfaces near intersections. J. Chem. Phys. 1991, 95, 1862–1876; (b) Yarkony, D. R. Diabolical Conical Intersections. Rev. Mod. Phys. 1996, 68, 985–1013; (c) Bernardi, F.; Olivucci, M.; Robb, M. A. Potential energy surface crossings in organic photochemistry. Chem. Soc. Rev. 1996, 25, 321–328. (2) (a) Gozem, S.; Melaccio, F.; Valentini, A.; Filatov, M.; Huix-Rotllant, M.; Ferré, N.; Frutos, L. M.; Angeli, C.; Krylov, A. I.; Granovsky, A. A.; Lindh, R.; Olivucci, M. Shape of Multireference, Equation-of-Motion Coupled-Cluster, and Density Functional Theory Potential Energy Surfaces at a Conical Intersection. J. Chem. Theory ´ Comput. 2014, 10, 3074–3084; (b) Tuna, D.; Lefrancois, D.; Wolanski, Ł.; Gozem, S.; ´ T.; Dreuw, A.; Olivucci, M. Assessment of Approximate Schapiro, I.; Andruniow, Coupled-Cluster and Algebraic-Diagrammatic-Construction Methods for Groundand Excited-State Reaction Paths and the Conical-Intersection Seam of a RetinalChromophore Model. J. Chem. Theory Comput. 2015, 11, 5758–5781. (3) (a) Siegbahn, P.; Heiberg, A.; Roos, B.; Levy, B. Comparison of the Super-CI and the Newton-Raphson Scheme in the Complete Active Space SCF Method. Phys. Scripta 1980, 21, 323–327; (b) Roos, B.; Taylor, P.; Siegbahn, P. A Complete Active Space SCF Method(CASSCF) Using a Density-matrix Formulated Super-CI Approach. Chem. ¨ J.; Heiberg, A.; Roos, B. The ComPhys. 1980, 48, 157–173; (c) Siegbahn, P.; Almlof, plete Active Space SCF (CASSCF) Method in a Newton-Raphson Formulation with Application to the HNO Molecule. J. Chem. Phys. 1981, 74, 2384–2396; (d) Roos, B.


Page 38 of 50

Page 39 of 50 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 Complete Active Space Self-Consistent Field Method and its Applications in Electronic Structure Calculations. Adv. Chem. Phys. 1987, 69, 399–445. ˚ Roos, B. O. Second-order perturbation theory with (4) Andersson, K.; Malmqvist, P.-A.; a complete active space self-consistent field reference function. J. Chem. Phys. 1992, 96, 1218–1226. (5) Diffenderfer, R. N.; Yarkony, D. R. Use of the State-Averaged MCSCF Procedure Application to Radiative Transitions in MgO. J. Phys. Chem. 1982, 86, 5098–5105. ¨ (6) Shiozaki, T.; Gyorffy, W.; Celani, P.; Werner, H.-J. Communication: Extended multistate complete active space second-order perturbation theory: Energy and nuclear gradients. J. Chem. Phys. 2011, 135, 081106. (7) (a) Hohenstein, E. G.; Luehr, N.; Ufimtsev, I. S.; Mart´ınez, T. J. An atomic orbitalbased formulation of the complete active space self-consistent field method on graphical processing units. J. Chem. Phys. 2015, 142, 224103; (b) Snyder, J., J. W.; Hohenstein, E. G.; Luehr, N.; Mart´ınez, T. J. An Atomic Orbital-Based Formulation of Analytical Gradients and Nonadiabatic Coupling Vector Elements for the StateAveraged Complete Active Space Self-Consistent Field Method on Graphical Processing Units. J. Chem. Phys. 2015, 143, 154107. (8) (a) Meyer, H.-D.; Manthe, U.; Cederbaum, L. S. The multi-configurational timedependent Hartree approach. Chem. Phys. Lett. 1990, 165, 73–78; (b) Meyer, H.-D.; Gatti, F.; Worth, G. A. Multidimensional quantum dynamics; John Wiley & Sons, 2009; (c) Worth, G.; Meyer, H.-D.; Cederbaum, L. The effect of a model environment on the S 2 absorption spectrum of pyrazine: A wave packet study treating all 24 vibrational ¨ modes. J. Chem. Phys. 1996, 105, 4412–4426; (d) Worth, G.; Meyer, H.-D.; Koppel, H.; 39 ACS Paragon Plus Environment


Cederbaum, L.; Burghardt, I. Using the MCTDH wavepacket propagation method to describe multimode non-adiabatic dynamics. Int. Rev. Phys. Chem. 2008, 27, 569–606. (9) Worth, G. A.; Burghardt, I. Full quantum mechanical molecular dynamics using Gaussian wavepackets. Chem. Phys. Lett. 2003, 368, 502–508. (10) (a) Mart´ınez, T. J.; Levine, R. D. First-principles molecular dynamics on multiple electronic states: A case study of NaI. J. Chem. Phys. 1996, 105, 6334–6341; (b) BenNun, M.; Mart´ınez, T. J. Ab initio molecular dynamics study of cis-trans photoisomerization in ethylene. Chem. Phys. Lett. 1998, 298, 57 – 65; (c) Virshup, A. M.; Punwong, C.; Pogorelov, T. V.; Lindquist, B. A.; Ko, C.; Mart´ınez, T. J. Photodynamics in Complex Environments: Ab Initio Multiple Spawning Quantum Mechanical/Molecular Mechanical Dynamics. J. Phys. Chem. B 2009, 113, 3280–3291. (11) Snyder, J. W.; Curchod, B. F. E.; Mart´ınez, T. J. GPU-Accelerated State-Averaged CASSCF Interfaced with Ab Initio Multiple Spawning Unravels the Photodynamics of Provitamin D3. J. Phys. Chem. Lett. 2016, 7, 2444–2449. (12) (a) Tully, J. C. Molecular dynamics with electronic transitions. J. Chem. Phys. 1990, 93, 1061; (b) Hammes-Schiffer, S.; Tully, J. C. Proton transfer in solution: Molecular dynamics with quantum transitions. J. Chem. Phys. 1994, 101, 4657–4667. (13) Wan, C.-C.; Schofield, J. Mixed quantum-classical molecular dynamics: Aspects of the multithreads algorithm. J. Chem. Phys. 2000, 113, 7047–7054. (14) (a) Wang, L.; Akimov, A.; Prezhdo, O. V. Recent Progress in Surface Hopping: 20112015. J. Phys. Chem. Lett. 2016, 7, 2100; (b) Subotnik, J. E.; Jain, A.; Landry, B.; Petit, A.; Ouyang, W.; Bellonzi, N. Understanding the Surface Hopping View of Elec40 ACS Paragon Plus Environment

Page 40 of 50

Page 41 of 50 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


tronic Transitions and Decoherence. Ann. Rev. Phys. Chem. 2016, 67, 387–417; (c) Jaeger, H. M.; Fischer, S.; Prezhdo, O. V. Decoherence-induced surface hopping. J. Chem. Phys. 2012, 137, 22A545; (d) Granucci, G.; Persico, M. Critical appraisal of the fewest switches algorithm for surface hopping. J. Chem. Phys. 2007, 126, 134114; (e) Gao, X.; Thiel, W. Non-Hermitian surface hopping. Phys. Rev. E 2017, 95, 013308; (f) Tapavicza, E.; Tavernelli, I.; Rothlisberger, U. Trajectory Surface Hopping within Linear Response Time-Dependent Density-Functional Theory. Phys. Rev. Lett. 2007, 98, 023001; (g) Jasper, A. W.; Nangia, S.; Zhu, C.; Truhlar, D. G. Non-Born-Oppenheimer Molecular Dynamics. Acc. Chem. Res. 2006, 39, 101. ¨ (15) (a) Prezhdo, O. V. Mean field approximation for the stochastic Schrodinger equation. J. Chem. Phys. 1999, 111, 8366; (b) Li, X.; Tully, J. C.; Schlegel, H. B.; Frisch, M. J. Ab initio Ehrenfest dynamics. J. Chem. Phys. 2005, 123, 084106; (c) Bedard-Hearn, M. J.; Larsen, R. E.; Schwartz, B. J. Mean-field dynamics with stochastic decoherence (MF-SD): A new algorithm for nonadiabatic mixed quantum/classical moleculardynamics simulations with nuclear-induced decoherence. J. Chem. Phys. 2005, 123, 234106; (d) Akimov, A. V.; Long, R.; Prezhdo, O. V. Coherence penalty functional: A simple method for adding decoherence in Ehrenfest dynamics. J. Chem. Phys. 2014, 140, 194107. (16) Ha, J.-K.; Lee, I. S.; Min, S. K. Surface Hopping Dynamics beyond Nonadiabatic Couplings for Quantum Coherence. J. Phys. Chem. Lett. 2018, 9, 1097–1104. (17) Kazaryan, A.; Heuver, J.; Filatov, M. Excitation Energies from Spin-Restricted Ensemble-Referenced Kohn-Sham Method: A State-Average Approach. J. Phys. Chem. A 2008, 112, 12980–12988.



(18) Filatov, M. Assessment of density functional methods for obtaining geometries at conical intersections in organic molecules. J. Chem. Theory Comput. 2013, 9, 4526–4541. (19) Filatov, M. Spin-restricted ensemble-referenced Kohn-Sham method: basic principles and application to strongly correlated ground and excited states of molecules. WIREs Comput. Mol. Sci. 2015, 5, 146–167. (20) Filatov, M. In Density-functional methods for excited states; Ferré, N., Filatov, M., HuixRotllant, M., Eds.; Top. Curr. Chem.; Springer: Heidelberg, 2016; Vol. 368; pp 97–124. (21) Filatov, M.; Liu, F.; Mart´ınez, T. J. Analytical derivatives of the individual state energies in ensemble density functional theory method. I. General formalism. J. Chem. Phys. 2017, 147, 034113. (22) Hunter, G. Conditional probability amplitudes in wave mechanics. Int. J. Quantum Chem. 1975, 9, 237–242. (23) Abedi, A.; Maitra, N. T.; Gross, E. K. U. Exact factorization of the time-dependent electron-nuclear wave function. Phys. Rev. Lett. 2010, 105, 123002. (24) (a) Abedi, A.; Maitra, N. T.; Gross, E. K. U. Correlated electron-nuclear dynamics: Exact factorization of the molecular wave-function. J. Chem. Phys. 2012, 137, 22A530; (b) Abedi, A.; Agostini, F.; Suzuki, Y.; Gross, E. K. U. Dynamical steps that bridge piecewise adiabatic shapes in the exact time-dependent potential energy surface. Phys. Rev. Lett. 2013, 110, 263001; (c) Agostini, F.; Abedi, A.; Suzuki, Y.; Min, S. K.; Maitra, N. T.; Gross, E. K. U. The exact electronic back-reaction on classical nuclei in non-adiabatic charge transfer. J. Chem. Phys. 2015, 142, 084303.


Page 42 of 50

Page 43 of 50 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


(25) (a) Min, S. K.; Agostini, F.; Gross, E. K. U. Coupled-trajectory quantum-classical approach to electronic decoherence in nonadiabatic processes. Phys. Rev. Lett. 2015, 115, 073001; (b) Agostini, F.; Min, S. K.; Abedi, A.; Gross, E. K. U. Quantum-Classical Nonadiabatic Dynamics: Coupled- vs Independent-Trajectory Methods. J. Chem. Theory Comput. 2016, 12, 2127–2143; (c) Min, S. K.; Agostini, F.; Tavernelli, I.; Gross, E. K. U. Ab Initio Nonadiabatic Dynamics with Coupled Trajectories: A Rigorous Approach to Quantum (De)Coherence. J. Phys. Chem. Lett. 2017, 8, 3048–3055. (26) (a) Valone, S. M. A one-to-one mapping between one-particle densities and some n-particle ensembles. J. Chem. Phys. 1980, 73, 4653–4655; (b) Lieb, E. H. Density functionals for Coulomb systems. Int. J. Quantum Chem. 1983, 24, 243–277; (c) Perdew, J. P.; Parr, R. G.; Levy, M.; Balduz Jr., J. L. Density-Functional Theory for Fractional Particle Number: Derivative Discontinuities of the Energy. Phys. Rev. Lett. 1982, 49, 1691–1694; (d) Englisch, H.; Englisch, R. Hohenberg-Kohn Theorem and Non-VRepresentable Densities. Physica 1983, A121, 253–268; (e) Englisch, H.; Englisch, R. Exact Density Functionals for Ground-State Energies. I. General Results. Phys. Stat. Sol. (b) 1984, 123, 711–721; (f) Englisch, H.; Englisch, R. Exact Density Functionals for Ground-State Energies II. Details and Remarks. Phys. Stat. Sol. (b) 1984, 124, 373–379. (27) (a) Gross, E. K. U.; Oliveira, L. N.; Kohn, W. Rayleigh-Ritz variational principle for ensembles of fractionally occupied states. Phys. Rev. A 1988, 37, 2805–2808; (b) Gross, E. K. U.; Oliveira, L. N.; Kohn, W. Density-functional theory for ensembles of fractionally occupied states. I. Basic formalism. Phys. Rev. A 1988, 37, 2809–2820; (c) Oliveira, L. N.; Gross, E. K. U.; Kohn, W. Density-functional theory for ensembles of fractionally occupied states. II. Application to the He atom. Phys. Rev. A 1988, 37,



2821–2833; (d) Oliveira, L. N.; Gross, E. K. U.; Kohn, W. Ensemble-Density Functional Theory. Int. J. Quantum Chem.: Quantum Chem. Symp. 1990, 24, 707–716. (28) Nikiforov, A.; Gamez, J. A.; Thiel, W.; Huix-Rotllant, M.; Filatov, M. Assessment of approximate computational methods for conical intersections and branching plane vectors in organic molecules. J. Chem. Phys. 2014, 141, 124122. (29) Huix-Rotllant, M.; Nikiforov, A.; Thiel, W.; Filatov, M. In Density-functional methods for excited states; Ferré, N., Filatov, M., Huix-Rotllant, M., Eds.; Top. Curr. Chem.; Springer: Heidelberg, 2016; Vol. 368; pp 445–476. (30) (a) Garavelli, M.; Bernardi, F.; Olivucci, M.; Vreven, T.; Klein, S.; Celani, P.; Robb, M. A. Potential-energy surfaces for ultrafast photochemistry Static and dynamic aspects. Faraday Discuss. 1998, 110, 51–70; (b) Garavelli, M.; Bernardi, F.; Robb, M. A.; Olivucci, M. The short-chain acroleiniminium and pentadieniminium cations: towards a model for retinal photoisomerization. A CASSCF/PT2 study. J. Mol. Struct.: THEOCHEM 1999, 463, 59–64; (c) Page, C. S.; Olivucci, M. Ground and excited state CASPT2 geometry optimizations of small organic molecules. J. Comput. Chem. 2003, 24, 298–309; (d) Sinicropi, A.; Migani, A.; De Vico, L.; Olivucci, M. Photoisomerization acceleration in retinal protonated Schiff-base models. Photochem. Photobiol. Sci. 2003, 2, 1250–1255; (e) Fantacci, S.; Migani, A.; Olivucci, M. CASPT2//CASSCF and TDDFT//CASSCF Mapping of the Excited State Isomerization Path of a Minimal Model of the Retinal Chromophore. J. Phys. Chem. A 2004, 108, 1208–1213; (f) Valsson, O.; Filippi, C. Photoisomerization of Model Retinal Chromophores: Insight from Quantum Monte Carlo and Multiconfigurational Perturbation Theory. J. Chem. Theory Comput. 2010, 6, 1275–1292; (g)


Page 44 of 50

Page 45 of 50 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


Gozem, S.; Huntress, M.; Schapiro, I.; Lindh, R.; Granovsky, A. A.; Angeli, C.; Olivucci, M. Dynamic Electron Correlation Effects on the Ground State Potential Energy Surface of a Retinal Chromophore Model. J. Chem. Theory Comput. 2012, 8, 4069–4080; (h) Gozem, S.; Melaccio, F.; Lindh, R.; Krylov, A. I.; Granovsky, A. A.; Angeli, C.; Olivucci, M. Mapping the Excited State Potential Energy Surface of a Retinal Chromophore Model with Multireference and Equation-of-Motion Coupled-Cluster Methods. J. Chem. Theory Comput. 2013, 9, 4495–4506. (31) (a) Weingart, O.; Schapiro, I.; Buss, V. Bond torsion affects the product distribution in the photoreaction of retinal model chromophores. J. Mol. Model. 2006, 12, 713–721; (b) Weingart, O.; Buss, V.; Robb, M. A. Excited state molecular dynamics of retinal model chromophores. Phase Transitions 2005, 78, 17–24. (32) Barbatti, M.; Ruckenbauer, M.; Szymczak, J. J.; Aquino, A. J. A.; Lischka, H. Nonadiabatic excited-state dynamics of polar π-systems and related model compounds of biological relevance. Phys. Chem. Chem. Phys. 2008, 10, 482–494. (33) 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. (34) Weingart, O. The role of HOOP-modes in the ultrafast photo-isomerization of retinal models. Chem. Phys. 2008, 349, 348–355. (35) Schapiro, I.; Ryazantsev, M. N.; Frutos, L. M.; Ferre, N.; Lindh, R.; Olivucci, M. The Ultrafast Photoisomerizations of Rhodopsin and Bathorhodopsin Are Modulated by



Bond Length Alternation and HOOP Driven Electronic Effects. J. Am. Chem. Soc. 2011, 133, 3354–3364. (36) Johnson, P. J. M.; Halpin, A.; Morizumi, T.; Prokhorenko, V. I.; Ernst, O. P.; Miller, R. J. D. Local vibrational coherences drive the primary photochemistry of vision. Nat. Chem. 2015, 7, 980–986. ¨ (37) Szalay, P. G.; Muller, T.; Gidofalvi, G.; Lischka, H.; Shepard, R. Multiconfiguration Self-Consistent Field and Multireference Configuration Interaction Methods and Applications. Chem. Rev. 2012, 112, 108–181. (38) Liu, L.; Liu, J.; Mart´ınez, T. J. Dynamical Correlation Effects on Photoisomerization: Ab Initio Multiple Spawning Dynamics with MS-CASPT2 for a Model transProtonated Schiff Base. J. Phys. Chem. B 2016, 120, 1940–1949. (39) Huix-Rotllant, M.; Filatov, M.; Gozem, S.; Schapiro, I.; Olivucci, M.; Ferré, N. Assessment of Density Functional Theory for Describing the Correlation Effects on the Ground and Excited State Potential Energy Surfaces of a Retinal Chromophore Model. J. Chem. Theory Comput. 2013, 9, 3917–3932. (40) Zen, A.; Coccia, E.; Gozem, S.; Olivucci, M.; Guidoni, L. Quantum Monte Carlo Treatment of the Charge Transfer and Diradical Electronic Character in a Retinal Chromophore Minimal Model. J. Chem. Theory Comput. 2015, 11, 992–1005. (41) Granucci, G.; Persico, M.; Zoccante, A. Including quantum decoherence in surface hopping. J. Chem. Phys. 2010, 133, 134111. (42) (a) Filatov, M.; Shaik, S. A spin-restricted ensemble-referenced Kohn-Sham method and its application to diradicaloid situations. Chem. Phys. Lett. 1999, 304, 429–437; 46 ACS Paragon Plus Environment

Page 46 of 50

Page 47 of 50 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


(b) Moreira, I. d. P. R.; Costa, R.; Filatov, M.; Illas, F. Restricted ensemble-referenced Kohn-Sham versus broken symmetry approaches in density functional theory: Magnetic coupling in Cu binuclear complexes. J. Chem. Theory Comput. 2007, 3, 764–774. (43) Filatov, M.; Shaik, S. Spin-restricted density functional approach to the open-shell problem. Chem. Phys. Lett. 1998, 288, 689–697. (44) (a) Ufimtsev, I. S.; Mart´ınez, T. J. Quantum Chemistry on Graphical Processing Units. 1. Strategies for two-electron integral evaluation. J. Chem. Theory Comput. 2008, 4, 222–231; (b) Ufimtsev, I. S.; Mart´ınez, T. J. Quantum Chemistry on Graphical Processing Units. 2. Direct Self-Consistent-Field (SCF) Implementation. J. Chem. Theory Comput. 2009, 5, 1004–1015; (c) Ufimtsev, I. S.; Mart´ınez, T. J. Quantum Chemistry on Graphical Processing Units. 3. Analytical Energy Gradients and First Principles Molecular Dynamics. J. Chem. Theory Comput. 2009, 5, 2619–2628; (d) Ufimtsev, I. S.; Mart´ınez, T. J. Graphical Processing Units for Quantum Chemistry. Comput. Sci. Eng. 2008, 10, 26–34; (e) Titov, A. V.; Ufimtsev, I. S.; Luehr, N.; Mart´ınez, T. J. Generating Efficient Quantum Chemistry Codes for Novel Architectures. J. Chem. Theory Comput. 2013, 9, 213–221; (f) Song, C.; Wang, L.-P.; Mart´ınez, T. J. Automated Code Engine for Graphical Processing Units: Application to the Effective Core Potential Integrals and Gradients. J. Chem. Theory Comput. 2016, 12, 92–106. (45) Rohrdanz, M. A.; Martins, K. M.; Herbert, J. M. A long-range-corrected density functional that performs well for both ground-state properties and time-dependent density functional theory excitation energies, including charge-transfer excited states. J. Chem. Phys. 2009, 130, 054112. (46) Krishnan, R.; Binkley, J. S.; Seeger, R.; Pople, J. A. Self-consistent molecular orbital 47 ACS Paragon Plus Environment


methods. XX. A basis set for correlated wave functions. J. Chem. Phys. 1980, 72, 650– 654. (47) Kästner, J.; Carr, J. M.; Keal, T. W.; Thiel, W.; Wander, A.; Sherwood, P. DL-FIND: An Open-Source Geometry Optimizer for Atomistic Simulations. J. Phys. Chem. A 2009, 113, 11856–11865. (48) Levine, B.; Coe, J. D.; Mart´ınez, T. J. Optimizing conical intersections without derivative coupling vectors: application to multistate multireference second-order perturbation theory (MS-CASPT2). J. Phys. Chem. B 2008, 112, 405–413. (49) Valsson, O.; Angeli, C.; Filippi, C. Excitation energies of retinal chromophores: critical role of the structural model. Phys. Chem. Chem. Phys. 2012, 14, 11015–11020. (50) Huix-Rotllant, M.; Filatov, M.; Gozem, S.; Schapiro, I.; Olivucci, M.; Ferré, N. Assessment of density functional theory for describing the correlation effects on the ground and excited state potential energy surfaces of a retinal chromophore model. J. Chem. Theory Comput. 2013, 9, 3917–3932. (51) Gozem, S.; Huntress, M.; Schapiro, I.; Lindh, R.; Granovsky, A. A.; Angeli, C.; Olivucci, M. Dynamic Electron Correlation Effects on the Ground State Potential Energy Surface of a Retinal Chromophore Model. J. Chem. Theory Comput. 2012, 8, 4069–4080. (52) Gozem, S.; Krylov, A. I.; Olivucci, M. Conical Intersection and Potential Energy Surface Features of a Model Retinal Chromophore: Comparison of EOM-CC and Multireference Methods. J. Chem. Theory Comput. 2013, 9, 284–292.


Page 48 of 50

Page 49 of 50 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


(53) Kästner, J.; Carr, J. M.; Keal, T. W.; Thiel, W.; Wander, A.; Sherwood, P. DL-FIND: An Open-Source Geometry Optimizer for Atomistic Simulations. J. Phys. Chem. A 2009, 113, 11856–11865. (54) Oppenheim, I.; Ross, J. Temperature Dependence of Distribution Functions in Quantum Statistical Mechanics. Phys. Rev. 1957, 107, 28–32. (55) Davies, R. W.; Davies, K. T. R. On the Wigner Distribution Function for an Oscillator. Ann. Phys. 1975, 89, 261–273. (56) Nikiforov, A.; Gamez, J. A.; Thiel, W.; Filatov, M. Computational Design of a Family of Light-Driven Rotary Molecular Motors with Improved Quantum Efficiency. J. Phys. Chem. Lett. 2016, 7, 105–110. (57) Weingart, O.; Lan, Z.; Koslowski, A.; Thiel, W. Chiral Pathways and Periodic Decay in cis -Azobenzene Photodynamics. J. Phys. Chem. Lett. 2011, 2, 1506–1509. (58) Kube, S.; Lasser, C.; Weber, M. Monte Carlo sampling of Wigner functions and surface hopping quantum dynamics. J. Comp. Phys. 2009, 228, 1947 – 1962. (59) Bahar, I.; Lezon, T. R.; Bakan, A.; Shrivastava, I. H. Normal Mode Analysis of Biomolecular Structures: Functional Mechanisms of Membrane Proteins. Chem. Rev. 2010, 110, 1463–1497.



Graphical Abstract

trajectory surface hopping from exact factorization

+

18%

ensemble DFT for excited states

cis:trans = 63:37

hν

82%

τ(S1) = 99 ± 51 fs

Figure 13


Page 50 of 50

Direct Non-Adiabatic Dynamics by Mixed Quantum-Classical

Recommend Documents