Gaussians for Electronic and Rovibrational Quantum Dynamics

The assumptions underpinning the adiabatic Born–Oppenheimer (BO) approximation are broken for molecules interacting with attosecond laser pulses, which generate complicated coupled electronic-nuclear wave packets that generally will have components of electronic and dissociation continua as well as bound-state contributions. The conceptually most straightforward way to overcome this challenge is to treat the electronic and nuclear degrees of freedom on equal quantum-mechanical footing by not invoking the BO approximation at all. Explicitly correlated Gaussian (ECG) basis functions have proved successful for non-BO calculations of stationary molecular states and energies, reproducing rovibrational absorption spectra with very high accuracy. In this Article, we present a proof-of-principle study of the ability of fully flexible ECGs (FFECGs) to capture the intricate electronic and rovibrational dynamics generated by short, high-intensity laser pulses. By fitting linear combinations of FFECGs to accurate wave function histories obtained on a large real-space grid for a regularized 2D model of the hydrogen atom and for the 2D Morse potential, we demonstrate that FFECGs provide a very compact description of laser-driven electronic and rovibrational dynamics.


Introduction
With the advent of new technology for manipulating atoms and molecules with intense ultrashort (atto-and femtosecond) laser pulses, there is an urgent need for further development of accurate and reliable quantum-dynamics (QD) tools for simulations of events involved in such manipulations.Such simulations are needed both to guide the design of experiments and to ensure correct interpretation of observations.Atomic and molecular QD simulations involving interaction of these systems with ultrashort intense laser pulses can be carried out by integrating the time-dependent Schrödinger equation (TDSE) on a real-space grid or by using an expansion of the wave function of the system in terms of basis functions whose lin-ear and non-linear parameters are adjusted along with the number of basis functions during the propagation.The focus of this work is an analysis of the necessary features such basis functions must possess to be effective for simulating laser-induced dynamics of an atomic or a molecular system.
The dynamics induced by the interaction of intense electromagnetic radiation with nuclei and electrons of a molecule necessitate an accurate account of the coupling of these two types of particles.The broad intensity distribution in the frequency domain of ultrashort laser pulses implies that a large number of bound and continuum states are involved in the dynamics, including rovibrational, electronic, and collective states where the motion of both nuclei and electrons are simultaneously excited.To describe such an intricate situation, the separation of the nuclear and electronic motions-a hallmark of the adiabatic Born-Oppenheimer (BO) approximation 1,2 -must not be assumed and, ideally, all particles forming the system should be treated on an equal footing.Such an approach is tested in the present work using two simple 2D models.Also, an extension of the approach to simulate the laser-induced dynamics of attosecond atomic and molecular events involved in attosecond experiments [3][4][5] is discussed.
The coupled nuclear-electronic motion is highly correlated, as the electrons, particularly the core electrons, generally follow the nuclei very closely and the nuclei stay apart from each other due to their strong Coulomb repulsion and large masses.The situation is markedly different for the electrons whose wave functions, due to much lower masses, more significantly overlap.To best describe these effects using a basis-set approach, one needs to expand the wave function in terms of functions that explicitly depend on nucleus-nucleus, nucleuselectron, and electron-electron distances, i.e., the explicitly-correlated functions (ECFs).In the first part of this work, we review the ECFs used in the atomic and molecular calculations of stationary bound states and we discuss the features of these functions that may make them particularly useful in QM calculations of atomic and molecular systems.][8][9][10][11][12][13][14][15][16][17][18][19][20] In the second part, two-dimensional time-propagation calculations are performed for two model systems involving Coulomb and Morse potentials using a grid approach.Next, the time-dependent grid wave functions are fitted with ECGs that are chosen to best represent the key features that appear in the wave function due to the interaction of the system with ultra-short intense laser pulses.
3][24] In electronic-structure theory the Gaussians are real-valued functions centered at the atomic nuclei and contracted to form atomic orbitals which, in turn, form a non-orthogonal basis for the expansion of molecular orbitalssee, e.g., Ref. 21 for a detailed account.Such Gaussians have also been used for the study of many-electron dynamics, [25][26][27] although important highly nonlinear phenomena such as ionization processes and high harmonic generation cannot be properly accounted for.The latter can to a certain extent be ameliorated by augmenting the standard basis with Gaussians fitted to continuum (Bessel, Coulomb, or Slater-type) functions; see the recent review by Coccia and Luppi and references therein for more details. 28e core idea of Heller's approach [22][23][24] to vibrational dynamics is to use complex-valued Gaussians (Gaussian wave packets), which are exact solutions for harmonic potentials.Unlike in electronic-structure theory, the Gaussian parameters are now time-dependent variational parameters.For anharmonic potentials, however, the equations of motion for the Gaussian parameters quickly become ill-conditioned and one resorts to locally harmonic approximations of the potential, which is a reasonable approach as long as the wave packet is sufficiently localized-see, e.g., Ref. 29 for a recent review of the so-called thawed Gaussian approach.
1][32][33] In all cases, however, the ill-conditioned equations of motion is a serious obstacle.In this work, we investigate the ability of complexvalued Gaussians to represent complicated dynamics of electrons and nuclei by fitting to accurate grid-based solutions of the time-dependent Schrödinger equation.

Non-BO Hamiltonian
The total non-relativistic all-particle non-BO molecular Hamiltonian describing the interaction of a neutral molecule with a uniform, time-dependent electric field defining the x-axis of a laboratory-fixed coordinate frame can be rigorously separated into a center-of-mass (COM) kinetic-energy operator and the internal Hamiltonian, 6,8,34 Ĥ(t).The separation is accomplished by transforming the total Hamiltonian from Cartesian laboratory coordinates, R i , i = 1, . . ., N (N is the total number of particles in the molecule) to a new Cartesian coordinate system, parallel to the laboratory frame, where the first three coordinates are the coordinates of the COM and the remaining coordinates are internal coordinates.The origin of the internal frame is chosen at a reference particle, typically the heaviest nucleus, which is taken to be particle number 1 such that r The internal Hamiltonian then takes the following form (using atomic units throughout): 6 where F(t) is the time-dependent electric-field strength, M 1 is the mass of the reference particle (particle 1), |r i |, and the prime denotes vector/matrix transposition.The Hamiltonian (1) represents n interacting particles with masses equal to the reduced masses moving in the central Coulomb potential of the reference particle.We refer to these particles as pseudo-particles because, while they have the same charges as the original particles, their masses are different.For F(t) = 0 the Hamiltonian (1) is fully symmetric (isotropic or atom-like) with respect to all rotations around the center of the internal coordinate system and its eigenfunctions transform as irreducible representations of the fully symmetric group of rotations.When F(t) ̸ = 0, however, the symmetry is reduced to cylindrical about the field direction, here the x-axis.
For a diatomic system, after separation of the center of mass motion, the internal Hamil-tonian used in the non-BO time-evolution calculations represents a motion of the second nucleus and the electrons (with their masses replaced by the respective reduced masses) in the field of the charge of the first nucleus (the reference nucleus) located in the center of the internal coordinate system.The potential acting on the second nucleus that results from the interaction of this nucleus with the charge of the reference nucleus and the electrons can effectively be represented by a Morse-like potential.An important effect that also determines the electronic-nuclear dynamics of the system is the electrostatic attractive interaction of each of the electrons with the reference nucleus located at the center of the coordinate system.Thus, the Coulombic and Morse interactions investigated in this work are central to understanding the molecular dynamics.The interactions which are present in the internal Hamiltonian, but are not investigated in this work, are two particle interactions involving the second nucleus and the electrons, and the inter-electron interactions.
To represent these types of interactions, models involving more than two dimensions would be needed and, thus, they are not investigated in this work.However, based on our ECG stationary-state calculations, we expect that FFECGs should perform very well in describing these interactions.

ECGs used in very accurate non-BO calculations of stationary states of small atoms and molecules
To achieve high accuracy needed in the computations we will use an approach that is both adaptive in space and time.Mesh-free complex explicitly-correlated Gaussian (CECG) functions, that are free to warp and roam in space, will be the main tool.The Adamowicz group has used various types of ECGs and CECGs in very accurate non-BO atomic and molecular calculations of stationary bound states for over two decades.7][8][9] The simplest ECG with real non-linear parameters used to calculate an S state of an n-electron atom (S-ECG) is: where r is vector of 3n internal Cartesian coordinates of the electrons and A is a 3n × 3n real symmetric positive-definite matrix of the non-linear parameters.A has the following block structure: A = A ⊗ I 3 , where A is a n × n real dense symmetric positive-definite matrix and I 3 is a 3 × 3 identity matrix, while symbol ⊗ denotes the Kronecker product.
Such representation of matrix A ensures that the exponential part of the basis function is invariant with respect to 3D rotations.
S-ECGs can alternatively be represented as: where the first factor is a product of n orbitals and the second factor is a product of n(n+1)/2 pair functions explicitly dependent on the squares of all inter-electron distances, r 2 ij .The methods allow for very accurate calculations of the spectra of small atoms and molecules when the leading relativistic and QED corrections are also included in the calculations.
The non-BO ECG calculations for S, P , D, and F states of atomic systems with 2-5 electrons [35][36][37][38][39][40][41][42][43][44][45][46][47] are among the most accurate in the literature.As the ECGs explicitly depend on the distances between the particles (electrons and nuclei), they very efficiently represent the coupled nucleus-electron motions and allow very accurate accounting of the inter-particle correlation effects.These effects are indispensable in non-BO calculations, because, as mentioned, the Coulomb interactions make particles with alike charges avoid each other and particles with opposite charges to follow each other.This effect can also be very effectively described with the ECGs.
A challenge in non-BO ECG calculations of stationary ground and excited states of small atoms and molecules is to accurately describe radial and angular oscillations of the non-BO wave functions of highly excited states.Three types of ECGs have been used to describe these features.These are: (a) molecular ECGs with pre-exponential multipliers in the form of powers of the internuclear distances.The functions are called "power" ECGs (PECGs) and they have the following form: where m i and m ij are even non-negative integers, and A, as defined before, is a real symmetric positive-definite 3n × 3n matrix.PECGs have been used in molecular non-BO calculations for small molecules; 44,[48][49][50][51][52][53] (b) complex single-center ECGs.The works that are particularly relevant to this project concern implementation of algorithms for performing very accurate calculations on bound states of small molecules that employ complex ECGs (CECGs): where A, as defined before, is a real symmetric matrix and B is also a real symmetric matrix with the same structure as (c) real ECGs with shifted centers (SECGs) of the form: where q is 3n real vector of the Gaussian shifts and A, as defined before, is a real symmetric matrix of the non-linear parameters.9][60][61][62][63][64][65] Including real shifts in the Gaussians enables to describe radial and angular polarization of the wave function.These types of deformations can be also described by linear combinations of spherical-harmonics factors, though the shifts may be a more effective way for the task.
The purpose of this work is to develop and test an ECG basis set to be employed in quantum-dynamics simulations of molecular systems exposed to an ultrashort laser pulse within the semiclassical electric-dipole approximation.The proposed basis is tested by fitting wave packets obtained as solutions of the time-dependent Schrödinger equation with a grid-based method for two two-dimensional (2D) model systems.The models represent two main features that need to be described in a QD simulation of a molecule, i.e., the electrostatic interaction represented by a Coulombic potential and the rovibrational interaction represented for a diatomic molecule by a Morse potential.In the next section, before the ECG basis functions for QD molecular simulations are introduced, the grid-based calculations of the two models are described and discussed.

2D Model calculations using a grid approach
Our ultimate goal is to solve the TDSE, i ψ(t) − Ĥ(t)ψ(t) = 0 for the non-BO internal Hamiltonian (1) representing a molecule interacting with a short intense laser pulse.For a diatomic molecule, there are two major interactions that need to be described.The first is the repulsive interaction of the pseudo-nucleus with the charge of the reference nucleus located in the center of the coordinate system, and the second is the attractive Coulomb interaction of a pseudo-electron with the charge of the reference nucleus.Due to the screening effect of the former interaction by the electrons, the interaction potential is not Coulombic but is more appropriately represented by a Morse-type potential.Thus, at the very minimum, in selecting an ECG basis set for solving the non-BO TDSE, one should verify if the chosen basis is capable of describing the laser-induced dynamics of a single particle interacting with the charge of the reference particle with the Morse and Coulomb potentials.For the verification we use an elementary 2D model Hamiltonian of the form: For the electron we use the soft Coulomb potential, V (x, y) = −(x 2 + y 2 + 1/2) −1/2 , which mimics the nuclear potential of a hydrogen-like atom, and the Morse potential is given by 2 , with D e = 0.17449, r e = 1.4011, and α = 1.4556.The charge and (reduced) mass are set to q = −1, µ = 1 for the Coulomb model, and q = 1, µ = 1605.587for the Morse model.The electric-field strength is nonzero only in the time interval t 0 < t < t 1 , where it is equal to: where E 0 denotes the maximum electric field amplitude.In our calculations for both models we set t 0 = 0.For the Coulomb model we set ω = 0.25 a.u., t 1 = 60 a.u. and E 0 = 0.4 a.u., which corresponds to a laser pulse of wavelength ≈ 182 nm, consisting of 2.5 optical cycle.
For the Morse model we set ω = 0.0 a.u., t 1 = 20 a.u. and E 0 = 2.0 a.u., corresponding to a short, delta-like pulse (relative to the time-scale of the nuclear motion).Our laser pulse parameters are chosen not for their significance in relation to any experiment, but rather such that it generates complicated ionization and dissociation dynamics.We consider the around the center of the internal coordinate system in the xy plane.The ground-state wave function is obtained using inverse iterations using the conjugate-gradient method for the solution of large sparse linear systems.In the case of the Coulomb potential the groundstate wave function approximates a 2D 1s orbital and, in the case of the Morse potential, the ground-state wave function has a "torus" shape and is practically zero at the coordinate center, peaks at r e , and then again goes to zero at larger distances.In Fig. 2 and Fig. 3, some snapshots of the time evolution of the wave functions are shown for the Coulomb and Morse simulations, respectively.For each case, the real and imaginary parts of the wave function, as well as the wave-function absolute value, are plotted.In both cases the respective wave functions become increasingly more complicated with many features, more deformed and oscillatory, and more diffused.
It is interesting to know to what extent the time-evolving wave packets for the two considered models involve contributions from higher angular momenta.For both systems, at t = 0, the initial wave packets, i.e. the corresponding wave functions, are fully rotationally symmetric, i.e., for both Lz ψ(0) = 0, where L z = −i ∂ ∂θ is the angular momentum operator in 2d and θ is the angular coordinate in planar polar coordinate system.Moreover, [ Ĥ(t), Lz ] = 0 before and after the pulse begins and after it ends.The non-vanishing of the commutator for t 0 < t < t 1 implies that ψ(t) is not, in general, an eigenfunction of Lz .
In Figure 4, the angular momentum probability distributions are shown for both model systems.Clearly, the pulse induces high angular momenta in the wave function.
Fitting of the wave functions obtained in the Coulomb and Morse simulations with ECGs is described in the next section.Based on the analysing of the simulation results, it is clear that Gaussians used in the fitting need to be capable of describing the ground-state wave function, the cylindrical deformation and oscillation of the function due to the interaction with the field, the diffusion of the function associated with possible ionization or/and dissociation of the system, and the coupling of the motion of all particles forming the system including the non-BO coupling of the motions of the electrons and the nuclei.Most of these features appear in the calculations of molecular static ground and excited states and an excellent performance of the ECGs have been well documented in those calculations.The features, that do not appear in the static calculations are ionization and dissociation.It necessitates that the wave function is allowed to acquire some plane-wave character.This can   be achieved by allowing the shift vector and parameter matrix A in Gaussian (5) to become complex.Such Gaussians, named by us the fully flexible explicitly correlated Gaussians (FFECGs), have the following form: or alternatively the following form: where A and B, as defined before, are real symmetric 3n × 3n matrices, and p and q are 3n vectors.The FFECGs are fully flexible complex multi-particle Gaussians that can provide a basis set for expanding a time-evolving non-BO wave packet of a molecular system interacting with a short fast-varying intense laser pulse.As shown in the next section, a linear combination of FFECGs can be used to very accurately represent the ground-state wave functions of the two models considered in this work.It can also describe the timedependent oscillations of the time-evolving wave packet.Also, due to making the Gaussian shift vectors to be complex, the ionization and dissociation processes can be described.And finally, FFECGs can very effectively represent the coupled and highly correlated motions of the electrons and nuclei forming the system.They allow the electrons and the nuclei to be treated on an equal footing in the calculation.

Fitting the grid-based wave packet with FFECGs
The 2D wave packet obtained in the time-dependent grid calculation at each particular time step for each of the two models is fitted with a linear combination of FFECGs.In this case FFECG has the following form: where p = (p x , p y ) ′ and q = (q x , q y ) ′ .Thus, in the 2D case ϕ(r), is dependent on six real non-linear parameters, a, b, p x , p y , q x , and q y .As the time-dependent calculations for both models are initiated with the corresponding ground-state wave functions, which are real and spherically symmetric, appropriate FFECGs need to be used.For the hydrogen model, these FFECGs are simple spherical Gaussians centered at x = y = 0 with a ̸ = 0 and b = p x = p y = q x = q y = 0.A linear combination of FFECGs of this kind can fit the ground-state wave function obtained in the grid calculation with very good accuracy.For the Morse model, one needs to only use FFECGs with a ̸ = 0, b ̸ = 0, and p x = p y = q x = q y = 0 to generate a good fit to the ground-state wave function obtained in the grid calculation.
Additionally, for each FFECG pair used, the two Gaussians need to have the same a, but their b's should be b and −b, so they can be "contracted" for form the following real function: where K is a constant.A linear combination of such contracted FFECGs provides a very accurate fit to the ground-state wave function for the Morse model.
The fitting of the grid ground-state wave functions with FFECGs are done using the standard least-squares implementation from the SciPy Python library. 67The method is also used to fit the wave functions obtained in consecutive time steps obtained in the grid simulation.First, the number of FFECGs used is the same as for the ground state, but the parameters frozen at zero for the ground-state wave function are now unfrozen and optimized.For both models we set the threshold for the least squares cost function equal to 10 −5 .When the assumed accuracy of the fitting cannot be achieved, additional FFECGs are added to the basis set.Each addition includes a group of FEECGs with parameters: The fitting that involves optimization of all linear and nonlinear parameters of the enlarged FFECG basis set continues for some number of the following time steps until it is determined that the fitting process is no longer successful.At that point new FFECGs are added to the basis set and the wave packet is refitted.The fitting for the Coulomb and Morse models is shown in Figs. 5 and 6.In the figures, the wave packets obtained in the grid calculations are compared with the corresponding FFECG fits for some selected time points.
Also, in Fig. 7 the changes of the cost function representing the accuracy of the fitting, as well as the time-resolved size of the basis set, are plotted for the Coulomb and Morse models.Both Figs. 5 and 6, and Fig. 7 show that the fitting accuracy is satisfactory for both models despite the wave function becoming increasingly more complicated.Naturally, as this happens, the number of FFECGs has to be constantly increased.
In general, the Morse model seems to be somewhat more difficult to represent using the variational Rayleigh-Ritz functional or Rothe variational functional) usually are more accurate in some spatial domains than the others.It is difficult to a priori determine which imperfections of the representation of the wave function will be amplified and which will be suppressed when a particular observable is calculated.
To better elucidate the accuracy of the FFECG fits obtained in this work, we compare time-resolved observables obtained for the grid wave packet and the fitted wave packet.The calculated observables are the expectation value of the field-free Hamiltonian, the ground-   We note that the ground-state survival probability of the final state is small or zero for both models, and that the angular-momentum expectation value indicates involvement of highly excited rotational states.The final energy for the Coulomb model is positive, indicating an unbound (i.e., ionized) state.For both models, the increase of the dipole moment during the dynamics indicate the large, asymmetric spreading of the wave packet.
It is quite remarkable that so relatively few FFECGs are required to accurately reproduce such complicated dynamics.

Conclusion
The grid approach is used to obtain solution of the time- Finally, this work represents a preliminary step in the application of FFECGs to describe the coupled electronic-nuclear dynamics in atomic and molecular systems.In the future work involving FFECGs and the non-BO nuclear-electronic quantum dynamics, the Rothe method 68-70 will be employed to propagate the wave packet.The approach, also known as the adaptive method of time layers, 71 relies on reformulating the time-dependent variational principle into a series of minimizations of the Rothe functional at consecutive time steps.This is an alternative to the standard real-time propagation techniques based on the Dirac-Frenkel variational principle and propagated using, e.g., Runge-Kutta methods.The ECG optimization protocol developed in the present work to fit a linear combination of ECG to the grid-based wave packet will be applied to minimize the Rothe functional with respect to the linear and non-linear parameters of the ECGs.In our works on the variational calculations of molecular stationary states with real and complex ECGs we have developed procedures for the variational minimization of the energy functional which employs the analytical energy gradient determined with respect to the ECG nonlinear parameters.The use of the gradient has significantly expedited the functional minimization and enabled to obtain non-BO energies and the corresponding wave functions whose accuracy by far exceeds the results obtained by others.The gradient-based approach will also be used in the optimization of the ECGs parameters carried out through the minimization of the Rothe functional.The high efficiency of the computer code for the optimization of the ECGs in the stationary-state calculations has been also achieved by deriving the algorithms for calculating the necessary N-particle matrix elements (i.e. the overlap, Hamiltonian, and gradient matrix elements) using the matrix differential calculus and by coding them using highly parallel, vectorized, and GPU-enabled strategies.These strategies will be used in the Rothe time-propagation calculations.
where b is a real dense symmetric n × n matrix).It was shown that CECGs can very effectively describe high-frequency radial oscillations of the wave function of highly vibrationally excited states.The angular oscillations can be described by adding Cartesian spherical harmonics as pre-exponential multipliers to the Gaussians.CECGs have been used in non-BO calculations of molecular Σ, Π, and ∆ bound rovibrational states[54][55][56][57] It has been shown that CECGs are equally, if not more, efficient as PECGs in describing radially and angularly oscillating wave functions of states located near the dissociation threshold;

Figure 1 :
Figure 1: Shape of the laser pulse used in the simulation of the Coulomb model (left) and of the Morse model (right)

Figure 2 :
Figure 2: Snapshots illustrating the time evolution of the Coulomb model wave function during the grid-based simulation in the time interval from t = 0.0 to t = 54.0 a.u.As the time advances, the wave function becomes progressively more complicated, with nonlinear phase, amplitude oscillations, and localized features.To facilitate visualization, the wave function values have been rescaled so that the maximum value remains constant across all timeframes, matching the maximum value at t = 0. 14

Figure 3 :Figure 4 :
Figure 3: Same as in Fig. 2 but for the Morse model in the time interval from t = 0.0 a.u. to t = 270.0a.u.

FFECGs
than the Coulomb model.The reason for this behavior can be attributed to the maximum of the density of the second nucleus in the Morse model being shifted away from the reference nucleus by some distance (in the ground-state this distance is approximately equal to the equilibrium internuclear distance).This type of shifting does not happen in the Coulomb model.The shifting of the density maximum away from the reference nucleus required including more ECGs in the wave function.However, if this is done, the ECG expansion of the wave packet in the Morse model should be equally accurate as it is for the Coulomb model.The least-square fitting of a linear combination of ECGs to a grid wave function produces a wave packet that provides a representation of the grid function with uniform spatial quality.ECG representations obtained by minimization of energy-based functionals (e.g.

Figure 5 :
Figure 5: FFECG fits to the selected time frames extracted from the grid-based simulation trajectory of the Coulomb model.The error represents the difference in either the absolute value, the real part, or the imaginary part between the grid wave function and the optimized linear combination of Gaussian functions 22

Figure 6 :Figure 7 :
Figure 6: Same as Fig. 5 but for the Morse model 24

Figure 8 :
Figure 8: Left: Time-resolved observables computed for the Coulomb model, with the reference grid wave function (dashed lines), and with the fitted Gaussian wave function (solid lines).Right: Differences between grid and Gaussian observables.From top to bottom: expectation value of the field-free Hamiltonian, projection of the wave packet onto the groundstate wave function (the initial state), dipole moment expectation value, expectation value of the squared z-component of the angular momentum.

Figure 9 :
Figure 9: Same as Fig. 8 but for the Morse model.insets display evolution of the observables over the duration of the laser pulse (0-20 a.u).
Schrödinger equation for two 2D model systems that represent features which appear in quantum-dynamics timepropagation of the wave packet representing a diatomic neutral molecule interacting with a short intense laser pulse and performed without assuming the Born-Oppenheimer approximation.The grid wave functions obtained in consecutive time steps are fitted with a combination Gaussian functions that are 2D versions of more the general fully-flexible explicitly correlated Gaussians (FFECGs) with complex exponential parameters and complex shifts of the Gaussian centers.The fitting procedure employs the least-square procedure and involves growing the basis set of the Gaussians to provide a uniformly good fit for a representative set of time points obtained from the grid time propagation.The two models considered in the calculations involve a single-particle in the central potential represented by an attractive Coulomb interaction and a Morse potential.Based on the results obtained in the calculations we can expect that FFECGs will provide a good basis set for laser-induced non-BO dynamics of a diatomic molecule.Work on implementing FFECGs in molecular QD simulations is in progress.