Challenges
Datasets
Workspaces
Papers
Discussions
Leaderboard
Log inSign up
Challenges
Datasets
Workspaces
Papers
Discussions
Leaderboard
Blog
Job Board
Q3AS

© 2026 Aqora Quantum S.A.S.

TermsPrivacyLegal Notice

arXiv:2605.07440

On the single-Hessian Gaussian wavepacket dynamics

Davide Barbiero, Jiří J. L. Vaníček

physics.chem-phmath-phphysics.comp-phquant-phSubmitted May 8, 2026• Announced May 11, 2026
arXivHTMLPDF

Abstract

Single-Hessian Gaussian wavepacket dynamics (GWD) significantly reduces the computational burden of Heller's local harmonic GWD, while maintaining comparable accuracy in approximating vibronic spectra. Here, we provide a new, symplectic derivation of the equations of motion of single-Hessian GWD and show that, unlike the local harmonic version, this method conserves the non-canonical symplectic structure on the manifold of Gaussian wavepackets and$-$for bounded dynamics in smooth potentials$-$avoids the drift of energy. Our numerical results suggest that, despite being much more efficient than the local harmonic variant, the single-Hessian GWD exhibits the same $\mathcal{O}(\hbar)$ asymptotic error in averages of observables. To further accelerate numerical simulations, we implement high-order time-stepping geometric integrators that are time-reversible and conserve the norm and symplectic structure exactly, regardless of the time step. In addition, we present explicit expressions for the exact evolution of the width of a single-Hessian Gaussian wavepacket in a general potential, as well as for the exact evolution of the whole wavepacket in a global harmonic potential. Using on-the-fly ab initio Gaussian wavepacket dynamics on the first excited-state surface of ammonia, we numerically confirm the conservation of geometric properties by these integrators and demonstrate that high-order integrators can enhance both accuracy and computational efficiency. We also compute the photoelectron spectrum of the difluorocarbene anion and the absorption spectrum of methylamine, and find that, in comparison with experiment, single-Hessian GWD outperforms global harmonic models and matches the accuracy of local harmonic GWD. Finally, we identify which spectral features are sensitive to the choice of reference Hessian.

Journal reference
J. Chem. Phys. 164, 234114 (2026)
DOI
10.1063/5.0327624

Related papers

Nonzero-temperature vibronic spectra of polyatomic molecules from a zero-temperature classical trajectory

Davide Barbiero, Jiří J. L. Vaníček

By combining coherence thermofield dynamics with the single-Hessian approximation, we enable simulations of low- to medium-resolution vibronic spectra of weakly anharmonic systems at nonzero temperatures, at negligible additional cost relative to zero-temperature calculations. Single-Hessian coherence thermofield Gaussian wavepacket dynamics is exact in any harmonic potential, provided that the reference Hessian is that of the final surface. When applied to Morse systems of increasing anharmonicity and varying temperature, this method successfully captures excited-state anharmonicity and key temperature-dependent spectral features, including hot bands and broadening. By combining the method with on-the-fly ab initio dynamics, we demonstrate its utility by computing the absorption spectra of naphthalene, aminocoumarin C450, and phenyl radical, and the photoelectron spectrum of SeO$_{2}^{-}$ . Within the ab initio single-Hessian approximation, after the zero-temperature spectrum is obtained at the cost of classical molecular dynamics (on the order of hours), all nonzero-temperature spectra are computed in seconds.

physics.chem-phCross-listedarXiv:2608.20075•Aug 21, 2026

A comparison of local and global single Gaussian approximations to time dynamics: One-dimensional systems

Ralph Eric Turner, Robert F. Snider

The detailed calculation of the dynamics of a chemical system is usually not considered due to the size and cost of the computation. It is thus useful to examine various approximation methods. Such methods first need to be tried out on simple systems, like one-dimensional motion. Here two approaches to approximating the solutions of the Schrödinger and von Neumann equations by single time-dependent Gaussians are explored and contrasted, explicitly for one-dimensional barrier penetration. The first approach, in which no tunneling occurs, is local in nature and characterized by an expansion of the equations of motion to second order about the average position of the Gaussian wave packet or about the average position and momentum of the Gaussian Wigner function. This approach was first introduced by Heller [E. J. Heller, J. Chem. Phys. 62, 1544 (1975)]. Here both Heller’s approach and a more general truncation method are considered. Indeed tunneling can be incorporated if second-order terms in the quantal von Neumann equation are included. However, the resulting dynamics is unstable for kinetic energies where the exact wavepacket (and Wigner function) splits into nonnegligible parts that represent reflection and transmission. In contrast, the second approach is a global method which is obtained by applying appropriate closure approximations to the equations of motion for the first- and second-order position and momentum expectation values. This method allows tunneling and is stable at all kinetic energies. It is also possible to approximate the global equations in such a way that the local equations are obtained.

The Journal of Chemical Physics•Jul 15, 1987•9 citations

Symplectic semiclassical wave packet dynamics

Tomoki Ohsawa, Melvin Boon-tiong Leok

The paper gives a symplectic-geometric account of semiclassical Gaussian wave packet dynamics. We employ geometric techniques to ‘strip away’ the symplectic structure behind the time-dependent Schrödinger equation and incorporate it into semiclassical wave packet dynamics. We show that the Gaussian wave packet dynamics is a Hamiltonian system with respect to the symplectic structure, apply the theory of symplectic reduction and reconstruction to the dynamics, and discuss dynamic and geometric phases in semiclassical mechanics. A simple harmonic oscillator example is worked out to illustrate the results. We show that the reduced semiclassical harmonic oscillator dynamics is completely integrable by finding the action–angle coordinates for the system, and calculate the associated dynamic and geometric phases explicitly. We also propose an asymptotic approximation of the potential term that provides a practical semiclassical correction term to the approximation by Heller. Numerical results for a simple one-dimensional example show that the semiclassical correction term realizes a semiclassical tunneling.

arXiv:1302.1139•Journal of Physics A Mathematical and Theoretical•Sep 19, 2013•31 citations

Efficient geometric integrators for nonadiabatic quantum dynamics. I. The adiabatic representation

Seonghoon Choi, Jiří J. L. Vaníček

Geometric integrators of the Schrödinger equation conserve exactly many invariants of the exact solution. Among these integrators, the split-operator algorithm is explicit and easy to implement but, unfortunately, is restricted to systems whose Hamiltonian is separable into kinetic and potential terms. Here, we describe several implicit geometric integrators applicable to both separable and nonseparable Hamiltonians and, in particular, to the nonadiabatic molecular Hamiltonian in the adiabatic representation. These integrators combine the dynamic Fourier method with the recursive symmetric composition of the trapezoidal rule or implicit midpoint method, which results in an arbitrary order of accuracy in the time step. Moreover, these integrators are exactly unitary, symplectic, symmetric, time-reversible, and stable and, in contrast to the split-operator algorithm, conserve energy exactly, regardless of the accuracy of the solution. The order of convergence and conservation of geometric properties are proven analytically and demonstrated numerically on a two-surface NaI model in the adiabatic representation. Although each step of the higher order integrators is more costly, these algorithms become the most efficient ones if higher accuracy is desired; a thousand-fold speedup compared to the second-order trapezoidal rule (the Crank-Nicolson method) was observed for a wavefunction convergence error of 10−10. In a companion paper [J. Roulet, S. Choi, and J. Vaníček, J. Chem. Phys. 150, 204113 (2019)], we discuss analogous, arbitrary-order compositions of the split-operator algorithm and apply both types of geometric integrators to a higher-dimensional system in the diabatic representation.

arXiv:1902.04661•The Journal of Chemical Physics•May 28, 2019•24 citations

Resonance Raman spectroscopy from ab initio Hagedorn wavepacket dynamics

Davide Barbiero, Léa Zupan, Jiří J. L. Vaníček

We present a practical, ab initio time-dependent method using Hagedorn wavepackets to simulate resonance Raman (RR) spectra of polyatomic molecules. Hagedorn functions---Gaussians multiplied by specific polynomials---are used to represent RR initial and final states because these functions are exact solutions to the time-dependent Schrödinger equation for at-most-quadratic potentials and can be propagated at zero cost beyond that of propagating the guiding Gaussian. Using efficient recursive formulae to compute overlaps between Hagedorn wavepackets, we can evaluate RR excitation profiles for arbitrary spectral signals, such as fundamental, overtone, combination, and hot bands. We then construct the Stokes and anti-Stokes RR spectra from these profiles. We first validate the method in a two-dimensional displaced, distorted, and Duschinsky-rotated harmonic model against numerically exact split-operator calculations. Then, we apply the method to compute RR spectra of anthracene by performing dynamics on a 66-dimensional harmonic potential energy surface constructed from density functional theory calculations.

physics.chem-phCross-listedarXiv:2608.19985•Aug 21, 2026

On Hagedorn wavepackets associated with different Gaussians

Jiří J. L. Vaníček, Zhan Tong Zhang

Abstract Hagedorn functions are carefully constructed generalizations of Hermite functions to the setting of many-dimensional squeezed and coupled harmonic systems. Wavepackets formed by superpositions of Hagedorn functions have been successfully used to solve the time-dependent Schrödinger equation exactly in harmonic systems and variationally in anharmonic systems. To evaluate typical observables, such as position or kinetic energy, it is sufficient to consider orthonormal Hagedorn functions with a single Gaussian center. Instead, we derive various relations between Hagedorn bases associated with different Gaussians, including their overlaps, which are necessary for evaluating quantities nonlocal in time, such as the time correlation functions needed for computing spectra. First, we use the Bogoliubov transformation to obtain the commutation relations between the ladder operators associated with different Gaussians. Then, instead of using numerical quadrature, we employ these commutation relations to derive exact recurrence relations for the overlap integrals between Hagedorn functions with different Gaussian centers. Finally, we present numerical experiments that demonstrate the accuracy and efficiency of our algebraic method as well as its suitability for treating problems in spectroscopy and chemical dynamics.

arXiv:2405.07880•Journal of Physics A Mathematical and Theoretical•Jan 30, 2025•6 citations
Single-Hessian Gaussian wavepacket dynamics (GWD) significantly reduces the computational burden of Heller's local harmonic GWD, while maintaining comparable accuracy in approximating vibronic spectra. Here, we provide a new, symplectic derivation of the equations of motion of single-Hessian GWD and show that, unlike the local harmonic version, this method conserves the non-canonical symplectic structure on the manifold of Gaussian wavepackets and−-−for bounded dynamics in smooth potentials−-−avoids the drift of energy. Our numerical results suggest that, despite being much more efficient than the local harmonic variant, the single-Hessian GWD exhibits the same O(ℏ)\mathcal{O}(\hbar)O(ℏ) asymptotic error in averages of observables. To further accelerate numerical simulations, we implement high-order time-stepping geometric integrators that are time-reversible and conserve the norm and symplectic structure exactly, regardless of the time step. In addition, we present explicit expressions for the exact evolution of the width of a single-Hessian Gaussian wavepacket in a general potential, as well as for the exact evolution of the whole wavepacket in a global harmonic potential. Using on-the-fly ab initio Gaussian wavepacket dynamics on the first excited-state surface of ammonia, we numerically confirm the conservation of geometric properties by these integrators and demonstrate that high-order integrators can enhance both accuracy and computational efficiency. We also compute the photoelectron spectrum of the difluorocarbene anion and the absorption spectrum of methylamine, and find that, in comparison with experiment, single-Hessian GWD outperforms global harmonic models and matches the accuracy of local harmonic GWD. Finally, we identify which spectral features are sensitive to the choice of reference Hessian.
By combining coherence thermofield dynamics with the single-Hessian approximation, we enable simulations of low- to medium-resolution vibronic spectra of weakly anharmonic systems at nonzero temperatures, at negligible additional cost relative to zero-temperature calculations. Single-Hessian coherence thermofield Gaussian wavepacket dynamics is exact in any harmonic potential, provided that the reference Hessian is that of the final surface. When applied to Morse systems of increasing anharmonicity and varying temperature, this method successfully captures excited-state anharmonicity and key temperature-dependent spectral features, including hot bands and broadening. By combining the method with on-the-fly ab initio dynamics, we demonstrate its utility by computing the absorption spectra of naphthalene, aminocoumarin C450, and phenyl radical, and the photoelectron spectrum of SeO2−_{2}^{-}2−​ . Within the ab initio single-Hessian approximation, after the zero-temperature spectrum is obtained at the cost of classical molecular dynamics (on the order of hours), all nonzero-temperature spectra are computed in seconds.
The detailed calculation of the dynamics of a chemical system is usually not considered due to the size and cost of the computation. It is thus useful to examine various approximation methods. Such methods first need to be tried out on simple systems, like one-dimensional motion. Here two approaches to approximating the solutions of the Schrödinger and von Neumann equations by single time-dependent Gaussians are explored and contrasted, explicitly for one-dimensional barrier penetration. The first approach, in which no tunneling occurs, is local in nature and characterized by an expansion of the equations of motion to second order about the average position of the Gaussian wave packet or about the average position and momentum of the Gaussian Wigner function. This approach was first introduced by Heller [E. J. Heller, J. Chem. Phys. 62, 1544 (1975)]. Here both Heller’s approach and a more general truncation method are considered. Indeed tunneling can be incorporated if second-order terms in the quantal von Neumann equation are included. However, the resulting dynamics is unstable for kinetic energies where the exact wavepacket (and Wigner function) splits into nonnegligible parts that represent reflection and transmission. In contrast, the second approach is a global method which is obtained by applying appropriate closure approximations to the equations of motion for the first- and second-order position and momentum expectation values. This method allows tunneling and is stable at all kinetic energies. It is also possible to approximate the global equations in such a way that the local equations are obtained.
The paper gives a symplectic-geometric account of semiclassical Gaussian wave packet dynamics. We employ geometric techniques to ‘strip away’ the symplectic structure behind the time-dependent Schrödinger equation and incorporate it into semiclassical wave packet dynamics. We show that the Gaussian wave packet dynamics is a Hamiltonian system with respect to the symplectic structure, apply the theory of symplectic reduction and reconstruction to the dynamics, and discuss dynamic and geometric phases in semiclassical mechanics. A simple harmonic oscillator example is worked out to illustrate the results. We show that the reduced semiclassical harmonic oscillator dynamics is completely integrable by finding the action–angle coordinates for the system, and calculate the associated dynamic and geometric phases explicitly. We also propose an asymptotic approximation of the potential term that provides a practical semiclassical correction term to the approximation by Heller. Numerical results for a simple one-dimensional example show that the semiclassical correction term realizes a semiclassical tunneling.
Geometric integrators of the Schrödinger equation conserve exactly many invariants of the exact solution. Among these integrators, the split-operator algorithm is explicit and easy to implement but, unfortunately, is restricted to systems whose Hamiltonian is separable into kinetic and potential terms. Here, we describe several implicit geometric integrators applicable to both separable and nonseparable Hamiltonians and, in particular, to the nonadiabatic molecular Hamiltonian in the adiabatic representation. These integrators combine the dynamic Fourier method with the recursive symmetric composition of the trapezoidal rule or implicit midpoint method, which results in an arbitrary order of accuracy in the time step. Moreover, these integrators are exactly unitary, symplectic, symmetric, time-reversible, and stable and, in contrast to the split-operator algorithm, conserve energy exactly, regardless of the accuracy of the solution. The order of convergence and conservation of geometric properties are proven analytically and demonstrated numerically on a two-surface NaI model in the adiabatic representation. Although each step of the higher order integrators is more costly, these algorithms become the most efficient ones if higher accuracy is desired; a thousand-fold speedup compared to the second-order trapezoidal rule (the Crank-Nicolson method) was observed for a wavefunction convergence error of 10−10. In a companion paper [J. Roulet, S. Choi, and J. Vaníček, J. Chem. Phys. 150, 204113 (2019)], we discuss analogous, arbitrary-order compositions of the split-operator algorithm and apply both types of geometric integrators to a higher-dimensional system in the diabatic representation.
We present a practical, ab initio time-dependent method using Hagedorn wavepackets to simulate resonance Raman (RR) spectra of polyatomic molecules. Hagedorn functions---Gaussians multiplied by specific polynomials---are used to represent RR initial and final states because these functions are exact solutions to the time-dependent Schrödinger equation for at-most-quadratic potentials and can be propagated at zero cost beyond that of propagating the guiding Gaussian. Using efficient recursive formulae to compute overlaps between Hagedorn wavepackets, we can evaluate RR excitation profiles for arbitrary spectral signals, such as fundamental, overtone, combination, and hot bands. We then construct the Stokes and anti-Stokes RR spectra from these profiles. We first validate the method in a two-dimensional displaced, distorted, and Duschinsky-rotated harmonic model against numerically exact split-operator calculations. Then, we apply the method to compute RR spectra of anthracene by performing dynamics on a 66-dimensional harmonic potential energy surface constructed from density functional theory calculations.
Abstract Hagedorn functions are carefully constructed generalizations of Hermite functions to the setting of many-dimensional squeezed and coupled harmonic systems. Wavepackets formed by superpositions of Hagedorn functions have been successfully used to solve the time-dependent Schrödinger equation exactly in harmonic systems and variationally in anharmonic systems. To evaluate typical observables, such as position or kinetic energy, it is sufficient to consider orthonormal Hagedorn functions with a single Gaussian center. Instead, we derive various relations between Hagedorn bases associated with different Gaussians, including their overlaps, which are necessary for evaluating quantities nonlocal in time, such as the time correlation functions needed for computing spectra. First, we use the Bogoliubov transformation to obtain the commutation relations between the ladder operators associated with different Gaussians. Then, instead of using numerical quadrature, we employ these commutation relations to derive exact recurrence relations for the overlap integrals between Hagedorn functions with different Gaussian centers. Finally, we present numerical experiments that demonstrate the accuracy and efficiency of our algebraic method as well as its suitability for treating problems in spectroscopy and chemical dynamics.