New coherent state representation for the imaginary time propagator with applications to forward-backward semiclassical initial value representations of correlation functions

There have been quite a few attempts in recent years to provide an initial value coherent state representation for the imaginary time propagator exp(-betaH). The most notable is the recent time evolving Gaussian approximation of Frantsuzov and Mandelshtam [J. Chem. Phys. 121, 9247 (2004)] which may be considered as an expansion of the imaginary time propagator in terms of coherent states whose momentum is zero. In this paper, a similar but different expression is developed in which exp(-betaH) is represented in a series whose terms are weighted phase space averages of coherent states. Such a representation allows for the formulation of a new and simplified forward-backward semiclassical initial value representation expression for thermal correlation functions.     

Imaginary time Gaussian dynamics of the Ar3 cluster

Semiclassical Gaussian approximations to the Boltzmann operator have become an important tool for the investigation of thermodynamic properties of clusters of atoms at low temperatures. Usually, numerically expensive thawed Gaussian variants are applied. In this article, we introduce a numerically much cheaper frozen Gaussian approximation to the imaginary time propagator with a width matrix especially suited for the dynamics of clusters. The quality of the results is comparable to that of thawed Gaussian methods based on the single-particle ansatz. We apply the method to the argon trimer and investigate the dissociation process of the cluster. The results clearly show a classical-like transition from a bounded moiety to three free particles at a temperature T ≈ 20 K, whereas previous studies of the system were not able to resolve this transition. Quantum effects, i.e., differences with the purely classical case manifest themselves in the low-temperature behavior of the mean energy and specific heat as well as in a slight shift of the transition temperature. We also discuss the influence of an artificial confinement of the atoms usually introduced to converge numerical computations. The results show that restrictive confinements often implemented in studies of clusters can influence the thermodynamic properties drastically. This finding may have implications on other studies of atomic clusters.     

Path integral ground state with a fourth-order propagator: application to condensed helium

Ground state properties of condensed helium are calculated using the path integral ground state (PIGS) method. A fourth-order approximation is used as short (imaginary) time propagator. We compare our results with those obtained with other quantum Monte Carlo (QMC) techniques and different propagators. For this particular application, we find that the fourth-order propagator performs comparably to the pair product approximation, and is far superior to the primitive approximation. Results obtained for the equation of state of condensed helium show that PIGS compares favorably to other QMC methods traditionally utilized for this type of calculation.     

A prefactor free semiclassical initial value series representation of the propagator

A new class of prefactor free semiclassical initial value representations (SCIVR) of the quantum propagator is presented. The derivation is based on the physically motivated demand, that on the average in phase space and in time, the propagator obey the exact quantum equation of motion. The resulting SCIVR series representation of the exact quantum propagator is also free of prefactors. When using a constant width parameter, the prefactor free SCIVR propagator is identical to the frozen Gaussian propagator of Heller [J. Chem. Phys. 75, 2923 (1981)]. A numerical study of the prefactor free SCIVR series is presented for scattering through a double slit potential, a system studied extensively previously by Gelabert et al. [J. Chem. Phys. 114, 2572 (2001)]. As a basis for comparison, the SCIVR series is also computed using the optimized Herman-Kluk SCIVR. We find that the sum of the zeroth order and the first order terms in the series suffice for an accurate determination of the diffraction pattern. The same exercise, but using the prefactor free propagator series needs also the second order term in the series, however the numerical effort is not greater than that needed for the Herman-Kluk propagator, since one does not need to compute the monodromy matrix elements at each point in time. The numerical advantage of the prefactor free propagator grows with increasing dimensionality of the problem.     

Renormalization of the frozen Gaussian approximation to the quantum propagator

The frozen Gaussian approximation to the quantum propagator may be a viable method for obtaining "on the fly" quantum dynamical information on systems with many degrees of freedom. However, it has two severe limitations, it rapidly loses normalization and one needs to know the Gaussian averaged potential, hence it is not a purely local theory in the force field. These limitations are in principle remedied by using the Herman-Kluk (HK) form for the semiclassical propagator. The HK propagator approximately conserves unitarity for relatively long times and depends only locally on the bare potential and its second derivatives. However, the HK propagator involves a much more expensive computation due to the need for evaluating the monodromy matrix elements. In this paper, we (a) derive a new formula for the normalization integral based on a prefactor free HK propagator which is amenable to "on the fly" computations; (b) show that a frozen Gaussian version of the normalization integral is not readily computable "on the fly"; (c) provide a new insight into how the HK prefactor leads to approximate unitarity; and (d) how one may construct a prefactor free approximation which combines the advantages of the frozen Gaussian and the HK propagators. The theoretical developments are backed by numerical examples on a Morse oscillator and a quartic double well potential.     

A justification for a nonadiabatic surface hopping Herman-Kluk semiclassical initial value representation of the time evolution operator

A justification is given for the validity of a nonadiabatic surface hopping Herman-Kluk (HK) semiclassical initial value representation (SC-IVR) method. The method is based on a propagator that combines the single surface HK SC-IVR method [J. Chem. Phys. 84, 326 (1986)] and Herman's nonadiabatic semiclassical surface hopping theory [J. Chem. Phys. 103, 8081 (1995)], which was originally developed using the primitive semiclassical Van Vleck propagator. We show that the nonadiabatic HK SC-IVR propagator satisfies the time-dependent Schrodinger equation to the first order of variant Planck's over 2pi and the error is O(variant Planck's over 2pi(2)). As a required lemma, we show that the stationary phase approximation, under current assumptions, has an error term variant Planck's over 2pi(1) order higher than the leading term. Our derivation suggests some changes to the previous development, and it is shown that the numerical accuracy in applications to Tully's three model systems in low energies is improved.     

A semiclassical initial-value representation for quantum propagator and boltzmann operator

Starting from the position-momentum integral representation, we apply the correction operator method to the derivation of a uniform semiclassical approximation for the quantum propagator and then extend it to approximate the Boltzmann operator. In this approach, the involved classical dynamics is determined by the method itself instead of given beforehand. For the approximate Boltzmann operator, the corresponding classical dynamics is governed by a complex Hamiltonian, which can be described as a pair of real Hamiltonian systems. It is demonstrated that the semiclassical Boltzmann operator is exact for linear systems. A quantum propagator in the complex time is thus proposed and preliminary numerical results show that it is a reasonable approximation for calculating thermal correlation functions of general systems. © 2018 Wiley Periodicals, Inc.     

Time-dependent self-consistent field approximation for semiclassical dynamics using gaussian wavepackets in phase space

Gaussian wavepackets in phase space, which are constructed to have the exact first and second moments with respect to the coordinates and momenta, are used to develop self-consistent equations of motion for the semiclassical time evolution of interacting anharmonic systems. The equations apply to pure as well as to mixed states and may, therefore, be particularly useful for molecular dynamics in condensed phases.     

Quadratic response functions in a second-order polarization propagator framework

The linear and quadratic response functions have been derived for an exact state, based on an exponential parametrization of the time evolution consisting of products of exponentials for orbital rotations and for higher-order excitations. Truncating the linear response function such that the response function itself and its pole structure is correct to second order in M?ller-Plesset perturbation theory, we arrive at the second-order polarization propagator approximation (SOPPA). Previous derivations of SOPPA have used the superoperator formalism, making the extension of SOPPA to quadratic and higher order response functions difficult. The derivation of the quadratic response function is described in detail, allowing molecular properties such as hyperpolarizabilities, two-photon cross sections, and excited-state properties to be calculated using the SOPPA model.     

Long-time and unitary properties of semiclassical initial value representations

We numerically compare the semiclassical "frozen Gaussian" Herman-Kluk propagator [Chem. Phys. 91, 27 (1984)] and the "thawed Gaussian" propagator put forward recently by Baranger et al. [J. Phys. A 34, 7227 (2001)] by studying the quantum dynamics in some nonlinear one-dimensional potentials. The reasons for the lack of long-time accuracy and norm conservation in the latter method are uncovered. We amend the thawed Gaussian propagator with a global harmonic approximation for the stability of the trajectories and demonstrate that this revised propagator is a true alternative to the Herman-Kluk propagator with similar accuracy.     

Quantum mechanical correlation functions, maximum entropy analytic continuation, and ring polymer molecular dynamics

The maximum entropy analytic continuation (MEAC) and ring polymer molecular dynamics (RPMD) methods provide complementary approaches to the calculation of real time quantum correlation functions. RPMD becomes exact in the high temperature limit, where the thermal time betavariant Planck's over 2pi tends to zero and the ring polymer collapses to a single classical bead. MEAC becomes most reliable at low temperatures, where betavariant Planck's over 2pi exceeds the correlation time of interest and the numerical imaginary time correlation function contains essentially all of the information that is needed to recover the real time dynamics. We show here that this situation can be exploited by combining the two methods to give an improved approximation that is better than either of its parts. In particular, the MEAC method provides an ideal way to impose exact moment (or sum rule) constraints on a prior RPMD spectrum. The resulting scheme is shown to provide a practical solution to the "nonlinear operator problem" of RPMD, and to give good agreement with recent exact results for the short-time velocity autocorrelation function of liquid parahydrogen. Moreover these improvements are obtained with little extra effort, because the imaginary time correlation function that is used in the MEAC procedure can be computed at the same time as the RPMD approximation to the real time correlation function. However, there are still some problems involving long-time dynamics for which the RPMD+MEAC combination is inadequate, as we illustrate with an example application to the collective density fluctuations in liquid orthodeuterium.     

Finite temperature application of the corrected propagator method to reactive dynamics in a condensed-phase environment

The recently proposed mixed quantum-classical method is extended to applications at finite temperatures. The method is designed to treat complex systems consisting of a low-dimensional quantum part (the primary system) coupled to a dissipative bath described classically. The method is based on a formalism showing how to systematically correct the approximate zeroth-order evolution rule. The corrections are defined in terms of the total quantum Hamiltonian and are taken to the classical limit by introducing the frozen Gaussian approximation for the bath degrees of freedom. The evolution of the primary system is governed by the corrected propagator yielding the exact quantum dynamics. The method has been tested on a standard model system describing proton transfer in a condensed-phase environment: a symmetric double-well potential bilinearly coupled to a bath of harmonic oscillators. Flux correlation functions and thermal rate constants have been calculated at two different temperatures for a range of coupling strengths. The results have been compared to the fully quantum simulations of Topaler and Makri [J. Chem. Phys. 101, 7500 (1994)] with the real path integral method.     

Error analysis of molecular dynamics and fractal time approximants from a combinatorial perspective

Trotter's theorem forms the theoretical basis of most modern molecular dynamics. In essence this theorem states that a time displacement operator (a Lie operator) constructed by exponentiating a sum of noncommuting operators can be approximated by a product of single operators provided the time interval is "very small." In theory "very small" implies infinitesimally small (at which point the approximate product becomes exact), while in practical analysis a finite time interval is divided into several small subintervals or steps. It follows, therefore, that the larger the number of steps the better the approximation to the exact time displacement operator. The question therefore arises: How many steps are sufficient? For bounded operators, standard theorems are available to provide the answer. In this paper we show that a very simple combinatorial formula can be derived which allows the computation of the global differences (as a function of the number of steps) between the Taylor coefficients of the exact time displacement operator and an approximate one constructed by using a finite number of steps. The formula holds for both bounded and nonbounded operators and shows, quantitatively, what is qualitatively expected-that the error decreases with increasing number of steps. Furthermore, the formula applies irrespective of the complexity of the system, boundary conditions, or the thermodynamic ensemble employed for averaging the initial conditions. The analysis yields explicit expressions for the Taylor coefficients which are then used to compute the errors. In the case of the algorithmically based practical numerical simulations in which fixed, albeit small, steps are repeatedly applied, the rise in the number of steps does not reduce the size of the steps but increases the total time of interest. The combinatorial formula shows that, here, the errors diverge. Furthermore, this work can be used to supplement other efforts such as the use of shadow Hamiltonians where the truncation of the series expansion of the latter will produce errors in the higher order propagator moments.     

Extended hydrodynamic approach to quantum-classical nonequilibrium evolution. I. Theory

A mixed quantum-classical formulation is developed for a quantum subsystem in strong interaction with an N-particle environment, to be treated as classical in the framework of a hydrodynamic representation. Starting from the quantum Liouville equation for the N-particle distribution and the corresponding reduced single-particle distribution, exact quantum hydrodynamic equations are obtained for the momentum moments of the single-particle distribution coupled to a discretized quantum subsystem. The quantum-classical limit is subsequently taken and the resulting hierarchy of equations is further approximated by various closure schemes. These include, in particular, (i) a Grad-Hermite-type closure, (ii) a Gaussian closure at the level of a quantum-classical local Maxwellian distribution, and (iii) a dynamical density functional theory approximation by which the hydrodynamic pressure term is replaced by a free energy functional derivative. The latter limit yields a mixed quantum-classical formulation which has previously been introduced by I. Burghardt and B. Bagchi, Chem. Phys. 134, 343 (2006).     

Study of simulation method of time evolution of atomic and molecular systems by quantum electrodynamics

We discuss a method to follow step‐by‐step time evolution of atomic and molecular systems based on quantum electrodynamics. Our strategy includes expanding the electron field operator by localized wavepackets to define creation and annihilation operators and following the time evolution using the equations of motion of the field operator in the Heisenberg picture. We first derive a time evolution equation for the excitation operator, the product of two creation or annihilation operators, which is necessary for constructing operators of physical quantities such as the electronic charge density operator. We, then, describe our approximation methods to obtain time differential equations of the electronic density matrix, which is defined as the expectation value of the excitation operator. By solving the equations numerically, we show “electron‐positron oscillations,” the fluctuations originated from virtual electron‐positron pair creations and annihilations, appear in the charge density of a hydrogen atom and molecule. We also show that the period of the electron‐positron oscillations becomes shorter by including the self‐energy process, in which the electron emits a photon and then absorbs it again, and it can be interpreted as the increase in the electron mass due to the self‐energy. © 2014 Wiley Periodicals, Inc.     

The Kramers-Moyal expansion for electrochemical stochastic diffusion

It is proved that there is a general stochastic equation, according to which any random process in the transient mode can be presented by spatially homogeneous Kramers-Moyal expansion. In the electrochemical stochastic diffusion, an integral of the fluctuation component of electrode potential over the time plays the role of spatial coordinate. Based on these two facts, we derived a spatially homogeneous Kramers-Moyal expansion for the propagator of electrochemical stochastic diffusion. By using the limiting transition to long observation times, we obtained a time and spatially homogeneous asymptotic Kramers-Moyal expansion for the propagator of asymmetric non-Gaussian electrochemical stochastic diffusion. Under the conditions of Gaussian electrochemical noise, the asymptotic Kramers-Moyal expansion turns into the Einstein stochastic diffusion equation. The method of determining time and spatially homogeneous asymptotic Kramers-Moyal expansion for the propagator of asymmetric non-Gaussian electrochemical stochastic diffusion may be useful in the stochastic theory of slow electrochemical discharge and in the electrochemical noise diagnostics.     

A semiclasical justification for the use of non-spreading wavepackets in dynamics calculations

A justification is given for the use of non-spreading or frozen gaussian packets in dynamics calculations. In this work an initial wavefunction or quantum density operator is expanded in a complete set of grussian wavepackets. It is demonstrated that the time evolution of this wavepacket expansion for the quantum wavefunction or density is correctly given within the approximations employed by the classical propagation of the avarage position and momentum of each gaussian packet, holding the shape of these individual gaussians fixed. The semiclassical approximation is employed for the quantum propagator and the stationary phase approximation for certain integrals is utilized in this derivation. This analysis demonstrates that the divergence of the classical trajectories associated with the individual gaussian packets accounts for the changes in shape of the quantum wavefunction or density, as has been suggested on intuitive grounds by Heller. The method should be exact for quadratic potentials and this is verified by explicitly applying it for the harmonic oscillator example.     

Continuum limit frozen Gaussian approximation for the reduced thermal density matrix of dissipative systems

A continuum limit frozen Gaussian approximation is formulated for the reduced thermal density matrix for dissipative systems. The imaginary time dynamics is obtained from a novel generalized Langevin equation for the system coordinates. The method is applied to study the thermal density in a double well potential in the presence of Ohmic-like friction. We find that the approximation describes correctly the delocalization of the density due to quantization of the vibrations in the well. It also accounts for the friction induced reduction of the tunneling density in the barrier region.     

Real time correlation function in a single phase space integral beyond the linearized semiclassical initial value representation

It is shown how quantum mechanical time correlation functions [defined, e.g., in Eq. (1.1)] can be expressed, without approximation, in the same form as the linearized approximation of the semiclassical initial value representation (LSC-IVR), or classical Wigner model, for the correlation function [cf. Eq. (2.1)], i.e., as a phase space average (over initial conditions for trajectories) of the Wigner functions corresponding to the two operators. The difference is that the trajectories involved in the LSC-IVR evolve classically, i.e., according to the classical equations of motion, while in the exact theory they evolve according to generalized equations of motion that are derived here. Approximations to the exact equations of motion are then introduced to achieve practical methods that are applicable to complex (i.e., large) molecular systems. Four such methods are proposed in the paper--the full Wigner dynamics (full WD) and the second order WD based on "Wigner trajectories" [H. W. Lee and M. D. Scully, J. Chem. Phys. 77, 4604 (1982)] and the full Donoso-Martens dynamics (full DMD) and the second order DMD based on "Donoso-Martens trajectories" [A. Donoso and C. C. Martens, Phys. Rev. Lett. 8722, 223202 (2001)]--all of which can be viewed as generalizations of the original LSC-IVR method. Numerical tests of the four versions of this new approach are made for two anharmonic model problems, and for each the momentum autocorrelation function (i.e., operators linear in coordinate or momentum operators) and the force autocorrelation function (nonlinear operators) have been calculated. These four new approximate treatments are indeed seen to be significant improvements to the original LSC-IVR approximation.     

Time-dependent approaches for the calculation of intersystem crossing rates

We present three formulas for calculating intersystem crossing rates in the Condon approximation to the golden rule by means of a time-dependent approach: an expression using the full time correlation function which is exact for harmonic oscillators, a second-order cumulant expansion, and a short-time approximation of this expression. While the exact expression and the cumulant expansion require numerical integration of the time correlation function, the integration of the short-time expansion can be performed analytically. To ensure convergence in the presence of large oscillations of the correlation function, we use a Gaussian damping function. The strengths and weaknesses of these approaches as well as the dependence of the results on the choice of the technical parameters of the time integration are assessed on four test examples, i.e., the nonradiative S(1) ? T(1) transitions in thymine, phenalenone, flavone, and porphyrin. The obtained rate constants are compared with previous results of a time-independent approach. Very good agreement between the literature values and the integrals over the full time correlation functions are observed. Furthermore, the comparison suggests that the cumulant expansion approximates the exact expression very well while allowing the interval of the time integration to be significantly shorter. In cases with sufficiently high vibrational density of states also the short-time approximation yields rates in good agreement with the results of the exact formula. A great advantage of the time-dependent approach over the time-independent approach is its excellent computational efficiency making it the method of choice in cases of large energy gaps, large numbers of normal modes, and high densities of final vibrational states.