Scholarly article on topic 'Ab initio multiple cloning algorithm for quantum nonadiabatic molecular dynamics'

Ab initio multiple cloning algorithm for quantum nonadiabatic molecular dynamics Academic research paper on "Chemical sciences"

Share paper
Academic journal
The Journal of Chemical Physics
OECD Field of science

Academic research paper on topic "Ab initio multiple cloning algorithm for quantum nonadiabatic molecular dynamics"

The Journal of Chemical Physics

Ab initio multiple cloning algorithm for quantum nonadiabatic molecular dynamics

Dmitry V. Makhov, William J. Glover, Todd J. Martinez, and Dmitrii V. Shalashilin

Citation: The Journal of Chemical Physics 141, 054110 (2014); doi: 10.1063/1.4891530 View online:

View Table of Contents: Published by the AIP Publishing

Articles you may be interested in

On-the-fly ab initio molecular dynamics with multiconfigurational Ehrenfest method J. Chem. Phys. 137, 22A506 (2012); 10.1063/1.4734313

Nonadiabatic ab initio molecular dynamics of photoisomerization in bridged azobenzene J. Chem. Phys. 137, 204305 (2012); 10.1063/1.4767459

Ab initio quantum dynamical study of the multi-state nonadiabatic photodissociation of pyrrole J. Chem. Phys. 135, 154310 (2011); 10.1063/1.3651536

Ab initio nonadiabatic quantum dynamics of cyclohexadiene/hexatriene ultrafast photoisomerization J. Chem. Phys. 124, 084313 (2006); 10.1063/1.2171688

Molecular dynamics simulation with an ab initio potential energy function and finite element interpolation: The photoisomerization of cis-stilbene in solution J. Chem. Phys. 108, 8773 (1998); 10.1063/1.475397

A |p| The Journal of /Al I Chemical PhysicsÈ

Meet The New Deputy Editors

David E. Manolopoulos

James L.M\j¥A Skinnèt

/8\ CrossMark

m^^Htm £dick for updates

Ab initio multiple cloning algorithm for quantum nonadiabatic molecular dynamics

Dmitry V. Makhov,1 William J. Glover,2 Todd J. Martinez,2 and Dmitrii V. Shalashilin1

1 Department of Chemistry, University of Leeds, Leeds LS2 9JT, United Kingdom 2Department of Chemistry and The PULSE Institute, Stanford University, Stanford, California 94305, USA and SLAC National Accelerator Laboratory, Menlo Park, California 94025, USA

(Received 16 April 2014; accepted 6 June 2014; published online 6 August 2014)

We present a new algorithm for ab initio quantum nonadiabatic molecular dynamics that combines the best features of ab initio Multiple Spawning (AIMS) and Multiconfigurational Ehrenfest (MCE) methods. In this new method, ab initio multiple cloning (AIMC), the individual trajectory basis functions (TBFs) follow Ehrenfest equations of motion (as in MCE). However, the basis set is expanded (as in AIMS) when these TBFs become sufficiently mixed, preventing prolonged evolution on an averaged potential energy surface. We refer to the expansion of the basis set as "cloning," in analogy to the "spawning" procedure in AIMS. This synthesis of AIMS and MCE allows us to leverage the benefits of mean-field evolution during periods of strong nonadiabatic coupling while simultaneously avoiding mean-field artifacts in Ehrenfest dynamics. We explore the use of time-displaced basis sets, "trains," as a means of expanding the basis set for little cost. We also introduce a new bra-ket averaged Taylor expansion (BAT) to approximate the necessary potential energy and nonadiabatic coupling matrix elements. The BAT approximation avoids the necessity of computing electronic structure information at intermediate points between TBFs, as is usually done in saddle-point approximations used in AIMS. The efficiency of AIMC is demonstrated on the nonradiative decay of the first excited state of ethylene. The AIMC method has been implemented within the AIMS-MOLPRO package, which was extended to include Ehrenfest basis functions. © 2014 Author(s). All article content, except where otherwise noted, is licensed under a Creative Commons Attribution 3.0 Unported License. []


The dynamics of electronically nonadiabatic transitions is fundamental to many photoinduced processes in chemistry and biochemistry. Examples include photosynthetic light harvesting, the first events in visual reception in rhodopsin, and fluorescent sensors such as the Green Fluorescent Protein (GFP). As the dynamical events in these problems are initiated by electronic excitation, multiple electronic states and nonadiabatic transitions between these states are invariably involved. Furthermore, both nuclear and electronic coherence effects can be crucial, as recently shown in the context of photosynthetic light harvesting.1-5 An acceptable model for nonadiabatic processes must be able to describe coherent phenomena, including also the loss of coherence, i.e., decoherence. Ideally, this is achieved by a quantum mechanical description of the nuclear degrees of freedom. At the same time, such models must be able to describe the potential energy surfaces and their nonadiabatic couplings near the conical intersections6,7 (true degeneracies of two or more electronic states) that often promote electronic transitions. This can be difficult using traditional analytic force fields and provides impetus (even more so than for ground state processes) for ab initio molecular dynamics (AIMD) methods8-10 where the electronic wavefunctions are solved simultaneously with the nuclear dynamics.

Trajectory-based mixed quantum-classical methods are perhaps the simplest approach for modeling nonadiabatic

phenomena in the context of ab initio molecular dynamics. There are two extremes here, according to the form of the potential energy surface used to guide the trajectories. This can be chosen as an average over the populated electronic states (Ehrenfest mean-field method),11-13 in which case a single trajectory represents evolution on multiple electronic states. Mean-field methods can often be quite accurate for very short times, but fail when the wavefunction splits into multiple parts (one or more on each of the involved electronic states). This wavepacket splitting cannot be described because of the chosen ansatz (which demands that the trajectory describes electronic superposition while maintaining a single phase space center). Indeed, this was well recognized by Meyer and Miller, who proposed incorporating a binning procedure to resolve the resulting dynamics into pure electronic states.11 The other extreme is the surface hopping method,14,15 where each trajectory evolves on a pure electronic state and nona-diabatic transitions are described as stochastic hops between different electronic states. The widely used "fewest-switches" surface-hopping15 leverages an Ehrenfest-like dynamics for the quantum amplitudes (but not the phase space centers of the trajectory) in order to determine when to stochastically switch between electronic states. The description of electronic coherence with a single trajectory (even one which always evolves on a pure electronic state) leads to difficulties in describing these coherences correctly. Indeed, the method sometimes seems to enhance coherence (introducing interferences where none should be present) and other times seems to do the


141, 054110-1

opposite (averaging over interferences that should be present), and these characteristics of the method were recognized at the outset.15 On the other hand, surface-hopping provides a much better (but still incorrect) description of the expected equilibration behavior compared to the Ehrenfest method that often leads to an unphysical heating of the quantum mechanical subsystem.16-18

A major difficulty of the mixed quantum-classical trajectory methods is that there is no simple way to improve their description of quantum coherence. Neither the Ehrenfest nor surface-hopping schemes contain any limit where they are guaranteed to converge to the correct quantum mechanical wavepacket evolution. This can be quite unsettling when considering nonadiabatic AIMD methods, which can be quite computationally time-consuming. Ideally, it would be possible to improve the description of quantum mechanical coherence and verify the robustness of the results. Indeed, this is exactly why the first implementations19-22 of nonadi-abatic AIMD eschewed Ehrenfest or surface-hopping methods for the dynamics. Instead, the ab initio multiple spawning (AIMS) method was developed,10,23,24 exploiting the physical intuition underlying the surface-hopping method15 but augmenting it with a basis set expansion perspective. The basis set was expressed in terms of traveling Gaussian wavepackets (trajectory basis functions or TBFs), due to their well-known connections to classical mechanics.25,26 Because of its formulation as a basis set expansion, the AIMS method is guaranteed to converge to the exact solution of the Schrodinger equation for a large enough basis set and this has been demonstrated for model problems.27,28

The development of AIMS occurred during a revival of interest in Gaussian wavepacket methods (motivated largely by their suitability for AIMD), and was quickly followed by the development of the coupled coherent states29,30 (CCS) and variational multiconfigurational Gaussian wavepacket31,32 (vMCG) methods. The vMCG method is perhaps the most flexible of all of these, using a variational procedure to determine the equations of motion for the phase space centers of the TBFs. Unfortunately, this also introduces strong couplings between the equations of motion for the phase space centers and the quantum amplitudes, leading to stiff differential equations which are numerically difficult to integrate. Another issue with vMCG in the context of AIMD is that it is not easy to employ adiabatic representations of the electronic states (which are the natural representation provided by electronic structure methods) and a "diabatization" procedure has to be developed. Nevertheless an ab initio on-the-fly direct dynamics version of vMCG does exist.33

The AIMS10,24,34-36 method uses a much simpler choice for the evolution of the TBFs: their phase space centers evolve classically on a specific electronic state and the basis set is expanded adaptively as needed, in a process called spawning. The spawning becomes particularly important near the intersection with another electronic state, where more Gaussians are spawned on the second PES. The growing adaptive basis follows the dynamics more accurately than a basis of simple classical trajectories.

Recently the Multiconfigurational Ehrenfest (MCE)37-40 dynamics was introduced. This method also employs pre-

determined trajectories to guide the basis set, but the MCE recipe is different from that of AIMS: the MCE technique uses Ehrenfest mean field trajectories with the Hamiltonian averaged (with proper weights) over two or more electronic PESs. The MCE approach is a fully coupled basis set expansion of TBFs (like AIMS) and thus is also guaranteed to converge to the correct solution of the Schrodinger equation for a sufficiently large number of TBFs. It has the advantage that TBFs do not separate as quickly as they would in AIMS, which could aid in convergence at short time. Unfortunately, the average potential guiding the Ehrenfest trajectories is likely to hamper convergence at longer times. An intermediate representation would be desirable and might lead to faster convergence of the underlying dynamics.

In this paper we propose a new algorithm that combines the advantages of AIMS and MCE methods. The new ab initio multiple cloning (AIMC) algorithm introduces Ehrenfest Configuration Cloning (ECC), which expands the basis set of TBFs representing the wavefunction. The cloning procedure is introduced when any of the TBFs becomes strongly mixed (corresponding to a coherent superposition of multiple electronic states with similar amplitudes) and the forces on these electronic states are considerably different (such that representing the evolution of the wavepacket by its average will be a poor approximation). This will typically happen after a TBF passes through/near a conical intersection or region of strong nonadiabatic coupling. Cloning reexpresses the "mixed state" Ehrenfest configuration as a superposition of distinct TBFs, each dominated by a single electronic state. By cloning, we sidestep the prolonged mixed state propagation which is known to lead to erroneous retention of coherence41,42 and violations of detailed balance.15-17 The cloning procedure is analogous to spawning in AIMS, and the AIMS and AIMC methods can be viewed as complementary approaches.

In AIMD methods, the computational bottleneck is invariably the solution of the electronic structure problem to obtain potential energy surfaces, forces, and nonadiabatic couplings. Here, we use the idea of time-displaced basis functions43 or "trains"44 in order to increase the size of the TBF basis set without increasing the number of electronic structure calculations. Time-displaced basis functions follow each other along the same trajectory, i.e., are displaced from each other by a time shift, allowing multiple reuse of electronic structure results. This means of increasing the basis set is particularly efficient in AIMD, where the trains greatly increase the basis set size at almost no extra cost.

Finally, we introduce a bra-ket averaged Taylor expansion method (BAT) for calculating the matrix elements of the potential energy surface and nonadiabatic couplings. The BAT approximation can be computed using only the electronic structure information along each TBF and does not require any further electronic structure calculations at intermediate molecular geometries. A form of this approximation has been previously used for AIMD with MCE dynamics,37 but here we show that it can easily take advantage of both potential energies and gradients along the trajectories. Our AIMC algorithm has been implemented within the AIMS-MOLPRO package.36


A. Working equations

We first provide the equations of the MCE approach. The Hamiltonian of a system of electrons and nuclei is given as

H = Hel(r, R) - 2V«M-1Vr,

MCE dynamics, the wave-function is represented as

(R, r,i)) = £c„(t)|^„(R, r,i)>, (2)

where f n are trajectory basis functions (TBFs) composed of nuclear and electronic parts,

\fn(R, r,i)) = |x„(R,f)) a1(t )\< (r; R)> ,

where Hel is the fixed-nucleus electronic Hamiltonian and where < are electronic wavefunctions (assumed to be M-1 is a diagonal matrix of inverse masses of the atoms. In orthonormal) and xn are Gaussian nuclear basis functions,

x„(R, R„(t), P„(t)) =

exp - a(R - R„)2 + -P„(R - R„) + -y„(t)

where Rn and Pn are the phase space centers of the nth basis function and Ndof is the number of degrees of freedom, i.e., the length of the vector Rn. The parameter a determines the width of the Gaussian in each degree of freedom, and we choose it to be atom dependent and constant, i.e., these are frozen Gaussian basis functions.26 Note that there are three formally redundant sets of variables in Eq. (2): the TBF amplitudes (cn), the Ehrenfest amplitudes for a given TBF (an), and the phase of the frozen Gaussian basis function (yn). The redundancy of the phase is resolved by choosing it semiclas-sically such that

dYn dt

As in AIMS,10,23,24,27 the equations of motion for the phase space centers are not chosen variationally,45 but rather prescribed. Unlike AIMS, we do not use Hamilton's equations for these basis function parameters, but instead choose to use Ehrenfest equations of motion. The electronic representation is taken to be adiabatic in this work (although this is not necessary and a diabatic representation is also possible). The Ehrenfest equations of motion for the phase space centers of a trajectory in the adiabatic representation are given as46 (here shown for the nth TBF)

Rn = M-1Pn,

P = F ,

Fn = E |an|2V«yi(Rn)

+ £ £ (an)*a"dj(Rn) [V;(Rn) - Vj(Rn)], i j=i

where V/R) is the adiabatic electronic energy of the Ith electronic state for molecular coordinates R, and dI7 is the nona-

diabatic coupling vector connecting the Ith and Jth electronic states,

du(Rn) = <i (r; R n)\V«\<j (r; R n))r.

As is commonly done, we assume here and throughout that the second derivative of the electronic wavefunction can be neglected, i.e., that

< |Vr M-1 V« <j >«q.

The equations of motion for the Ehrenfest amplitudes of each TBF are

an = 4EHeJ(RRn, Pn)a"j,

where the electronic Hamiltonian matrix is defined as ' V; (Rn), I = J

HJ(R„, PJ =

. (1Q)

-ih PnM-1d;j(Rn), I = J

Both V;(Rn) and d;j(Rn) can be obtained from standard electronic structure codes for many methods, including complete active space self-consistent field (CASSCF)47,48 and perturbation corrected variants (CASPT2).49-51

With the above prescribed equations for the evolution of the phase and Ehrenfest amplitudes, we can now obtain the equations of motion for the TBF amplitudes cn by inserting the ansatz into the time-dependent Schrödinger equation, giving

J2(fJt Mn(t )>Cn (t)

= - £ (Hmn - ihifjt )\d\fn (t )) ) Cn(t), (11)

(tm\tn> = J2 KxJ^Xn) = Smn

is an overlap matrix and

Hmn = Y! [tfXmVl \ H \anjXnVj) iJ

are the matrix elements of the full (electronic and nuclear) Hamiltonian fromEq. (1). Theterm ((t)| d \ (t)), which originates from the time dependence of the basis functions, can easily be calculated analytically,

dfn dt

+ (Xm\Xn)J2 (amaj, (14)

dxn\ — d — d

"TT = Rn<Xm^^rlXn)+Pn<Xm^lXn

dt I dRn dPn

+ J. Yn <Xm\Xn)'

Detailed formulas for the integrals in Eqs. (14) and (15) are given in Appendices A and B.

It remains to deal with the integrals in Eq. (13). In general, these cannot be evaluated completely analytically, at least when the electronic structure is being solved "on the fly." Neglecting the second derivative coupling term (see above),

the matrix element becomes <Xm(PllH lxnVj )

= &iJ <Xml - 2 VR M VR lXn)+SIJ <XmlVj (R)lXn

-h2 <Xmldij (R)M-1 Vr lXn).

The first term (kinetic energy) is easily calculated analytically (see Appendix A). Approximations are needed for the two remaining terms. A saddle point approximation (SPA) has been used previously in AIMS, expanding V(R) or d/7(R) in a Taylor series centered about the maximum of the xmXn product (the "centroid" of the two basis functions). This has the drawback of requiring electronic structure calculations at the centroids, implying O(N^BF) electronic structure calculations at each time step. In practice, many of these can be (and are) neglected because the basis functions are far apart. The practical scaling in terms of number of electronic structure calculations per time step is thus intermediate between O(NTBF) and O(NTbf). Here, we show that a slightly different (bra-ket averaged Taylor expansion or BAT) approximation can be much more efficient. Instead of a single Taylor expansion, we average the result of two Taylor expansions, each one centered about the maximum of one of the basis functions involved in the matrix element. For the second term in Eq. (16), expansion to first order yields

<XmlVI (R)|Xj ) ^ <XmlXn,

Vi (Rn) + Vi (Rmr

<Xm l (R - Rn)lXn)VRVi (Rj) + <Xm l (R - Rm)lXn)VfiVi (Rm)

Note that the values of the potential energy and its gradient at these two positions (Rn and Rm) were already required in order to propagate the phase space centers in Eq. (6). Thus, there are no extra electronic structure calculations associated with the BAT approximation (in contrast to the SPA). A zeroth order BAT approximation is used for the third term in Eq. (16):

<Xmidi/(R)M1 Vr |Xn> « Yh XmlXnXPmM-Mj (Rm)

+ Pn M-1di7 (Rn )). (18)

An additional advantage of the BAT approximation for the matrix elements in Eq. (13) is that one can separate the evolution of TBFs from the solution of the TBF amplitudes cn Thus, it is possible to independently evolve all of the desired TBFs using Eqs. (5), (6), and (9). Provided all of the electronic structure information at each time step for each of the TBFs was saved, one then has all the needed information to solve Eq. (11) which determines the TBF amplitudes cn and

thus provides the time-evolved nuclear wavefunction. This is difficult with the SPA used previously in AIMS.

We have previously used a zeroth order BAT approximation for the potential energy matrix elements in MCE.37 The first-order BAT approximation should be better and, as noted above, incurs no extra cost. Just as in the SPA,10 it is of course possible to go to second order, but this requires the evaluation of second derivatives of the potential energy surfaces, which is usually quite costly.

A rather subtle point is worthy of some discussion here. The electronic wave function in Eq. (3) depends parametrically on the nuclear coordinates R. However, because of the BAT approximation for the matrix elements, it is only ever calculated at the center of the wave packet, i.e., for R = Rn (t). One can thus reformulate the above MCE approach with the alternative ansatz (compare Eqs. (2) and (3)):

(R, r,t)) = £Cj (t) |fn(R, r,t)),

\1r„(R, r,t)) = \X„(R,t))

V al (t) \^(r; R (t)))

In this case, the nonadiabatic coupling no longer originates from the kinetic energy operator, since the electronic wave-function does not depend on the nuclear coordinates (but only on the basis function parameters). Instead, nonadiabatic coupling terms will appear from the time derivative of i.e., the analog of Eq. (14) will have new terms involving the nonadiabatic coupling. Fortunately, the final equations are unchanged as long as the electronic wave function does not depend too strongly on the nuclear coordinates. A detailed derivation is provided in Appendix B.

B. Basis improvement 1. Coherent state trains

Apart from the BAT approximation for matrix elements, the only approximation in MCE is the finite size of the basis set, i.e., we only use a limited number of TBFs. Thus, it would be useful to expand the basis set, if this were possible at little or no computational expense. Since most of the expense in ab initio molecular dynamics will be in the computation of the electronic structure (energies, gradients, and nonadiabatic couplings), one might thus consider adding time-displaced variants of each TBF into the simulation. This has been suggested previously, and is known as time-displaced basis sets43 or trains.44 The basic idea is illustrated in Figure 1. The upper panel (Figure 1(a)) depicts an evolving set of TBFs, with each TBF shown as a circle evolving in time along a dashed line. We now place additional TBFs along the same trajectory (dashed line) corresponding to each of the existing TBFs, as shown in Figure 1(b). The basis set has been expanded, but no new electronic structure information is needed since the new TBFs all evolve exactly the same as the TBFs in the smaller basis set.

This idea is straightforwardly implemented within the context of a sequential procedure where all TBF evolution is calculated prior to the solution of Eq. (11) for the TBF amplitudes. As an example, consider that we want each "train"

to consist of Ntrain+1 time-displaced TOF^ with Ninitial distinct initial conditions. Our calculation will thus comprise Ninitiai(N¡rain+1) TBFs, while the calculation without time-displacement would have a much smaller basis set of only Ninitiai TBFs. We define a time shift Attrain corresponding to the temporal spacing between the "cars" in the "train" (this should be an integer multiple of the time step used to integrate the trajectories). We calculate the complete time evolution

up to tSim + AttrainNtrain/2 of the Ninitial TBFs and then back-

propagate each of the initial conditions by AttrainNtrain/2. For each initial condition, we thus now have all data required to specify any of the time-displaced TBFs for the entire simulation time tsim. Each of the "cars" in the "train" is now specified by the initial condition TBF it corresponds to and also a time

shift. We choose Attrain to be small enough that neighboring basis functions have significant overlap (and thus significant coupling). However, one should not choose it so small that the

FIG. 1. Graphical illustration of the propagation on a small Ehrenfest trajectory guided basis (a), a basis of Ehrenfest trajectory guided trains (b), and trains with cloning (c).

expansion of the basis is meaningless or that unnecessary linear dependence is introduced (numerically complicating the solution of Eq. (11)). A basis set of several trains, as shown in 1(b), can be made quite large and dense at almost no additional cost (in the context of ab initio molecular dynamics, where the electronic structure is much more costly than the solution of Eq. (11)) because electronic structure information is reused within each train many times.

2. Basis-function cloning

As in any mean field approach, the fundamental difficulty with Ehrenfest propagation occurs when the average is no longer a faithful representation of the whole. In the context of nonadiabatic dynamics, this occurs when an Ehrenfest trajectory has significant population on multiple electronic states with widely varying forces. The force that governs the Ehrenfest trajectory becomes an unphysical average of the different forces on each electronic state and the Ehrenfest trajectory can follow a path that does not resemble the evolution predicted by the Schrodinger equation. In MCE, there are multiple

interacting Ehrenfest trajectories. Thus, the failure of any of these Ehrenfest trajectories could be ameliorated in principle by the solution of Eq. (11), if there were a large enough number of TBFs. Unfortunately, experience suggests that practical calculations are far from the number of TBFs needed to remedy this defect. The solution we propose here follows the spawning idea introduced in AIMS, i.e., adaptive expansion of the basis set of TBFs. We expand the basis set by "cloning" a new TBF when the Ehrenfest force differs significantly from the force that would be felt on one of the electronic states that the Ehrenfest amplitudes predict to have significant population. The difference between the force on the Ith state and the averaged force is given by

= VRVi -

{an) *aj VRVj.

We also want to consider the population on the Ith state, since there is no point in cloning to the Ith state if the corresponding Ehrenfest amplitudes indicate it is insignificantly populated. Thus, we define the "breaking force" which will trigger cloning as

F^ = (an)* ai AFi n

Multiplication by the inverse mass matrix gives an acceleration, and we use the norm of this acceleration as the trigger which signals the need to clone, i.e., when |M-1Fjrn| > Sclone, a new TBF will be cloned from the nth TBF onto the Ith electronic state. The cloning threshold Sclone needs to be determined empirically. Choosing a value which is too small will lead to a large rate of basis set expansion with no benefit, while a value which is too large will lead back to the original MCE method.

As in spawning, we choose the basis set expansion such that it does not alter the wavefunction at the time of spawning/cloning. In AIMS, this simply means the newly spawned TBF enters the calculation with a vanishing amplitude. For AIMC, the procedure is slightly different. The Ehrenfest TBF corresponds to a superposition of multiple electronic states, evolving on an averaged potential energy surface. When we clone, we create two copies of the same TBF, but we then adjust the Ehrenfest amplitudes such that one of the clones corresponds to a pure electronic state and the other contains all the remaining electronic states (in the case of two electronic states which we consider here, both clones correspond to pure electronic states). The sum of the two cloned TBFs is then exactly equivalent to the original "parent" TBF. A schematic description of the cloning procedure is shown in Figure 1(c).

Prior to cloning, the TBF is represented as an Ehrenfest configuration f n as in Eq. (3). After cloning, this TBF is replaced by two TBFs, fn and f:

lfn) = lXn

x Vpi ) + ^ 0 x VP.

lfU) = lXn) 0 x \pi) +

V1 - lan l 2

J2a" 1 Pj>

The TBF amplitudes for the two new TBFs are set to

c'n = Cj l an | , (25)

= cJ 1 -

As can be easily verified, the sum of the two clone TBFs is precisely the same as the parent TBF:

cnlfn) = c'nlfj) + cjlfn

However, the ensuing dynamics will be different—the positions and momenta of the fn clone will evolve on the Ith pure electronic state and those of the fH clone will evolve on a mixed electronic state that has no contribution from the Ith electronic state. When the nonadiabatic coupling increases, both of these TBFs will become mixed states again, and they may be subject to further cloning.

In the presently coded version of the AIMC algorithm, we also demand that the nonadiabatic coupling be small before cloning can occur, i.e., |dI7| < Snac. This ensures that cloning happens outside the region of strong transitions, limiting the rate of basis set growth. This is not strictly necessary and future work will investigate the effects of lifting this restriction. In the specific case of two electronic states, the breaking force that triggers cloning (Eq. (22)) can be written as

Fbr 1,n

Fbr 2,n

= \ananI VV - V2).

This expression makes it clear that IFf' | is largest when both states have equal amplitudes, vanishes when only one state is populated, and increases as the difference between the forces on the two states increases. These are precisely the cases where the Ehrenfest dynamics of the TBF becomes suspect. Of course, the algorithm for cloning can be further refined and this will be investigated in future work.


In order to benchmark the AIMC method, we performed calculations of the photodynamics of ethylene after n * excitation. As ethylene is the simplest molecule with a C=C double bond, it has been extensively studied both experimentally and computationally.22,50-58 The AIMC and AIMS simulations were carried out with a modified version of AIMS-MOLPRO,36 which was extended to include Ehrenfest dynamics for the TBFs. Electronic structure calculations were performed with the CASSCF method using an active space of two electrons in two orbitals and equally weighted state-averaging over the three lowest singlet states, i.e., SA-3-CAS(2/2). The electronic wavefunction is represented with the 6-31G** basis set. The width a of the Gaussian trajectory functions was taken to be 4.7 bohr-2 for hydrogen and 22.7 bohr-2 for carbon atoms, as used in the AIMS-MOLPRO package.59 The simulations are restricted to two electronic states—S0 and Sj.

We compare the standard MCE method (without cloning) to the AIMC and AIMS methods. For MCE and AIMC, we also compare the influence of trains on the results. Matrix elements for AIMC were evaluated using the BAT

approximation as discussed above, while the AIMS calculations use the SPA for this purpose. For the AIMC simulations, we set the cloning thresholds 5done and 5nac to 5 x 10-6 and 2 x 10-3 atomic units, respectively. The similar "spawning threshold" in AIMS was set to 4.0 atomic units. When the norm of the nonadiabatic coupling vector exceeds the spawning threshold, the AIMS calculations will spawn a new TBF (see previous work10 for details). In the AIMS calculations, spawning is prohibited for TBFs that are insignificantly populated (norm of TBF amplitude less than 0.1). A similar procedure controls the growth of the basis set in the AIMC calculations, where we allow each TBF to clone at most 3 times.

In all cases, the simulations are initiated with 200 distinct initial TBFs, sampled randomly from a Wigner distribution corresponding to the ground vibrational state in the harmonic approximation, and then projected onto the S1 state. For the MCE and AIMC calculations with train basis sets, we use a time-shift Attrain of 25 atomic time units («0.6 fs), which corresponds to a nearest neighbor overlap of «0.6. The number of basis functions in each train, Ntrain, is set to 100 for MCE/AIMC calculations. Initially each train is independent from all other trains. For calculations without cloning, each of the 200 trains is propagated independently so that the size of the basis set remains equal to 100 and does not change during the simulations. In contrast, the size of basis set increases every time cloning occurs in AIMC. Moreover, each cloning event increases the number of trains, which are coupled with each other. Figure 2 shows the growth of the average number of TBFs in the AIMC simulations. Each of the initial TBFs gives rise to an average of «4 further TBFs through the cloning process. In the AIMS calculations, each initial condition spawns an average of 4.1 new TBFs, giving a similar basis set growth rate compared to AIMC.

An example of wavefunction spreading over a train basis set with Ntrain = 200 for ethylene is shown in Figure 3. The initial population is localized on the middle "car" of the set of TBFs in a "train" corresponding to the same initial condition. As time progresses, the window of TBFs along the trajectory that are being included in the calculation also shifts forward.

■g 400

Ï? 200 CO

100 Time / fs

FIG. 2. Growth of the average number of TBFs in an AIMC calculation (with 200 initial conditions) using trains with Ntmin = 100. On average, each TBF produces an additional 4 TBFs through cloning.

50 100 150

Basis function number n

FIG. 3. Time evolution of ethylene's wave function in a train basis set. The basis is moving along a quasi-classical trajectory so that the maximum of amplitude remains in the middle of the train.

Thus, the corresponding population of the "cars" stays localized in the middle TBF, but spreads over the cars in the train. As Figure 3 shows, for this case, the wavefunction remains largely localized on the 100 TBFs in the middle of the train. On the basis of this and similar preliminary studies, we concluded that Ntrain = 100 would be sufficient for the MCE and AIMC calculations in this case.

Figure 4 compares the average ground-state population as a function of time given by MCE and AIMC calculations, with and without basis set expansion using trains. The population was averaged over 200 non-interacting trajectories with initial conditions randomly generated from the ground-state vibrational Wigner distribution. The ground-state population evolution predicted by calculations, with and without train basis functions, is very similar. The results are somewhat smoother when the train basis set is used, but the difference

100 Time / fs

FIG. 4. Ground-state population dynamics of ethylene following excitation to the spectroscopic 1Bu (n*) state. We compare calculations using Ehrenfest dynamics with a basis of one Gaussian function (MCE), coherent state trains (MCE/train, Ntmin = 100), one Gaussian function with basis function cloning (AIMC), and coherent state trains with basis function cloning (AIMC/train, Ntrain = 100). Results are averaged over 200 initial conditions and error bars (95% confidence interval) are shown.

100 Time / fs

FIG. 5. Comparison of time evolution of ground state population after n* photoexcitation from AIMS (red) and AIMC (using trains, Ntmin = 100, blue) calculations with SA-3-CAS(2/2)/6-31G**. Results are averaged over 200 initial conditions and error bars (95% confidence interval) are shown.

is within the error bars. One can conclude that, at least for this problem where the conical intersection is highly peaked (and thus nonadiabatic transitions are both ultrafast and ultraefficient), the benefits of the train basis set expansion are not pronounced. Of course, one should also recall that this basis set expansion carries no cost (in the context of ab initio molecular dynamics). Future work will investigate the degree to which the representation of the nuclear wavefunction is improved by the use of time-displaced basis set expansion.

The initial population dynamics predicted by MCE and AIMC are similar. However, as the ground state population increases the predictions begin to deviate, and by 100 fs they are quite different. The relaxation rate predicted by AIMC is significantly faster. This behavior is as expected - as the population on the ground and excited states becomes nearly equal - the Ehrenfest dynamics of the TBFs will become that of the average of the two electronic states. This will tend to keep the TBFs in the region near the conical intersection longer and population transfer in both directions (to the upper state and to the lower state) will be equally probable. This is a manifestation of the violations of detailed balance which are a well known difficulty in Ehrenfest dynamics.16,17 The cloning procedure in AIMC solves this problem by allowing the TBFs to separate and evolve on adiabatic states.

In Figure 5, we compare the results from AIMS and AIMC calculations. The AIMS results for ethylene have been previously verified by direct comparison to experiment using CASPT2 for the electronic structure.52 Gratifyingly, the AIMC and AIMS methods agree almost quantitatively on the population dynamics in this case, and are certainly indistinguishable within the error bars shown.


We propose a new algorithm for quantum molecular dynamics (AIMC) which combines the best features of two existing methods: MCE and AIMS. The algorithm uses an ensemble of Ehrenfest trajectories to guide the basis of nuclear

Gaussian coherent states in place of the classical equations of motion used in AIMS. Thus, the trajectories move on an averaged potential energy surface. As is well-known,15 this can lead to problematic behavior when the electronic states in the average are very different. We alleviate this by expanding the basis set, similar to spawning, in a procedure we call cloning. As for spawning in AIMS, cloning does not change the wavefunction—it merely reexpresses it in an expanded representation. However, it does lift the restrictions of the ansatz, i.e., the newly cloned trajectory basis functions now feel the force appropriate to an adiabatic state and are free to separate.

We have also introduced a new approximation (BAT) for the required matrix elements that is an alternative to the SPA used in AIMS. The BAT approximation averages truncated Taylor expansions centered on the two TBFs involved in the matrix element. By comparison, the SPA evaluates the matrix element by a truncated Taylor series expansion about a single point intermediate between the position centers of the two trajectory basis functions. The advantage of the BAT approximation for the matrix element is that it can be computed without any additional information beyond the energies, gradients, and nonadiabatic couplings at the position centers of each TBF. There is no need for evaluation of these quantities at intermediate positions. Hence the number of electronic structure calculations needed scales linearly with the number of TBFs, as compared to a formally quadratic scaling when the SPA is used.

We have implemented the AIMC algorithm in the AIMS module of the MOLPRO package and tested it on the benchmark case of photodynamics in ethylene. We showed that the AIMS and AIMC algorithms give comparable results for the excited state lifetime in this molecule. In this test case, there does not appear to be any benefit from the use of time-displacement, i.e., trains, to increase the size of the nuclear basis set. However, there may be advantages to time-displacement in other cases or for more sensitive properties than the electronic population (e.g., the absorption spectrum60,61 or multi-dimensional spectroscopies).


This work was supported by EPSRG Grant No. EP/J001481/1 International Collaboration in Chemistry: New First Principles Methods for Nonadiabatic Dynamics and NSF CHE-11-24515. We also gratefully acknowledge Rick Heller who suggested combining features of the AIMS and MCE methods during HellerFest in October 2010.


The equations for integrals involving complex Gaussian basis functions, which appear in our formalism, are as follows:

( a n 1 0 i — —--

<Xm|Xn> = exp i - 2 AR2 - — AP2 + -(PmRm - PnRn

+ RAP) + - Ay),

(Xm\R\Xn) =

AP + R

(Xm\Xn ), (A2)

(Xm\^R\Xn) =

--P - aAR

(Xm \ Xn

(Xm\n \Xn) =

-2aARP - a + a2AR2 - Ph fi2

(Xm \ Xn

(Xm \ ,tt \ Xn )

(Xm\ \Xn)

- P + aAR fi

(Xm \ Xn

l i --AP--AR

4afi2 2h

(Xm \ Xn

where AR = Rn - Rm, AP = Pn - Pm, Ay = Y„ - Ym,

■¡7 R +R =; P +P

IJ _ n m IJ _ n m


Let us consider the ansatz of Eqs. (19) and (20) for the full wave function

№ (t )) = £ Cn (t ) \fn(t ))

(f m\H \f n

= £ «Yj Xm\{pm\ ( K (r, R) - 2m vR) \PJ )X)

= £ «Y J {(Xm№\Hel (r, R )\pJ )\Xn) - Tmn{vm\<P"j )) ,

where Tmn = -h2 /2M{xm\^R\xn) are the kinetic energy matrix elements.

Comparing to Eq. (14), Eq. (B6) contains an extra term which includes {<pm\ d \W"); likewise, comparing to Eq. (16), Eq. (B7) does not have the term containing {<pN\VR\<pM) and {WN\V2R\wm). The difference is due to the fact that \wi) = W/(r; R)) is a function of R and, therefore, is acted upon by the nuclear kinetic energy operator. On the other hand, the function \wn) = W (r; Rn)) depends on the time dependent Rn and therefore is subject to time differentiation.

Approximations must again be introduced in order to evaluate the matrix elements. As in the previous case, we assume that on the scale of the Gaussian width, the electronic wave function depends weakly on Rt, the position of the nuclei. Then for pairs of basis functions with sufficient nuclear overlap we can write

Pm \ J * j

where the basis functions are the Ehrenfest configurations composed, as before, of nuclear and electronic parts,

\ fn(t)) = \Xn(t))J2 an (t)\pf

but unlike the standard ansatz in Eq. (2), the electronic part here is "attached" to the center of the nuclear wave packet:

\ pJ) = \ cPl (r; Rn(t ))).

Substituting the ansatz (B1) into the Schrodinger equation i h = H ^ one obtains

- (fm \ fn

= £ (fm \ H \ fn)Cn - ifiJ2 (fm \ d I fn)Cn. (B4)

The matrix elements are as follows:

(fm \ fn)=T, (am) ^npm \ Pn )(Xm \ X,

(fm \ Tt I fn) = £ (am) JTP) (Xm \ Xn d IJ

+ £ (am)*aJ}{pm \ Jt \PJ)(Xm \Xn) IJ

+ £ (amy-a^Pm jXm \ ^ \Xn), (B6)

m dt \J

pmwAPm)+p \-\PJ

Rmdij(Rm) + Rndij(Rn)

i.e., the nonadiabatic coupling matrix element between two wave functions localized at points Rm and Rn is an average between the two "local" nonadiabatic couplings. For m = n it is not an approximation.

We assume the following approximation for the Born-Oppenheimer potential energy matrix elements:

(xm \ {vm\Hel (r, R)\pJ) \ Xn)

« {Xm \ <^m(r; R) + (Rm - R)VR<(r; R) \ Hei(r, R)

\ pj (r; R) + (Rn - R)VrpJ (r; R)) \ Xn)

« <Xm \ <Pm(r, R) \ Kl(r, R) \Pj(r, R)) \Xj)

+<Xm \ <pm(r, R) \Kl(r, R) \ (Rj - R)VrpJ(r, R)) \Xn). + <Xm \ <(Rm - R)VRPm(r, R) \ VLel(r, R) \pj(r, R)) \Xn) = SIj<Xm \ VI(R) \Xn)

+<Xm\Vj(R)(Rj - R)<pm(r, R) \VrpJ(r, R))\ Xn) + <Xm \ Vi(R)(Rm - R)<VRPm(r, R) \pj(r, R)) \Xj).

The last two terms can be disregarded because far from intersection the nonadiabatic coupling terms are small and near the intersection where VI - Vu the two terms cancel out. Then using the BAT approximation of Eq. (17) for the first term in (B10) we get

^0^1 (Vr, r, R)|^n))|Xn> - Su^n^ V (Rm) + Vl (Rn ) )

/<Xm|(R - Rm)|Xn>VRVi(Rm) +QU(R - Rj|Xn>VRVi(R„ A _ ^

As before

\xJ = Rj [Xm

d d R„

IXn) + Pn (Xm

d P„

Altogether, under the same assumption of smoothness of the electronic wave function and its weak dependence on the nuclear coordinates, the equations for the amplitudes cn become the same as for those in the ansatz of Eqs. (2) and (3). The only difference is that the nonadiabatic coupling term

R d,,(R )+R d,,(R )

m ij V m /_j ij V j J

originates now from

HKXm \ a"X„

the time dependence of the electronic basis functions but not the kinetic energy operator. Although currently both forms of the total wave function given by Eqs. (2) and (3) or Eqs. (19) and (20) obey the same equation, it is not clear which of the two approaches will be more convenient for generalization in future, and for the sake of future reference we present here both of them.

1 Y.-C. Cheng and G. R. Fleming, "Dynamics of light harvesting in photosynthesis," Ann. Rev. Phys. Chem. 60, 241 (2009).

2V. Blanchet, M. Z. Zgierski, T. Seideman, and A. Stolow, "Discerning vibronic molecular dynamics using time-resolved photoelectron spec-troscopy," Nature (London) 401, 52 (1999).

3E. Collini, C. Y. Wong, K. E. Wilk, P. M. G. Curmi, P. Brumer, and G. D. Scholes, "Coherently wired light-harvesting in photosynthetic marine algae at ambient temperature," Nature (London) 463, 644 (2010).

4J. L. Herek, W. Wohlleben, R. J. Cogdell, D. Zeidler, and M. Motzkus, "Quantum control of energy flow in light harvesting," Nature (London) 417, 533 (2002).

5D. G. Kuroda, C. P. Singh, Z. Peng, and V. D. Kleiman, "Mapping excited-state dynamics by coherent control of a dendrimers photoemission efficiency," Science 326, 263 (2009).

6D. R. Yarkony, "Diabolical conical intersections," Rev. Mod. Phys. 68, 985 (1996).

7B. Lasorne, G. A. Worth, and M. A. Robb, "Excited-state dynamics," Wiley Interdiscip. Rev.-Comput. Mol. Sci. 1, 460 (2011).

8R. Car and M. Parrinello, "Unified approach for molecular dynamics and density-functional theory," Phys. Rev. Lett. 55, 2471 (1985).

9M. E. Tuckerman, P. J. Ungar, T. von Rosenvinge, and M. L. Klein, "Ab initio molecular dynamics simulations," J. Phys. Chem. 100, 12878 (1996).

10M. Ben-Nun and T. J. Martinez, "Ab initio quantum molecular dynamics," Adv. Chem. Phys. 121, 439 (2002).

11 H.-D. Meyer and W. H. Miller, "A classical analog for electronic degrees of freedom in nonadiabatic collision processes," J. Chem. Phys. 70, 3214 (1979).

12G. D. Billing, "On the use of Ehrenfest's theorem in molecular scattering," Chem. Phys. Lett. 100, 535 (1983).

13J. W. Negele, "The mean-field theory of nuclear structure and dynamics," Rev. Mod. Phys. 54, 913 (1982).

14J. C. Tully and R. K. Preston, "Trajectory surface hopping approach to nonadiabatic molecular collisions: The reaction of H+ with D2," J. Chem. Phys. 55, 562(1971).

15J. C. Tully, "Molecular dynamics with electronic transitions," J. Chem. Phys. 93, 1061 (1990).

16P. V. Parandekar and J. C. Tully, "Mixed quantum-classical equilibrium," J. Chem. Phys. 122, 094102 (2005).

17P. V. Parandekar and J. C. Tully, "Detailed balance in Ehrenfest mixed quantum-classical dynamics," J. Chem. Theor. Comput. 2, 229 (2006).

18J. R. Schmidt, P. V. Parandekar, and J. C. Tully, "Mixed quantum-classical equilibrium: Surface-hopping," J. Chem. Phys. 129, 044104 (2008).

19T. J. Martinez and R. D. Levine, "First-principles molecular dynamics on multiple electronic states: A case study of NaI," J. Chem. Phys. 105, 6334 (1996).

20T. J. Martinez and R. D. Levine, "Dynamics of the collisional electron transfer and femtosecond photodissociation of NaI on ab initio electronic energy curves," Chem. Phys. Lett. 259, 252 (1996).

21T. J. Martinez, "Ab initio molecular dynamics around a conical intersection: Li(2p)+H2," Chem. Phys. Lett. 272, 139 (1997).

22M. Ben-Nun and T. J. Martinez, "Ab initio molecular dynamics study of cis-trans photisomerization in ethylene," Chem. Phys. Lett. 298, 57 (1998).

23T. J. Martinez, M. Ben-Nun, and G. Ashkenazi, "Classical/quantal method for multistate dynamics: A computational study," J. Chem. Phys. 104, 2847 (1996).

24T. J. Martinez, M. Ben-Nun, and R. D. Levine, "Multi-electronic-state molecular dynamics: A wave function approach with applications," J. Phys. Chem. 100, 7884 (1996).

25E. J. Heller, "Time-dependent approach to semiclassical dynamics," J. Chem. Phys. 62, 1544 (1975).

26E. J. Heller, "Frozen Gaussians: A very simple semiclassical approximation," J. Chem. Phys. 75, 2923 (1981).

27M. Ben-Nun and T. J. Martinez, "Nonadiabatic molecular dynamics: Validation of the multiple spawning method for a multidimensional problem," J. Chem. Phys. 108, 7244 (1998).

28M. Ben-Nun and T. J. Martinez, "A continuous spawning method for nona-diabatic dynamics and validation for the zero-temperature spin-boson problem," Isr. J. Chem. 47, 75 (2007).

29D. V. Shalashilin and M. S. Child, "Time dependent propagation in phase space," J. Chem. Phys. 113, 10028 (2000).

30D. V. Shalashilin and M. S. Child, "The phase space CCS approach to quantum and semiclassical molecular dynamics for high-dimensional systems," Chem. Phys. 304, 103 (2004).

31 I. Burghardt, H. D. Meyer, and L. S. Cederbaum, "Approaches to the approximate treatment of complex molecular systems by the multiconfigura-tion time-dependent Hartree method," J. Chem. Phys. 111, 2927 (1999).

32G. A. Worth and I. Burghardt, "Full quantum mechanical molecular dynamics using Gaussian wavepackets," Chem. Phys. Lett. 368, 502 (2003).

33G. A. Worth, M. A. Robb, and B. Lasorne, "Solving the time-dependent Schrodinger equation for nuclear motion in one step: Direct dynamics of non-adiabatic systems," Mol. Phys. 106, 2077 (2008).

34M. Ben-Nun and T. J. Martinez, "Photodynamics of ethylene: Ab initio studies of conical intersections," Chem. Phys. 259, 237 (2000).

35S. Yang, J. D. Coe, B. Kaduk, and T. J. Martinez, "An 'optimal' spawning algorithm for adaptive basis set expansion in nonadiabatic dynamics," J. Chem. Phys. 130, 134113 (2009).

36B. G. Levine, J. D. Coe, A. M. Virshup, and T. J. Martiinez, "Implementation of ab initio multiple spawning in the Molpro quantum chemistry package," Chem. Phys. 347, 3 (2008).

37K. Saita and D. V. Shalashilin, "On-the-fly ab initio molecular dynamics with multiconfigurational Ehrenfest method," J. Chem. Phys. 137, 22A506 (2012).

38 D. V. Shalashilin, "Multiconfigurational Ehrenfest approach to quantum coherent dynamics in large molecular systems," Faraday Discuss. 153, 105 (2011).

39 D. V. Shalashilin, "Nonadiabatic dynamics with the help of multiconfigurational Ehrenfest method: Improved theory and fully quantum 24D simulation of pyrazine," J. Chem. Phys. 132, 244111 (2010).

40D. V. Shalashilin, "Quantum mechanics with the basis set guided by Ehrenfest trajectories: Theory and application to spin-boson model," J. Chem. Phys. 130, 244101 (2009).

41E. R. Bittner and P. J. Rossky, "Quantum decoherence in mixed quantum-classical systems: Nonadiabatic processes," J. Chem. Phys. 103, 8130 (1995).

42 O. V. Prezhdo and P. J. Rossky, "Mean-field molecular dynamics with surface hopping," J. Chem. Phys. 107, 825 (1997).

43 M. Ben-Nun and T. J. Martinez, "Exploiting temporal nonlocality to remove scaling bottlenecks in nonadiabatic quantum dynamics," J. Chem. Phys. 110,4134 (1999).

44D. V. Shalashilin and M. S. Child, "Basis set sampling in the method of coupled coherent states: Coherent state swarms, trains and pancakes," J. Chem. Phys. 128, 054102 (2008).

45M. Ronto and D. V. Shalashilin, "Numerical implementation and test of the modified variational multiconfigurational Gaussian method for high-dimensional quantum dynamics," J. Phys. Chem. A 117(32), 6948 (2013).

46N. L. Doltsinis and D. Marx, "First principles molecular dynamics involving excited states and nonadiabatic transitions," J. Theor. Comput. Chem. 01,319(2002).

47B. O. Roos, "The complete active space self-consistent field method and its applications in electronic structure calculations," Adv. Chem. Phys. 69, 399 (1987).

48B. H. Lengsfield III and D. R. Yarkony, "Nonadiabatic interactions between potential energy surfaces: Theory and applications," Adv. Chem. Phys. 82, 1 (1992).

49B. O. Roos, "Theoretical studies of electronically excited states using mul-ticonfigurational perturbation theory," Acc. Chem. Res. 32, 137 (1999).

50T. Mori, W. J. Glover, M. S. Schuurman, and T. J. Martinez, "Role of Ryd-berg states in the photochemical dynamics of ethylene," J. Phys. Chem. A 116, 2808 (2012).

51H. Tao, B. G. Levine, and T. J. Martinez, "Ab initio multiple spawning dynamics using multi-state second-order perturbation theory," J. Phys. Chem. A 113, 13656 (2009).

52H. Tao, T. K. Allison, T. W. Wright, A. M. Stooke, C. Khurmi, J. van Tilborg, Y. Liu, R. W. Falcone, A. Belkacem, and T. J. Martinez, "Ultra-fast internal conversion in ethylene. I. The excited state lifetime," J. Chem. Phys. 134,244306(2011).

53R. J. Sension and B. S. Hudson, "Vacuum ultraviolet resonance Raman studies of the excited electronic states of ethylene," J. Chem. Phys. 90, 1377 (1989).

54E. F. Cromwell, A. Stolow, M. J. J. Vrakking, and Y. T. Lee, "Dynamics of ethylene photodissociation from rovibrational and transla-tional energy distributions of H2 products," J. Chem. Phys. 97, 4029 (1992). 2

55V. Stert, H. Lippert, H. H. Ritze, and W. Radloff, "Femtosecond time-resolved dynamics of the electronically excited ethylene molecule," Chem. Phys. Lett. 388, 144 (2004).

56A. Viel, R. P. Krawczyk, U. Manthe, and W. Domcke, "Photoinduced dynamics of ethene in the N, V, and Z valence states: A six-dimensional nonadiabatic quantum dynamics investigation," J. Chem. Phys. 120, 11000 (2004).

57A. Viel, R. P. Krawczyk, U. Manthe, and W. Domcke, "The sudden-polarization effect and its role in the ultrafast photochemistry of ethene," Ang. Chem. Int. Ed. 42, 3434 (2003).

58K. Kosma, S. A. Trushin, W. Fuss, and W. E. Schmid, "Ultrafast dynamics and coherent oscillations in ethylene and ethylene-d(4) excited at 162nm," J. Phys. Chem. A 112, 7514 (2008).

59A. L. Thompson, C. Punwong, and T. J. Martinez, "Optimization of width parameters for quantum dynamics with frozen Gaussian basis sets," Chem. Phys. 370, 70 (2010).

60M. Ben-Nun and T. J. Martinez, "Electronic absorption and resonance Raman spectroscopy from ab initio quantum molecular dynamics," J. Phys. Chem. A 103, 10517(1999).

61 M. Sulc, H. Hernandez, T. J. Martinez, and J. Vanicek, "Relation of exact Gaussian basis methods to the dephasing representation: Theory and application to time-resolved electronic spectra," J. Chem. Phys. 139, 034112 (2013).