Direct methods for non-adiabatic dynamics: connecting the single-set variational multi-configuration Gaussian (vMCG) and Ehrenfest perspectives

In this article, we outline the current state-of-the-art “on-the-fly” methods for non-adiabatic dynamics, highlighting the similarities and differences between them. We derive the equations of motion for both the Ehrenfest and variational multi-configuration Gaussian (vMCG) methods from the Dirac–Frenkel variational principle. We explore the connections between these two methods by presenting an alternative derivation of the vMCG method, which gives the Ehrenfest equations of motion when taking the appropriate limits.


Introduction
Chemistry is ultimately about the reorganisation of electrons, which causes the breaking and formation of bonds i.e. the motion of atoms resulting in new molecules. The time evolution of a (non-relativistic) molecular system is determined by the time-dependent Schrödinger equation [1]: where Φ(r, R, t) is the total molecular wave function, and r and R are the electronic and nuclear coordinates respectively. The Hamiltonian operator reads:Ĥ(r, R) =T n (R) +T e (r) + V n−e (r, R) =T n (R) +Ĥ e (r, R), witĥ T n (R) the kinetic energy operator of the nuclei,T e (r) the kinetic energy operator of the electrons, V n−e (r, R) which includes all inter-particles interactions (electron-electron, nucleus-nucleus and electron-nucleus) andĤ e (r, R) = T e (r)+V n−e (r, R) the electronic Hamiltonian for fixed nuclei R. Unfortunately, this equation is impossible to solve analytically for more than two particles, i.e. the hydrogen atom. The field of theoretical chemistry is thus dominated by developments of numerical and approximate methods, that can be used to treat other atoms and molecules.
Conventional Born-Oppenheimer molecular dynamics only allows one to simulate nuclear motion on a single potential energy surface and therefore does not describe non-adiabatic processes involving non-radiative electronic transitions. Standard grid-based quantum mechanical simulations are expensive computationally because of their exponential scaling with the system size. A separate bottleneck is the computation and fitting of the potential energy surfaces prior to any dynamics calculation. Reducing the number of nuclear degrees of freedom of the system is sometimes done to make the calculation feasible but the validity of this approximation is limited [2][3][4]. "On-the-fly" dynamics methods have been developed to address these issues. They are referred to as direct dynamics methods since the potential energy surfaces are calculated as needed along nuclear trajectories, thus avoiding the pre-computation of globally fitted surfaces, and sampling only the relevant regions of the potential energy surfaces. These nuclear trajectories are used to describe the nuclear wave packet motion, i.e. the nuclear wave packet is expanded in the basis of nuclear trajectories via expansion coefficients.
The several on-the-fly methods which are able to describe non-adiabatic dynamics treat the electrons quantum mechanically. The major feature that differentiates between them is the treatment of the nuclear motion itself through the basis nuclear trajectories. Do the basis trajectories obey quantum or classical mechanics? Are the basis trajectories coupled or independent? Does each basis trajectory evolve on a single potential energy surface (at a time) or does it follow the gradient of a superposition of electronic states and therefore evolve on an effective potential energy surface? This last point raises an important distinction. In the former case, a different set of basis trajectories is used for each electronic state. In technical terms, a multi-set formalism is used. In the latter case, one set of basis trajectories is used to treat the dynamics in all electronic states: a single-set formalism is used, where each basis trajectory has a time-dependent coefficient for every electronic states considered. Figure 1 attempts to summarise the relationships among the several theoretical strategies and methods that have been developed and applied by different groups: 1 -The direct dynamics variational multi-configuration Gaussian (DD-vMCG) method of Worth, Burghardt and Lasorne [5][6][7][8][9] describes the nuclear wave packet using a basis set of time-dependent Gaussian functions that evolve quantum mechanically and are variationally coupled. This means that the  evolution of not only the wave packet expansion coefficients, but also of the parameters (eg. mean position and momentum) of every Gaussian basis function (GBF) is determined by the time-dependent Schrödinger equation. Two different formalisms exist: (i) the multi-set formalism in which a different set of GBF is used for each electronic state (there can be a different number of GBF for each state) and (ii) the single-set formalism in which the same set of GBF is used for all electronic states (the electronic states are included as an extra degree of freedom described by a finite basis labelling the states).
-Full multiple spawning (FMS) and ab initio multiple spawning (AIMS) methods, developed by Martínez, Ben-Nun and Levine [10,11], are similar to the multi-set formalism of the DD-vMCG method: the nuclear wave packet is described by a set of GBF, each belonging to a given electronic state. The main difference is that it uses classical mechanics to generate the GBF trajectories and thus, only the evolution of the wave packet expansion coefficients is determined by solving the time-dependent Schrödinger equation. This implies that certain quantum mechanical phenomena (i.e. non-adiabatic effects and tunnelling) may not be well-described unless a large number of trajectories are computed. However, another important feature of the AIMS method is that the basis set is expanded adaptively in regions of strong non-adiabatic coupling where the wave function bifurcates, allowing an efficient description of state crossing processes. This feature has been extended to allow the description of tunnelling effects as well [12].
-The surface hopping method, developed by Tully [13,14], is also a classicalbased direct dynamics method but where the nuclear wave packet is represented by a swarm of independent trajectories, each evolving on a single adiabatic potential energy surface at any given time and able to "hop" from one surface to another.
-The Ehrenfest method [15][16][17][18][19] is the equivalent of the surface hopping method in the single-set formalism: each independent classical trajectory follows the gradient of a superposition of electronic adiabatic states, i.e. evolves on an effective potential energy surface.
-The coupled-coherent states (CCS) and multi-configurational Ehrenfest (MCE) methods of Shalashilin and co-workers [20][21][22][23] build on Ehrenfest trajectories and couple them together. Note that there exists a variant of the MCE method where the basis set of Ehrenfest trajectories is expanded adaptively as in the AIMS method; this method is called ab initio multiple cloning (AIMC) [24]. They differ from the DD-vMCG single-set method since the basis trajectories obey classical and not quantum mechanics.
A new challenge for on-the-fly methods able to describe non-adiabatic dynamics has risen with the emergent field of "attosecond measurement". (It now becomes possible to experimentally study of the dynamics of electrons and nuclei in molecules, on the attosecond (1 as = 10 −18 s) time scale.) The synthesis of attosecond pulses relies on the time-energy uncertainty principle ∆E∆t ≥h: collecting together coherent light sources with an energy bandwidth ∆E of several eV gives a pulse with duration τ ≈h/∆E of a few attoseconds. The broad spectral bandwidth of such short pulses leads to the (initial) coherent population of several electronic states, thus breaking the Born-Oppenheimer approximation. The system is not confined any more to be in a single stationary electronic state but is now a non-stationary superposition of electronic states, called an electronic wave packet. With the population of a superposition of electronic states, the convenient picture of a single potential energy surface is lost: the nuclei "feel" the multiple coupled potential energy surfaces. Such coupled electron-nuclear dynamics is often called charge-directed reactivity in the literature [25,26].
We will argue that the "single-set" class of methods are more natural for the simulation of such non-adiabatic dynamics since each nuclear trajectory "feels" all populated electronic states. We choose to present in more detail the simplest approach -the Ehrenfest method -and the most accurate one -the DD-vMCG method. In Section 2, we derive the Dirac-Frenkel variational principle which is then applied to derive the equations of motion of the Ehrenfest method in Section 3 and of the vMCG method in Section 4. In Section 5, we suggest an alternative derivation of the vMCG method which gives equivalent equations of motion that result into the Ehrenfest equations of motion when taking the one-configuration and classical limits.

Dirac-Frenkel variational principle
A way to obtain approximate solutions to the time-dependent Schrödinger equation (1), in methods such as the Ehrenfest and DD-vMCG methods, is to apply the variational principle to a guess form of the molecular wave function Φ(t). The resulting equations of motion for the time-dependent parameters of the wave function assure that the evolution of these parameters optimally represents the true evolution of the wave function.
In the so-called Lagrangian form, the time-dependent variational principle of quantum mechanics [27] is formulated from the action-type integral [28] for the time evolution of a wave function in a given time interval [t 0 , t 1 ]. The integrand in the action integral equation (2) can be interpreted as the expectation value of the deviation with respect to the exact time-dependent Schrödinger equation (1). The idea is that the desired wave function Φ(t) should minimise this deviation within the given time interval. The aim is thus to determine an "optimal" wave function Φ(t) such that the action integral is stationary with respect to (small) variations of the form Φ(t) + δΦ(t). Here the variations are required to vanish at the boundaries of the time interval: δΦ(t 0 ) = δΦ(t 1 ) = 0. The variational principle then reads A more convenient form of the time-dependent variational principle may be obtained by rewriting δS[Φ] in such a way that the variational principle applies directly to the integrand in the time-integration rather than to the action integral. Starting from one can replace the term involving δΦ(t) using integration by parts to arrive at Taking into account that the integration boundaries t 0 and t 1 are arbitrary, one may conclude that the integral (5) vanishes if and only if [29] Re δΦ(t)|ih ∂ ∂t −Ĥ|Φ(t) = 0 This establishes an instantaneous form of a time-dependent variational principle. If together with the variations δΦ(t) the variations iδΦ(t) are also permitted, one arrives at the related form Both forms are contained in the so-called Dirac-Frenkel variational principle (DFVP) [30,31] and can be written as where δΦ denotes all possible variations of Φ with respect to the parameters.
In the next sections, we apply the DFVP to derive the Ehrenfest and DD-vMCG methods.

Wave function ansatz
The simplest possible guess form of the total molecular wave function is a product ansatz, which separates the electronic and nuclear variables: where Ψ (r, t) and χ(R, t) are the electronic and nuclear wave functions respectively. This is the first approximation made in the Ehrenfest method. Note that this approximation is called a single-configuration ansatz for the total wave function. It is mentioned in passing that this product ansatz (9) is different from the Born-Oppenheimer ansatz [32] for separating the electronic and nuclear variables even in its one-determinant limit, where only a single electronic eigenstate ψ s is included in the expansion. 2 One deficiency of the ansatz (9) is the fact that the wave function does not have the possibility to decohere: the populated electronic states in Ψ (r, t) must share the same nuclear wave packet χ(R, t) by definition of the total wave function. 3 The neglect of electronic decoherence could lead to non-physical asymptotic behaviours in the case of bifurcating paths.
In order to simplify the appearance of the expressions at a later stage of the derivation [36], a phase factor is introduced for the total wave function and also some internal phase factors for the two individual wave functions with E(t) = χΨ |Ĥ e |χΨ R,r and E tot = χΨ |Ĥ|χΨ R,r . The indices R and r indicate the coordinates of integration.

Equations of motion
The equations of motion of the nuclear and electronic parts are obtained by applying the DFVP (8) to the ansatz (11). Here,

Equation of motion for the electronic part
Let us first apply the DFVP to the variation of the electronic part δΨ : and ih δΦ| ∂Φ ∂t We obtain: More explicitly, it reads: The index i refers to the electrons; m e is used to denote the mass of an electron.

Equation of motion for the nuclear part
Let us now apply the DFVP to the variation of the nuclear part δχ: and ih δΦ| ∂Φ ∂t We obtain: More explicitly, it reads: The index I refers to the nuclei; M I is used to denote the mass of the nucleus I.
The set of coupled equations (18) and (22) are the basis of the timedependent self-consistent field (TD-SCF) method [30,37], also referred to as time-dependent Hartree (TDH) when the nuclear wave function χ(R, t) is written as a simple product of (time-dependent) one-dimensional functions. 4 The mean-field origin of the TD-SCF approach imposes limitations, as discussed already above. To understand the consequence of using the ansatz (9), let us for instance look closer at the second term on the right hand side of equation (18). The interaction between electrons at points r in space and nuclei at points R is weighted by the probability that the nuclei are at these particular points R. This is the effective potential experienced by the electrons due to the nuclei. The corresponding remark can be made about the second term on the right hand side of equation (22).
According to the set of coupled equations (18) and (22), the feedback between electronic and nuclear degrees of freedom is described in a mean-field manner, in both directions. In other words, both electrons and nuclei move in time-dependent effective potentials obtained from appropriate expectation values of the nuclear and electronic wave functions respectively. More exactly, the nuclear wave packet on a particular potential energy surface would evolve on its electronic state under coupling with the other electronic states. When simulating the nuclear motion induced by a superposition of electronic states (i.e. several electronic states being initially populated), we expect the Ehrenfest method to be valid at short times before the nuclear wave packets belonging to different electronic states move too far apart from each other.

Classical limit for nuclear motion
The Ehrenfest method is the classical analogue to the TD-SCF method [38] and therefore inherits the same limitation. It is obtained by taking the classical limit of equations (18) and (22). To do that in equation (22), the nuclear wave function is (exactly) rewritten as allowing for a complex phase S [39]. After inserting ansatz (23) in (22), we obtain: Equation (24) is entirely equivalent to the original equation (22). The righthand side term (proportional toh) may be thought of as a time-dependent "quantum correction". The classical Hamilton-Jacobi equation is obtained when taking the limith → 0: 5

∂S ∂t
The resulting equation (25) is thus equivalent to Newton's equation of motion, where P I = − → ∇ I S is the classical momentum of nucleus I: In equation (18), we can replace χ(R, t) by a delta function at the classical trajectory R(t): =Ĥ e (r; R(t)) · Ψ (r, t; R) 5 In many textbooks, the nuclear wave function is rather rewritten in a polar coordinate system in terms of an amplitude A and a phase S which are both considered to be real [40]: χ(R, t) = A(R, t) · exp ī h S(R, t) . This ansatz leads to two equations for the amplitude and phase functions, often called the hydrodynamic formulation. The equation of motion for the phase function is called the quantum Hamilton-Jacobi equation. It has an extra term relative to the classical Hamilton-Jacobi equation which involves the second derivative of the amplitude (not that of the phase) and may be thought of as a time-dependent quantum potential. This is the Bohm interpretation of quantum dynamics. It is tempting to think of the quantum Hamilton-Jacobi equation as a natural way to pass to the classical limit. However, in neglecting the quantum potential, ∇ 2 I A must be small, which is a condition that the wave packet be broad [39]! The fact that narrow wave packets would signify a breakdown in classical-quantum correspondence must be an artefact.
Note that now the electronic wave function Ψ depends parametrically on R(t) through V n−e (r, R(t)) and thusĤ e (r; R(t)). By treating the nuclear motion classically, we lose the spatial delocalisation of the nuclei and their motion is now described by a classical trajectory. To obtain a realistic description of the dynamics of the system, one mimics the initial nuclear wave packet distribution by propagating a swarm of independent trajectories starting with sampled classical positions R and momentum P of the nuclei.
Equations (26) and (27) define the Ehrenfest method. 6 They allow transfer of energy between quantum and classical degrees of freedom such that the total energy is conserved [41]. The a priori construction of the potential energy surfaces is avoided from the outset by solving numerically the coupled set of equations simultaneously "on-the-fly" for each nuclear geometry R(t) generated along the trajectory.

DD-vMCG method
The DD-vMCG method comes from the multi-configurational time-dependent Hartree (MCTDH) method. The latter was originally devised as an efficient grid-based quantum dynamics solution by using the variational principle to derive the equations of motion for a flexible wave function ansatz. The resulting evolution of the time-dependent basis functions means that the basis set remains optimally small, i.e. the molecular wave function is very compact.
To remove the restrictions of the grid, the G-MCTDH method was introduced [5], where some of the multi-dimensional basis functions are replaced by parametrised Gaussian functions. If only Gaussian basis functions are included in the G-MCTDH method, one naturally arrives at the vMCG method. Its direct dynamics implementation is known as DD-vMCG [9].

Wave function ansatz
The vMCG method uses a set of multidimensional Gaussian basis functions (GBF) to represent the guess form of the total molecular wave function. The ansatz spreads over a set of coupled electronic states {φ s } and reads Here, we used the single-set formalism as described before, i.e. the GBF χ j is the same for all electronic states. A consequence is that a particular GBF χ j is constrained to move identically on every electronic state considered.
One can use the Heller expression for a multidimensional separable GBF [42]: a real-valued Gaussian function (spatial amplitude envelope) multiplied by a Fourier function (plane wave giving a group velocity), where along each degree of freedom α (nuclear coordinate R α ): σ α is the width (spatial standard deviation) kept fixed during the simulation for numerical stability and taken as equal for all GBF, R jα and p jα are the mean position and momentum, defining the trajectory followed by the centre of the function in the phase space. Note that here, we set the real part of the complex factor to zero and we use the imaginary part to keep the GBF normalised.

Equations of motion
The equations of motion of the expansion coefficients and the GBF parameters are obtained by applying the DFVP (8). Here, with {λ jα } the parameters of the GBF χ j , i.e. the mean position and momentum {R jα , p jα }.

Equations of motion for the expansion coefficients
Let us first apply the DFVP to the variation of the coefficients δA (s) and withĤ (ss ) = φ s |Ĥ|φ s , H (ss ) jl = χ j |Ĥ (ss ) |χ l , S jl = χ j |χ l and τ jl = χ j |χ l the elements of the Hamiltonian matrix, overlap matrix and timederivative matrix. We obtain: In matrix notation, it reads: or ihȦ Equation (37) differs from the standard MCTDH equation only by the term including the time-derivative matrix, which accounts for the non-orthogonality of the Gaussian basis set.

Equations of motion for the GBF parameters
Let us now apply the DFVP to the variation of the GBF parameters δλ jα : and with ρ Using equation (37) to expandȦ (s) l in the second term in the right hand side gives (with matrix notation): Let us define the vector Λ collecting all the time-dependent parameters defining each GBF λ jα and the matrices C and Y as follows: We then obtain: ih lβ C jα,lβλlβ = Y jα (44) Equations (37) and (46) are the equations of motion of the vMCG method. Again, the a priori construction of the potential energy surfaces is avoided from the outset by solving numerically the coupled set of equations simultaneously "on-the-fly" in the DD-vMCG method. Note that the GBF are coupled both directly and through the expansion coefficients. The vMCG method in principle allows the nuclear wave packets on the different electronic states to move apart from each other. This gives a better description of the wave function in the case of bifurcating paths and of electronic decoherence.

Link between Ehrenfest and DD-vMCG methods
In this Section, we attempt to explicitly show the link between the Ehrenfest and DD-vMCG methods.

Alternative wave function ansatz for DD-vMCG
The relation between the guess forms for the total molecular wave function used in the derivation of the Ehrenfest and DD-vMCG methods is better highlighted if one rewrites the DD-vMCG ansatz (28) as where the summation over the electronic states is now implicit: It is important to note that here the electronic wave functions are normalised but not orthogonal: Ψ j |Ψ l r = 0.
The DD-vMCG ansatz now appears as the sum of several products of electronic and nuclear parts, i.e. configurations, B j (t) being the time-dependent expansion coefficients of these configurations. It is clear that the Ehrenfest ansatz (9) is obtained by taking the one-configuration limit, i.e. keeping only one term of the sum.
As in the Ehrenfest derivation, we introduce a phase factor to the total wave function and some internal phase factors for the two individual wave functions with E j (t) = χ j Ψ j |Ĥ e |χ j Ψ j R,r = H jj (t) and E tot j = χ j Ψ j |Ĥ|χ j Ψ j R,r = H tot jj .

Alternative equations of motion for DD-vMCG
Now let us apply the DFVP (8) to the new (equivalent) DD-vMCG ansatz (49).

Equation of motion for the expansion coefficients B j
Let us first apply the DFVP to the variation of the expansion coefficients δB j : and with H tot jl = χ j Ψ j |Ĥ|χ l Ψ l R,r , S el jl = Ψ j |Ψ l r and τ el jl = Ψ j |Ψ l r . Thus, we obtain ih l S jl S el jlḂl = l H tot jl − ih(τ jl S el jl + S jl τ el jl ) + E l (t)S jl S el jl B l (55)

Equation of motion for the electronic part Ψ j
Let us now apply the DFVP to the variation of the electronic part δΨ j : and Note that so far, no approximation has been made (except for the molecular wave function guess form). These "new" equations of motion are thus equivalent to equations (37) and (46).

One-configuration and classical limits
When taking the one-configuration limit, the wave function ansatz becomes: It corresponds to the DD-vMCG ansatz when only one GBF is used in the single-set formalism. The nuclear wave packet is thus constrained to be and remain Gaussian.
Let us now take the one-configuration limit in the equations of motion. The equation of motion for the expansion coefficients (55) becomes: The equation of motion for the nuclear part (61) becomes: As expected, equations (64) and (65) are identical to the equations (18) and (22) defining the TDH method, with the additional constraint that the nuclear wave packet has a Gaussian shape. The Ehrenfest equations (26) and (27) are obtained by taking the classical limit as shown above in Section 3.3.
The Ehrenfest method is thus obtained by taking: (i) the one-configuration limit of the vMCG method constraining the nuclear wave packet to remain Gaussian and be the same for all electronic states; and (ii) the classical limit for the nuclei making the Gaussian infinitely narrow.
We note in passing that, in the same spirit as CCS and MCE methods of Shalashilin and co-workers [21], one could post-process independent Ehrenfest trajectories and form a nuclear wave packet by putting floating GBF on top of the trajectories. Equation (55) would govern the time evolution of the amplitudes of the different trajectories.

Conclusion
In this article, we have presented an overview of the current state-of-the-art "on-the-fly" methods for non-adiabatic dynamics, introducing the similarities and differences between them (Figure 1). We have also derived two of the methods: the Ehrenfest method and the DD-vMCG method. In summary, both methods treat the electronic degrees of freedom quantum mechanically. To describe the nuclear degrees of freedom, the DD-vMCG method uses a basis of coupled quantum trajectories while the Ehrenfest method uses a swarm of independent classical trajectories.
We have presented derivations of the two methods from the DFVP with unified notations and we have proposed an alternative derivation for the DD-vMCG method, which gives the Ehrenfest equations of motion when taking the appropriate limits. By doing so, we mention the possibility of a hybrid method in which one could post-process independent Ehrenfest trajectories to form a nuclear wave packet with coupled amplitudes of the different trajectories (in the same spirit as CCS and MCE methods).