A new class of ensemble conserving algorithms for approximate quantum dynamics: Theoretical formulation and model problems

Finite-Temperature Correlation Functions Obtained from Combined Real- and Imaginary-Time Propagation of Variational Thawed Gaussian Wavepackets

We apply the so-called variational Gaussian wavepacket approximation (VGA) for conducting both real- and imaginary-time dynamics to calculate thermal correlation functions. By considering strongly anharmonic systems, such as a quartic potential and a double-well potential at high and low temperatures, it is shown that this method is partially able to account for tunneling. This is contrary to other popular many-body methods, such as ring polymer molecular dynamics and the classical Wigner method, which fail in this respect. It is a historical peculiarity that no one has considered the VGA method for representing both the Boltzmann operator and the real-time propagation. This method should be well suited for molecular systems containing many atoms.

Incorporating Nuclear Quantum Effects in Molecular Dynamics with a Constrained Minimized Energy Surface

The accurate incorporation of nuclear quantum effects in large-scale molecular dynamics (MD) simulations remains a significant challenge. Recently, we combined constrained nuclear-electronic orbital (CNEO) theory with classical MD and obtained a new approach (CNEO-MD) that can accurately and efficiently incorporate nuclear quantum effects into classical simulations. In this Letter, we provide the theoretical foundation for CNEO-MD by developing an alternative formulation of the equations of motion for MD. In this new formulation, the expectation values of quantum nuclear positions evolve classically on an effective energy surface that is obtained from a constrained energy minimization procedure when solving for the quantum nuclear wave function, thus enabling the incorporation of nuclear quantum effects in classical MD simulations. For comparison with other existing approaches, we examined a series of model systems and found that this new MD approach is significantly more accurate than the conventional way of performing classical MD and generally outperforms centroid MD and ring-polymer MD in describing vibrations in these model systems.

Universal Matsubara time decay of quantum autocorrelations for Boltzmann particles

The general properties of time dependent autocorrelations in many-body quantum systems are here analyzed at thermodynamic equilibrium in the Boltzmann canonical ensemble at temperature T, by means of the exponential expansion theory (EET). It is shown that the Kubo-Martin-Schwinger (KMS) symmetry applied to the exponential expansion of the correlation leads to the existence of two different sets of decay modes (channels) here indicated as "Matsubara modes" and "system modes," respectively. The Matsubara modes are a series of pure decay channels with time constants representing a direct action of the thermostat upon the correlation, with a characteristic principal decay time τ_{1}=ℏ/(2πk_{B}T), where ℏ and k_{B} are the Planck and Boltzmann constants, and T is the temperature. Moreover, the KMS condition implies that the amplitudes pertaining to the even and odd contribution of the system modes to the quantum correlation are not independent. These two properties are quantum mechanical in nature and "universal," in the sense that they are present for any autocorrelation of a quantum system at equilibrium at a temperature T. The Matsubara modes' contribution to the time behavior of a quantum correlation is limited to times of the order of τ_{1}, which however can be comparable with some of the characteristic decay times of the system modes. In addition, since the parameters representing the overall time behavior of the quantum correlation can be given in terms of the parameters of its Kubo transform, the EET representation turns out to be useful in calculations exploiting the outputs of some widespread quantum simulation methods. A discussion of the properties of these relations is described in detail with numerical examples. The case of the velocity autocorrelation function of para hydrogen at low temperature is also reported as a final example for a real system.

Microscopic collective dynamics in liquid neon-deuterium mixtures: Inelastic neutron scattering and quantum simulations

In this paper a combined neutron scattering and quantum simulation study of the collective dynamics in liquid Ne-D_{2} mixtures, at a temperature of T=30K and in the wave-vector transfer range 4 nm^{-1}<q<51nm^{-1}, is presented. Two D_{2} concentrations are investigated, one close to 25% molar and the other close to 50% molar, together with pure Ne. The dynamic structure factor for the centers of mass of the two molecular species is extracted from the neutron scattering data and subsequently compared with that obtained from three different quantum simulation methods, such as ring polymer molecular dynamics and two slightly different versions of the Feynman-Kleinert approach. A general agreement is found, even though some discrepancies both among simulations, and between simulations and experimental data, can be observed. In order to clarify the physical meaning of the present spectroscopic results, an analysis of the longitudinal current spectral maxima is carried out showing the peculiarities of the D_{2} center-of-mass dynamics in these mixtures. A comparison with the centroid molecular dynamics results obtained for the D_{2} center-of-mass self-dynamics in the same liquid mixtures is finally proposed.

A Divergence-Free Wigner Transform of the Boltzmann Operator Based on an Effective Frequency Theory

The centroid effective frequency representation of path integrals as developed by Feynman and Kleinert was originally aimed at calculating partition functions and related quantities in the canonical ensemble. In its path integral formulation, only closed paths were relevant. This formulation has been used by the present authors in order to calculate the many-body Wigner function of the Boltzmann operator, which includes also open paths. This usage of the theory outside of the original intention can lead to mathematical divergence issues for potentials with barriers, particularly at low temperature. In the present paper, we modify the effective frequency theory of Feynman and Kleinert by also including open paths in its variational equations. In this way, a divergence-free approximation to the Boltzmann operator matrix elements is derived. This generalized version of Feynman and Kleinert’s formulation is thus more robust and can be applied to all types of barriers at all temperatures. This new version is used to calculate the Wigner functions of the Boltzmann operator for a quartic oscillator and for a double well potential and both static and dynamic properties are studied at several temperatures. The new theory is found to be essentially as precise as the original one. Its advantage is that it will always deliver a well-defined, even if approximate, Wigner function, which can, for instance, be used for sampling initial conditions for molecular dynamics simulations. As will be discussed, the theory can be systematically improved by including higher-order Fourier modes into the nonquadratic part of the trial action.

Machine learning phase space quantum dynamics approaches

Derived from phase space expressions of the quantum Liouville theorem, equilibrium continuity dynamics is a category of trajectory-based phase space dynamics methods, which satisfies the two critical fundamental criteria: conservation of the quantum Boltzmann distribution for the thermal equilibrium system and being exact for any thermal correlation functions (even of nonlinear operators) in the classical and harmonic limits. The effective force and effective mass matrix are important elements in the equations of motion of equilibrium continuity dynamics, where only the zeroth term of an exact series expansion of the phase space propagator is involved. We introduce a machine learning approach for fitting these elements in quantum phase space, leading to a much more efficient integration of the equations of motion. Proof-of-concept applications to realistic molecules demonstrate that machine learning phase space dynamics approaches are possible as well as competent in producing reasonably accurate results with a modest computation effort.

Ring-polymer, centroid, and mean-field approximations to multi-time Matsubara dynamics

Based on a recently developed generalization of Matsubara dynamics to the multi-time realm, we present a formal derivation of multi-time generalizations of ring-polymer molecular dynamics, thermostatted ring-polymer molecular dynamics (TRPMD), centroid molecular dynamics (CMD), and mean-field Matsubara dynamics. Additionally, we analyze the short-time accuracy of each methodology. We find that for multi-time correlation functions of linear operators, (T)RPMD is accurate up to order t³, while CMD is only correct up to t, indicating a degradation in the accuracy of these methodologies with respect to the single-time counterparts. The present work provides a firm justification for the use of path-integral-based approximations for the calculation of multi-time correlation functions.

Time dependence of quantum correlation functions

In the past few years, the exponential expansion analysis of time autocorrelation functions has provided profound insight into the leading microscopic processes driving the atomic-scale dynamics and has made it possible to highlight the presence and the role of various relaxation channels through which the fundamental correlation functions decay with time. Here we apply this method to the determination of the full time dependence of a correlation function c(t) in a quantum system at nonzero temperature, by making explicit its relationship with its Kubo transform c_{K}(t), which in some cases can be approximately computed with the presently available quantum simulation techniques. We obtain an exact expression for c(t) in terms of the exponential modes that describe the time behavior of c_{K}(t). The relative importance of the various modes in determining the overall shape of c(t) can then be studied in detail. This work extends to the full time domain the results of a previous paper [Guarini et al., Phys. Rev. Lett. 123, 135301 (2019)PRLTAO0031-900710.1103/PhysRevLett.123.135301], in which we employed the same method to calculate the zero time value of the velocity autocorrelation function, to obtain a microscopic description of the quantum mean kinetic energy in a fluid. In particular, we show that the decay constants and the frequencies of the dominant microscopic modes of c(t) are the same as those of c_{K}(t), but the dynamics of the quantum system also contains an additional term decaying on a time scale determined solely by temperature of the system.

Accurate calculation of zero point energy from molecular dynamics simulations of liquids and their mixtures

The two-phase thermodynamic (2PT) method is used to compute the zero point energy (ZPE) of several liquids and their mixtures. The 2PT method uses the density of states (DoS), which is computed from the velocity autocorrelation (VAC) function obtained from a short classical molecular dynamics trajectory. By partitioning the VAC and the DoS of a fluid into solid and gaslike components, quantum mechanical corrections to thermodynamical properties can be computed. The ZPE is obtained by combining the partition function of the quantum harmonic oscillator with the vibrational part of the solidlike DoS. The resulting ZPE is found to be in excellent agreement with both experimental and ab initio results. Solvent effects such as hydrogen bonding and polarization can be included by the utilization of ab initio density functional theory based molecular dynamics simulations. It is found that these effects significantly influence the DoS of water molecules. The obtained results demonstrate that the 2PT model is a powerful method for efficient ZPE calculations, in particular, to account for solvent effects and polarization.

Dynamical Origin of the Total and Zero-Point Kinetic Energy in a Quantum Fluid

By applying an exponential mode analysis to ring polymer molecular dynamics simulations of dense fluid parahydrogen, we find that the dynamical processes establishing the time behavior of the Kubo velocity autocorrelation function display the same nature as those already observed in high-density classical fluids. This result permits us to demonstrate that the exponential mode decomposition is a unique tool to identify which dynamical processes lead to one of the most notable properties of quantum fluids: the large value of the mean kinetic energy per particle and the importance of the zero-temperature quantum effects in determining it.

Sampling the thermal Wigner density via a generalized Langevin dynamics

The Wigner thermal density is a function of considerable interest in the area of approximate (linearized or semiclassical) quantum dynamics where it is employed to generate initial conditions for the propagation of appropriate sets of classical trajectories. In this paper, we propose an original approach to compute the Wigner density based on a generalized Langevin equation. The stochastic dynamics is nontrivial in that it contains a coordinate-dependent friction coefficient and a generalized force that couples momenta and coordinates. These quantities are, in general, not known analytically and have to be estimated via auxiliary calculations. The performance of the new sampling scheme is tested on standard model systems with highly nonclassical features such as relevant zero point energy effects, correlation between momenta and coordinates, and negative parts of the Wigner density. In its current brute force implementation, the algorithm, whose convergence can be systematically checked, is accurate and has only limited overhead compared to schemes with similar characteristics. We briefly discuss potential ways to further improve its numerical efficiency.

Multi-time formulation of Matsubara dynamics

Matsubara dynamics has recently emerged as the most general form of a quantum-Boltzmann-conserving classical dynamics theory for the calculation of single-time correlation functions. Here, we present a generalization of Matsubara dynamics for the evaluation of multitime correlation functions. We first show that the Matsubara approximation can also be used to approximate the two-time symmetrized double Kubo transformed correlation function. By a straightforward extension of these ideas to the multitime realm, a multitime Matsubara dynamics approximation can be obtained for the multitime fully symmetrized Kubo transformed correlation function. Although not a practical method, due to the presence of a phase-term, this multitime formulation of Matsubara dynamics represents a benchmark theory for future development of Boltzmann preserving semiclassical approximations to general higher order multitime correlation functions. It also reveals a connection between imaginary time-ordering in the path integral and the classical dynamics of multitime correlation functions.

Inclusion of nuclear quantum effects for simulations of nonlinear spectroscopy

The computation and interpretation of nonlinear vibrational spectroscopy is of vital importance for understanding a wide range of dynamical processes in molecular systems. Here, we introduce an approach to evaluate multi-time response functions in terms of multi-time double symmetrized Kubo transformed thermal correlation functions. Furthermore, we introduce a multi-time extension of ring polymer molecular dynamics to evaluate these Kubo transforms. Benchmark calculations show that the approximations are useful for short times even for nonlinear operators, providing a consistent improvement over classical simulations of multi-time correlation functions. The introduced methodology thus provides a practical way of including nuclear quantum effects in multi-time response functions of non-linear optical spectroscopy.

Approximating Matsubara dynamics using the planetary model: Tests on liquid water and ice

Matsubara dynamics is the quantum-Boltzmann-conserving classical dynamics which remains when real-time coherences are taken out of the exact quantum Liouvillian [T. J. H. Hele et al., J. Chem. Phys. 142, 134103 (2015)]; because of a phase-term, it cannot be used as a practical method without further approximation. Recently, Smith et al. [J. Chem. Phys. 142, 244112 (2015)] developed a "planetary" model dynamics which conserves the Feynman-Kleinert (FK) approximation to the quantum-Boltzmann distribution. Here, we show that for moderately anharmonic potentials, the planetary dynamics gives a good approximation to Matsubara trajectories on the FK potential surface by decoupling the centroid trajectory from the locally harmonic Matsubara fluctuations, which reduce to a single phase-less fluctuation particle (the "planet"). We also show that the FK effective frequency can be approximated by a direct integral over these fluctuations, obviating the need to solve iterative equations. This modification, together with use of thermostatted ring-polymer molecular dynamics, allows us to test the planetary model on water (gas-phase, liquid, and ice) using the q-TIP4P/F potential surface. The "planetary" fluctuations give a poor approximation to the rotational/librational bands in the infrared spectrum, but a good approximation to the bend and stretch bands, where the fluctuation lineshape is found to be motionally narrowed by the vibrations of the centroid.

Formulation of state projected centroid molecular dynamics: Microcanonical ensemble and connection to the Wigner distribution

A derivation of quantum statistical mechanics based on the concept of a Feynman path centroid is presented for the case of generalized density operators using the projected density operator formalism of Blinov and Roy [J. Chem. Phys. 115, 7822-7831 (2001)]. The resulting centroid densities, centroid symbols, and centroid correlation functions are formulated and analyzed in the context of the canonical equilibrium picture of Jang and Voth [J. Chem. Phys. 111, 2357-2370 (1999)]. The case where the density operator projects onto a particular energy eigenstate of the system is discussed, and it is shown that one can extract microcanonical dynamical information from double Kubo transformed correlation functions. It is also shown that the proposed projection operator approach can be used to formally connect the centroid and Wigner phase-space distributions in the zero reciprocal temperature β limit. A Centroid Molecular Dynamics (CMD) approximation to the state-projected exact quantum dynamics is proposed and proven to be exact in the harmonic limit. The state projected CMD method is also tested numerically for a quartic oscillator and a double-well potential and found to be more accurate than canonical CMD. In the case of a ground state projection, this method can resolve tunnelling splittings of the double well problem in the higher barrier regime where canonical CMD fails. Finally, the state-projected CMD framework is cast in a path integral form.