TDDFT and Quantum-Classical Dynamics: A Universal Tool Describing the Dynamics of Matter

TDDFT and Quantum-Classical Dynamics: a Universal Tool Describing the Dynamics of Matter Federica Agostini1, Basile F. E. Curchod2, Rodolphe Vuilleumier3a,3b, Ivano Tavernelli4, E. K. U. Gross5 Abstract Time-dependent density functional theory (TDDFT) is currently the most efficient approach allowing to describe electronic dynamics in complex systems, from isolated molecules to the condensed-phase. TDDFT has been employed to investigate an extremely wide range of time-dependent phenomena, as spin dynamics in solids, charge and energy transport in nanoscale devices, and photoinduced exciton transfer in molecular aggregates. It is therefore nearly impossible to give a general account of all developments and applications of TDDFT in material science, as well as in physics and chemistry. A large variety of aspects are covered throughout these volumes, see e.g. Chapters X, Y (to be indicated). In the present chapter, we will limit our presentation to the description of TDDFT developments and applications in the field of quantum molecular dynamics simulations in combination with trajectory-based approaches for the study of nonadiabatic excited-state phenomena. We will present different quantum-classical strategies used to describe the coupled dynamics of electrons and nuclei underlying nonadiabatic processes. In addition, we will give an account of the most recent applications with the aim of illustrating the nature of the problems that can be addressed with the help of these approaches. The potential, as well as the limitations, of the presented methods are discussed, 1Laboratoire de Chimie Physique, UMR 8000 CNRS/University Paris-Sud, University ParisSaclay, 91405 Orsay, France, e-mail: [email protected] 2Department of Chemistry, Durham University, South Road, Durham DH1 3LE, United Kingdom, e-mail: [email protected] 3aPASTEUR, D´ epartement de chimie, ´ Ecole normale sup´ erieure, UPMC Univ. Paris 06, CNRS, PSL Research University, 75005 Paris, France. 3bSorbonne Universit´ es, UPMC Univ. Paris 06, ´ Ecole normale suprieure, CNRS, Processus d’activation s´ elective par transfert d’´ energie uni-´ electronique ou radiatif (PASTEUR), 75005 Paris, France, e-mail: [email protected] 4Zurich Research Laboratory, IBM Research GmbH, 8803 R¨ uschlikon, Switzerland, e-mail: [email protected] 5Max-Planck-Institut f¨ ur Mikrostrukturphysik, Weinberg 2, Halle 06120, Germany, e-mail: [email protected] 1 2 E. K. U. Gross et al. along with possible avenues for future developments in TDDFT and nonadiabatic dynamics. 1 Introduction Photoinduced isomerization processes, photosynthetic and photovoltaic energy conversion phenomena, charge and energy transport through molecular junctions, are all typical examples of, so-called, nonadiabatic processes. Nonadiabatic processes are characterized by a strong coupling between electronic and nuclear motion; in fact, nuclear motion is responsible for inducing electronic (nonadiabatic) transitions, and in turn, the time evolution of the electronic states also affects the nuclear dynamics at very short timescales (down to a few tens of fs). In this nonadiabatic regime, thus when the Born-Oppenheimer approximation breaks down, performing (quantum) molecular dynamics simulations is tremendously challenging. Accurate electronic structure properties are required to describe electronic dynamics, and to correctly drive the nuclear evolution. Identifying regions of nuclear configuration space where the electronic states are coupled, as avoided crossings and conical intersections, is essential to predict quantum yields. Efficient evolution techniques have to be employed to describe nuclear motion in order, for instance, to determine final molecular structures, or to account for possible quantum effects. Therefore, theoretical and numerical developments need to address the problem from the perspective of both electronic structure theory and nuclear quantum dynamics. Perhaps the most celebrated method to investigate excited electronic states is time-dependent density functional theory (TDDFT). TDDFT offers an in principle exact formalism for propagating the time-dependent electronic density and, within linear response theory, for calculating excitation energies as well as critical excited-state properties. It is therefore without any surprise that TDDFT became the electronic-structure method of choice to be coupled with nonadiabatic dynamics. Particularly successful has been the combination of TDDFT, employed to solve the electronic problem, with the description of nuclear motion in terms of trajectories that evolve or hop between coupled (electronic) potential energy surfaces. The most well-known method is Tully’s “fewest switches” trajectory surface hopping (Tully (1990)), which has evolved into a widely used and successful technique. The mean-field Ehrenfest dynamics is often employed to investigate explicitly the electronic dynamics, combined for example with real-time TDDFT Tavernelli et al (2005); Tavernelli (2006). Full multiple spawning (Mart´ ınez et al (1996); Mart´ ınez and Levine (1997); Ben-Nun and Mart´ ınez (1998); Ben-Nun et al (2000); Hack et al (2001); Ben-Nun and Mart´ ınez (2002); Virshup et al (2008)) propagates coupled Gaussian functions along classical trajectories, whereas the coupled-trajectory mixed quantum-classical (CT-MQC) scheme (Min et al (2015)) derived from the Exact Factorization (Abedi et al (2010)) is based on the propagation of trajectories along a time-dependent potential energy surface (Abedi et al (2013a)). Other techniques like the quantum-classical Liouville equation (Kapral and Ciccotti (1999); TDDFT and Quantum-Classical Dynamics 3 Nielsen et al (2000); Kapral (2006)), Bohmian dynamics (Wyatt et al (2001); Lopreore and Wyatt (2002); Rassolov and Garashchuk (2005); Curchod and Tavernelli (2013a)), variational multiconfiguration Gaussians (Worth et al (2004); Lasorne et al (2006, 2007); Worth et al (2008); Mendive-Tapia et al (2012); Richings et al (2015)), multiconfigurational Ehrenfest (Shalashilin (2010); Saita and Shalashilin (2012); Makhov et al (2017)) or linearization approaches to compute timecorrelations functions (Bonella and Coker (2005); Huo and Coker (2012); Dunkel et al (2008)) have been also proposed for nonadiabatic dynamics. Despite their differences, all the methods mentioned above are rooted in the Born-Huang representation of the total molecular wavefunction, i.e., an expansion in an infinite sum over the correlated Born-Oppenheimer electronic states. In contrast, the recently introduced Exact Factorization of the time-dependent molecular wavefunction offers a paradigm shift in our perception of nonadiabatic dynamics, away from the Born-Huang picture, and blaze a trail for the development of nonadiabatic techniques away from Born-Oppenheimer concepts. It is important to mention here that ensembles of trajectories, when properly constructed via the method of characteristics (Agostini et al (2018)), can represent, in principle arbitrarily closely, the solution of the underlying partial differential equation. Practical implementations, however, involve further-going approximations where, for instance, interference and tunnelling effects are neglected or only approximately taken into account. The advantage of trajectory-based method is that they circumvent the enormous numerical effort associated with quantum wavepackets propagation techniques, such as Multi Configuration Time Dependent Hartree approach (MCTDH) (Meyer et al (1990); Burghardt et al (1999); Wang and Thoss (2003); Meyer and Worth (2003)). By its very nature, this approach requires the computation of the relevant potential energy surfaces (PESs) and corresponding couplings before the actual propagation of nuclear wavepackets. This clearly implies an important computational effort that limits the applicability of this method to a small number of degrees of freedom (up to ∼10). In addition, the determination of the relevant degrees of freedom to include in the dynamics can also become a challenging problem, which requires some a priori knowledge of the “active” vibrational modes involved in the dynamics. Such wavefucntion-based nonadiabatic approaches are beyond the scope of this chapter and will not be discussed further. The goal of this chapter is to present in a self-contained manner the key theoretical concepts and equations of the most important methods cited above, starting from the electronic structure problem and (LR-)TDDFT, up to nuclear dynamics methods like Surface Hopping, Ehrenfest dynamics, and Ab Initio Multiple Spawning. To contrast with these Born-Huang-based methods, we also present the formalism of the Exact Factorization and introduce the reader to the first mixed-quantum classical algorithm derived from this formalism, coined coupled-trajectory mixed quantum classical (CT-MQC) dynamics. 4 E. K. U. Gross et al. 2 Coupled electron-nuclear dynamics in molecules In molecules and condensed phase systems, the time evolution of interacting electrons and nuclei is described by the time-dependent Schr¨ odinger equation ˆ HΨ(r,R,t) = i¯ h∂tΨ(r,R,t),(1) where the electron-nuclear wavefunction Ψ(r,R,t)describes the state of the system over time, and ˆ His the molecular Hamiltonian, i.e., ˆ H(r,R) = Nn ∑ ν=1 −¯ h2 2Mν ∇2 ν+ˆ Te(r)+Vee(r)+Vnn(R)+Ven(r,R) = Nn ∑ ν=1 −¯ h2 2Mν ∇2 ν+ˆ HBO(r,R).(2) Here, r= (r1,...,rNel ),R= (R1,...,RNn),Nel is the number of electrons and Nn the number of nuclei. The first term on the right-hand side of Eq. (2) is the nuclear kinetic energy, with ∇νindicating the spatial derivative with respect to the position of the nucleus ν, and Mνits mass, whereas ˆ HBO is the so-called BornOppenheimer (BO), or electronic, Hamiltonian. ˆ HBO is defined as the sum of the electronic kinetic energy, ˆ Te, the electron-electron, ˆ Vee, the nucleus-nucleus Vnn, and the electron-nucleus, Ven, interactions. Usually, the problem is reformulated adopting the Born-Huang expansion of the molecular wavefunction in the adiabatic basis. The adiabatic, or BO, states, ϕ(k) R(r), are defined as the eigenfunctions of the BO Hamiltonian, ˆ HBO(r,R)ϕ(k) R(r) = ε(k) BO(R)ϕ(k) R(r),(3) with eigenvalues ε(k) BO(R). The electronic time-independent problem is diagonalized at each nuclear position R, thus the eigenfunctions and eigenvalues depend on R. Nuclear positions are interpreted here as parameters, that label both the electronic states and the electronic energies. If Eq. (3) is solved for all nuclear configurations, ε(k) BO(R)identify the so-called BO potential energy surfaces (PESs). In the BornHuang expansion of the electron-nuclear wavefunction, Ψ(r,R,t) = ∑ k χk(R,t)ϕ(k) R(r),(4) the coefficients χk(R,t)clearly depend on nuclear positions and on time. These coefficients can be interpreted as the nuclear contributions corresponding to the electronic states included in the sum, and can be also referred to as nuclear wavepackets. In fact, it can be easily proven that the nuclear density, defined as the integral of |Ψ(r,R,t)|2over electronic coordinates, TDDFT and Quantum-Classical Dynamics 5 Zdr|Ψ(r,R,t)|2=∑ k|χk(R,t)|2,(5) can be written as the sum of adiabatic contributions, |χk(R,t)|2. Here, the orthogonality of the BO states Zdrϕ(l) R∗(r)ϕ(k) R(r) = Dϕ(l) Rϕ(k) REr=δlk (6) has been used. The Born-Huang expansion (4) is inserted in Eq. (1), that is then projected on ϕ(k) R∗(r)and integrated over r. A set of partial differential equations are derived for the expansion coefficients "Nn ∑ ν −¯ h2 2Mν ∇2 ν+ε(k) BO(R)#χk(R,t)+∑ l Fkl (R)χl(R,t) = i¯ h∂tχk(R,t).(7) The last term on the right-hand side is responsible for coupling the evolution of the k-th coefficient to all other coefficients, via the nonadiabatic couplings Fkl (R) =Zdrϕ(k) R∗(r)"Nn ∑ ν −¯ h2 2Mν ∇2 ν#ϕ(l) R(r) + Nn ∑ ν 1 MνZdrϕ(k) R∗(r)h−i¯ h∇νϕ(l) R(r)i·[−i¯ h∇ν],(8) arising from the effect of the nuclear kinetic energy operator on the parametric dependence of the BO states on R. In the most general case, the non-diagonal elements of Fkl (R)are non-zero and induce a coupling between different electronic states due to the motion of the nuclei. The nonadiabatic coupling term is responsible for exchanging “nuclear contributions” between the electronic adiabatic states k and l. The BO framework presented so far is widely adopted by a large community of physicists and chemists to interpret the coupled electron-nuclear problem under nonadiabatic conditions. However, such framework is not the only one, as we will discuss below. An alternative perspective on the coupled electron-nuclear problem has been recently proposed, the Exact Factorization of the electron-nuclear wavefunction (Abedi et al (2010, 2012)). In this framework, we make an Ansatz different from the BornHuang representation of the molecular wavefunction, namely Ψ(r,R,t) = ΦR(r,t)χ(R,t).(9) Here χ(R,t)is the nuclear wavefunction and ΦR(r,t)is the electronic wavefunction which parametrically depends on the nuclear positions and satisfies the partial normalization condition (PNC) 6 E. K. U. Gross et al. Zdr|ΦR(r,t)|2=1∀R,t.(10) The theorems introduced in (Abedi et al (2010, 2012)) prove the existence and uniqueness of Eq. (9), up to within a (R,t)-dependent gauge transformation. The PNC guarantees the interpretation of |χ(R,t)|2as the probability of finding the nuclear configuration Rat time t, and of |ΦR(r,t)|2itself as the conditional probability of finding the electronic configuration rat time tgiven the nuclear configuration R. The stationary variations (Frenkel (1934)) of the quantum mechanical action with respect to ΦR(r,t)and χ(R,t)lead to the derivation of the following equations of motion ˆ HBO(r,R)+ ˆ Ucoup en [ΦR,χ]−ε(R,t)ΦR(r,t) = i¯ h∂tΦR(r,t)(11) "Nn ∑ ν=1 [−i¯ h∇ν+Aν(R,t)]2 2Mν +ε(R,t)#χ(R,t) = i¯ h∂tχ(R,t),(12) where the PNC is enforced by means of Lagrange multipliers (Alonso et al (2013); Abedi et al (2013b)). The electron-nuclear coupling operator (Agostini et al (2015b)), ˆ Ucoup en [ΦR,χ] = Nn ∑ ν=1 1 Mν"[−i¯ h∇ν−Aν(R,t)]2 2 +−i¯ h∇νχ χ+Aν(R,t)−i¯ h∇ν−Aν(R,t),(13) the time-dependent potential energy surface (TDPES) (Abedi et al (2013a); Agostini et al (2013); Suzuki et al (2015); Agostini et al (2015a); Curchod et al (2016a); Suzuki and Watanabe (2016)), ε(R,t) = hΦR(t)|ˆ HBO +ˆ Ucoup en −i¯ h∂t|ΦR(t)ir,(14) and the time-dependent vector potential (Curchod and Agostini (2017)), Aν(R,t) = hΦR(t)|−i¯ h∇νΦR(t)ir(15) are responsible for the coupling between electrons and nuclei in a formally exact way. It is worth noting that the electron-nuclear coupling operator, ˆ Ucoup en [ΦR,χ], in the electronic equation (11), depends on the nuclear wavefunction and acts on the parametric dependence of ΦR(r,t)as a differential operator. This “pseudo-operator” includes the coupling to the nuclear subsystem beyond the parametric dependence in the BO Hamiltonian ˆ HBO(r,R). The symbol h·irindicates an integration over electronic coordinates only. The nuclear equation (12) has the particularly appealing form of a Schr¨ odinger equation that contains a time-dependent vector potential (15) and a time-dependent scalar potential (14) that govern the nuclear dynamics and yield the nuclear wavefunction. χ(R,t)is interpreted as the nuclear wavefunction since it leads to an N-body nuclear density, and an N-body current-density, which TDDFT and Quantum-Classical Dynamics 7 reproduce the true nuclear N-body density and current-density (Abedi et al (2012)) obtained from the full wavefunction Ψ(r,R,t). In order to connect the Born-Huang representation to the Exact Factorization, the electronic wavefunction ΦR(R,t)is expanded in terms of the BO states, similarly to what is done for the molecular wavefunction of Eq. (4), namely ΦR(r,t) = ∑ k Ck(R,t)ϕ(k) R(r).(16) The expansion coefficients in Eqs. (4) and (16) are related, χk(R,t) = Ck(R,t)χ(R,t),(17) by virtue of the factorization (9). Additionally, the PNC can be rewritten as ∑ k|Ck(R,t)|2=1∀R,t.(18) We point out that even in the case where the nuclear wavepacket splits into more than one BO PESs the full wavefunction is still a single product: the nuclear wavefunction has contributions (projections) on different BO PESs while the electronic wavefunction is a linear combination of the adiabatic states, but still we may write Ψ(r,R,t) = ei ¯ hS(R,t)r∑ l|χl(R,t)|2! ∑ k Ck(R,t)ϕ(k) R(r)!(19) where the first term in parenthesis is χ(R,t), with a phase S(R,t)determined by fixing the gauge freedom, and the second term in parenthesis is ΦR(r,t), using Eq. (16). In the absence of nonadiabatic couplings in Eq. (7), the evolution equations for the coefficients χk(R,t)decouple, and each nuclear contribution now evolves adiabatically according to the TDSE "Nn ∑ ν −¯ h2 2Mν ∇2 ν+ε(k) BO(R)#χk(R,t) = i¯ h∂tχk(R,t),(20) under the effect of a potential produced only by the electrons in the adiabatic state k. This is the essence of the BO approximation. Analogously, in the limit of infinite nuclear masses, Eqs. (11) and (12) (Scherrer et al (2015); Schild et al (2016); Eich and Agostini (2016); Scherrer et al (2017)) reduce to the fundamental equations of the BO approximation, namely the static electronic equation (3) and a nuclear evolution equation identical to Eq. (20) with χk(R,t)replaced by the nuclear wavefunction χ(R,t)of the Exact Factorization. If the nonadiabatic couplings cannot be neglected, the fully coupled electronnuclear problem, summarized in Eq. (7) or Eqs. (11) and (12), has to be solved. 8 E. K. U. Gross et al. Electronic dynamics is simulated at a quantum-mechanical level employing quantum-chemistry approaches, either based on the electronic wavefunction or on the electronic density. If either the adiabatic or the diabatic basis are used to characterize the electronic subsystem, electronic dynamics is implied in the time evolution of the expansion coefficients (see for instance Eqs. (4) or (16)), since the basis functions are time-independent. On the other hand, (real-time) TDDFT yields an explicit evolution of the electronic subsystem, as the electrons are represented via their time-dependent one-body density. As we will show below, real-time TDDFT can be combined with a mean-field solution of the coupled electron-nuclear dynamics. The LR formulation of TDDFT, instead, is able to provide information about the time-independent electronic properties, such as adiabatic forces and nonadiabatic couplings, needed for approaches based on the Born-Huang expansion. Possible extensions of TDDFT to solve the electronic equation of the Exact Factorization are currently under investigation (Requist and Gross (2016)). Section 3 is devoted to a thorough review of the basis of TDDFT and of LR-TDDFT Nuclear dynamics can be treated exactly or approximated at different levels, depending on the complexity of the system of interest. Simulation methods that retain the quantum character of nuclear dynamics are indeed very expensive, as the numerical cost for solving the quantum-mechanical problem scales exponentially with the number of degrees of freedom. Therefore, different strategies have been proposed over the years to make the problem numerically tractable. Quantum wavepacket propagation techniques aim at solving Eq. (7) either on grids (Lauvergnat and Nauts (2010, 2014); Sadri et al (2012)), or by expanding the nuclear wavepackets χk(R,t) on a basis where calculations are computationally cheaper (Meyer et al (1990); Burghardt et al (1999); Wang and Thoss (2003); Meyer and Worth (2003); Sadri et al (2014)). The major bottleneck of these approaches is the “pre-calculation” of the electronic properties, i.e., BO PESs and of the nonadiabatic couplings, needed to solve the nuclear equations. Attempts at solving exactly the coupled equations at the basis of the Exact Factorization are currently under investigations. On-the-fly calculations of electronic properties are instead possible, if only local nuclear information is necessary to solve (in an approximate way) Eq. (7). Full and ab initio multiple spawning methods (Mart´ ınez et al (1996); Mart´ ınez and Levine (1997); Ben-Nun and Mart´ ınez (1998); Ben-Nun et al (2000); Hack et al (2001); BenNun and Mart´ ınez (2002); Virshup et al (2008)), similarly to direct-dynamics techniques (Worth et al (2004); Lasorne et al (2006, 2007); Worth et al (2008); MendiveTapia et al (2012); Richings et al (2015)), employ a representation of the nuclear wavepackets in terms of moving Gaussian functions, that evolve along trajectories determined either variationally or classically. Trajectory-based quantum-classical schemes adopt a representation of nuclear dynamics in terms of purely classical trajectories, as in the Ehrenfest and surface-hopping methods. They are indeed numerically cheaper than the methods above, but the price to pay is sometimes the neglect of important quantum-mechanical features both in the nuclear dynamics and in the coupling between electronic and nuclear motion. Similarly to direct dynamics and full multiple spawning, evolving the nuclei along trajectories enables us to exploit the locality of classical dynamics for on-the-fly simulations, where electronic infor- TDDFT and Quantum-Classical Dynamics 9 mation is needed, and thus computed, only for the visited nuclear configurations. Trajectory-based solutions of Eqs. (11) and (12) have been proposed (Agostini et al (2014); Abedi et al (2014)), and the most recent developments (Min et al (2017)) will be reviewed in Section 4, along with Ehrenfest dynamics (Tully (1998)), trajectory surface hopping (Tully (1990); Doltsinis and Marx (2002); B¨ ockmann et al (2010); Jasper et al (2004, 2006); Subotnik et al (2013); Curchod and Tavernelli (2013b); Jaeger et al (2012); Fang and Hammes-Schiffer (1999); Tapavicza et al (2007a); Craig et al (2005); Akimov and Prezhdo (2014)) and full/ab-initio multiple spawning (Mart´ ınez et al (1996); Mart´ ınez and Levine (1997); Ben-Nun and Mart´ ınez (1998); Ben-Nun et al (2000); Hack et al (2001); Ben-Nun and Mart´ ınez (2002); Virshup et al (2008)). 3 Electronic dynamics: Time-dependent density functional theory 3.1 Time-dependent density functional theory The Hohenberg-Kohn (HK) theorem (Hohenberg and Kohn (1964)) of ground-state DFT states that knowledge of the ground-state density uniquely determines the external potential of the system (up to within a trivial constant) and thus the entire electronic Hamiltonian and the associated total ground-state energy. It is important to realize that ground-state DFT in nearly all applications is intimately tied to the BO approximation: the electronic density one calculates is the one produced by clamped nuclei. Then, by varying the positions of the clamped nuclei, ground-state DFT provides an efficient approach to map out the lowest BO PES and to calculate physical observables associated with the lowest BO PES, such as vibrational spectra, cohesive energies, barrier heights, etc. Higher BO PESs and the time evolution of systems strongly driven by external fields are not accessible wit ground-state DFT. In their seminal paper, Runge and Gross (Runge and Gross (1984)) proved a theorem that established a 1-1 correspondence between the time-dependent density and the time-dependent external potential for systems evolving from a given initial many-electron state, Φ0. The time evolution of the many-electron wavefunction is governed by the time-dependent Schr¨ odinger equation ˆ Hel(t)Φ(r,t) = i¯ h∂ ∂tΦ(r,t)(21) Φ(r,t0) = Φ0(r) with Hamiltonian ˆ Hel(t) = ˆ Te(r)+Vee(r)+vext (r,t).(22) 16 E. K. U. Gross et al. within LR-TDDFT, where the states |ϕ(0) Riand |ϕ(n) Ridescribe the ground state and nth electronic excited state wavefunctions, respectively. To achieve this goal, we will proceed by a direct comparison with the same quantity derived using manybody perturbation theory (MBPT). Therefore, we start with a short outline of the main linear-response equations in MBPT. From the definition of the retarded density-density response function χ(r,t,r0,t0) = ΠR(r,t,r0,t0) = −iθ(t−t0)hϕ(0) R|[ˆ ρ(r,t),ˆ ρ(r0,t0)]|ϕ(0) Ri hϕ(0) R|ϕ(0) Ri,(48) the change of an observable O, under the influence of a perturbation vext (r0,t0)in the linear-response regime is given by δO(t) = Z∞ 0 dt0ZdrZdr0o(r)vext (r0,t0)χ(r,t,r0,t0)(49) (here we consider an interaction of the form δvext (r0,t0) = v0(r0)E(t0)). If χdepends only on the difference (t−t0), the Fourier transform in time gives δO(ω) = ZdrZdr0o(r)v0(r0)E(ω)χ(r,r0,ω).(50) This expression can be rewritten, after a bit of algebra (Curchod et al (2013)), as a sum-over-states (SOS) formula δO(ω) = −2∑ n ωnhϕ(0) R|ˆ O|ϕ(n) Rihϕ(n) R|ˆv0E(ω)|ϕ(0) Ri ω2 n−ω2,(51) where |ϕ(n) Riand ωnare the true excitation energies and wavefunctions. Meanwhile, if we use the KS representation of LR-TDDFT as above, the change of observable is in matrix representation δO(ω) = ∑ i jσ,klτ oi jσχi jσ,klτ(ω)v0 klτE(ω),(52) where oi jσ=hφiσ|O(ω)|φjσiand v0 klτ=hφlτ|v0(r)|φkτi. Similarly, a SOS formula can also be derived for LR-TDDFT (see Refs. (Curchod et al (2013)) for a derivation), and reads δO(ω) = −2∑ n o†(A−B)1/2ZnZ† n(A−B)1/2 ω2 n−ω2v0E(ω).(53) with Znis related to the eigenvectors of Eq. (43) according to (Casida (2009)) Zn= (A−B)−1/2(Xn+Yn). Comparing the residues of LR-TDDFT response function Eq. (53) with the residues of the MBPT response function Eq. (51) at equal energy ωn, we obtain TDDFT and Quantum-Classical Dynamics 17 the following identity 1 hϕ(0) R|ˆ O|ϕ(n) Ri= (fiσ−fjσ)>0 ∑ i jσ 1 √ωn oi jσ(A−B)1/2Zni jσ.(54) This equation was derived by Casida (Casida (1995)) and then applied by Tavernelli et al. and Hu et al. for the calculation of the nonadiabatic coupling vectors between the ground state and an excited state. A similar equation was also given in Ref. (Chernyak and Mukamel (1996)). 3.3.2 The concept of auxiliary many-electron wavefunction It may be useful at this point to investigate the possibility to further simplify the definition and the calculation of matrix elements within LR-TDDFT by means of the definition of a set of ”auxiliary” multideterminantal many-electron wavefunctions based on KS orbitals. This route was first explored by Casida (Casida (1995)) to solve the assignment problem of the LR-TDDFT excited state transitions and then further developed by Tavernelli et. al (Tapavicza et al (2007b)) in relation to the calculation of matrix elements in the linear and second-order response regimes (Tavernelli et al (2009b,a, 2010)). In Ref. (Tavernelli et al (2009a)), we showed that defining the ground state manyelectron wavefunction hr1,r2,r3,...,rNel |˜ ϕ(0) Rias a Slater determinant of all occupied KS orbitals {φi}Nel i=1and the excited state wavefunction corresponding to the excitation energy ωnas hr1,r2,r3,...,rNel |˜ ϕ(n) Ri=∑ iaσrεa−εi ωn (Zn)iaσˆa† aσˆaiσhr1,r2,r3,...,rNel |˜ ϕ(0) Ri =∑ iaσ Cn iaσhr1,r2,r3,...,rNel |˜ ϕaσ R,iσi,(55) we obtain for any one-body operator of the form ˆ O=∑pqσopqσˆa† pσˆaqσ(where p,q are general indices) the correct linear-response expression for the matrix element hϕ(0) R|ˆ O|ϕ(n) Ri. Eq. (55) is derived from Eq. (54) where now the index iruns over all occupied and aover the unoccupied (virtual) KS orbitals and |˜ ϕaσ R,iσidenotes a singly-excited Slater determinant defined by the transition iσ→aσ. This theory was then successfully extended to the case of the calculation of matrix elements between two excited state wavefunctions, hϕ(n) R|ˆ O|ϕ(m) Rias will be briefly discussed in the next section on the calculation of nonadiabatic coupling vectors. 1As stated before, with no Hartree-Fock exchange contribution in the functional, (A−B)is diagonal and becomes (Casida (2009)): (A−B)iaσ,jbτ=δi,jδa,bδσ,τ(εaτ−εiτ). 18 E. K. U. Gross et al. It is important to further stress the fact that both auxiliary functions introduced above have a physical meaning only when used within LR-TDDFT for the calculation of matrix elements of the type h˜ ϕ(0) R|ˆ O|˜ ϕ(n) Riand eventually h˜ ϕ(n) R|ˆ O|˜ ϕ(m) Ri. The use of this representations of the ground state and excited state KS many-electron wavefunctions in other contexts is not justified. 3.3.3 Nonadiabatic coupling vectors within LR-TDDFT Using the concept of the auxiliary many-electron wavefunction approach described above, we can now propose an approach for the calculation of nonadiabatic vectors within LR-TDDFT. Couplings between ground and excited states We start from an alternative definition of the NACV (Epstein (1954)) (see also Ch. 5 of Ref. (Baer (2006)) for a complete discussion) between the ground (0) state and the nth excited state for a molecular system characterized by nuclear coordinates R in the configuration space (R3Nn) dγ 0n=hϕ(0) R|∂γˆ HBO|ϕ(n) Ri ε(n) BO(R)−ε(0) BO(R) (56) where γis an atomic label, ˆ HBO is the electronic Hamiltonian, and ∂γˆ HBO = ∂ˆ HBO/∂Rγ. Applying the results of the above sections on the evaluation of matrix elements of the form hϕ(0) R|ˆ O|ϕ(n) Riin LR-TDDFT to the NACV gives directly the desired expression dγ 0n= (fiσ−fjσ)>0 ∑ i jσ 1 (ωn)3/2hγ i jσ(A−B)1/2Zni jσ(57) where hγ i jσ=Rdr∂γˆ HBO φ∗ iσ(r)φjσ(r). This formula for the NACVs within LR-TDDFT was derived several times in the literature using slightly different formalisms. The first derivation was given by Chernyak and Mukamel (Chernyak and Mukamel (2000)) using a classical Liouville dynamics for the single-electron density matrix, followed by Baer (Baer (2002)). Later, Tavernelli et al. (Tapavicza et al (2007b); Tavernelli et al (2009b)) and Hu et al. (Hu et al (2007, 2008)) arrived to the same result (Eq. (57)) using the most widely used formulation based on Casida’s LR-TDDFT equations (Casida (1995)). Concerning the numerical implementation of Eq. (56) several approaches have also been proposed that differ mainly in the choice of the basis set and in the way the implicit dependence of the pseudopotentials on the nuclear positions is treated. Due to the technical nature of this subject, we will not go through the numerical TDDFT and Quantum-Classical Dynamics 19 details but better refer to the literature, which is very rich on this subject (Tavernelli et al (2009b); Hu et al (2007); Send and Furche (2010); Hu et al (2010, 2012)). Couplings between excited states LR-TDDFT only gives access strictly speaking to the couplings between ground and excited state. However, the concept of LR-TDDFT auxiliary many-electron wavefunctions can also be used as a good approximation, exact within the TammDancoff approximation, to compute couplings between excited states (Tavernelli et al (2010)), dkn. An exact derivation of these coupling terms beyond the linear response formalism of TDDFT was also proposed in the literature (Li and Liu (2014); Li et al (2014); Ou et al (2015)). However, this formalism implies the calculation of an exchange-correlation hyperkernel and leads to the critical appearance of divergences in the couplings as a result of the adiabatic approximations (Parker et al (2016)). 3.4 Nuclear forces within LR-TDDFT Excited-state dynamics using LR-TDDFT will also require the calculation of nuclear forces. Among the different approaches developed for the calculation of analytic derivatives, the Lagrangian method (Helgaker and Jørgensen (1989)) is of particular interest because of its numerical efficiency. However, the derivation of LR-TDDFT forces is technically involved and goes beyond the scope of this Chapter. We refer the interested reader to the abundant literature on the subject (Pulay (1987); Hutter (2003); Deglmann et al (2002); Rappoport and Furche (2005); Marx and Hutter (2009)). 4 Nuclear dynamics: Trajectory-based quantum-classical dynamics In this section, different approaches to nonadiabatic electron-nuclear dynamics will be presented, namely the Ehrenfest scheme (Tully (1998)), surface hopping (Tully (1990)), the coupled-trajectory mixed quantum-classical (CT-MQC) method derived from the Exact Factorization (Min et al (2015)), and full multiple spawning (Mart´ ınez et al (1996)). Their common feature is the use of trajectories to explore the nuclear configuration space, which are subject to the time-dependent effect of the electrons in the ground state as well as in the excited states. The electronic properties needed in the calculations can be determined on-the-fly based on ab initio electronic structure methods. For the purpose of this work, TDDFT and its LR formulation will be employed. Other approaches based on Bohmian trajectories are 20 E. K. U. Gross et al. also possible (Curchod et al (2011); Curchod and Tavernelli (2013a); Tavernelli (2013)) but they will not be discussed in this book chapter. Ehrenfest, surface hopping and CT-MQC are based on a purely classical description of nuclear motion, that is coupled to the quantum-mechanical evolution of the electrons. In the three approaches, a hypothesis is made to decompose the full TDSE into two coupled equations, one describing the evolution of the electronic subsystem, and the other describing the evolution of the nuclear subsystem. The main difference among them lies in the procedure followed for such decomposition. In particular, only the Exact Factorization starts from an Ansatz for the molecular wavefunction, which translates into exact coupled electronic and nuclear equations. Only in a second step, the nuclear evolution is modelled using classical trajectories. The full multiple spawning scheme, on the other hand, introduces an expansion in terms of Gaussian wavepackets to represent each nuclear coefficients χk(R,t)of the Born-Huang expansion (4). The parameters of the Gaussians are evolved classically, under the assumption that classical dynamics samples correctly the nuclear configuration space. Indeed, in the limit of an infinite number of Gaussians, full multiple spawning converges to an exact description of the electron-nuclear problem. 4.1 Ehrenfest dynamics To derive Ehrenfest decomposition, one makes the assumption that the full wavefunction can be written as a single product of a purely electronic Φ(r,t)and a purely nuclear χ(R,t)wavefunction, Ψ(r,R,t) = ei ¯ hRt 0dt0EBO(t0)Φ(r,t)χ(R,t).(58) Here, the time-dependent phase on the right-hand side is inserted to simplify the following equations derived from such an Ansatz, thus the energy EBO(t)is chosen as EBO(t) = ZdrΦ∗(r,t)i¯ h∂tΦ(r,t).(59) The product form of the molecular wavefunction in Eq. (58) is clearly uncorrelated, and in this initial Ansatz lies the fundamental approximation of the Ehrenfest scheme. When Eq. (58) in inserted into the molecular TDSE (1), the coupled equations TDDFT and Quantum-Classical Dynamics 21 ˆ Te(r)+ ˆ Vee(r)+ZdRχ∗(R,t)ˆ Vnn(R)+ ˆ Ven(r,R)χ(R,t)Φ(r,t) = i¯ h∂tΦ(r,t) (60) ∑ ν −¯ h2 2Mν ∇2 ν+ZdrΦ∗(r,t)ˆ HBO(r,R)Φ(r,t)χ(R,t) = i¯ h∂tχ(R,t) (61) are derived, by averaging over the instantaneous nuclear, in Eq. (60), and electronic, in Eq. (61), state. In both equations, the wavefunctions Φ(r,t)and χ(R,t)are supposed to be normalized. Therefore, Eq. (60) describes the evolution of the electrons in the mean field created by the nuclei, whereas the nuclei move according to Eq. (61) in the mean field of the electrons. A quantum-classical algorithm can be derived from Eqs. (60) and (61) by approximating classically the nuclear equation, that is by determining the force to propagate the nuclei as trajectories. A standard procedure can be followed, by introducing a complex-phase representation of χ(R,t), and by only considering terms O(¯ h0)in the asymptotic expansion of the complex phase in powers of ¯ h(Van Vleck (1928)). The equation for the zeroth order term S(R,t)of this expansion is thus obtained, namely ∂tS(R,t) = −∑ ν [∇νS(R,t)]2 2Mν +ZdrΦ∗(r,t)ˆ HBO(r,R)Φ(r,t).(62) This Hamilton-Jacobi-like equation can be solved via characteristics, thus yielding the expression of the classical (Ehrenfest) force as Fν(t) = −∇νZdrΦ∗(r,t)ˆ HBO(r,R)Φ(r,t).(63) The classical approximation is also introduced in the electronic evolution equation (60). Here, the nuclear density |χ(R,t)|2is approximated as a product of δfunctions centered at each time at the position of the classical nuclei, that is at the position of the classical trajectory R(I)(t). Therefore, the TDSE describing the evolution of Φ(r,t)becomes ˆ HBO r,R(I)(t)Φr,R(I)(t),t=i¯ h∂tΦr,R(I)(t),t.(64) The electronic wavefunction acquires an implicit dependence on the nuclear positions, expressed as the classical trajectory, via the dependence of the BO Hamiltonian on R(I)(t). The trajectory Iof the nucleus νis determined by solving Newton’s equation with force F(I) ν(t) = ZdrΦ∗(r,R(I)(t),t)h−∇νˆ HBO(r,R(I)(t))iΦ(r,R(I)(t),t),(65) 22 E. K. U. Gross et al. where now a label (I)has been introduced to show that along a trajectory, Eqs. (64) and (65) have to be evolved consistently. Multiple trajectories can also be employed, to “wash out” some of the details of the coherent evolution along a single trajectory. Nuclear and electronic observables can thus be determined as averages over this ensemble of trajectories. As it provides the true time-dependent electronic density, TDDFT can be used within an Ehrenfest dynamics scheme to perform nonadiabatic molecular dynamics. The mapping of the nuclear equation (Eq. (65)) into the DFT formalism is straightforward and only requires the description of the forces h−∇νˆ HBO(r,R(I)(t))ias a functional of the time-dependent density ρ(r,t). If we replace the expectation value of the electronic Hamiltonian with the DFT energy evaluated with the exchangecorrelation potential vxc[ρ]|ρ(r)←ρ(r,t), the gradient with respect to the nuclear coordinates can be performed analytically as in the case of the adiabatic BO dynamics and the Car-Parrinello (Car and Parrinello (1985)) molecular dynamics schemes (Marx and Hutter (2009)). 4.1.1 Application of Ehrenfest dynamics combined with TDDFT As Ehrenfest dynamics gives a direct access to electronic dynamics, it is a method of choice to investigate the dynamics of the electronic density and subsequent nuclear dynamics after a strong perturbation. Such perturbation can be induced by the action of an external light pulse or through the collision with a highly-charged particle, generating either an electronic excitation or, in the some other cases, electron abstraction (Tavernelli et al (2005); Tavernelli (2006); Castro et al (2004); Li et al (2005); Yagi and Takatsuka (2005); Andrade et al (2009); Moss et al (2009); Liang et al (2010); Gaigeot et al (2010); Lopez-Tarifa et al (2011); Elliott and Maitra (2012); Tavernelli (2015)). The latter takes place when an XUV attosecond pulse interacts with a molecule and leads to a core ionization. In a Born-Huang picture, such an ultrafast ionization leads to the generation of an electronic wavepacket, i.e., the generation of a coherent superposition of different nuclear contributions on a large number of electronic states (the number of electronic states being considered depends on the bandwidth of the ionizing pulse). Ehrenfest dynamics combined with TDDFT offers an alternative to the Born-Huang picture by only requiring the generation of an initial electronic density to represent the initial ionized state. Mart´ ın et al. employed this strategy to study the role of nuclear motion in the electronic dynamics upon XUV ultrafast ionization of a small amino acid, glycine (Lara-Astiaso et al (2017)). The one-electron ionization generated by the sub-300-as XUV pulse generates an electronic wavepacket that can be described by a coherent superposition of more than ten one-hole states, in an energy domain ranging from 17 to 35 eV (the pulse bandwidth). The electronic density corresponding to this electronic wavepacket, ρ(r,t0), is used as initial condition for two simulations: (i) real-time TDDFT with frozen nuclei and (ii) real-time TDDFT combined with Ehrenfest dynamics. As a result of the nature of the electronic wavepacket, the time-evolution of the unpaired electron (with respect to the initial density) shows that the electron TDDFT and Quantum-Classical Dynamics 23 migrates over the entire molecular scaffold with a dynamics that is characterized by only few, system dependent, frequencies (Fig. 1). This is the behavior that one would expect in the case the ionized electron is removed from one given localized orbital. Comparing the two panels of Fig. 1, we observe that nuclear motion starts altering the electronic dynamics already after the first 10 fs of dynamics, emphasizing the importance of including nuclear dynamics in such simulations. However, it is important to note that the mean-field character of Ehrenfest dynamics might hamper a more detailed study of the electronic wavepacket dynamics, in particular due to the underestimation of decoherence and dephasing effects at longer time scales (Vacher et al (2017)). The ease of the Ehrenfest formalism combined with the efficiency of TDDFT offer nevertheless a valid tool for the study of the short-time electronic wavepacket dynamics in molecular systems. Fig. 1: Spin density differences at different times after interaction of a XUV attosecond pulse with the glycine molecule. The initial conditions correspond to a geometry obtained after thermalization at 100 K. Adapted from Chemical Physics Letters, 683, M. Lara-Astiaso, A. Palacios, P. Decleva, I. Tavernelli, F. Mart´ ın, Role of electron-nuclear coupled dynamics on charge migration induced by attosecond pulses in glycine, 357, Copyright (2017), with permission from Elsevier. 4.2 Surface hopping Surface-hopping decomposition is derived under the preliminary assumption that the nuclei evolve along classical trajectories. Therefore, the BO Hamiltonian acquires an implicit time dependence via its dependence on the nuclear coordinates. A TDSE is proposed in this way, namely ˆ HBO(R(I)(t))ΦR(I)(t)(r,t) = i¯ h∂tΦR(I)(t)(r,t),(66) for the electronic wavefunction, which is, itself, dependent on the classical trajectories. As done above, a classical trajectory is labeled by the index (I). An expansion in the adiabatic basis is introduced for ΦR(I)(t)(r,t), 24 E. K. U. Gross et al. ΦR(I)(t)(r,t) = ∑ k Ck(t)ϕ(k) R(I)(t)(r),(67) and Eq. (66) yields ˙ C(I) k(t) = −i ¯ hε(k) BO R(I)(t)C(I) l(t)−∑ l C(I) l(t) Nn ∑ ν=1 P(I) ν(t) Mν·dν,kl R(I)(t)(68) Here, the BO PES ε(k) BO (R)and the nonadiabatic coupling vectors, i.e., dν,kl (R) = Zdrϕ(k) R∗(r)∇νϕ(l) R(r) = Dϕ(k) R∇νϕ(l) REr,(69) which are functions of the nuclear coordinates, are evaluated at the instantaneous positions along the trajectories, they thus become functions of the trajectory itself. Henceforth, a superscript (I)will be introduced to indicate this dependence on the trajectory. The surface-hopping scheme takes its name from the idea suggested for the evolution of the classical nuclear trajectories, namely that a trajectory evolves according to one adiabatic BO force, determined as (minus) the gradient of the BO potential energy surface (PES), until a stochastic hop occurs onto another BO PES. The classical (surface-hopping) force can then be written as F(I) ν(t) = −∇νε∗ BO,(70) with the symbol ∗indicating that the force-state is selected stochastically at each time step. The discontinuity in the force, and thus in the potential energy, for a given trajectory, is compensated by a discontinuity in the velocity, and thus in the kinetic energy, that guarantees energy conservation. The hopping scheme fewestswitches (Tully (1990)) prescribes that the trajectory Ihops from surface kto surface laccording to the probability Pk→l=max  0,−2dt C(I) k(t) 2ℜhC(I) k∗(t)C(I) l(t)i∑ ν P(I) ν(t) Mν·d(I) lk,ν  ,(71) with dt the integration time step. The major drawback of the surface-hopping scheme is the (over)coherent evolution of the electronic coefficients coupled to the classical (independent) trajectories. The issue has been well-documented in the literature (Subotnik et al (2013); Bittner and Rossky (1995); Curchod and Tavernelli (2013b); Gao and Thiel (2017)), and several schemes have been proposed (Shenvi et al (2011b,a); Shenvi and Yang (2012); Subotnik and Shenvi (2011b,a); Jaeger et al (2012); Jasper and Truhlar (2007); Granucci and Persico (2007)) to cure or alleviate this shortcoming. TDDFT and Quantum-Classical Dynamics 25 4.2.1 Application of Surface Hopping combined with LR-TDDFT Surface hopping has been used to study a large number of excited-state mechanisms, and we refer the interested reader to specialized reviews Barbatti (2011); Curchod et al (2013); Persico and Granucci (2014) for a list of applications. The application presented here highlights the combination of surface hopping with LRTDDFT (using the concepts developed in Sec. 3.2), including implicitly and explicitly the role of spin-orbit coupling as well as explicit solvent effects. Ruthenium (II) trisbipyridine, [Ru(bpy)3]2+, is an inorganic molecule recognized for its extremely efficient intersystem crossing process, i.e., when the molecule changes, in this particular case, from a singlet electronic state to a triplet electronic state (Cannizzo et al (2006); Gawelda et al (2006)). [Ru(bpy)3]2+ is initially photoexcited in a singlet metal-to-ligand-charge-transfer (1MLCT) state before it rapidly relaxes among other 1MLCT or 3MLCT, as a result of the high density of states; the overall dynamics to the triplet states has been observed experimentally in water within a ∼50 fs timescale (Fig. 2). In the first theoretical study (Tavernelli et al (2011)), the excited-state dynamics of the [Ru(bpy)3]2+ in water was studied by employing surface hopping with LRTDDFT, in a QM/MM formalism where water molecules were treated classically. Intersystem-crossing events were analyzed a posteriori, monitoring the crossings between singlet and triplet states and evaluating spin-orbit coupling from qualitative rules. Owing to the cost of the overall dynamics, this study was limited to only two trajectories. Nevertheless, both trajectories indicated an ultrafast decay of the molecule towards triplet states in less than 50 fs, in good correlation with experimental evidences. The MLCT character of the different excited states implies that an electron moves from the central metal to one (or two, depending on the state) solvent-exposed bipyridine ligands. Hence, the simulation showed that water molecules in the first solvation shell can rapidly rearrange in a non-diffusive rotation around the hydrogen bond axis to stabilize an extra charge located on a close ligand. An explicit treatment of solvent molecules is central to capture such effects as well as a proper ordering of the different electronic states. In a more recent study (Atkins and Gonz´ alez (2017)), surface hopping combined within LR-TDDFT was used to simulate the excited-state dynamics of [Ru(bpy)3]2+ in gas phase, but with the explicit treatment of spin-orbit coupling in a perturbative ZORA formalism Wang and Ziegler (2005) and a larger number (101) of trajectories. This study confirmed the ultrafast decay of the original 1MLCT population towards triplet states, already at the early time of the dynamics. Horizontal intersystem crossing processes were observed, followed by ultrafast nonadiabatic dynamics among the triplet states (Fig. 3). Thanks to normal-mode and principal component analysis, the authors could identify that the motion of both the ruthenium and the coordinated nitrogens is activated, even within such short timescale, leading potentially to the intersystem crossing events. 32 E. K. U. Gross et al. to the equivalent left-open ring structure (ii). In Fig. 5 (upper panels) we show the electronic populations |C(I) k(t)|2as functions of time and for the selected trajectories, along with the energy profiles (lower panels) of the three adiabatic states considered here and the gauge-invariant (GI) part of the TDPES (the first two terms on the right-hand side of Eq. (14)). The TDPES provides information about the “active” electronic state: if one wants to connect the interpretation of the dynamics based on the Exact Factorization to the standard perspective in terms of wavepackets evolving “on” different adiabatic surfaces, it is instructive to compare the TDPES with the BO PES, as done in Fig. 5 (lower panels). The upper panels of Fig. 5 confirm that the coupling region is encountered by trajectories of type (iii) (right panels) at later times if compared to trajectories of type (i) or (ii) (left panels). Additionally, the S2/S1population exchange is very sharp for type (i), and smooth for (iii). Observing the TDPES, we can argue that after about 5 fs, trajectories (i) encounter a steep S2potential, that directs them towards the conical intersection. Trajectories (iii) are trapped in a region of flat potential, that prevents them from a fast de-excitation to S1. Subsequently, very different paths are undertaken, and thus different region of the S1PES are explored. Towards the end of the simulated dynamics, only trajectories of type (i) are expected to relax to the ground state S0, as confirmed by the closing of the energy gap between S0and S1at 25 fs (lower left panel of Fig. 5). At this time, the trajectory of type (iii) is evolving on a portion of the BO PES S1that is located at about 4 eV from the ground state. 4.4 Full and Ab Initio Multiple Spawning 4.4.1 Full Multiple Spawning Full multiple spawning (FMS) (Mart´ ınez et al (1996); Mart´ ınez and Levine (1997); Ben-Nun and Mart´ ınez (1998); Ben-Nun et al (2000); Hack et al (2001); Ben-Nun and Mart´ ınez (2002); Virshup et al (2008)) proposes to expand the nuclear amplitudes in the Born-Huang expansion in a linear combination of frozen multidimensional Gaussian functions. (Heller (1981)) But these Gaussians functions do not form a static grid; on the contrary, they are evolving over time in both position and momentum space to better adapt to the evolution of the nuclear wavefunctions, forming a set of Trajectory Basis Functions (TBFs). Ψ(r,R,t) = ∑ k NT BFs,k ∑ I C(k) I(t)˜ χ(k) IR;R(k) I(t),P(k) I(t),γ(k) I(t),αϕ(k) R(r),(80) where C(k) I(t)is the complex coefficient for the TBF Ievolving on electronic state (k)(used here as a label) and ˜ χ(k) IR;R(k) I(t),P(k) I(t),γ(k) I(t),αis the travelling multidimensional Gaussian Ion state (k)with mean position R(k) I(t), momenta TDDFT and Quantum-Classical Dynamics 33 P(k) I(t), phase γ(k) I(t), and frozen width α. In FMS, the TBFs follow classical trajectories, i.e., R(t)and P(t)are propagated according to Hamilton’s equation of motion (and the phase is integrated semi-classically). We note that this classical propagation of the TBFs does not imply that the method is semiclassical in itself, as the TBFs are only a support for the propagation of the nuclear wavefunctions. Indeed, in the limit of a large number of TBFs (NTBFs), the FMS expansion would be exact. In fact, in the limit of an infinite number of Gaussian functions, their dynamics is redundant and we have a (infinitely fine) grid. We note that other methods were proposed where the TBFs follow Ehrenfest trajectories (Shalashilin (2009, 2010); Saita and Shalashilin (2012); Makhov et al (2017)) (multiconfiguration Ehrenfest, MCE) or quantum trajectories (Worth et al (2004); Lasorne et al (2006, 2007); Worth et al (2008); Mendive-Tapia et al (2012); Richings et al (2015)) (variational Multiconfiguration Gaussian, vMCG), as well as mixed strategies (Makhov et al (2014); Meek and Levine (2016); Izmaylov and Joubert-Doriol (2017); Joubert-Doriol et al (2017)). One can express the time-dependent Schrdinger equation (1) in the basis of TBFs by inserting Eq. (80) in the former, left multiplying by h˜ χ(k) JR;R(l) J(t),P(l) J(t),γ(l) J(t),αϕ(l) R(r)i∗ and integrating over both nuclear and electronic coordinates, leading, in atomic units, to (Ben-Nun and Mart´ ınez (2002)): d dt Cl(t) = −i(S−1 ll )"Hll −i˙ SllCl+∑ k6=l HlkCk#.(81) for each electronic state lconsidered. The nonorthonormality of the Gaussian basis result in overlap matrices Sll and ˙ Sll, with elements (Sll)JI =h˜ χ(l) J|˜ χ(l) IiRand ˙ SllJI =h˜ χ(l) J|∂ ∂t|˜ χ(l) IiR. We note that in Eq. (81) Slk =˙ Slk =0, ∀l6=k, due to the orthonormality of the electronic basis. As mentioned earlier, the trajectories in FSSH are uncoupled. This is not the case in FMS and TBFs are coupled thanks to the Hamiltonian matrix Hin Eq. (81). Let us consider the Hamitonian matrix element between two TBFs Jand Ievolving in adiabatic electronic states: HIJ kl =h˜ χ(k) I|ˆ Tnuc|˜ χ(l) JiRδkl +h˜ χ(k) I|ε(k) BO|˜ χ(l) JiRδkl −h˜ χ(k) I| 3N ∑ ρ=1 1 Mρhϕ(k) R|∂ ∂Rρ|ϕ(l) Rir ∂ ∂Rρ|˜ χ(l) JiR −h˜ χ(k) I| 3N ∑ ρ=1 1 2Mρhϕ(k) R|∂2 ∂R2 ρ|ϕ(l) Rir|˜ χ(l) JiR(82) 34 E. K. U. Gross et al. If the two TBFs are in the same electronic state k, they will be coupled via the first two terms on the r.h.s of Eq. (82): the nuclear kinetic energy operator and the (adiabatic) electronic energy. If the two TBFs belong to different electronic states kand l, then a term containing the nonadiabatic coupling vectors (third term on the r.h.s) and the second-order nonadiabatic couplings (fourth term on the r.h.s of Eq. (82)) will ensure their coupling. We note that the diagonal second-order nonadiabatic couplings (also known as “Diagonal Born-Oppenheimer Correction”) will contribute an additional intrastate coupling. FMS can easily incorporate additional source of couplings between TBFs like spin-orbit coupling (Curchod et al (2016b)), tunneling effects (Ben-Nun and Mart´ ınez (2000)), or an external electromagnetic field (Mignolet et al (2016)), for example. 4.4.2 Spawning algorithm Up to this point, we have discussed the formal equations of FMS when a large number of TBFs is considered. However, FMS, as its name indicates, proposes to replace the large number of trajectory basis functions by an algorithm, coined the spawning algorithm, that will adapt the size of the basis set dynamically during the simulation. In other words, the spawning algorithm proposes the following alteration: NT BFs →NT BFs(t). In short, a TBF is evolving on a given PES and, as soon as a sizable coupling with a different electronic state is recorded, a spawning mode is triggered and will determine if, where, and when a new TBF should be created (“spawned”) on the coupled state to maximize the coupling with the existing TBF and, therefore, the description of nonadiabatic transitions. Different versions of the spawning algorithm have been proposed and the interested reader is encouraged to read Refs. (Ben-Nun and Mart´ ınez (2002); Yang et al (2009)) for more information. It is perhaps important to note at this stage that the time-dependence of the number of TBFs, NT BFs(t), implies that the size of matrices and vectors in Eq. (81) will vary during the dynamics. 4.4.3 Ab Initio Multiple Spawning FMS is in principle exact. However, the coupling between TBFs given in Eq. (82) implies integration, and therefore knowledge, of electronic structure properties like PESs or nonadiabatic couplings over the full nuclear configuration space. Approximations are, therefore, needed to treat molecular systems in their full dimensionality. Let us start by performing a Taylor expansion of the electronic quantity of interest (here the electronic energy for example) around the centroid position of two TBFs Jand Ievolving in electronic state k:R(kk) JI =R(k) J+R(k) I 2: TDDFT and Quantum-Classical Dynamics 35 ε(k) BO(R) = ε(k) BO(R(kk) JI )+ 3N ∑ ρ (Rρ−R(kk) ρ,JI)∂ε(k) BO(R) ∂RρRρ=R(kk) ρ,JI +1 2 3N ∑ ρρ0 (Rρ−R(kk) ρ,JI)∂2ε(k) BO(R) ∂Rρ∂Rρ0Rρ=R(kk) ρ,JI ,Rρ0=R(kk) ρ0,JI (Rρ0−R(kk) ρ0,JI)+... Owing to the locality of Gaussian functions, one will consider that, to a good approximation, only the term of order zero can be retained; in other words ε(k) BO(R)≈ ε(k) BO(R(kk) JI ). (Ben-Nun and Mart´ ınez (2002)) This approximation, called the saddlepoint approximation of order zero (SPA0), strongly simplifies the coupling between TBFs as the Hamiltonian matrix elements becomes HIJ kl =h˜ χ(k) I|ˆ Tnuc|˜ χ(l) JiRδkl +ε(k) BO(R(kk) IJ )h˜ χ(k) I|˜ χ(l) JiRδkl − 3N ∑ ρ=1 1 Mρhϕ(k) R|∂ ∂Rρ|ϕ(l) Rir|Rρ=R(kl) ρ,IJ h˜ χ(k) I|∂ ∂Rρ|˜ χ(l) JiR(83) (note that we also dropped the terms depending on second-order couplings). Hence, calculating the coupling between TBF Jand Ionly require the extra calculation of electronic structure quantities at their mutual centroid position, which can easily be achieved in an on-the-fly dynamics scheme. The SPA0 allows to port the in principle exact framework of FMS to the nonadiabatic dynamics of molecules. The resulting nonadiabatic method within the SPA0 is often called Ab Initio Multiple Spawning (AIMS)2. AIMS has been coupled with different level of electronic structure like SA-CASSCF (Levine et al (2008)), MS-CASPT2 (Tao et al (2009)), FOMO-CASCI (Pijeau et al (2017)), or LR-TDDFT (Curchod et al (2017)). 4.4.4 Applications AIMS has recently (Curchod et al (2017)) been interfaced with an implementation of LR-TDDFT accelerated by graphical processing units (GPUs) (Isborn et al (2011)), offering an important speed-up for the calculation of electronic energies, analytical gradients, and nonadiabatic coupling terms. The combination of AIMS and GPUaccelerated LR-TDDFT was employed to shed light on the excited-state dynamics of 4-N,N’-dimethylaminobenzonitrile (DMABN). DMABN is a molecule known to exhibit dual fluorescence depending on its environment, and it was proposed that the nature of the two emitting states is correlated with a twist of the dimethylamino (DMA) group Grabowski et al (2003) (see molecular structure in the inset of Fig. 6). However, DMABN is photoexcited into its second excited state, S2, and relaxes to the first excited state S1where emission will occur at a later time. One question is 2We note that an additional approximation is commonly employed within AIMS – the independent first generation approximation – that approximates the initial nuclear wavepacket at time t=0 by a set of independent parent TBFs, i.e., coupling between parent TBFs is neglected. Couplings between all the TBFs spawned by each parent TBF are preserved.(Ben-Nun and Mart´ ınez (2002)) 36 E. K. U. Gross et al. Fig. 6: Left panel: Population traces (with standard errors) of the different electronic states considered in the AIMS/LR-TDDFT dynamics. The grey dashed line indicates the number of TBFs during the dynamics. Right panel: Twist angle of the DMA group for the entire swarm of TBFs. The line width is proportional to the TBF population. Inset: representation of the DMABN molecule, with the DMA group highlighted in red. Adapted with permission from Curchod BFE , Sisto A, Mart´ ınez TJ (2016) Ab initio multiple spawning photochemical dynamics of DMABN using GPUs. The Journal of Physical Chemistry A 121(1):265-276. Copyright (2016) American Chemical Society. then: does the S2/S1nonadiabatic dynamics imply a twist of the DMA group, which could potentially have an influence on the fate of the molecule at later time in the S1state. Early static calculations have predicted that the S2/S1transfer should be fast, and that a twist of the DMA was not required for the nonadiabatic transition to occur (G´ omez et al (2005)). AIMS combined with GPU-accelerated LR-TDDFT confirmed that the S2population decays to S1in less than <50 fs (blue line, left panel of Fig. 6). Also, the transfer of the nuclear wavepacket to S1is not correlated with the torsion of the DMA group, as showed by the projection of the TBFs on the twist coordinate (right panel of Fig. 6). The DMA group in fact starts its rotation only after the molecule reached S1(t>50 fs, see right panel of Fig. 6), and such a transfer implies a change in the electronic character of the wavepacket from charge transfer in S2around the Franck-Condon region to locally excited after the nonadiabatic transfer to S1. On a more technical note, the grey dashed line on the left panel of Fig. 6 indicates the number of TBFs during the AIMS dynamics. Starting from 21 parent TBFs, the spawning algorithm rapidly creates a large amount of child TBFs to ensure a good basis for the propagation in the different nonadiabatic regions encountered. TDDFT and Quantum-Classical Dynamics 37 5 Conclusions The aim of this chapter has been to provide a broad overview of quantum molecular dynamics methods to simulate nonadiabatic phenomena in isolated and condensedphase systems, adopting a quantum-classical perspective. The combination of classical and quantum-mechanical approaches is, in fact, capable to describe, accurately and efficiently, dynamical processes involving electronic excited states. The assumption that a good description of the nuclear dynamics can be achieved based on classical mechanics is based on the fact that at the molecular scale atoms (usually heavier than a proton) are associated to a relatively short de Broglie wavelength, e.g., shorter than the typical scale of variation of the potential. Electrons, on the other hand, require a purely quantum-mechanically description, which in this work has been achieved using the framework of time-dependent density functional theory. The challenges faced in the development of theoretical and numerical approaches for nonadiabatic dynamics are copious, justifying the rise over the years of a multitude of strategies to tackle these type of problems, some of which have been described in this review chapter. TDDFT has emerged as a very powerful method to describe electronic excitedstate dynamics, both for molecular systems and in condensed phase. TDDFT describes the evolution of the electronic density in a time-dependent external potential and within the linear response regime, LR-TDDFT allows the calculation of excitedstate properties, e.g., energies, forces and nonadiabatic couplings at a modest computational cost. The straightforward combination of the fundamental TDDFT theorem (leading to the time-dependent Kohn-Sham equations) with the classical motion of the nuclei gives rise to the conceptually simple Ehrenfest dynamics scheme. We described some of its applications to molecular processes involving explicit timedependent electronic wavepacket dynamics, pointing out its limitations related to the mean-field character of the underlying approximations. The strong interplay of electronic and nuclear motion is, however, highly nontrivial, and effects beyond the mean-field approximations are difficult to capture with trajectory-based approaches. In particular, the accurate description of the electronic “nonadiabatic effects” on the nuclear dynamics requires the derivation of a suitable theoretical framework that enables the separation of electronic and nuclear motions. Building on the TDDFT or LR-TDDFT description of the electronic dynamics within either the Born-Huang or the Exact-Factorization framework we have then reviewed three quantum-classical nonadiabatic simulation techniques that extend beyond the mean-field Ehrenfest approach, namely surface hopping, coupledtrajectory mixed quantum-classical dynamics and multiple spawning, along with some of their recent applications. The common denominator of all these techniques is that they combine “on-the-fly” trajectory-based dynamics with the computational advantages of TDDFT. In the past years, large progresses have thus been accomplished in mixed quantumclassical nonadiabatic dynamics and we hope that this chapter will help newcomers to engage in this exciting field of research. We believe that TDDFT-based nonadiabatic dynamics can become the method of choice for treating photochemical 38 E. K. U. Gross et al. and photophysical processes of complex molecular systems in the gas and condensed phases. While it is hard to predict the outcome of current and future developments, the improvement of the exchange-correlation functionals, for instance going beyond the adiabatic approximation, is one of the key issues that should be addressed, aiming to get more reliable excited-state energies, nonadiabatic couplings, nuclear forces and electronic dissipation effects.Finally, the development of coupling schemes to describe the interaction with electromagnetic fields described quantum mechanically are underway (Flick et al (2017b); Dimitrov et al (2017); Flick et al (2017a)) and will open new research areas in the field of cavity quantum electrodynamics. References Abedi A, Maitra NT, Gross EKU (2010) Exact factorization of the time-dependent electron-nuclear wave function. Phys Rev Lett 105(12):123,002 Abedi A, Maitra NT, Gross EKU (2012) Correlated electron-nuclear dynamics: Exact factorization of the molecular wave-function. J Chem Phys 137(22):22A530 Abedi A, Agostini F, Suzuki Y, Gross EKU (2013a) Dynamical steps that bridge piecewise adiabatic shapes in the exact time-dependent potential energy surface. Phys Rev Lett 110(26):263,001 Abedi A, Maitra NT, Gross EKU (2013b) Reply to comment on “Correlated electron-nuclear dynamics: Exact factorization of the molecular wave-function”. J Chem Phys 139(8):087,102 Abedi A, Agostini F, Gross EKU (2014) Mixed quantum-classical dynamics from the exact decomposition of electron-nuclear motion. Europhys Lett 106(3):33,001 Adamo C, Jacquemin D (2013) The calculations of excited-state properties with time-dependent density functional theory. Chemical Society Reviews 42(3):845–856 Agostini F, Abedi A, Suzuki Y, Gross EKU (2013) Mixed quantum-classical dynamics on the exact time-dependent potential energy surfaces: A novel perspective on non-adiabatic processes. Mol Phys 111(22-23):3625 Agostini F, Abedi A, Gross EKU (2014) Classical nuclear motion coupled to electronic nonadiabatic transitions. J Chem Phys 141(21):214,101 Agostini F, Abedi A, Suzuki Y, Min SK, Maitra NT, Gross EKU (2015a) The exact forces on classical nuclei in non-adiabatic charge transfer. J Chem Phys 142(8):084,303 Agostini F, Min SK, Gross EKU (2015b) Semiclassical analysis of the electron-nuclear coupling in electronic non-adiabatic processes. Ann Phys 527(9-10):546–555 Agostini F, Min SK, Abedi A, Gross EKU (2016) Quantum-classical non-adiabatic dynamics: Coupledvs. independent-trajectory methods. J Chem Theory Comput 12(5):2127–2143 Agostini F, Tavernelli I, Ciccotti G (2018) Nuclear quantum effects in electronic (non)adiabatic dynamics. submitted x:xxx Akimov AV, Prezhdo OV (2014) Advanced capabilities of the PYXAID program: Integration schemes, decoherence effects, multiexcitonic states, and field-matter interaction. Journal of Chemical Theory and Computation 10:789 Alonso JL, Clemente-Gallardo J, Echeniche-Robba P, Jover-Galtier JA (2013) Comment on “Correlated electron-nuclear dynamics: Exact factorization of the molecular wave-function”. J Chem Phys 139:087,101 Andrade X, Castro A, Zueco D, Alonso J, Echenique P, Falceto F, Rubio A (2009) Modified Ehrenfest formalism for efficient large-scale ab initio molecular dynamics. J Chem Theory Comput 5(4):728–742 TDDFT and Quantum-Classical Dynamics 39 Atkins AJ, Gonz´ alez L (2017) Trajectory surface-hopping dynamics including intersystem crossing in [ru (bpy) 3] 2+. The Journal of Physical Chemistry Letters 8(16):3840–3845 Baer M (2006) Beyond Born-Oppenheimer: Electronic Nonadiabatic Coupling Terms and Conical Intersections. John Wiley & Sons, Inc. Baer R (2002) Non-adiabatic couplings by time-dependent density functional theory. Chem Phys Lett 364:75–79 Barbatti M (2011) Nonadiabatic dynamics with trajectory surface hopping method. WIREs Comput Mol Sci 1:620–633 Ben-Nun M, Mart´ ınez TJ (1998) Nonadiabatic molecular dynamics: Validation of the multiple spawning method for a multidimensional problem. J Chem Phys 108:7244–7257 Ben-Nun M, Mart´ ınez TJ (2002) Ab initio quantum molecular dynamics. Advances in Chemical Physics 121:439–512 Ben-Nun M, Mart´ ınez TJ (2000) A multiple spawning approach to tunneling dynamics. The Journal of Chemical Physics 112(14):6113–6121 Ben-Nun M, Quenneville J, Mart´ ınez TJ (2000) Ab initio multiple spawning: Photochemistry from first principles quantum molecular dynamics. J Phys Chem A 104:5161–5175 Bittner ER, Rossky PJ (1995) Quantum decoherence in mixed quantum-classical systems: Nonadiabatic processes. J Chem Phys 103:8130 B¨ ockmann M, Doltsinis N, Marx D (2010) Unraveling a chemically enhanced photoswitch: Bridged azobenzene. Angew Chemie Int Ed 49:3382 Bonella S, Coker DF (2005) LAND-map, a linearized approach to nonadiabatic dynamics using the mapping formalism. J Chem Phys 122:194,102–13 Burghardt I, Meyer HD, Cederbaum LS (1999) Approaches to the approximate treatment of complex molecular systems by the multiconfiguration time-dependent Hartree method. J Chem Phys 111:2927 Cannizzo A, van Mourik F, Gawelda W, Zgrablic G, Bressler C, Chergui M (2006) Broadband femtosecond fluorescence spectroscopy of [Ru(bpy)3]2+. Angew Chem Int Ed 45:3174–3176 Car R, Parrinello M (1985) Unified approach for molecular dynamics and density-functional theory. Phys Rev Lett 55:2471 Casida M, Huix-Rotllant M (2012) Progress in time-dependent density-functional theory. Annu Rev Phys Chem 63(1):287–323, http://www.annualreviews.org/doi/pdf/10. 1146/annurev-physchem-032511-143803 Casida ME (1995) Time-dependent density-functional response theory for molecules. In: Chong DP (ed) Recent Advances in Density Functional Methods, Singapore, World Scientific, p 155 Casida ME (2009) Time-dependent density-functional theory for molecules and molecular solids. J Mol Struc (Theochem) 914(1-3):3–18 Casida ME, Gutierrez F, Guan J, Gadea FX, Salahub D, Daudey JP (2000) Charge-transfer correction for improved time-dependent local density approximation excited-state potential energy curves: Analysis within the two-level model with illustration for H2O and LiH. J Chem Phys 113:7062 Castro A, Marques MAL, Rubio A (2004) Propagators for the time-dependent Kohn-Sham equations. J Chem Phys 121(8):3425–3433, DOI 10.1063/1.1774980 Cave R, Zhang F, Maitra N, Burke K (2004) A dressed TDDFT treatment of the 21Agstates of butadiene and hexatriene. Chem Phys Lett 389(1):39–42 Chernyak V, Mukamel S (1996) Size-consistent quasiparticle representation of nonlinear optical susceptibilities in many-electron systems. J Chem Phys 104(2):444–459, DOI 10.1063/1. 470843, URL http://link.aip.org/link/?JCP/104/444/1 Chernyak V, Mukamel S (2000) Density-matrix representation of nonadiabatic couplings in timedependent density functional (TDDFT) theories. J Chem Phys 112:3572–3579 Cordova F, Doriol LJ, Ipatov A, Casida ME, Filippi C, Vela A (2007) Troubleshooting timedependent density-functional theory for photochemical applications: Oxirane. J Chem Phys 127:164,111 Craig CF, Duncan WR, Prezhdo OV (2005) Trajectory surface hopping in the time-dependent kohn-sham approach for electron-nuclear dynamics. Phys Rev Lett 95:163,001 40 E. K. U. Gross et al. Curchod BFE, Agostini F (2017) On the dynamics through a conical intersection. J Phys Chem Lett 8:831 Curchod BFE, Tavernelli I (2013a) On trajectory-based nonadiabatic dynamics: Bohmian dynamics versus trajectory surface hopping. J Chem Phys 138:184,112 Curchod BFE, Tavernelli I (2013b) On trajectory-based nonadiabatic dynamics: Bohmian dynamics versus trajectory surface hopping. J Chem Phys 138:184,112 Curchod BFE, Tavernelli I, Rothlisberger U (2011) Trajectory-based solution of the nonadiabatic quantum dynamics equations: an on-the-fly approach for molecular dynamics simulations. Phys Chem Chem Phys 13:3231–3236 Curchod BFE, Rothlisberger U, Tavernelli I (2013) Trajectory-based nonadiabatic dynamics with time-dependent density functional theory. ChemPhysChem 14(7):1314–1340 Curchod BFE, Agostini F, Gross EKU (2016a) An exact factorization perspective on quantum interferences in nonadiabatic dynamics. J Chem Phys 145:034,103 Curchod BFE, Rauer C, Marquetand P, Gonz´ alez L, Mart´ ınez T (2016b) Communication: Gaims– generalized ab initio multiple spawning for both internal conversion and intersystem crossing processes. The Journal of Chemical Physics 144(10):101,102 Curchod BFE, Sisto A, Mart´ ınez TJ (2017) Ab initio multiple spawning photochemical dynamics of dmabn using gpus. The Journal of Physical Chemistry A 121(1):265–276 Dancoff SM (1950) Phys Rev 78:382 Deglmann P, Furche F, Ahlrichs R (2002) An efficient implementation of second analytical derivatives for density functional methods. Chem Phys Lett 362(5-6):511 – 518, DOI 10. 1016/S0009-2614(02)01084-9, URL http://www.sciencedirect.com/science/ article/pii/S0009261402010849 Dimitrov T, Flick J, Ruggenthaler M, Rubio A (2017) Exact functionals for correlated electronphoton systems. New J Phys 19:113,036 Dobson JF, B¨ unner MJ, Gross EKU (1997) Time-dependent density functional theory beyond linear response: An exchange-correlation potential with memory. Phys Rev Lett 79(10):1905 Doltsinis NL, Marx D (2002) Nonadiabatic Car-Parrinello molecular dynamics. Phys Rev Lett 88:166,402 Dreuw A, Head-Gordon M (2004) Failure of time-dependent density functional theory for long-range charge-transfer excited states: The zincbacteriochlorin-bacteriochlorin and bacteriochlorophyll-spheroidene complexes. J Am Chem Soc 126:4007–4016 Dreuw A, Head-Gordon M (2005) Single-reference ab initio methods for the calculation of excited states of large molecules. Chem Rev 105:4009 Dreuw A, Weisman J, Head-Gordon M (2003) Long-range charge-transfer excited states in timedependent density functional theory require non-local exchange. J Chem Phys 119:2943 Dunkel ER, Bonella S, Coker DF (2008) Iterative linearized approach to nonadiabatic dynamics. J Chem Phys 129:114,106 Eich FG, Agostini F (2016) The adiabatic limit of the exact factorization of the electron-nuclear wave function. J Chem Phys 145:054,110 Elliott P, Maitra NT (2012) Propagation of initially excited states in time-dependent densityfunctional theory. Phys Rev A 85:052,510 Elliott P, Furche F, Burke K (2009) 3 excited states from time-dependent density functional theory. Reviews in computational chemistry 26:91 Elliott P, Goldson S, Canahui C, Maitra NT (2011) Perspectives on double-excitations in TDDFT. Chem Phys 391(1):110 – 119 Epstein S (1954) Note on perturbation theory. Am J Phys 22:613 Fang JY, Hammes-Schiffer S (1999) Improvement of the internal consistency in trajectory surface hopping. J Phys Chem A 103:9399–9407 Flick J, Appel H, Ruggenthaler M, Rubio A (2017a) Cavity Born-Oppenheimer approximation for correlated electron-nuclear-photon systems. J Chem Theory Comput 13:1616–1625 Flick J, Ruggenthaler M, Appel H, Rubio A (2017b) Atoms and molecules in cavities, from weak to strong coupling in quantum-electrodynamics (QED) chemistry. Proc Nat Ac Sci 114:3026– 3034 TDDFT and Quantum-Classical Dynamics 41 Frenkel J (1934) Wave mechanics, Clarendon, Oxford edn Furche F (2001) On the density matrix based approach to time-dependent density functional response theory. J Chem Phys 114:5982–5992 Gaigeot MP, Lopez-Tarifa P, Martin F, Alcami M, Vuilleumier R, Tavernelli I, Herv´ edu Penhoat MA, Politis MF (2010) Theoretical investigation of the ultrafast dissociation of ionised biomolecules immersed in water: Direct and indirect effects. Mutat Res-Rev Mutat 704(1-3):45–53, URL http://www.sciencedirect.com/science/article/ pii/S1383574210000086 Gao X, Thiel W (2017) Non-hermitian surface hopping. Phys Rev E 95:013,308 Garashchuk S, Rassolov VA (2003) Quantum dynamics with Bohmian trajectories: energy conserving approximation to the quantum potential. Chem Phys Lett 376:358 Gawelda W, Johnson M, de Groot FMF, Abela R, Bressler C, Chergui M (2006) Electronic and molecular structure of photoexcited [Ru(II)(bpy)3]2+probed by picosecond x-ray absorption spectroscopy. J Am Chem Soc 128:5001–5009 G´ omez I, Reguero M, Boggio-Pasqua M, Robb MA (2005) Intramolecular charge transfer in 4aminobenzonitriles does not necessarily need the twist. Journal of the American Chemical Society 127(19):7119–7129 Grabo T, Petersilka M, Gross EKU (2000) Molecular excitation energies from time-dependent density functional theory. J Mol Struc (Theochem) 501-502:353–367 Grabowski ZR, Rotkiewicz K, Rettig W (2003) Structural changes accompanying intramolecular electron transfer: focus on twisted intramolecular charge-transfer states and structures. Chemical reviews 103(10):3899–4032 Granucci G, Persico M (2007) Critical appraisal of the fewest switches algorithm for surface hopping. J Chem Phys 126:134,114 Gritsenko O, Baerends E (2004) Asymptotic correction of the exchange–correlation kernel of timedependent density functional theory for long-range charge-transfer excitations. J Chem Phys 121:655 Gross E, Kohn W (1990) Time-dependent density-functional theory. Advances in quantum chemistry 21:255–291 Gross EKU, Kohn W (1985) Local density-functional theory of frequency-dependent linear response. Phys Rev Lett 55:2850–2852 Gross EKU, Ullrich CA, Gossmann UJ (1994) Density functional theory of time-dependent systems. In: Gross EKU, Dreizler RM (eds) Density Functional Theory, Plenum, New York, pp 149–171 Gross EKU, Dobson J, Petersilka M (1996) Density functional theory of time-dependent phenomena. In: Nalewajski RF (ed) Density Functional Theory II, Topics in Current Chemistry, vol 181, Springer, Berlin, pp 81–172 Hack MD, Wensmann AM, Truhlar DG, Ben-Nun M, Mart´ ınez TJ (2001) Comparison of full multiple spawning, trajectory surface hopping, and converged quantum mechanics for electronically nonadiabatic dynamics. J Chem Phys 115:1172 Helgaker T, Jørgensen P (1989) Configuration-interaction energy derivatives in a fully variational formulation. Theor Chem Acc 75:111–127 Heller EJ (1981) Frozen gaussians: A very simple semiclassical approximation. The Journal of Chemical Physics 75(6):2923–2931 Hellgren M, Gross EKU (2012) Discontinuities of the exchange-correlation kernel and chargetransfer excitations in time-dependent density-functional theory. Phys Rev A 85:022,514 Hirata S, Head-Gordon M (1999) Time-dependent density functional theory within the TammDancoff approximation. Chem Phys Lett 314:291–299 Hohenberg P, Kohn W (1964) Inhomogeneous electron gas. Phys Rev B 136:B864 Hsu C, Hirata S, Head-Gordon M (2001) Excitation energies from time-dependent density functional theory for linear polyene oligomers: butadiene to decapentaene. J Phys Chem A 105(2):451–458