Research

Phase space with coordinate-momentum variables bridges a rigorous formulation of classical mechanics and an exact interpretation of quantum mechanics. Conventional infinite coordinate-momentum phase space had been developed for systems with continuous degrees of freedom. We have proposed a generic framework to construct rigorous representations on constraint coordinate-momentum phase space (CPS) for F-state quantum systems, where exact equations of motion on phase space are simply linear and that the sum of state population is one is satisfied for each phase point.

In the CPS formulation, creation and annihilation operators are defined such that the space is complete for any combined excitation. Commutation and anti-commutation relations are then naturally derived, which show that the underlying degrees of freedom are neither bosons nor fermions. It leads to a generalized commutation relation with a commutator matrix. The CPS formulation with commutator variables is diffeomorphic to the complex Stiefel manifolds. By establishing CPS for discrete degrees of freedom and employing infinite phase space for continuous DOFs, we have developed the exact generalized coordinate-momentum phase space formulations for composite systems in chemistry, physics, information, etc.

It lays the solid foundation of developing consistent trajectory-based nonadiabatic dynamics methods. It offers a novel means to derive and understand isomorphic models where electronic and nuclear DOFs are treated on the same footing as pioneered by the Meyer-Miller-Stock-Thoss mapping Hamiltonian model. We have used a complex Stiefel manifold U(F)/U(F-1) and the Abel integral equation to construct a new class of exact phase space representations of two-state systems and proved that the heuristic Cotton-Miller triangle window function approach is a special case of this class.

We have further introduced the isomorphism between the multi-state Hamiltonian and the second-quantized many-electron Hamiltonian (with only 1-electron interactions) and demonstrated that all methods developed for the former can be used for the latter, and vice versa. It indicates that phase space mapping approaches for multi-state nonadiabatic systems can be directly used for studying nonequilibrium quantum transport processes, e.g., the electronic current of the resonant level (Landauer) model for a quantum dot state coupled to two electrodes, the carrier mobility of the Holstein model for organic semiconductors, and so on.

We have proposed nonadiabatic field (NaF), a conceptually novel approach for nonadiabatic dynamics with independent trajectories on quantum coordinate-momentum phase space. NaF involves the nonadiabatic force arising from the nonadiabatic coupling between different electronic states, in addition to the adiabatic force contributed by a single adiabatic electronic state. NaF is capable of faithfully describing the main features of the quantum mechanical behavior of both electrons and nuclei in nonadiabatic transition processes in a broad regime, from systems where the relevant electronic states keep coupled in a wide range and the nonadiabatic force plays an important role, for which it is difficult for surface hopping approaches to describe relatively-long time behavior, to systems where the bifurcation characteristic of nuclear motion is essential, which Ehrenfest(-like) dynamics methods fail to capture. NaF has been used to study population dynamics, electron transfer reaction rates, carrier mobilities, electronic spectra, cavity-modified dynamics for both gas phase and condensed phase systems.

We were the first to develop a practical quantum dynamics method satisfying the two critical fundamental criteria: conservation of the quantum Boltzmann distribution for the thermal equilibrium system and being exact for any thermal correlation functions (of even nonlinear operators, i.e., nonlinear functions of position or momentum operators) in the classical and harmonic limits. The so-called path integral Liouville dynamics (PILD) [J. Chem. Phys. 140, 224107 (2014)] was derived from the elaborate combination of the imaginary time path integral and the three new theoretical frameworks that we formulated for generating trajectory-based dynamics methods in the quantum phase space [J. Chem. Phys. 134, 104101 (2011); 134, 104102 (2011); 134, 194110 (2011). Scientia Sinica Chimica, 46, 27 (2016)], of which the equilibrium continuity dynamics (ECD) category is the most fundamental.

PILD was applied to vibrational dynamics of several simple but representative real molecular systems (OH, water, ammonia, and methane) [J. Chem. Phys. 144, 034307 (2016)]. The dipole-derivative autocorrelation function was employed to obtain the infrared spectrum as a function of temperature and isotopic substitution. Comparison to the exact vibrational frequency shows that PILD produces a reasonably accurate peak position with a relatively small full width at half maximum. PILD offers a potentially useful trajectory-based quantum dynamics approach to compute vibrational spectra of molecular systems.

More recently, we have shown that the equations of motion of ECD involves only the zeroth-order term of an exact series expansion of the phase space propagator and that higher order corrections can consistently improve the accuracy [J. Chem. Phys. 154, 184104 (2021)]. The effective force and effective mass matrix are important elements in ECD. 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.

Path integral molecular dynamics (PIMD) is a prevailing computational tool to study nuclear quantum statistics for real molecular systems at thermal equilibrium. We have employed the phase space evolution operator to develop a unified framework that includes both stochastic and deterministic thermostatting algorithms, including such as the Andersen thermostat, Langevin dynamics, Nosé-Hoover chain (including Nosé-Hoover thermostat), etc. [J. Chem. Phys. 145, 024103 (2016); 147, 034109 (2017)] and covers various Monte Carlo or MD barostatting algorithms [J. Chem. Theory Comput. (2025)]. Compared to most conventional PIMD algorithms that can be unified in the “side”/“end” schemes of the framework, the "middle" thermostat scheme introduced by us significantly improves the sampling efficiency as well as accuracy for general molecular systems of canonical ensembles (with constant volume and temperature) as well as of isobaric-isothermal ensembles (with constant pressure and temperature). In addition to numerical examples for real systems, we have presented the analytical analysis to prove that the “middle” scheme leads to exact results for PIMD in the harmonic limit regardless of the thermostat parameter and of the finite time interval (as long as the propagation is stable).

In the unified framework we have revealed that, in addition to the usual density evolution, there exists another type of discrete evolution that may not correspond to a continuous, real dynamical counterpart of the Langevin equation [J. Chem. Phys. 147, 184104 (2017)]. This virtual dynamics case is also amenable to the desired Boltzmann distribution for various stochastic thermostatting methods [Chin. J. Chem. Phys. 30, 735 (2017)]. By analyzing the asymptotic distribution and characteristic correlation time, we have shown that, for good numerical performance in efficiency as well as accuracy, one may choose the thermostat parameter (of the “middle” scheme) in a wide range from around the optimal value to the high value limit.

We have further used the unified framework to develop a novel practical PIMD methodology in either of the diabatic and adiabatic representations for studying exact quantum statistics of large multi-electronic-state systems in thermal equilibrium when the Born-Oppenheimer approximation, Condon approximation, and harmonic bath approximation are no longer valid [J. Chem. Phys. 148, 102319 (2018)].

More recently, we have employed the "middle" scheme for molecular dynamics and path integral molecular dynamics simulations for efficient and accurate evaluation of density, isobaric heat capacity, isothermal compressibility, thermal expansion coefficient, and other physical properties of general molecular systems under conditions of controlled pressure and temperature.

Understanding vibrational spectra (infrared (IR) and Raman spectra) at the atomistic level is important for the elucidation of dynamical processes in complex/large molecular systems. The correlation approach offers a bridge to connect the microscopic motion to the experimental spectrum. Because the dipole moment/polarizability (or its time derivative) is often a highly nonlinear function of coordinates or/and momenta, we have employed the linearized semiclassical initial value representation for quantum dynamical simulations of liquid water (and heavy water) under ambient conditions based on an ab initio based, flexible, polarizable model. Besides applying our new PIMD algorithms in the “middle” scheme, we have proposed a new technique (by using the finite difference along the direction of a related vector) to significantly reduce the computational cost for preparing the initial condition for the correlation function.

It is shown that quantum dynamical effects play a critical role in reproducing the peaks in the intermediate region between the librational and bending bands, those between the bending and stretching bands, and the double-peak in the stretching band in the experimental isotropic Raman spectrum. In contrast, quantum dynamical effects are important but less decisive in the anisotropic Raman spectrum. By selectively freezing either the intramolecular O-H stretching or H-O-H bending mode, we demonstrate that the peak in the intermediate region (2000-2400 cm-1) of the isotropic Raman spectrum arises from the interplay of the stretching and bending motions while a substantial part of the peak in the same intermediate region of the anisotropic Raman spectrum may be attributed to the combined motion of the bending and librational modes.