Simulating Vibronic Spectra without Born-Oppenheimer Surfaces
Kevin Lively, Guillermo Albareda, Shunsuke A. Sato, Aaron Kelly, Angel Rubio
SSimulating Vibronic Spectra withoutBorn-Oppenheimer Surfaces
Kevin Lively, † Guillermo Albareda, † , ‡ , ¶ Shunsuke A. Sato, † , § Aaron Kelly, ∗ , † , (cid:107) andAngel Rubio ∗ , † , ¶ , ⊥ † Max Planck Institute for the Structure and Dynamics of Matter and Center forFree-Electron Laser Science, Luruper Chaussee 149, 22761 Hamburg, Germany ‡ Institute of Theoretical and Computational Chemistry, University of Barcelona, Martí iFranquès 1-11, 08028 Barcelona, Spain ¶ Nano-Bio Spectroscopy Group and ETSF, Universidad del País Vasco, 20018 SanSebastían, Spain § Center for Computational Sciences, University of Tsukuba, Tsukuba 305-8577, Japan (cid:107)
Department of Chemistry, Dalhousie University, Halifax B3H 4R2, Canada ⊥ Center for Computational Quantum Physics (CCQ), Flatiron Institute, 162 Fifth avenue,New York NY 10010, USA
E-mail: [email protected]; [email protected] a r X i v : . [ phy s i c s . c h e m - ph ] J a n bstract We show how vibronic spectra in molecular systems can be simulated in an efficientand accurate way using first principles approaches without relying on the explicit useof multiple Born-Oppenheimer potential energy surfaces. We demonstrate and analysethe performance of mean field and beyond mean field dynamics techniques for the H molecule in one-dimension, in the later case capturing the vibronic structure quite ac-curately, including quantum Franck-Condon effects. In a practical application of thismethodology we simulate the absorption spectrum of benzene in full dimensionality us-ing time-dependent density functional theory at the multi-trajectory mean-field level,finding good qualitative agreement with experiment. These results show promise for fu-ture applications of this methodology in capturing phenomena associated with vibroniccoupling in more complex molecular, and potentially condensed phase systems. Graphical TOC Entry a priori knowledge of the BO states that are involved,along with the associated potential energy surfaces and nonadiabatic couplings.An alternative strategy to the sum over states in the BO basis is to take a coordinatespace perspective, and construct the response function for the system of interest from directtime-propagation of the system in that picture.
This invariably requires some level of ap-proximation in the representation dynamics of the electronic and nuclear degrees of freedom,with different consequences for their coupling depending on the method chosen. The mixedquantum-classical Ehrenfest approach is a practical approximation to the fully quantum me-chanical dynamics of the system, and despite it’s approximate dynamics , a formally exactrepresentation of the quantum equilibrium structure of the correlated electronic and vibra-tional degrees of freedom can be included in a multi-trajectory Ehrenfest (MTEF) simulationthrough the use of the Wigner representation.
In this case, the Wigner transform mapsthe vibrational quantum states onto phase space distributions of continuous position andmomentum coordinates which can be sampled by an appropriate Monte Carlo procedure tocapture the quantum equilibrium structure of the problem. The limitations of the Ehrenfest3pproach, and other independent trajectory semiclassical methods, are well known andwhile there have been many attempts to ameliorate these shortcomings, with some excep-tions, most rely on the BO framework in their implementation.
In this work wetake a different approach to go beyond mean field theory based on the recently introducedinteracting conditional wave function (ICWF) formalism, which is able to capture corre-lated electronic and nuclear dynamics.
We apply MTEF and ICWF dynamics to anexactly solvable one dimensional H model, and show that these methods are able to recoverelectron-nuclear correlations in linear vibronic spectra without the need to calculate multipleBO surfaces . In addition, we show that the MTEF method can be easily extended to abinitio non-adiabatic molecular dynamics simulations by calculating the vibronic spectra forbenzene, where we find good agreement with experimental results.The linear spectrum of a system is given by the Fourier transform of time correlation func-tion (TCF) C AB ( t ) = (cid:104) [ ˆ A ( t ) , ˆ B ] (cid:105) of the transition dipole operator, ˆ µ , C µµ ( t ) = (cid:104) ˆ µ ( t )ˆ µ (0) (cid:105) , (unless otherwise stated all expressions are in atomic units): I ( ω ) = 4 πω c (cid:90) ∞−∞ dt e iωt (cid:104) [ˆ µ ( t ) , ˆ µ ] (cid:105) = 8 πω c (cid:60) (cid:90) ∞ dt e iωt T r (cid:16) ˆ µ ( t ) ˆ ρ eq ˆ µ ( t = 0) (cid:17) , (1)where the trace occurs over nuclear and electronic degrees of freedom, ˆ ρ eq is the equilibriumdensity matrix for the coupled system, and we evolve ˆ µ ( t ) in the Hilbert representation.Traditionally vibronic spectra are explained by invoking the Frank-Condon approximationin the BO picture, where the electronic system is instantly excited, thus promoting theunperturbed ground state nuclear system to a different electronic surface. If one has accessto the electronic states involved in a particular spectral range then the contributions tothe spectrum due to each electronic transition can be identified by resolving the transitiondipole operator in the basis of the electronic states of interest, and the vibronic side peaks ofthat transition can be calculated by propagating the initial state’s nuclear subsystem under4he effect of the non-equilibrium electronic occupation. When it is feasible to resolve thenuclear wavefunction dynamics, this can be one of the most accurate methods of calculatingmolecular vibronic spectra. Although resolving eq. (1) in the BO framework is a powerful analysis tool, it is compu-tationally impractical for systems with many nuclear degrees of freedom, particularly whenone desires spectra over multiple surfaces. One can bypass this computational bottleneck byrepresenting the system in a real space basis and using the “ δ -kick" method, which captureselectronic transitions to all dipole-transition allowed states (resolved on the grid) within asingle calculation by utilising the dipole response to a perturbative, but impulsive externalfield H field = E ( t )ˆ µ , i.e. with E ( t ) = κδ ( t ) and κ << . Using first order perturbationtheory, the dipole response (cid:104) ∆ µ ( t ) (cid:105) = (cid:104) µ ( t ) (cid:105) − (cid:104) µ (0) (cid:105) can be written in powers of the field, (cid:104) ∆ µ ( t ) (cid:105) = i T r (cid:16) (cid:2) ˆ µ I ( t ) , ˆ µ I (0) (cid:3) ˆ ρ eq (cid:17) κ + O ( κ ) , (2)where ˆ µ I ( t ) is evolved in the interaction representation. Hence, the linear response spectramay also be obtained via the following relation, C µµ ( t ) = − iκ (cid:104) ∆ µ ( t ) (cid:105) , (3)provided the strength of the perturbing field κ is sufficiently small. This δ -kick approachonly requires the initial state of the full system as input, followed by time propagation for asufficient interval so as to obtain the desired energy resolution. Importantly, this techniquecan also serve as a foundation for calculating non-linear optical response spectra. While the methods described above are formally equivalent, difference between the calcu-lated spectra can arise when approximations are made. Here we briefly describe two methodsfor performing coupled electron nuclear dynamics simulations, the quantum-classical mean-field MTEF method, and the ICWF formalism which was designed to go beyond the meanfield limit. 5 typical approach to Ehrenfest theory is to assume a separable electronic-nuclear wavefunction ansatz, take the classical limit of the nuclear portion, and initialise the nuclei atthe equilibrium position with zero nuclear momentum.
This single trajectory Ehrenfest(STEF) method is often employed when a mixed quantum-classical method is needed tocouple electronic and nuclear dynamics, in some cases providing a stark difference in elec-tronic dynamics compared to fixed nuclei. Although attempts at capturing quantizedvibrational effects in STEF with the δ -kick method have been made, they can containunphysical spectral features (see SI).An alternative route to Ehrenfest is also possible in the density matrix picture, andproceeds via the quantum-classical Liouville equation. The major difference is that thisrepresentation results in a multi-trajectory Ehrenfest picture of the dynamics, where the ini-tial quantum statistics of the correlated system can, in principle, be captured exactly. Here,we outline the evolution equations, and offer more details in the supporting information.The time evolution of the reduced electronic density is ddt ˆ ρ e ( t ) = − i (cid:104) ˆ H Effe,W ( X ( t )) , ˆ ρ e ( t ) (cid:105) , (4)where the subscript W refers to the partial Wigner transform over the nuclei, X = ( R , P ) isa collective variable for the nuclear position R and momentum P , and the effective electronicmean-field Hamiltonian is ˆ H Effe,W ( X ( t )) = ˆ H e + ˆ H en,W ( X ( t )) , where ˆ H e refers to the electronicportion of the hamiltonian, and ˆ H en to the electron nuclear coupling. The nuclear dynamics isrepresented as an ensemble of N independent Wigner phase-space trajectories, ρ n,W ( X , t ) = N (cid:80) Ni δ ( X i − X i ( t )) , that evolve according to Hamilton’s equations of motion generatedfrom the effective nuclear mean-field Hamiltonian, ∂ R i ∂t = ∂H Effn,W ∂ P i , ∂ P i ∂t = − ∂H Effn,W ∂ R i H Effn,W = H n,W ( R i ( t )) + T r e (cid:16) ˆ ρ e ( t ) ˆ H en,W ( X i ( t )) (cid:17) . (5)6he average value of any observable, (cid:104) O ( t ) (cid:105) , can then be written as follows, (cid:104) O ( t ) (cid:105) = T r e (cid:90) d X ˆ O W ( X , t ) ˆ ρ W ( X , (6)which can be evaluated by sampling initial conditions from ˆ ρ W ( X , , and evolving theexpectation value of the observable with according to the above equations of motion. Usingthis dynamics method in conjunction with the BO basis representation to evaluate eq.s (1)and (4 - 6) ultimately leads to the following equations of motion, with sums over BO statesdenoted by a , (see SI for details) ∂ t ρ aa (cid:48) e = − iρ aa (cid:48) e ( t )( (cid:15) a ( R i ( t )) − (cid:15) a (cid:48) ( R i ( t )))+ (cid:88) a (cid:48)(cid:48) P i ( t ) M (cid:16) ρ aa (cid:48)(cid:48) e ( t ) d ia (cid:48)(cid:48) a (cid:48) ( t ) − d iaa (cid:48)(cid:48) ( t ) ρ a (cid:48)(cid:48) a (cid:48) e ( t ) (cid:17) ∂ t R i ( t ) = P i ( t ) /M∂ t P i ( t ) = (cid:88) a − ∂ R (cid:15) a ( R i ( t )) ρ aae ( t )+ (cid:88) aa (cid:48) (cid:60) (cid:104)(cid:0) (cid:15) a ( R i ( t )) d iaa (cid:48) ( t ) − (cid:15) a (cid:48) ( R i ( t )) d ia (cid:48) a ( t ) (cid:1) ρ a (cid:48) ae ( t ) (cid:105) ∂ t µ aa (cid:48) W ( R i ( t )) = iµ aa (cid:48) W ( R i ( t ))( (cid:15) a ( R i ( t )) − (cid:15) a (cid:48) ( R i ( t ))) . (7)Where (cid:15) a ( R ) are the BO surfaces and d aa (cid:48) are the non-adiabatic coupling vectors (NACVs)between states a and a (cid:48) .In contrast to the previous expression, utilising MTEF in the real space δ -kick approachrequires initialising the electronic wave function as the BO eigenstate for each initially sam-pled nuclear geometry. The δ − kick is applied and the electronic wave function is propagatedusing the time dependent Schrödinger equation equivalent to eq. (4) alongside the nucleiaccording to eq. (5). Calculating the spectrum via MTEF dynamics in the BO picture isfrom here on referred to as MTEF-BO, and calculating it via the δ − kick method is referredto as MTEF-kick. 7oving beyond semi-classical dynamics, the formally exact CWF method and it’s prac-tical ICWF implementation are recently developed methods which have shown to be able tocapture non-equilibrium correlated nuclear-nuclear and electron-nuclear phenomena beyondthe mean field limit. This approach is based on taking single-particle slices (the CWFs)of the time-dependent wave function of full system, and approximating the equations of mo-tion for these CWFs by the Hermitian components of the sliced Hamiltonian, and finally, inthe ICWF extension, utilising these electron-nuclear CWFs as a basis of Hartree productsin a wave function ansatz.Here we describe an implementation of this approach utilising the static and time-dependent variational principles for the expansion coefficients in a static
CWF basis. Thebasis is chosen via sampling electronic and nuclear positions ( r α , R α ) , α ∈ { , . . . , N c } ,where r and R are understood to be collective position variables, from initial guesses tothe electronic and nuclear densities. These are used to construct the Hermitian limit of theCWF propagators h αe ( r i ) = − ∇ r i + N e (cid:88) j (cid:54) = i V ee ( r i , r αj ) + N n (cid:88) l V en ( r i , R αl ) h αn ( R l ) = − M l ∇ R l + N e (cid:88) j V en ( R l , r αj ) + N n (cid:88) m (cid:54) = l V nn ( R l , R αm ) (8)for a system with N e electrons and N n nuclear degrees of freedom. Taking eigenstates of h αe ( r i ) and h αn ( R l ) , denoted φ α ( r i ) and χ α ( R l ) respectively, as our CWF basis we write thefollowing wave function ansatz: Ψ( r , R , t ) = N c (cid:88) α C α ( t ) N e (cid:89) i φ α ( r i ) N n (cid:89) l χ α ( R l ) , (9)where we have taken a Hartree product of electronic and nuclear CWFs for each degree offreedom. While the Hartree product over electronic degrees of freedom has been sufficientfor accuracy in applications of ICWF so far, this ansatz can in principle be trivially extended8o have fermionic anti-symmetry via inclusion of Slater determinants. We then utilise theDirac-Frenkel variational procedure to develop equations of motion for (cid:126)C ( t ) , which leadsto the following evolution equation for the expansion coefficients, ddt (cid:126)C = − i S − H (cid:126)C, (10)where S αβ = N e (cid:89) i (cid:90) d r i ( φ α ( r i )) ∗ φ β ( r i ) N n (cid:89) l (cid:90) d R l ( χ α ( R l )) ∗ χ β ( R l ) , H αβ = N e (cid:89) i N n (cid:89) l (cid:90) d R l d r i ( φ α ( r i ) χ α ( R l )) ∗ ˆ H ( r , R ) φ β ( r i ) χ β ( R l ) . In practice S may be nearly singular, but its inverse can be approximated by the Moore-Penrose pseudoinverse. The ground state wave function is obtained from this approachusing imaginary time evolution, and the δ − kick spectra (ICWF-kick) is calculated byapplying the perturbative field to the CWFs at time zero and recalculating the S and H matrices, equivalent to propagating in the interaction representation. This "closed-loop" ofinitial state preparation and time-propagation ensures that our ICWF approach is a fullyself-consistent method that increases in accuracy with increasing N c , and requires no BOstate information.To investigate the performance of the MTEF and ICWF approaches to vibronic spectrallineshapes we studied the vibronic transitions in an exactly solvable one dimensional modelsystem for molecular Hydrogen. The total Hamiltonian can be written in the center ofmass frame in atomic units as ˆ H ( r , r , R ) = − ∂ R µ n − (cid:88) i =1 ∂ r i µ e + 1 (cid:112) ( r − r ) + 1 + 1 R − (cid:88) i =1 (cid:113) ( r i + R ) + 1 + 1 (cid:113) ( r i − R ) + 1 (11)9here µ n = m p / and µ e = 2 m p / (2 m p + 1) are the reduced nuclear and electronic masses, R is the internuclear separation, and r i are the electronic positions. We take the protonmass to be m p = 1836 . The electronic and nuclear degrees of freedom were each resolvedon grids for the numerically exact solution and ICWF-kick approaches, while the MTEF-kick electronic wave functions were time evolved on the ( r , r ) grid, and the MTEF-BOinformation was calculated by solving the electronic subsystem across the nuclear grid; seethe computational methods section for more details. A kick strength of κ = 10 − a.u. wassufficient to generate the kick spectra within the linear response regime and, unless otherwisestated, a total propagation time of , a.u. ≈ f s was used to generate the spectra. I ()[ a ] Exact010 I ()[ a ] MTEF-BO8.5 9.0 9.5 10.0 10.5 11.0 11.5Energy [eV]01020 I ()[ a ] (a)(b)(c) STEF-kickMTEF-kick
Figure 1: 1D H , S ← S spectra calculated via the MTEF-TCF, MTEF-kick and STEF-kick approaches, with the exact peak placements overlaid as dashed vertical lines. Spectralcross sections are reported in square Bohr radii a . For clarity the STEF-kick spectrum hasbeen multiplied by a factor of 0.175 to match the scale of the MTEF-kick results.In Fig. 1 we show mean field spectra calculated both with (MTEF-BO), and without(MTEF-kick) the use of multiple BO surfaces for the absorption from S to S in comparison10ith the numerically exact results. We see that in the BO picture the MTEF methodrecovers the vibronic absorbtion peak placement quite accurately for the first five peaks,with a broadening occurring for the higher energy peaks that leads to a loss of structure.This broadening of the spectral signal is due to the well-known fact that the MTEF dynamicsdoes not preserve the correct quantum statistics and thus cannot fully capture the electron-nuclear correlation in the problem (see the SI for a detailed discussion of this issue). Thepre-peak features in Fig. 1b are also unphysical artefacts of MTEF. The MTEF-BO spectrawere converged to within graphical accuracy using N = 50 , trajectories.Focusing on the MTEF-kick results in Fig. 1c, we see that this approach recovers vi-bronic side peak structures again without any BO surface information, albeit with inaccuratespacing, while STEF-kick captures only the vertical electronic transition from the minimumof the S surface. The MTEF-kick spectra converged to within graphical accuracy using N = 30 , trajectories. The average peak spacing in the MTEF-kick spectra is approxi-mately . eV; this corresponds remarkably well with the natural frequency of the harmonicapproximation to the ground state surface expanded around the equilibrium geometry, whichis also . eV in this case. This result is unsurprising as the electronic kick induces a verysmall population transfer to the upper surface proportional to the square of the kick strength,which results in the mean forces on the nuclei in MTEF-kick essentially corresponding tothose of the initial state. The influence of initial state on the MTEF-kick spectra is furtherdemonstrated by analysing the emission spectra in Fig. 2. The initial state here was cho-sen by hand as the lowest lying nuclear state on the S surface. Once again we see thatMTEF-BO recovers the peak placement quite well, while the MTEF-kick data has a less ac-curate vibronic spacing. Fitting the MTEF-kick peaks, we find an excellent correspondencebetween mean spacing of the five lowest energy MTEF peaks and the excited surface naturalfrequency of . eV.For ICWF-kick, we found that N c = 4096 and mixing the three lowest energy CWFeigenstates was sufficient to obtain quite accurate results. In Fig. 3 we demonstrate that the11 I ()[ a ] Exact020 I ()[ a ] MTEF-BO8.2 8.4 8.6 8.8 9.0 9.2 9.4Energy [eV]02040 I ()[ a ] (a)(b)(c) STEF-kickMTEF-kick
Figure 2: S ← S spectra compared between the MTEF-TCF, MTEF-kick and STEF-kickapproaches, with exact peak placement overlaid as dashed vertical lines. MTEF nuclearinitial conditions are sampled from the lowest lying vibrational state on S . The sign ofall spectra here is inverted for ease of comparison to other figures, and for legibility theSTEF-kick spectrum was multiplied by a factor of . to match the MTEF-kick spectramaximum. 12 I ()[ a ] Exact01020 I ()[ a ] ICWF-kick8.5 9.0 9.5 10.0 10.5 11.0 11.5Energy [eV]01020 I ()[ a ] (a)(b)(c) MTEF-kick
Figure 3: S ← S spectra of the ICWF-kick and MTEF-kick methods, with the exact peaksplacement overlaid as dashed lines. 13CWF ansatz used in a variational context achieves a much more accurate vibronic spacingthan the MTEF-kick approach, without the failing of peak broadening or unphysical spectralnegativity apparent in the MTEF-BO results. The accuracy of these results underscoresthat the ICWF ansatz is a robust framework to capture the electronic and vibronic quantumdynamics, being accurate for not only the electron-nuclear correlation inherent to vibronicspectra, but also the electronic subsystem itself, which in the MTEF results was solvedexactly either on a grid or using explicit BO state information. The deviation from theexact results does grow with increasing energy, although this is ameliorated with increasing N c , and can in principle be eliminated at large enough values of N c (see SI).Finally we demonstrate the application of MTEF-kick to real 3D molecular systems usingthe ab initio Octopus real-space time dependent density functional theory (TDDFT) package to calculate the linear vibronic MTEF-kick spectra of Benzene. The initial conditionsfor the nuclear subsystem were obtained by calculating the normal mode frequencies anddynamical matrix of the molecule, and sampling Wigner transforms of the ground statewave functions in the harmonic approximation; see SI for more details. The adiabatic-LDAfunctional was used, along with norm-conserving Troullier-Martins pseudo-potentials, andthe trajectories were evolved for (cid:126) eV ≈ f s . A kick strength of κ = 5 x − a.u. was usedto generate the kick spectra within the linear response regime in this case, and the graphicalconvergence of the MTEF results was found to be achieved with N = 500 trajectories.In Fig. 4, we compare the MTEF-TDDFT-kick results to its STEF-TDDFT-kick coun-terpart, each scaled to match the peak intensity of an experimental data set for the opticalabsorption of benzene digitized from Ref. We see that there is remarkably good agreementacross the wide energy range available from experiment, before molecular dissociation path-ways become available around 13.8eV. Again, the full spectrum is resolved without resortingto calculations of individual transitions as would be required in a BO state calculation. Thethree STEF peaks in the 7eV region correspond to the energy range of the doubly degen-erate, dipole allowed E u ← A g , π ∗ ← π transition, with the energy degeneracy14 E x p e r i m e n t a l A r b . U n i t s MTEF-TDDFT-kickSTEF-TDDFT-kickExperiment
Figure 4: Experimental vibronic spectra for the lowest lying optical transitions of ben-zene compared to the MTEF and STEF kick spectra calculated with TDDFT. The spectralweights of the MTEF and STEF data sets have been scaled to the arbitrary units of theexperimental data set. 15rtificially lifted by the discrete grid. The experimental band preceding the central peak,in the range 6eV to 6.5eV is commonly ascribed to the dipole forbidden, but vibronicallyallowed B u ← A g transition, and in the MTEF-TDDFT-kick results we see a lowenergy tail extending through the 5eV-6.25eV range, well away from the STEF results, even-tually transitioning into the broad peak centered around 7eV. It’s reasonable to expect thatthe broadening of the MTEF signal relative to the experimental signal is due to the effectsdiscussed above that arise due to the mean field treatment.We have demonstrated that semi-classical MTEF simulations can capture vibronic struc-ture with the correct spectral sign in the region of the transition. Moreover, we have shownhow this can be achieved without using multiple BO surfaces via the δ -kick method, andthat the vibronic spacing predicted with the MTEF-kick approach matches the profile ofthe initial state. We have addressed these shortcomings by combining the ICWF formalismwith the δ -kick method, which provides more accurate vibronic spectra in a computationallyefficient and systematically improvable fashion. Finally, we demonstrated that MTEF-kickis easily applied to ab initio molecular systems by simulating the vibronic spectra of benzeneand finding good agreement to experimental results.These linear response results establish a solid basis for further investigations into non-linear response of field driven molecular systems utilising the practical and efficient MTEFand ICWF techniques along with ab initio electronic structure methods. Work in preparationby the present authors also explores the utility of ICWF with electron-electron and electron-nuclear correlated systems, and explores the response of these systems under nonperturbativeelectric fields. Furthermore we expect that MTEF-kick will improve in accuracy for periodicsystems, as changes in the electronic configuration are often to likely produce smaller changesin the nuclear forces than in molecular hydrogen. This makes this method interesting topursue in periodic systems in particular, where there is a dearth of theoretical frameworks for ab initio , nonpertubrative electron-nuclear coupling. Work in this direction is in progress,as is the implementation of the ICWF method within an ab initio framework for molecular16nd periodic systems.
Computational Methods
In the 1D H model, the electronic coordinates are each resolved on a a wide intervalwith spacing . a , while the nuclear grid extends to R max = 6 . a with . a spacing.Quadratic complex absorbing potentials were also added to the Hamiltonian to preventreflection from the simulation box edge (see SI). To generate the exact results we evolvedthe full wave function under the δ − Kick on the three dimensional electron-nuclear grid, whilefor MTEF-kick, the electronic subsystem’s Schrödinger equation, dependent on R i ( t ) , wassolved exactly on the two dimensional electronic grid for each trajectory. All wave functionswere time-propagated using a fourth-order Runge-Kutta integration scheme with a time-step size of ∆ t = 0 . a.u.. For the MTEF trajectories, the nuclear degree of freedom waspropagated via a veloctiy-Verlet type scheme with the same time-step size. An exponentialdamping mask function exp ( − γt ) was applied to all time dependent signals in the Fouriertransform, with the damping factor was set to damp the signal to 0.1% it’s strength at thefinal time.For the 1D H MTEF-BO results, the potential energy surfaces (cid:15) a ( R ) and µ aa (cid:48) W ( R ) werecalculated on a nuclear grid with ∆ R = 0 . a up to R max = 8 a , fit to a cubic splinefunction, and interpolated every . R . The NACV between S and S in this modelis numerically zero. These quantities were resolved for the first allowed dipole transition,between the ground state S and the second excited state S , and the results were found tobe well converged within about × trajectories.For the MTEF-TDDFT-kick simulations we used a real space grid formed from overlap-ping spheres of radius Å centered on the initial positions of the nuclei, with an isotropicgrid spacing of . Å, which was found to be sufficient to converge the energies of the lowestlying absorption lines. 17 cknowledgement
This work was supported by the European Research Council (ERC-2015-AdG694097), theCluster of Excellence Advanced Imaging of Matter’ (AIM), JSPS KAKENHI Grant Num-ber 20K14382, Grupos Consolidados (IT1249-19) and SFB925. The Flatiron Institute is adivision of the Simons Foundation.
Supporting Information Available
MTEF Equations of Motion
Starting from a density matrix representation of the full system, ˆ ρ , we Wigner trans-form over the nuclear subsystem, producing a unique mapping onto a nuclear position R and momentum P phase space X = ( R , P ) , where R and P are collective variables R = ( R , . . . , R N n ) , P = ( P , . . . , P N n ) , with R i , P i ∈ R d . The partial wigner transform isdefined for any operator as ˆ ρ W ( R , P ) = 1(2 π ) dN n (cid:90) d X e i P · X (cid:104) R − X | ˆ ρ | R + X (cid:105) , (12)leaving a Hilbert space operator character over the electronic degrees of freedom, dependenton the continuous nuclear phase space parameters. In general, developing equations of motionfor ˆ ρ W ( R , P ) , (or any operator), requires taking the partial Wigner transformation of theLiouville von-Neumann equation of motion for ρ : ∂ ˆ ρ W ∂t = − i (cid:16) ( ˆ H ˆ ρ ) W − ( ˆ ρ ˆ H ) W (cid:17) ( ˆ H ˆ ρ ) W = ˆ H W exp (cid:16) i Λ (cid:17) ˆ ρ W Λ = ←−∇ P · −→∇ R − ←−∇ R · −→∇ P g exp (cid:16) κ Λ (cid:17) f = ∞ (cid:88) s =0 κ s s ! s (cid:88) t =0 ( − t (cid:18) st (cid:19) (cid:2) ∂ s − t R ∂ t P f (cid:3) (cid:2) ∂ t R ∂ s − t P g (cid:3) . (13)18here the final line defines the “Moyal product” also known as the “star product”. Byexpressing the Poisson braket operator Λ , in terms of the ratio of masses between the nucleiand the electrons Λ = ( m/M ) Λ (cid:48) , and truncating the Moyal product of e ( m/M ) Λ (cid:48) at firstorder, one can arrive at the Quantum-Classical Liouville Equation (QCLE): i ∂∂t ˆ ρ W ( R , P ) = − i [ ˆ H W , ˆ ρ W ] + 12 (cid:16) ˆ H W , ˆ ρ W − ˆ ρ W , ˆ H W (cid:17) , (14)where { A ( R , P ) , B ( R , P ) } refers to the normal Poisson bracket.To derive MTEF equations of motion from the QCLE, one takes the mean field ap-proximation by assuming that the full system can be written as a sum of correlated anduncorrelated parts, ˆ ρ W ( X , t ) = ˆ ρ e ( t ) ρ n,W ( X , t ) + ˆ ρ corr,W ( X , t ) , (15)and then neglecting the contribution of the correlated part in the dynamics . Note that whilethe ensuing dynamics do not explicitly treat the effect of subsystem correlation, the initialstate generally is correlated, and therefore is implicitly included in the dynamics.Under this approximation, the electronic density matrix is ˆ ρ e ( t ) = T r n (cid:16) ˆ ρ ( t ) (cid:17) = (cid:90) d X ˆ ρ W ( X , t ) , (16)and the nuclear (quasi) probability phase space distribution is ρ n ( X , t ) = T r e ( ˆ ρ W ( X , t )) .In the equations of motion resulting from inserting this approximation into the QCLE, theevolution of the reduced Wigner density of the nuclear subsystem can be exactly represented,via the method of characteristics, by a sufficiently large ensemble of multiple independenttrajectories, ρ n,W ( X , t ) = N (cid:80) Ni δ ( X i − X ( t )) . Each trajectory evolves according to Hamil-19on’s equations of motion generated from the mean-field effective Hamiltonian, ∂ R i ∂t = ∂H Effn,W ∂ P i , ∂ P i ∂t = − ∂H Effn,W ∂ R i H Effn,W = H n,W ( X i ( t )) + T r e (cid:16) ˆ H en,W ( X i ( t )) ˆ ρ ie ( t ) (cid:17) . (17)Where H n,W and H en,W refer to the partially Wigner transformed nuclear and electron-nuclear coupling operators, respectively. The electronic density associated with each trajec-tory , ρ ie ( t ) , evolves according to the following commutator: ddt ˆ ρ ie ( t ) = − i (cid:104) ˆ H e + ˆ H en,W ( X i ( t )) , ˆ ρ ie ( t ) (cid:105) . (18)The exact expression for the average value of any observable, (cid:104) O ( t ) (cid:105) , can be written as (cid:104) O ( t ) (cid:105) = T r e (cid:90) d X ˆ O W ( X , t ) ˆ ρ W ( X ,
0) =
T r e (cid:90) d X ˆ O W ( X ) ˆ ρ W ( X , t )= N (cid:88) i T r e (cid:16) ˆ O W ( X i ( t )) ˆ ρ ie ( t ) (cid:17) (19)The mean field limit of this expression simple corresponds to evaluating the integral bysampling initial conditions for an ensemble of independent trajectories from ˆ ρ W ( X , , andthen generating the time evolution for each trajectory by approximating ˆ O W ( X , t ) by it’smean-field counterpart.Following the sampling of an initial nuclear condition, X i , from the Wigner distributionassociated to the nuclear subsystem wave function, the electronic system is initialised as: ( ˆ H e + ˆ H en,W ( R i )) φ a ( r ) = (cid:15) a ( R i ) φ a ( r ) , (20)i.e. implicitly as the BO electronic state at R i . Under this scheme, the electronic subsystem’sinitial conditions are implicitly correlated with the nuclear subsystem’s quantum statistics.In cases where the nuclear initial state is impractical to calculate exactly one may utilise20he normal modes of the molecular system, or phonon coordinates of a periodic system, totreat the full nuclear wavefunction as a Hartree product of N uncoupled harmonic oscillators,where N is the number of non-rotational and non-translational nuclear degrees of freedom: χ n ( R ) ≈ χ ( Q ) ⊗ . . . ⊗ χ N ( Q N ) χ i ( Q i ) = (cid:88) l c ( i ) l χ li ( Q i ) . (21)With c ( i ) l referring to the occupation of the l th excited state of normal mode i with wave-function, χ li , and Q i ( R ) the normal mode coordinate. Formally, this is exactly equivalent totaking a second order Taylor expansion approximation of the BO surface about the equilib-rium nuclear position R : H nuc ( R , P ) = (cid:88) l M l P l + (cid:88) lm
12 ( R l − R l ) ∂ V BO ∂ R l ∂ R m (cid:12)(cid:12)(cid:12)(cid:12) R ( R m − R m ) . (22)Defining the dynamical matrix, H lm = √ M l ∂ V∂ R l ∂ R m √ M m , and it’s diagonalizing unitary trans-form, D T H D = Ω , D T D = , where Ω ij = ω i δ ij , we construct the normal coordinatetransform for all non-rotational, non-translational (imaginary) ω i , (here we include (cid:126) forclarity): (cid:112) M l ( R l − R l ) = (cid:88) i D li q i , P l √ M l = (cid:88) i D li s i s i = (cid:112) (cid:126) ω i S i , q i = (cid:114) (cid:126) ω i Q i , (23)such that we obtain the nuclear Hamiltonian in dimensionless normal mode coodinates: H ( Q, S ) = (cid:88) i (cid:126) ω i S i + Q i ) . (24)Of course, the simple harmonic wave function solutions to the above Hamiltonian havewell known analytical expressions and are trivially Wigner transformed, the ground state21armonic oscillator wavefunction’s Wigner function for instance is: W ( Q, S ) = 1 π exp (cid:0) − S − Q (cid:1) . (25)We can therefore sample these transforms for ( Q, S ) and then use eq. (23) to back transformto from normal mode coordinates to cartesian coordinates. MTEF-BO Equations of Motion in the Born Oppenheimer Basis
In deriving the MTEF equations of motion in the BO basis, we start by writing the molecularhamiltonian in terms of position and momentum space operators for the electrons (lightparticles), ˆ r, ˆ p and nuclei (heavy particles) ˆ R, ˆ P . These are again understood to be collectivevariables. ˆ H (ˆ r, ˆ p, ˆ R, ˆ P ) = 12 M ˆ P + ˆ h e (ˆ r, ˆ p, ˆ R )ˆ h (ˆ r, ˆ p, ˆ R ) = 12 ˆ p + ˆ V (ˆ r, ˆ R )ˆ V (ˆ r, ˆ R ) = ˆ V ee (ˆ r ) + ˆ V en (ˆ r, ˆ R ) + ˆ V nn ( ˆ R ) . (26)We then utilise a position representation in the nuclear dof by expanding in the space ofnuclear position states R = (cid:82) d R | R (cid:105) (cid:104) R | , leading to ˆ H ( R ) = − M ∇ R + ˆ h e (ˆ r, ˆ p, R ) (27)For a transition between two electronic states g and e , we can expand in the adiabatic basis | φ a ( R ) (cid:105) , ( a = g, e ) which are dependent on the nuclear positions R defined by, ˆ h e ( R ) | φ a ( R ) (cid:105) = (cid:15) a ( R ) | φ a ( R ) (cid:105) . (28)22aking the partial Wigner transform of eq. (27) leads to ˆ H W ( R , P ) = 12 M P + ˆ h e,W (ˆ r, ˆ p, R ) (29)where ˆ h e,W ( R ) is the normal electronic hamiltonian operator, now dependent on R in theWigner nuclear phase space. Starting with the separability approximation for the densityoperator, and neglecting correlations, we have ˆ ρ W = ˆ ρ e ρ n ( R , P ) , with ∂ t ˆ ρ e = − i (cid:104) T r X (cid:104) ˆ h e,W ( R ) (cid:105) , ˆ ρ e (cid:105) (30)where T r X (cid:104) . . . (cid:105) = (cid:82) . . . d R d P , and P scalar terms are cancelled by the commutator. Weare of course interested in evaluating the dipole-dipole correlation function: C µµ ( t ) = (cid:90) d R d P T r e (cid:110) ˆ µ W ˆ σ ( t ) (cid:111) = (cid:90) d R d P T r e (cid:110) ˆ µ W ( t )ˆ σ (0) (cid:111) , (31)where ˆ σ = [ˆ µ W , ˆ ρ W ] , and we resolve the dipole operator as ˆ µ W ( R , t = 0) = − ˆ r + Z R R = | φ a (cid:105) (cid:104) φ a | − ˆ r | φ a (cid:48) (cid:105) (cid:104) φ a (cid:48) | + δ aa (cid:48) Z R R = R µ ge ( R ) µ eg ( R ) R (32)Where Z R refers to the ionic charge of each nuclei. In practice we can neglect the intra-state R term as we are focused entirely on the transition dipole moment.Taking the initial state as the ground state, ( | Ψ (cid:105) = | χ g φ g (cid:105) ) ˆ ρ W ( R , P ,
0) = ρ ng ( R , P ) , (33)23eads to ˆ σ (0) = [ˆ µ W , ˆ ρ W ( R , P , ρ ng ( R , P ) − µ ge ( R ) µ eg ( R ) 0 . (34)And therefore the correlation function becomes C µµ ( t ) = (cid:90) d R d P ( µ geW ( R , t ) σ eg (0) + µ egW ( R , t ) σ ge (0))= (cid:90) d R d P (cid:0) µ geW ( R , t ) µ egW ( R , − µ egW ( R , t ) µ geW ( R , (cid:1) ρ ng ( R , P ) . (35)We can construct an identical quantity from a different initial condition as a superpositionstate ( | Ψ (cid:105) = √ | χ g (cid:105) ( | φ g (cid:105) + i | φ e (cid:105) ) ) giving, ˆ˜ ρ W ( R , P,
0) = ρ ng ( R , P ) 12 − ii (36)For this different initial condition we propagate ˜ C µµ ( t ) = i (cid:90) dRdP (cid:16) µ geW ( R , t ) µ egW ( R , − µ egW ( R , t ) µ geW ( R , (cid:17) ρ ng ( R , P )= i C µµ ( t ) (37)With this different initial condition, we take the MTEF form of the nuclear density arisingfrom the Monte Carlo integration described above, ρ n ( R , P ) = 1 N (cid:88) i δ ( R − R i ( t )) δ ( P − P i ( t )) . (38)The subsequent equations of motion for the system are for the electronic density, needed for24he nuclear trajectories are: ∂ t ˜ ρ aa (cid:48) e = − i ˜ ρ aa (cid:48) e ( t )( (cid:15) a ( R i ( t )) − (cid:15) a (cid:48) ( R i ( t )))+ (cid:88) a (cid:48)(cid:48) P i ( t ) M (cid:16) ˜ ρ aa (cid:48)(cid:48) e ( t ) d ia (cid:48)(cid:48) a (cid:48) ( t ) − d iaa (cid:48)(cid:48) ( t ) ˜ ρ a (cid:48)(cid:48) a (cid:48) e ( t ) (cid:17) ∂ t R i ( t ) = P i ( t ) /M∂ t P i ( t ) = 12 (cid:88) aa (cid:48) (cid:16) F aa (cid:48) W ( t ) ˜ ρ a (cid:48) ae ( t ) + ˜ ρ aa (cid:48) e ( t ) F a (cid:48) aW ( t ) (cid:17) = (cid:88) aa (cid:48) (cid:60) (cid:104) F aa (cid:48) W ( t ) ˜ ρ a (cid:48) ae ( t ) (cid:105) = (cid:88) a − ∂ R (cid:15) a ( R i ( t )) ˜ ρ aae ( t )+ (cid:88) aa (cid:48) (cid:60) (cid:104)(cid:0) (cid:15) a ( R i ( t )) d iaa (cid:48) ( t ) − (cid:15) a (cid:48) ( R i ( t )) d ia (cid:48) a ( t ) (cid:1) ˜ ρ a (cid:48) ae ( t ) (cid:105) (39)Where in the last two equations we have used the identity d iaa (cid:48) ( t ) = (cid:104) φ a | ∂ R i φ a (cid:48) (cid:105) | R i ( t ) = − ( d ia (cid:48) a ( t )) ∗ , to manipulate F aa (cid:48) W ( t ) = − (cid:104) φ a ( R ) | ∂ R ˆ H W | φ a (cid:48) ( R ) (cid:105) | R i ( t ) . Note that for transitionslike the S /S transition 1D H focused on in the body of this paper, the non-adiabaticcoupling vector (NACV) d aa (cid:48) = 0 , means that the mean field force acting on the nuclei is atall times a superposition of the S and S surfaces.These are propagated alongside the dipole matrix element equations of motion, neededfor the correlation function: ∂ t µ aa (cid:48) W ( R i ( t )) = iµ aa (cid:48) W ( R i ( t ))( (cid:15) a ( R i ( t )) − (cid:15) a (cid:48) ( R i ( t ))) . (40) STEF Spectral Negativity
As mentioned in the main text, previous work by Goings et. al employed STEF-kick dy-namics simulations to calculate spectra in fully ab-initio 3D H by initialzing the nucleargeometry in non-equilibrium ‘compressed’ geometries. Geometries were selected correspond-ing to expected vibrational energies from Boltzmann distributions at arbitrary temperatures25nd the δ − Kick method was used to excite the electronic subsystem. Furthermore, only the magnitude of the spectral response was depicted, which does not show the spectral negativ-ity resulting from initialising the mean field simulations in a non-equilibrium state. Here weutilise the canonical initial conditions of the STEF-BO picture for the 1D H model. Theelectronic occupation is equal for each of the two surfaces ivolved in the transition, and thenuclear initial condition corresponds to the equillibrium geometry of the initial surface. In I ()[ a ] Exact8.5 9.0 9.5 10.0 10.5 11.0010 I ()[ a ] MTEF-BO8.5 9.0 9.5 10.0 10.5 11.050050 I ()[ a ] (a)(b)(c) STEF-BO
Figure 5: 1D H S ← S absorption spectra, comparing exact, MTEF-BO and STEF-BO.Fig. 5c we see the results of STEF-BO for the S ← S region of the spectrum, showingthat this only captures positive spectral intensities in the vicinity of the exact results, withaccurate peak placement only at the MTEF level. Furthermore the contributions to theunphysical pre-peak features of individual trajectories become apparent in the low energytail. For completeness we also feature the S ← S results in Fig. 6, which demonstrate thesame features of correct spectral sign only in the region of the exact results and alternatingsign elsewhere. 26 I ()[ a ] Exact8.5 9.0 9.5 10.0 10.5 11.0020 I ()[ a ] MTEF-BO8.5 9.0 9.5 10.0 10.5 11.00100 I ()[ a ] (a)(b)(c) STEF-BO
Figure 6: 1D H S ← S absorption spectra, comparing exact, MTEF-BO and STEF-BO. Application to Displaced Harmonic Oscillator Model
In order to investigate the limitations of MTEF, we can utilise a model which captures theessential physics of the S /S
1D H transition which was focused on in the first portionof the main text. Recall that for this transition, the NACV’s between the two electronicadiabatic states are zero, that is (cid:104) φ a ( R ) | ∂ R φ a (cid:48) ( R ) (cid:105) = 0 ∀ a, a (cid:48) in the BO basis, with a, a (cid:48) restricted to S /S This means that matrix elements for the partially Wigner transformedmolecular hamiltonian can be written as ˆ H W ( R, P ) = P M (cid:15) g ( R ) 00 (cid:15) e ( R ) . (41)As described in detail in the first section of this SI, MTEF is rooted in a mean field ap-proximation to the QCLE, which is itself the first order expansion of the partially Wigner27ransformed Liouville von-Neumann equation. Taking eq. (13) to second order provides, ∂ ˆ ρ W ∂t = − i (cid:104) ˆ H W , ˆ ρ W (cid:105) + 12 (cid:16) { ˆ H W , ˆ ρ W } − { ˆ ρ W , ˆ H W } (cid:17) − i (cid:16)(cid:104) ∂ P ˆ H W , ∂ R ˆ ρ W (cid:105) + (cid:104) ∂ R ˆ H W , ∂ P ˆ ρ W (cid:105)(cid:17) (42)Which in our model Hamiltonian eq. (41) becomes, ∂ρ aa (cid:48) W ∂t = − i ( (cid:15) a ( R ) − (cid:15) a (cid:48) ( R )) ρ aa (cid:48) W + (cid:20)
12 ( ∂ R (cid:15) a ( R ) + ∂ R (cid:15) a (cid:48) ( R )) ∂ P − PM ∂ R (cid:21) ρ aa (cid:48) W − i (cid:0) ∂ R (cid:15) a ( R ) − ∂ R (cid:15) a (cid:48) ( R ) (cid:1) ∂ P ρ aa (cid:48) W + O (cid:0) ( m/M ) (cid:1) (43)Such that the error in time propagation resultant from taking only the first order expansion,compared to the second, is proportional to the difference in energy surface curvature.If we take the analytically solvable Displaced Harmonic Oscillator (DHO) model byusing surfaces (cid:15) a ( R ) = ω a ( R − D a ) + E a , we see that for identical surfaces ω e = ω g thatthe nd order and higher terms in the Wigner transformed Liouville von-Neumann equationare zero, rendering the QCLE exact for this case.To demonstrate the effect of varying surface curvature, we took parameters similar toharmonic surface fits to the BO surfaces in 1D H , and for simplicity, took the FC approxima-tion alongside setting µ aa (cid:48) ( R ) = (1 − δ aa (cid:48) ) a.u.. We solve the exact and MTEF-TCF spectrafor the DHO with different values of ω e relative to ω g by propagating for T f = 2 · a.u..In Fig. (7) we see iin the left column that for identical upper and lower surfaces, mean fieldtheory is of course exact, and for varying surfaces, MTEF displays a peak broadening andprepeak features. The origin of this broadening is from an effective damping in the timedependent signal, shown in Fig. (8). Some More Detail on the ICWF Method .0 10.0 11.0 12.002040 I ()[ a ] e = 1 g Exact 9.0 10.0 11.0 12.002040 e = 0.9 g e = 0.8 g I ()[ a ] MTEF 9.0 10.0 11.0 12.001020 9.0 10.0 11.0 12.0010209.0 10.0 11.0 12.0Energy [eV]50050100 I ()[ a ] STEF 9.0 10.0 11.0 12.0Energy [eV]50050100 9.0 10.0 11.0 12.010050050100
Figure 7: Spectra for the DHO model with several excited and ground state surface fre-quencies in each column. Each row compares exact, MTEF-BO and STEF-BO results re-spectively, with the Exact peak placement for each column overlaid across each as verticaldashed lines. 29 .00.50.00.51.0 C [ a . u . ] MTEF [ C ] I ()[ a ] MTEF Spectra e = 1 g C [ a . u . ] I ()[ a ] e = 0.9 g C [ a . u . ] I ()[ a ] e = 0.8 g Figure 8: DHO MTEF time dependent dipole-dipole correlation signal in the left columnand the resulting spectra in the right column, with the relative surface curvature denotedin the right column legend, and exact spectral peaks overlaid as vertical black dashed lines.For clarity, the time dependent signal is curtailed at · a.u..30he conditional wave function (CWF) approach can be developed starting from the fullmolecular wave function for electrons and nuclei, Ψ( r , R , t ) , which can be formally decom-posed in terms of the CWFs of each subsystem: ψ αe ( r , t ) := (cid:90) d R δ ( R α ( t ) − R )Ψ( r , R , t ) , (44) ψ αn ( R , t ) := (cid:90) d r δ ( r α ( t ) − r )Ψ( r , R , t ) . (45)From these definitions one can show that the CWFs, ψ αe ( t ) and ψ αn ( t ) , obey non-Hermitianequations of motion involving complex potentials which are functionals of the full wavefunction and cause the time-evolution of the individual CWFs to be non-unitary. Therecently developed Interacting-CWF (ICWF) method avoids the direct calculation of thesenonlocal complex potentials by positing the following multiconfigurational CWF basis ansatzfor the full many-body wave function: Ψ( r, R , t ) = N c (cid:88) α =1 C α ( t ) ψ αe ( r , t ) ψ αn ( R , t ) . (46)The basis functions in this sum are chosen to be single particle CWFs that satisfy themean-field, or Hermitian, limit of the CWF equations in which the complex potentials triv-ially vanish. The upper limit of the sum, N c , refers to the total number of configurations,which can be stochastically sampled. Including interactions between the trajectories in theensemble through the coefficients C ( t ) = { C ( t ) , ..., C N c ( t ) } corrects the Hermitian-CWFevolution. The time evolution of these coefficients is obtained by inserting eq. (46) directlyinto the TDSE.As described in the text, for the kick spectra adapted ICWF algorithm, the CWFs areinstead selected as eigenstates of the Hermitian propagators, and used as a static basis.The imaginary and real time equations of motion for the expansion coefficient (cid:126)C are then31olved using the respective variational principles, allowing for a completely closed-loopalgorithm for wave function preparation and propagation.To generate the kick spectra, after preparing the ground state (cid:126)C (0) , the relevant degreeof freedom of the kick operator exp ( − iκ ˆ µ ) is applied to each CWF, the Hamiltonian andinverse overlap matrices are reconstructed, and (cid:126)C is propagated to the desired time. Thisprocedure is equivalent to propagating in the interaction representation, with ˆ V I ( t ) = κδ ( t )ˆ µ .Since these matrices are only constructed at time zero, this algorithm is extremely efficient,requiring only the propagation of a N c × vector by a N c × N c matrix. For comparison,the 1D H2 MTEF-kick results reported here required the propagation of , trajectorieseach consisting of × electronic wave functions. With a parallelized implementationand hardware allowing approximately traj/hr, this equates to roughly compute hours.The ICWF N c = 4096 results reported in the main body by contrast require computehours on the same hardware.With increasing non-redundant variational parameters, one is guaranteed to better cap-ture the initial state and minimize the error of time dependent propagation. As an exampleof the convergence properties of ICWF-kick, see Fig. 9. These spectra are the result of utilis-ing only lowest energy hermitian propagator eigenstates and propagating for T f = 1500 a.u.with a mask function W ( x ) = 1 − x + 2 x , for x = t/T f applied to the time signal in theFourier Transform. The more accurate N c = 4096 results in the main body are initialisedusing mixes of various excited eigenstates of the propagators. Theoretical and practicaldevelopments are underway to implement this method in arbitrary ab-initio settings. Complex Absorbing Potentials
Quadratic complex absorbing potentials of the following form were used in all simulationsof the one dimensional H model: W e ( r i ) = − iη (cid:2) ( r i − r l ) Θ( r l − r i ) + ( r i − r r ) Θ( r i − r r ) (cid:3) W n ( R ) = − iη ( R − R r ) Θ( R − R ) , (47)32 I ()[ a ] Exact010 I ()[ a ] N c : 4096010 I ()[ a ] N c : 2048010 I ()[ a ] N c : 10249.2 9.4 9.6 9.8 10.0 10.2 10.4 10.6 10.7 10.9Energy [eV]010 I ()[ a ] N c : 512 Figure 9: Convergence of the 1D H spectra for ICWF-kick with different numbers ( N c ) oflowest energy eigenstate CWF bases. 33here Θ is the Heaviside function, and η was set to 0.1Ha/ a for both subsystems.The electronic CAP cut offs, r l and r r , were placed a from the walls, while the nuclearCAP start was set at R = 5 . a . References (1) May, V.; Kühn, O.
Charge and Energy Transfer Dynamics in Molecular Systems: ThirdEdition ; 2011.(2) Ullrich, C. A. Time-Dependent Density-Functional Theory: Concepts and Applications.
Oxford Graduate Texts ,(3) Wigner, E. On the quantum correction for thermodynamic equilibrium.
Physical Review ,(4) Case, W. B. Wigner functions and Weyl transforms for pedestrians.
American Journalof Physics , , 937–946.(5) Grunwald, R.; Kelly, A.; Kapral, R. Quantum Dynamics in Almost Classical Environ-ments ; 2009.(6) Jasper, A. W.; Zhu, C.; Nangia, S.; Truhlar, D. G. Introductory lecture: Nonadiabaticeffects in chemical dynamics. Faraday Discussions. 2004.(7) Karsten, S.; Ivanov, S. D.; Bokarev, S. I.; Kühn, O. Quasi-classical approaches tovibronic spectra revisited.
Journal of Chemical Physics ,(8) Tully, J. C. Mixed quantum-classical dynamics.
Faraday Discussions ,(9) Kapral, R. Progress in the theory of mixed quantum-classical dynamics. Annual Reviewof Physical Chemistry. 2006. 3410) Lee, M. K.; Huo, P.; Coker, D. F. Semiclassical Path Integral Dynamics: PhotosyntheticEnergy Transfer with Realistic Environment Interactions.
Annual Review of PhysicalChemistry ,(11) Agostini, F.; Min, S. K.; Abedi, A.; Gross, E. K. U. Quantum-Classical NonadiabaticDynamics: Coupled- vs Independent-Trajectory Methods.
Journal of Chemical Theoryand Computation , , 2127–2143.(12) Talotta, F.; Agostini, F.; Ciccotti, G. Quantum Trajectories for the Dynamics in theExact Factorization Framework: A Proof-of-Principle Test. The Journal of PhysicalChemistry A , , 6764–6777.(13) Tully, J. C. Molecular dynamics with electronic transitions. The Journal of ChemicalPhysics ,(14) Donoso, A.; Martens, C. C. Simulation of Coherent Nonadiabatic Dynamics UsingClassical Trajectories.
The Journal of Physical Chemistry A , , 4291–4300.(15) Shalashilin, D. V. Multiconfigurational Ehrenfest approach to quantum coherent dy-namics in large molecular systems. Faraday Discussions ,(16) Mignolet, B.; Curchod, B. F. A walk through the approximations of ab initio multiplespawning.
Journal of Chemical Physics ,(17) Nijjar, P.; Jankowska, J.; Prezhdo, O. V. Ehrenfest and classical path dynamics withdecoherence and detailed balance.
Journal of Chemical Physics ,(18) Albareda, G.; Appel, H.; Franco, I.; Abedi, A.; Rubio, A. Correlated electron-nucleardynamics with conditional wave functions.
Physical Review Letters ,(19) Albareda, G.; Bofill, J. M.; Tavernelli, I.; Huarte-Larranaga, F.; Illas, F.; Rubio, A. Con-ditional born-oppenheimer dynamics: Quantum dynamics simulations for the modelporphine.
Journal of Physical Chemistry Letters ,3520) Albareda, G.; Abedi, A.; Tavernelli, I.; Rubio, A. Universal steps in quantum dynamicswith time-dependent potential-energy surfaces: Beyond the Born-Oppenheimer picture.
Physical Review A ,(21) Albareda, G.; Kelly, A.; Rubio, A. Nonadiabatic quantum dynamics without potentialenergy surfaces.
Physical Review Materials ,(22) Tokmakoff, A. Time-Dependent Quantum Mechanics and Spectroscopy.
Lecture ,(23) Raab, A.; Worth, G. A.; Meyer, H.-D.; Cederbaum, L. S. Molecular dynamics ofpyrazine after excitation to the S2 electronic state using a realistic 24-mode modelHamiltonian.
The Journal of Chemical Physics ,(24) Vendrell, O.; Meyer, H. D. Multilayer multiconfiguration time-dependent Hartreemethod: Implementation and applications to a Henon-Heiles Hamiltonian and topyrazine.
Journal of Chemical Physics ,(25) Yabana, K.; Bertsch, G. Time-dependent local-density approximation in real time.
Physical Review B - Condensed Matter and Materials Physics ,(26) De Giovannini, U.; Brunetto, G.; Castro, A.; Walkenhorst, J.; Rubio, A. Simulatingpump-probe photoelectron and absorption spectroscopy on the attosecond timescalewith time-dependent density functional theory.
ChemPhysChem ,(27) McLachlan, A. D. A variational solution of the time-dependent Schrodinger equation.
Molecular Physics ,(28) Vacher, M.; Bearpark, M. J.; Robb, M. A. Direct methods for non-adiabatic dynam-ics: connecting the single-set variational multi-configuration Gaussian (vMCG) andEhrenfest perspectives.
Theoretical Chemistry Accounts ,(29) Li, X.; Tully, J. C.; Schlegel, H. B.; Frisch, M. J. Ab initio Ehrenfest dynamics.
Journalof Chemical Physics , 3630) Andrea Rozzi, C.; Maria Falke, S.; Spallanzani, N.; Rubio, A.; Molinari, E.; Brida, D.;Maiuri, M.; Cerullo, G.; Schramm, H.; Christoffers, J. et al. Quantum coherence con-trols the charge separation in a prototypical artificial light-harvesting system.
NatureCommunications ,(31) Krumland, J.; Valencia, A. M.; Pittalis, S.; Rozzi, C. A.; Cocchi, C. Understanding real-time time-dependent density-functional theory simulations of ultrafast laser-induceddynamics in organic molecules.
The Journal of Chemical Physics , , 54106.(32) Goings, J. J.; Lingerfelt, D. B.; Li, X. Can Quantized Vibrational Effects Be Obtainedfrom Ehrenfest Mixed Quantum-Classical Dynamics? Journal of Physical ChemistryLetters ,(33) Kapral, R.; Ciccotti, G. Mixed quantum-classical dynamics.
Journal of ChemicalPhysics ,(34) Broeckhove, J.; Lathouwers, L.; Kesteloot, E.; Van Leuven, P. On the equivalence oftime-dependent variational principles.
Chemical Physics Letters ,(35) Lubich, C. On variational approximations in quantum molecular dynamics.
Mathemat-ics of Computation ,(36) Ohta, K. Time-dependent variational principle with constraints for parametrized wavefunctions.
Physical Review A - Atomic, Molecular, and Optical Physics ,(37)
Generalized Inverses: Theory and Applications ; Springer New York: New York, NY,2003; pp 201–256.(38) Kosloff, R.; Tal-Ezer, H. A direct relaxation method for calculating eigenfunctions andeigenvalues of the schrödinger equation on a grid.
Chemical Physics Letters ,(39) Shi, T.; Demler, E.; Ignacio Cirac, J. Variational study of fermionic and bosonic systemswith non-Gaussian states: Theory and applications.
Annals of Physics ,3740) Kreibich, T.; Lein, M.; Engel, V.; Gross, E. K. Even-harmonic generation due tobeyond-born-oppenheimer dynamics.
Physical Review Letters ,(41) Lein, M.; Kreibich, T.; Gross, E. K.; Engel, V. Strong-field ionization dynamics of amodel H2 molecule.
Physical Review A - Atomic, Molecular, and Optical Physics ,(42) Bandrauk, A. D.; Shon, N. H. Attosecond control of ionization and high-order harmonicgeneration in molecules.
Physical Review A - Atomic, Molecular, and Optical Physics ,(43) Tancogne-Dejean, N.; Oliveira, M. J.; Andrade, X.; Appel, H.; Borca, C. H.; Le Bre-ton, G.; Buchholz, F.; Castro, A.; Corni, S.; Correa, A. A. et al. Octopus, a com-putational framework for exploring light-driven phenomena and quantum dynamics inextended and finite systems.
Journal of Chemical Physics ,(44) Gross, E. K. U.; Maitra, N. T.
Introduction to TDDFT ; 2012.(45) Koch, E. E.; Otto, A. Optical absorption of benzene vapour for photon energies from6 eV to 35 eV.
Chemical Physics Letters , , 476–480.(46) Gingell, J. M.; Marston, G.; Mason, N. J.; Zhao, H.; Siggel, M. R. F. On the electronicspectroscopy of benzyl alcohol. Chemical Physics , , 443–449.(47) Borges, I.; Varandas, A. J.; Rocha, A. B.; Bielschowsky, C. E. Forbidden transitions inbenzene. Journal of Molecular Structure: THEOCHEM. 2003.(48) Ridolfi, E.; Trevisanutto, P. E.; Pereira, V. M. Expeditious computation of nonlinearoptical properties of arbitrary order with native electronic interactions in the timedomain. Phys. Rev. B , , 245110.(49) Verlet, L. Computer "experiments" on classical fluids. I. Thermodynamical propertiesof Lennard-Jones molecules. Physical Review ,3850) de Boor, C.
Springer-Verlag, New York ; 2001.(51) Fairlie, D. B. Moyal brackets, star products and the generalised Wigner function.
Chaos,Solitons & Fractals , , 365–371.(52) McKemmish, L. K.; McKenzie, R. H.; Hush, N. S.; Reimers, J. R. Quantum entangle-ment between electronic and vibrational degrees of freedom in molecules. Journal ofChemical Physics ,(53) Yabana, K.; Nakatsukasa, T.; Iwata, J.-I.; Bertsch, G. F. Real-time, real-space im-plementation of the linear response time-dependent density-functional theory. physicastatus solidi (b) ,243