Quantum rate dynamics for proton transfer reaction in a model system: Effect of the rate promoting vibrational mode

Hierarchical Equations of Motion for Quantum Chemical Dynamics: Recent Methodology Developments and Applications

ConspectusQuantum effects are critical to understanding many chemical dynamical processes in condensed phases, where interactions between molecules and their environment are usually strong and non-Markovian. In this Account, we review recent progress from our group in development and application of the hierarchical equations of motion (HEOM) method, highlighting its ability to address some challenging problems in quantum chemical dynamics.In the HEOM method, the bath degrees of freedom are represented using effective modes from exponential decomposition of the bath correlation function. Complex spectral densities and low temperature simulations often require a larger number of modes, making the simulations very expensive. Recent advances, such as the barycentric spectral decomposition (BSD) technique, can significantly reduce the number of effective modes, allowing to handle complex spectral densities and enabling simulations at very low temperatures, including near-zero temperature dynamics.Another key improvement in the computational efficiency is the use of tensor network methods like matrix product states and hierarchical tensor networks. These techniques allow for efficient HEOM propagation with thousands of effective modes, crucial for simulating large molecular systems interacting with multiple baths. This combination enables simulations of excitation energy transfer (EET) in systems like the Fenna-Matthews-Olson (FMO) complex and even larger systems with experimentally determined spectral densities.The versatility of the HEOM method is demonstrated through applications to a wide range of chemical dynamics problems. Simulations of EET and related ultrafast spectroscopy are first briefly covered. Applications of the HEOM to quantum tunneling effects in proton transfer reactions are then presented. Early works have studied the non-Kramers dependence of the rate constant as a function of bath friction due to deep tunneling and revealed vibrationally nonadiabatic dynamics within the so-called nontraditional view of proton transfer reactions. A recent work on the large kinetic isotope effects in soybean lipoxygenase also indicated that many quantum correction approximations to classical transition-state theory may fall short in describing deep tunneling effects.Charge transport and separation dynamics in organic semiconductors are another area where the HEOM method has been instrumental. We first demonstrate that the HEOM provides a unified description of both band-like and thermally assisted charge carrier transport in organic materials. The effect of non-nearest neighbor transitions is then investigated by combining generalized master equations with exact memory kernels. The HEOM method also enables simulation of charge separation in organic photovoltaics (OPVs) and reveals how factors such as external electric fields, entropy, and charge delocalization influence the charge separation barrier and dynamics.Moreover, HEOM has been applied to investigate hydrogen atom scattering on the Au(111) surface and vibrational energy relaxation at molecule-metal interfaces. These studies provide deeper insights into how electron-hole pair excitations and temporary charge transfer states influence the nuclear motion, offering a new framework for simulating nonadiabatic dynamics on metal surfaces.In summary, the HEOM method has developed into a robust tool for simulating quantum effects in condensed phases. Future developments in algorithm efficiency and computational power will likely expand its applicability to even more complex systems.

Quantum rates in dissipative systems with spatially varying friction

We investigate whether making the friction spatially dependent on the reaction coordinate introduces quantum effects into the thermal reaction rates for dissipative reactions. Quantum rates are calculated using the numerically exact multi-configuration time-dependent Hartree method, as well as the approximate ring-polymer molecular dynamics (RPMD), ring-polymer instanton methods, and classical molecular dynamics. By conducting simulations across a wide range of temperatures and friction strengths, we can identify the various regimes that govern the reactive dynamics. At high temperatures, in addition to the spatial-diffusion and energy-diffusion regimes predicted by Kramer's rate theory, a (coherent) tunneling-dominated regime is identified at low friction. At low temperatures, incoherent tunneling dominates most of Kramer's curve, except at very low friction, when coherent tunneling becomes dominant. Unlike in classical mechanics, the bath's influence changes the equilibrium time-independent properties of the system, leading to a complex interplay between spatially dependent friction and nuclear quantum effects even at high temperatures. More specifically, a realistic friction profile can lead to an increase (or decrease) of the quantum (classical) rates with friction within the spatial-diffusion regime, showing that classical and quantum rates display qualitatively different behaviors. Except at very low frictions, we find that RPMD captures most of the quantum effects in the thermal reaction rates.

Resonance Enhancement of Vibrational Polariton Chemistry Obtained from the Mixed Quantum-Classical Dynamics Simulations

We applied a variety of mixed quantum-classical (MQC) approaches to simulate the VSC-influenced reaction rate constant. All of these MQC simulations treat the key vibrational levels associated with the reaction coordinate in the quantum subsystem (as quantum states), whereas all other degrees of freedom (DOFs) are treated inside the classical subsystem. We find that, as long as we have the quantum state descriptions for the vibrational DOFs, one can correctly describe the VSC resonance condition when the cavity frequency matches the bond vibrational frequency. This correct resonance behavior can be obtained regardless of the detailed MQC methods that one uses. The results suggest that the MQC approaches can generate semiquantitative agreement with the exact results for rate constant changes when changing the cavity frequency, the light-matter coupling strength, or the cavity lifetime. The finding of this work suggests that one can use computationally economic MQC approaches to explore the collective coupling scenario when many molecules are collectively coupled to many cavity modes in the future.

Quantum dynamical effects of vibrational strong coupling in chemical reactivity

Recent experiments suggest that ground state chemical reactivity can be modified when placing molecular systems inside infrared cavities where molecular vibrations are strongly coupled to electromagnetic radiation. This phenomenon lacks a firm theoretical explanation. Here, we employ an exact quantum dynamics approach to investigate a model of cavity-modified chemical reactions in the condensed phase. The model contains the coupling of the reaction coordinate to a generic solvent, cavity coupling to either the reaction coordinate or a non-reactive mode, and the coupling of the cavity to lossy modes. Thus, many of the most important features needed for realistic modeling of the cavity modification of chemical reactions are included. We find that when a molecule is coupled to an optical cavity it is essential to treat the problem quantum mechanically to obtain a quantitative account of alterations to reactivity. We find sizable and sharp changes in the rate constant that are associated with quantum mechanical state splittings and resonances. The features that emerge from our simulations are closer to those observed in experiments than are previous calculations, even for realistically small values of coupling and cavity loss. This work highlights the importance of a fully quantum treatment of vibrational polariton chemistry.

Comparison of Matsubara dynamics with exact quantum dynamics for an oscillator coupled to a dissipative bath

We report the first numerical calculations in which converged Matsubara dynamics is compared directly with exact quantum dynamics with no artificial damping of the time-correlation functions (TCFs). The system treated is a Morse oscillator coupled to a harmonic bath. We show that, when the system-bath coupling is sufficiently strong, the Matsubara calculations can be converged by explicitly including up to M = 200 Matsubara modes, with the remaining modes included as a harmonic "tail" correction. The resulting Matsubara TCFs are in near-perfect agreement with the exact quantum TCFs, for non-linear as well as linear operators, at a temperature at which the TCFs are dominated by quantum thermal fluctuations. These results provide compelling evidence that incoherent classical dynamics can arise in the condensed phase at temperatures at which the statistics are dominated by quantum (Boltzmann) effects, as a result of smoothing of imaginary-time Feynman paths. The techniques developed here may also lead to efficient methods for benchmarking system-bath dynamics in the overdamped regime.

Application of the imaginary time hierarchical equations of motion method to calculate real time correlation functions

We investigate the application of the imaginary time hierarchical equations of motion method to calculate real time quantum correlation functions. By starting from the path integral expression for the correlated system–bath equilibrium state, we first derive a new set of equations that decouple the imaginary time propagation and the calculation of auxiliary density operators. The new equations, thus, greatly simplify the calculation of the equilibrium correlated initial state that is subsequently used in the real time propagation to obtain the quantum correlation functions. It is also shown that a periodic decomposition of the bath imaginary time correlation function is no longer necessary in the new equations such that different decomposition schemes can be explored. The applicability of the new method is demonstrated in several numerical examples, including the spin-Boson model, the Holstein model, and the double-well model for proton transfer reaction.

Nonequilibrium reaction rate theory: Formulation and implementation within the hierarchical equations of motion approach

The study of chemical reactions in environments under nonequilibrium conditions has been of interest recently in a variety of contexts, including current-induced reactions in molecular junctions and scanning tunneling microscopy experiments. In this work, we outline a fully quantum mechanical, numerically exact approach to describe chemical reaction rates in such nonequilibrium situations. The approach is based on an extension of the flux correlation function formalism to nonequilibrium conditions and uses a mixed real and imaginary time hierarchical equations of motion approach for the calculation of rate constants. As a specific example, we investigate current-induced intramolecular proton transfer reactions in a molecular junction for different applied bias voltages and molecule-lead coupling strengths.

Theory of vibrational polariton chemistry in the collective coupling regime

We theoretically demonstrate that the chemical reaction rate constant can be significantly suppressed by coupling molecular vibrations with an optical cavity, exhibiting both the collective coupling effect and the cavity frequency modification of the rate constant. When a reaction coordinate is strongly coupled to the solvent molecules, the reaction rate constant is reduced due to the dynamical caging effect. We demonstrate that collectively coupling the solvent to the cavity can further enhance this dynamical caging effect, leading to additional suppression of the chemical kinetics. This effect is further amplified when cavity loss is considered.

Comparative studies of IR spectra of deprotonated serine with classical and thermostated ring polymer molecular dynamics simulations

Here we report the vibrational spectra of deprotonated serine calculated from the classical molecular dynamics (MD) simulations and thermostated ring-polymer molecular dynamics (TRPMD) simulation with third-order density-functional tight-binding. In our earlier study [Inakollu and Yu, "A systematic benchmarking of computational vibrational spectroscopy with DFTB3: Normal mode analysis and fast Fourier transform dipole autocorrelation function," J. Comput. Chem. 39, 2067 (2018)] of deprotonated serine, we observed a significant difference in the vibrational spectra with the classical MD simulations compared to the infrared multiple photon dissociation spectra. It was postulated that this is due to neglecting the nuclear quantum effects (NQEs). In this work, NQEs are considered in spectral calculation using the TRPMD simulations. With the help of potential of mean force calculations, the conformational space of deprotonated serine is analyzed and used to understand the difference in the spectra of classical MD and TRPMD simulations at 298.15 and 100 K. The high-frequency vibrational bands in the spectra are characterized using Fourier transform localized vibrational mode (FT-νN AC) and interatomic distance histograms. At room temperature, the quantum effects are less significant, and the free energy profiles in the classical MD and the TRPMD simulations are very similar. However, the hydrogen bond between the hydroxyl-carboxyl bond is slightly stronger in TRPMD simulations. At 100 K, the quantum effects are more prominent, especially in the 2600-3600 cm-1, and the free energy profile slightly differs between the classical MD and TRPMD simulations. Using the FT-νN AC and the interatomic distance histograms, the high-frequency vibrational bands are discussed in detail.

Understanding the Large Kinetic Isotope Effect of Hydrogen Tunneling in Condensed Phases by Using Double-Well Model Systems

In recent years, many experiments have shown large kinetic isotope effects (KIEs) for hydrogen transfer reactions in condensed phases as evidence of strong quantum tunneling effects. Since accurate calculation of the tunneling dynamics in such systems still present significant challenges, previous studies have employed different types of approximations to estimate the tunneling effects and KIEs. In this work, by employing model systems consisting of a double-well coupled to a harmonic bath, we calculate the tunneling effects and KIEs using the numerically exact hierarchical equations of motion (HEOM) method. It is found that hydrogen and deuterium transfer reactions in the same system may show rather different behaviors, where hydrogen transfer is dominated by tunneling between the two lowest vibrational states and deuterium transfer is controlled by excited vibrational states close to the barrier top. The simulation results are also used to test the validity of various approximate methods. It is shown that the Wolynes theory of dissipative tunneling gives a good estimation of rate constants in the over-the-barrier regime, while the nonadiabatic reaction rate theory based on the Landau-Zener formula is more suitable for deep tunneling reactions.

Efficient propagation of the hierarchical equations of motion using the Tucker and hierarchical Tucker tensors

We develop new methods to efficiently propagate the hierarchical equations of motion (HEOM) by using the Tucker and hierarchical Tucker (HT) tensors to represent the reduced density operator and auxiliary density operators. We first show that by employing the split operator method, the specific structure of the HEOM allows a simple propagation scheme using the Tucker tensor. When the number of effective modes in the HEOM increases and the Tucker representation becomes intractable, the split operator method is extended to the binary tree structure of the HT representation. It is found that to update the binary tree nodes related to a specific effective mode, we only need to propagate a short matrix product state constructed from these nodes. Numerical results show that by further employing the mode combination technique commonly used in the multi-configuration time-dependent Hartree approaches, the binary tree representation can be applied to study excitation energy transfer dynamics in a fairly large system including over 10⁴ effective modes. The new methods may thus provide a promising tool in simulating quantum dynamics in condensed phases.

Modeling and Simulating the Excited-State Dynamics of a System with Condensed Phases: A Machine Learning Approach

Simulating the irreversible quantum dynamics of exciton- and electron-transfer problems poses a nontrivial challenge. Because the irreversibility of the system dynamics is a result of quantum thermal activation and dissipation caused by the surrounding environment, it is necessary to include infinite environmental degrees of freedom in the simulation. Because the capabilities of full quantum dynamics simulations that include the surrounding molecular degrees of freedom are limited, employing a system-bath model is a practical approach. In such a model, the dynamics of excitons or electrons are described by a system Hamiltonian, while the other degrees of freedom that arise from the environmental molecules are described by a harmonic oscillator bath (HOB) and system-bath interaction parameters. By extending on a previous study of molecular liquids [ J. Chem. Theory Comput. 2020, 16, 2099], here, we construct a system-bath model for exciton- and electron-transfer problems by means of a machine learning approach. We determine both the system and system-bath interaction parameters, including the spectral distribution of the bath, using the electronic excitation energies obtained from a quantum mechanics/molecular mechanics (QM/MM) simulation that is conducted as a function of time. Using the analytical expressions of optical response functions, we calculate linear and two-dimensional electronic spectra (2DES) for indocarbocyanine dimers in methanol. From these results, we demonstrate the capability of our approach to elucidate the nonequilibrium exciton dynamics of a quantum system in a nonintuitive manner.

Numerically "exact" approach to open quantum dynamics: The hierarchical equations of motion (HEOM)

An open quantum system refers to a system that is further coupled to a bath system consisting of surrounding radiation fields, atoms, molecules, or proteins. The bath system is typically modeled by an infinite number of harmonic oscillators. This system-bath model can describe the time-irreversible dynamics through which the system evolves toward a thermal equilibrium state at finite temperature. In nuclear magnetic resonance and atomic spectroscopy, dynamics can be studied easily by using simple quantum master equations under the assumption that the system-bath interaction is weak (perturbative approximation) and the bath fluctuations are very fast (Markovian approximation). However, such approximations cannot be applied in chemical physics and biochemical physics problems, where environmental materials are complex and strongly coupled with environments. The hierarchical equations of motion (HEOM) can describe the numerically "exact" dynamics of a reduced system under nonperturbative and non-Markovian system-bath interactions, which has been verified on the basis of exact analytical solutions (non-Markovian tests) with any desired numerical accuracy. The HEOM theory has been used to treat systems of practical interest, in particular, to account for various linear and nonlinear spectra in molecular and solid state materials, to evaluate charge and exciton transfer rates in biological systems, to simulate resonant tunneling and quantum ratchet processes in nanodevices, and to explore quantum entanglement states in quantum information theories. This article presents an overview of the HEOM theory, focusing on its theoretical background and applications, to help further the development of the study of open quantum dynamics.

Emergence of molecular chirality due to chiral interactions in a biological environment

We explore the interplay between tunneling process and chiral interactions in the discrimination of chiral states for an ensemble of molecules in a biological environment. Each molecule is described by an asymmetric double-well potential and the environment is modeled as a bath of harmonic oscillators. We carefully analyze different time-scales appearing in the resulting master equation at both weak- and strong-coupling limits. The corresponding results are accompanied by a set of coupled differential equations characterizing optical activity of the molecules. We show that, at the weak-coupling limit, chiral interactions prohibit the coherent racemization induced by decoherence effects and thus preserve the initial chiral state. At the strong-coupling limit, considering the memory effects of the environment, Markovian behavior is observed at long times.