Reactive Flux Correlation Function Research Articles

Transition rate theory, spectral analysis, and reactive paths.

The kinetics of a dynamical system dominated by two metastable states is examined from the perspective of the activated-dynamics reactive flux formalism, Markov state eigenvalue spectral decomposition, and committor-based transition path theory. Analysis shows that the different theoretical formulations are consistent, clarifying the significance of the inherent microscopic lag-times that are implicated, and that the most meaningful one-dimensional reaction coordinate in the region of the transition state is along the gradient of the committor in the multidimensional subspace of collective variables. It is shown that the familiar reactive flux activated dynamics formalism provides an effective route to calculate the transition rate in the case of a narrow sharp barrier but much less so in the case of a broad flat barrier. In this case, the standard reactive flux correlation function decays very slowly to the plateau value that corresponds to the transmission coefficient. Treating the committor function as a reaction coordinate does not alleviate all issues caused by the slow relaxation of the reactive flux correlation function. A more efficient activated dynamics simulation algorithm may be achieved from a modified reactive flux weighted by the committor. Simulation results on simple systems are used to illustrate the various conceptual points.

Hydrogen Bonds between Ions of Opposite and Like Charge in Hydroxyl-Functionalized Ionic Liquids: an Exhaustive Examination of the Interplay between Global and Local Motions and Intermolecular Hydrogen Bond Lifetimes and Kinetics.

Hydroxyl-functionalized ionic liquids (ILs) represent a new interesting class of ILs where hydrogen bonds (HBs) play an important role: here, "typical" HBs between cations and anions (ca) are competing with "atypical" HBs connecting pairs of cations (cc). We study the equilibrium and kinetics of (cc) and (ca) HBs in 1-(n-hydroxyalkyl)-pyridinium bis(trifluoromethlysulfonyl)imide [HOCnPy][NTf2] ILs by means of molecular dynamics simulations. (cc) HBs are found to be between 0.96 and 3.76 kJ mol-1 stronger than their (ca) counterparts, depending on the alkyl chain length. HB lifetimes and kinetics are analyzed by means of HB population and reactive flux correlation functions. Essentially, four different HB lifetimes have to be considered, spanning about 3 orders of magnitude, each valid in its own right and each associated with different aspects of HB breaking and HB reformation. The long-time limiting behavior of the HB population correlation function is controlled by diffusion of the ions and can be quantitatively described by analytical expressions. The short-time HB behavior is tied to the localized dynamics of the hydroxyl group exploring its local solvation environment. A minimalist kinetic two-domain model is introduced to realistically describe the time evolution of the HB population correlation function for both (ca) and (cc) HBs over 5 orders of magnitude. By employing the reactive flux method, we determine the kinetics of HB breaking, unaffected by diffusion processes. We determine both, the ultrafast upper boundary and the average rate of HB breaking, allowing recrossing-events during the transient relaxation time period. For sufficiently long alkyl chains, all those computed HB lifetimes indicate a higher kinetic stability of (cc) HBs over (ca) HBs; for short chains, it is vice-versa.

Dependence of the Rate of LiF Ion-Pairing on the Description of Molecular Interaction.

We present an analysis of the dynamics of ion-pairing of lithium fluoride (LiF) in aqueous solvent using both detailed molecular simulation as well as reduced models within a generalized Langevin equation (GLE) framework. We explored the sensitivity of the ion-pairing phenomena to the details of descriptions of molecular interaction, comparing two empirical potentials to explicit quantum based density functional theory. We find quantitative differences in the potentials of mean force for ion-pairing as well as time dependent frictions that lead to variations in the rate constant and reactive flux correlation functions. These details reflect differences in solvent response to ion-pairing between different representations of molecular interaction and influence anharmonicity of the dynamic response. We find that the short-time anharmonic response is recovered with a GLE parametrization. Recovery of the details of long time response may require extensions to the reduced model. We show that the utility of using a reduced model leads to a straightforward application of variational transition state theory concepts to the condensed phase system. The significance of this is reflected in the analysis of committor distributions and the variation of planar hypersurfaces, leading to an improved understanding of factors that determine the rate of LiF ion-pairing.

Insights on Hydrogen-Bond Lifetimes in Liquid and Supercooled Water

We study the temperature dependence of the lifetime of geometric and geometric/energetic water hydrogen-bonds (H-bonds), down to supercooled water, through molecular dynamics. The probability and lifetime of H-bonds that break either by translational or librational motions and those of energetic broken H-bonds, along with the effects of transient broken H-bonds and transient H-bonds, are considered. We show that the fraction of transiently broken energetic H-bonds increases at low temperatures and that this energetic breakdown is caused by oxygen-oxygen electrostatic repulsions upon too small amplitude librations to disrupt geometric H-bonds. Hence, differences between geometric and energetic continuous H-bond lifetimes are associated with large H-bond energy fluctuations, in opposition to moderate geometric fluctuations, within common energetic and geometric H-bond definition thresholds. Exclusion of transient broken H-bonds and transient H-bonds leads to H-bond definition-independent mean lifetimes and activation energies, ~11 kJ/mol, consistent with the reactive flux method and experimental scattering results. Further, we show that power law decay of specific temporal H-bond lifetime probability distributions is associated with librational and translational motions that occur on the time scale (~0.1 ps) of H-bond breaking /re-forming dynamics. While our analysis is diffusion-free, the effect of diffusion on H-bond probability distributions where H-bonds are allowed to break and re-form, switching acceptors in between, is shown to result in neither exponential nor power law decay, similar to the reactive flux correlation function.

Isomerization dynamics of a buckled nanobeam

We analyze the dynamics of a model of a nanobeam under compression. The model is a two-mode truncation of the Euler-Bernoulli beam equation subject to compressive stress applied at both ends. We consider parameter regimes where the first mode is unstable and the second mode can be either stable or unstable, and the remaining modes (neglected) are always stable. Material parameters used correspond to a silicon nanobeam. The two-mode model Hamiltonian is the sum of a (diagonal) kinetic energy term and a potential energy term. The form of the potential energy function suggests an analogy with isomerization reactions in chemistry, where "isomerization" here corresponds to a transition between two stable beam configurations. We therefore study the dynamics of the buckled beam using the conceptual framework established for the theory of isomerization reactions. When the second mode is stable the potential energy surface has an index one saddle, and when the second mode is unstable the potential energy surface has an index two saddle and two index one saddles. Symmetry of the system allows us to readily construct a phase space dividing surface between the two "isomers" (buckled states); we rigorously prove that, in a specific energy range, it is a normally hyperbolic invariant manifold. The energy range is sufficiently wide that we can treat the effects of the index one and index two saddles on the isomerization dynamics in a unified fashion. We have computed reactive fluxes, mean gap times, and reactant phase space volumes for three stress values at several different energies. In all cases the phase space volume swept out by isomerizing trajectories is considerably less than the reactant density of states, proving that the dynamics is highly nonergodic. The associated gap time distributions consist of one or more "pulses" of trajectories. Computation of the reactive flux correlation function shows no sign of a plateau region; rather, the flux exhibits oscillatory decay, indicating that, for the two-mode model in the physical regime considered, a rate constant for isomerization does not exist.

Transition-State Theory, Dynamics, and Narrow Time Scale Separation in the Rate-Promoting Vibrations Model of Enzyme Catalysis.

The power of transition-state theory (TST) for understanding enzymes is evidenced by its recent use in the design and synthesis of highly active de novo enzymes. However, dynamics can influence reaction kinetics, and some studies of rate-promoting vibrations even claim that dynamical theories instead of TST are needed to understand enzymatic reaction mechanisms. For the rate-promoting vibration (RPV) model of enzyme catalysis [Antoniou et al., J. Chem. Phys. 2004, 121, 6442], a reactive flux correlation function analysis shows that dynamical effects do slow the kinetics. However, the RPV model also shows extremely long-lived correlations because the RPV and the bath are not directly coupled. Additionally, earlier studies of the RPV model show a narrow time scale separation due to a small 5kT barrier. Thus earlier findings based on the RPV model may have little bearing on the properties of real enzymes. The intrinsic reaction coordinate (IRC) reveals that the RPV is an important component of the reaction coordinate at early and late stages of the pathway, but the RPV is not an important component of the reaction coordinate direction at the transition state. The unstable eigenmode from harmonic TST (which coincides with the IRC at the saddle point) gives a larger transmission coefficient than the coordinate used in the correlation functions of Antoniou et al. Thus while TST cannot predict the transmission coefficient, the RPV model suggests that TST can provide mechanistic insights on elementary steps in enzyme catalysis. Finally, we propose a method for using the transition-state ensemble as predicted from harmonic TST to distinguish promoting vibrations from other more mundane bath variables.

On the Effects of the External Surface on the Equilibrium Transport in Zeolite Crystals

With the aid of molecular simulation techniques (molecular dynamics, grand-canonical Monte Carlo, and reactive flux correlation function RFCF), the influence of the external surface on the equilibrium permeation of methane and ethane into and out of an AFI-type zeolite crystal has been studied. In particular, “extended dynamically corrected transition state theory”, which has been proven to describe the transport of tracers in periodic crystals correctly, has been applied to surface problems. The results suggest that the molecules follow paths that are close to the pore wall in the interior and also at the crystal surface. Moreover, the recrossing rate at the surface turns out to be non-negligible, yet, in contrast to the intracrystalline recrossing rate, remains almost constant over loading which gives indication to diffusive barrier crossing at the crystal surface. As a consequence of very different adsorption and desorption barriers, the corresponding permeabilities are shown to be not equal for one and the same condition (T and p). The critical crystal length, beyond which surface effects can be certainly neglected, is computed on basis of flux densities. Entrance/exit effects, in the present cases, are practically important solely for ethane at low pressures. The influence of the type of external surface on the surface flux is, hereby, rather small, because the transport at the surface is controlled by the slow supply from the gas phase. This has been evidenced by a simplified thermodynamic model that has been derived within this work and which is based on rapidly assessable simulation data. Finally, we propose a procedure for estimating the importance of different factors that have an impact on surface effects.

Quantum-classical Liouville dynamics of proton and deuteron transfer rates in a solvated hydrogen-bonded complex

Proton and deuteron transfer reactions in a hydrogen-bonded complex dissolved in a polar solution are studied using quantum-classical Liouville dynamics. Reactive-flux correlation functions that involve quantum-classical Liouville dynamics for species operators and quantum equilibrium sampling are used to calculate the rate constants. Adiabatic and nonadiabatic reaction rates are computed, compared, and analyzed. Large variations of the kinetic isotope effect (KIE) for this reaction have been observed in the literature, which depend on the nature of the approximate calculation used to estimate the proton and deuteron transfer rates. Our estimate of the KIE lies at the low end of the range of previously observed values, suggesting a rather small KIE for this reaction.

Quantum reaction rates and sampling of quantum equilibrium structure

The time dependent quantum reaction rate may be expressed in terms of the Kubo-transformed reactive flux correlation function involving the Kubo-transformed quantum spectral density function W¯A(X,0). The rate constant, which is obtained from the long time limit of this correlation function, can also be expressed in terms of a simple quantum canonical correlation function. In this context, two useful alternative expressions, WA(X,iℏβ/2) and ReWA(X,0), to the difficult-to-compute W¯A(X,0) function are compared. The spectral density WA(X,iℏβ/2) is a real and symmetric approximation to W¯A(X,0) while ReWA(X,0) is a much simpler form obtained from the quantum canonical correlation function. Computations on a simple parabolic barrier crossing problem and a two-level system reactive system confirm that the long time reaction rate constant is the same for all spectral density expressions. Our results show that the density function WA(X,iℏβ/2) may be more suitable for simple model systems while ReWA(X,0) may be more useful for complicated systems involving many degrees of freedom.

Quantum equilibrium structure for strong nonadiabatic coupling: Reaction rate enhancement

The quantum reactive flux correlation function is computed for a two-level system using an expression for the quantum equilibrium structure appropriate for strong nonadiabatic coupling, in conjunction with quantum–classical Liouville dynamics. The magnitude of the quantum mechanical enhancement of the reaction rate as a result of strong nonadiabatic coupling is studied. The reaction rate is found to increase strongly with an increase in the nonadiabatic coupling strength as well as with a decrease in the temperature. Equilibrium quantum effects increase the ground-state contribution to the rate constant but these effects decrease the excited-state contribution.

Calculation of reactive flux correlation functions for systems in a condensed phase environment: A multilayer multiconfiguration time-dependent Hartree approach

A numerically exact quantum mechanical approach is proposed to evaluate thermal rate constants for systems in a model condensed phase environment. Employing the reactive flux correlation function formalism, the approach efficiently combines the multilayer multiconfiguration time-dependent Hartree theory with an importance sampling scheme for thermal distribution of the initial states. The performance of the method is illustrated by applications to two models of condensed phase dynamics: the donor-acceptor electron transfer model also known as the spin-boson model and a model for proton transfer reactions in the condensed phase.

Analysis of kinetic isotope effects for nonadiabatic reactions

Factors influencing the rates of quantum mechanical particle transfer reactions in many-body systems are discussed. The investigations are carried out on a simple model for a proton transfer reaction that captures generic features seen in more realistic models of condensed phase systems. The model involves a bistable quantum oscillator coupled to a one-dimensional double-well reaction coordinate, which is in turn coupled to a bath of harmonic oscillators. Reactive-flux correlation functions that involve quantum-classical Liouville dynamics for chemical species operators and quantum equilibrium sampling are used to estimate the reaction rates. Approximate analytical expressions for the quantum equilibrium structure are derived. Reaction rates are shown to be influenced significantly by both the quantum equilibrium structure and nonadiabatic dynamics. Nonadiabatic dynamical effects are found to play the major role in determining the magnitude of the kinetic isotope effect for the model transfer reaction.

Path-integral virial estimator for reaction-rate calculation based on the quantum instanton approximation

The quantum instanton approximation is a type of quantum transition-state theory that calculates the chemical reaction rate using the reactive flux correlation function and its low-order derivatives at time zero. Here we present several path-integral estimators for the latter quantities, which characterize the initial decay profile of the flux correlation function. As with the internal energy or heat-capacity calculation, different estimators yield different variances (and therefore different convergence properties) in a Monte Carlo calculation. Here we obtain a virial (-type) estimator by using a coordinate scaling procedure rather than integration by parts, which allows more computational benefits. We also consider two different methods for treating the flux operator, i.e., local-path and global-path approaches, in which the latter achieves a smaller variance at the cost of using second-order potential derivatives. Numerical tests are performed for a one-dimensional Eckart barrier and a model proton transfer reaction in a polar solvent, which illustrates the reduced variance of the virial estimator over the corresponding thermodynamic estimator.

A molecular-dynamics study of a model SN1 dissociation reaction at the water liquid/vapor interface

The thermodynamics and dynamics of a model S(N)1 reaction: t-BuCl --> t-Bu+ + Cl- is studied at the water liquid/vapor interface using molecular-dynamics computer simulations. The empirical valence bond approach is used to couple two diabatic states, covalent and ionic, in the electronically adiabatic limit. Umbrella sampling calculations are used to calculate the potential of mean force along the reaction coordinate (defined as the t-Bu to Cl distance) in bulk water and in several locations at the interface. We find a significant increase of the dissociation barrier height and of the reaction free energy at the interface relative to the bulk. This is shown to be due to the reduced polarity of the interface. Reactive flux correlation function calculations show significant deviation of the rate constant from the transition-state theory: The transmission coefficients range from 0.49 in the bulk to 0.05 above the Gibbs surface. The low transmission coefficient at the interface despite the lower friction is shown to be due to slow vibrational relaxation.

Nonadiabatic chemical reactions

A mixed quantum-classical approach where the environment is treated classically and the reactive degrees of freedom are considered to be quantum mechanical can be used to describe many chemical reactions, such as proton and electron transfer processes. We present reactive flux correlation function expressions for the rate constants of nonadiabatic chemical reactions occurring in quantum-classical systems. By means of a two-state model coupled to a classical bath, we illustrate an efficient method of computation, based on a sequential short step propagation, which is applicable to quantum systems interacting with general classical condensed phase environments.

Quantum-classical dynamics of nonadiabatic chemical reactions

A reactive flux correlation function formalism for the calculation of rate constants for mixed quantum-classical systems undergoing nonadiabatic dynamics is presented. The linear response formalism accounts for the stationarity of the equilibrium density under quantum-classical dynamics and expresses the rate constant in terms of an ensemble of surface-hopping trajectories. Calculations are carried out on a model two-level system coupled to a nonlinear oscillator which is in turn coupled to a harmonic heat bath. Relevant microscopic species variables for this system include two stable states, corresponding to the ground state adiabatic surface, as well as another species corresponding to the excited state surface. The time-dependent rate constants for the model are evaluated in the adiabatic limit, where the dynamics is confined to the ground Born–Oppenheimer surface, and these results are compared with calculations that account for nonadiabatic transitions among the system states.

Combining semiclassical time evolution and quantum Boltzmann operator to evaluate reactive flux correlation function for thermal rate constants of complex systems

The semiclassical (SC) initial value representation (IVR) provides a way for including quantum effects into classical molecular dynamics simulations. Implementation of the SC-IVR to the thermal rate constant calculation, based on the reactive flux correlation function formalism, has two major obstacles: (1) the SC integrand may be highly oscillatory with respect to the initial phase space variables; and (2) matrix elements of the Boltzmannized flux operator, which are crucial in generating the initial (or final) distribution for the SC trajectories, are generally not available in analytic forms. In this paper, we present practical ways of overcoming these two barriers for the SC calculation of thermal rate constants. For the first problem, we show that use of a symmetric flux–flux correlation function, together with the generalized Filinov transformation technique, can significantly smooth the corresponding SC integrand and make the calculation practical for quite large systems. For the second problem, we propose a general method for evaluating matrix elements of the Boltzmannized flux operator “on-the-fly,” based on the combination of the imaginary-time path integral technique with the Metropolis random walk algorithm. Using these approaches, it is shown that thermal rate constants can be obtained for systems with more than 100 degrees of freedom, as well as for reactions in the deep tunneling regimes where quantum effects are significant.

Resolving the hydrogen bond dynamics conundrum

This paper analyzes dynamic properties of hydrogen bonds in liquid water. We use molecular dynamics simulation to calculate different probability densities that govern the time evolution of the formation and rupture of hydrogen bonds. We provide analytical connections between these functions. Excellent agreement with our simulation results is observed. We prove transition state theory rate constant to be identical to the inverse of the associated mean first passage time (hydrogen bond lifetime). Hence, the analysis establishes its Arrhenius temperature dependence. We give the explicit relation between reactive flux correlation function for the relaxation dynamics of hydrogen bonds, and their first passage time probability densities. All the different observations in the existing literature, associated with various estimates of hydrogen bonding times in liquid water that are affected (or not affected) by particular bond criteria, as well as by different definitions of hydrogen bond lifetimes applied in simulation, can be easily reconciled within the framework of reactive flux correlation function approach.

Forward–backward initial value representation for the calculation of thermal rate constants for reactions in complex molecular systems

The semiclassical (SC) initial value representation (IVR) provides a potentially practical way for including quantum effects into classical molecular dynamics simulations. The forward–backward (FB) version of the IVR provides an especially attractive way for calculating time correlation functions, in particular the reactive flux correlation function which determines chemical reaction rates. This paper presents a further analysis and development of the FB-IVR approach. Applications show that it is feasible and accurate for a reaction coordinate coupled to up to 40 degrees of freedom.

Generalization of the Linearized Approximation to the Semiclassical Initial Value Representation for Reactive Flux Correlation Functions

The semiclassical (SC) initial value representation (IVR) provides a general and practical approach for including quantum effects in classical molecular dynamics simulations. The linearized approximation (LA) to the SC-IVR simplifies the description much further, reducing it to the well-known classical Wigner model (i.e., a classical trajectory calculation with a Wigner distribution of initial conditions); the LA is able to describe quantum effects well for short times (t ≲ ℏβ) but not so for longer times. It is shown here how the full SC-IVR approach, which is able to describe quantum effects for long times, can be cast in a form very similar in structure to the LA, with specific application to flux correlation functions relevant to chemical reaction rates. This formulation may thus make it possible to carry out full SC-IVR calculations while still retaining much of the simplifying aspects of its linearized approximation.