A temperature control technique for nonequilibrium molecular simulation

We describe a dynamical approach to thermal regulation in molecular dynamics. Temperature is moderated by a control law and an additional variable, as in Nose dynamics, but whose influence on the system decreases as the system approaches equilibrium. This device enables approximation of microcanonical averages and autocorrelation functions consistent with a given target temperature. Moreover, we demonstrate that the suggested technique is effective for the control of heat dissipation in a nonequilibrium setting, first by showing that the temperature control correctly regulates heat introduced by a rapid change to the system, and then by studying the slow relaxation of vibrational degrees of freedom (e.g., due to bonded atoms) in a solvent bath.     

Reaction path potential for complex systems derived from combined ab initio quantum mechanical and molecular mechanical calculations

Combined ab initio quantum mechanical and molecular mechanical calculations have been widely used for modeling chemical reactions in complex systems such as enzymes, with most applications being based on the determination of a minimum energy path connecting the reactant through the transition state to the product in the enzyme environment. However, statistical mechanics sampling and reaction dynamics calculations with a combined ab initio quantum mechanical (QM) and molecular mechanical (MM) potential are still not feasible because of the computational costs associated mainly with the ab initio quantum mechanical calculations for the QM subsystem. To address this issue, a reaction path potential energy surface is developed here for statistical mechanics and dynamics simulation of chemical reactions in enzymes and other complex systems. The reaction path potential follows the ideas from the reaction path Hamiltonian of Miller, Handy and Adams for gas phase chemical reactions but is designed specifically for large systems that are described with combined ab initio quantum mechanical and molecular mechanical methods. The reaction path potential is an analytical energy expression of the combined quantum mechanical and molecular mechanical potential energy along the minimum energy path. An expansion around the minimum energy path is made in both the nuclear and the electronic degrees of freedom for the QM subsystem internal energy, while the energy of the subsystem described with MM remains unchanged from that in the combined quantum mechanical and molecular mechanical expression and the electrostatic interaction between the QM and MM subsystems is described as the interaction of the MM charges with the QM charges. The QM charges are polarizable in response to the changes in both the MM and the QM degrees of freedom through a new response kernel developed in the present work. The input data for constructing the reaction path potential are energies, vibrational frequencies, and electron density response properties of the QM subsystem along the minimum energy path, all of which can be obtained from the combined quantum mechanical and molecular mechanical calculations. Once constructed, it costs much less for its evaluation. Thus, the reaction path potential provides a potential energy surface for rigorous statistical mechanics and reaction dynamics calculations of complex systems. As an example, the method is applied to the statistical mechanical calculations for the potential of mean force of the chemical reaction in triosephosphate isomerase.     

The role of reagent internal excitation in collision experiments

Variations in the rates of endoergic reactions due to different reagent excitations at the same total energy are of a limited (positive or negative) range and reflect a dynamical bias. Not so for bulk experiments, where all non-selected degrees of freedom have a thermal distribution.     

Reversibility and diffusion in mandelythiamin decarboxylation. Searching dynamical effects in decarboxylation reactions

Decarboxylation of mandelylthiamin in aqueous solution is analyzed by means of quantum mechanics/molecular mechanics simulations including solvent effects. The free energy profile for the decarboxylation reaction was traced, assuming equilibrium solvation, while reaction trajectories allowed us to incorporate nonequilibrium effects due to the solvent degrees of freedom as well as to evaluate the rate of the diffusion process in competition with the backward reaction. Our calculations that reproduce the experimental rate constant show that decarboxylation takes place with a non-negligible free energy barrier for the backward reaction and that diffusion of carbon dioxide is very fast compared to the chemical step. According to these findings catalysts would not act by preventing the backward reaction.     

A model for energy transfer in inelastic molecular collisions applicable at steady state or non-steady state and for an arbitrary distribution of collision energies

A new model for energy exchange between translational and internal degrees of freedom in atom-molecule collisions has been developed. It is suitable for both steady state conditions (e.g., a large number of collisions with thermal kinetic energies) and non-steady state conditions with an arbitrary distribution of collision energies (e.g., single high-energy collisions). In particular, it does not require that the collision energies be characterized by a quasi-thermal distribution, but nevertheless it is capable of producing a Boltzmann distribution of internal energies with the correct internal temperature under quasi-thermal conditions. The energy exchange is described by a transfer probability density that depends on the initial relative kinetic energy, the internal energy of the molecule, and the amount of energy transferred. The probability density for collisions that lead to excitation is assumed to decrease exponentially with the amount of transferred energy. The probability density for de-excitation is obtained from microscopic reversibility. The model has been implemented in the ion trap simulation program ITSIM and coupled with an Rice-Rampsberger-Kassel-Marcus (RRKM) algorithm to describe the unimolecular dissociation of populations of ions. Monte Carlo simulations of collisional energy transfer are presented. The model is validated for non-steady state conditions and for steady state conditions, and the effect of the kinetic energy dependence of the collision cross-section on internal temperature is discussed. Applications of the model to the problem of chemical mass shifts in RF ion trap mass spectrometry are shown.     

Simulating quantum dynamics in classical environments

Methods for simulating the dynamics of composite systems, where part of the system is treated quantum mechanically and its environment is treated classically, are discussed. Such quantum–classical systems arise in many physical contexts where certain degrees of freedom have an essential quantum character while the other degrees of freedom to which they are coupled may be treated classically to a good approximation. The dynamics of these composite systems are governed by a quantum–classical Liouville equation for either the density matrix or the dynamical variables which are operators in the Hilbert space of the quantum subsystem and functions of the classical phase space variables of the classical environment. Solutions of the evolution equations may be formulated in terms of surface-hopping dynamics involving ensembles of trajectory segments interspersed with quantum transitions. The surface-hopping schemes incorporate quantum coherence and account for energy exchanges between the quantum and classical degrees of freedom. Various simulation algorithms are discussed and illustrated with calculations on simple spin-boson models but the methods described here are applicable to realistic many-body environments.     

From concepts to algorithms for the treatment of large amplitude internal motions and unimolecular reactions

The formalism of the large amplitude path Hamiltonian provides a theoretical framework for the study of dynamical problems ranging from anharmonic vibrations to unimolecular reactions. A hierarchy of models at different degrees of sophistication can be elaborated and some of them have been put into practical use through the development of the DiNa package. The simplest level requires just the characterization of all the stationary points encountered along a large amplitude path and their quadratic environments. Perturbation theory can next be used to analyze the role of anharmonicity in the vibrational modulation of physico-chemical observables for semirigid systems. An extension of the same approach to saddle points allows the computation of reliable reaction rates taking into account curvature and tunneling effects. Finally, full characterization of the harmonic valley surrounding the path allows study of vibrational modulation in flexible systems and energy flow between different degrees of freedom during chemical reactions.     

A local coherent-state approximation to system-bath quantum dynamics

A novel quantum method to deal with typical system-bath dynamical problems is introduced. Subsystem discrete variable representation and bath coherent-state sets are used to write down a multiconfigurational expansion of the wave function of the whole system. With the help of the Dirac-Frenkel variational principle, simple equations of motion--a kind of Schrodinger-Langevin equation for the subsystem coupled to (pseudo) classical equations for the bath--are derived. True dissipative dynamics at all times is obtained by coupling the bath to a secondary, classical Ohmic bath, which is modeled by adding a friction coefficient in the derived pseudoclassical bath equations. The resulting equations are then solved for a number of model problems, ranging from tunneling to vibrational relaxation dynamics. Comparison of the results with those of exact, multiconfiguration time-dependent Hartree calculations in systems with up to 80 bath oscillators shows that the proposed method can be very accurate and might be of help in studying realistic problems with very large baths. To this end, its linear scaling behavior with respect to the number of bath degrees of freedom is shown in practice with model calculations using tens of thousands of bath oscillators.     

Model system-bath Hamiltonian and nonadiabatic rate constants for proton-coupled electron transfer at electrode-solution interfaces

An extension of the Anderson-Newns-Schmickler model for electrochemical proton-coupled electron transfer (PCET) is presented. This model describes reactions in which electron transfer between a solute complex in solution and an electrode is coupled to proton transfer within the solute complex. The model Hamiltonian is derived in a basis of electron-proton vibronic states defined within a double adiabatic approximation for the electrons, transferring proton, and bath modes. The interaction term responsible for electronic transitions between the solute complex and the electrode depends on the proton donor-acceptor vibrational mode within the solute complex. This model Hamiltonian is used to derive the anodic and cathodic rate constants for nonadiabatic electrochemical PCET. The derivation is based on the master equations for the reduced density matrix of the electron-proton subsystem, which includes the electrons of the solute complex and the electrode, as well as the transferring proton. The rate constant expressions differ from analogous expressions for electrochemical electron transfer because of the summation over electron-proton vibronic states and the dependence of the couplings on the proton donor-acceptor vibrational motion. These differences lead to additional contributions to the total reorganization energy, an additional exponential temperature-dependent prefactor, and a temperature-dependent term in the effective activation energy that has different signs for the anodic and cathodic processes. This model can be generalized to describe both nonadiabatic and adiabatic electrochemical PCET reactions and provides the framework for the inclusion of additional effects, such as the breaking and forming of other chemical bonds.     

Product energy distribution for exothermic reactive collisions

A statistical-dynamic model for calculating product state distributions for reactive collisions is presented. For the vibrational-translational coupling Franck-Condon type approximations are adopted which lead to very good agreement with exact quantum calculations in the limit of the collinear arrangement. The rotational degrees of freedom are treated statistically. The energy distribution among vibrational translational and rotational degrees of freedom is discussed as a function of the masses and the attractive part of the potential. Extension from triatomic to polyatomic exchange reactions are considered. A detailed comparison with other simple models is made.     

Mixed Quantum Classical Reaction Rates based on the Phase Space Formulation of the Hierarchical Equations of Motion

In a previous work [J. Chem. Phys. 140 , 174105 (2014)], we have shown that a mixed quantum classical (MQC) rate theory can be derived to investigate the quantum tunneling effects in the proton transfer reactions. However, the method is based on the high temperature approximation of the hierarchical equation of motion (HEOM) with the Debye-Drude spectral density, and results in a multi-state Zusman type of equation. We now extend this theory to include quantum effects of the bath degrees of freedom. By writing the full HEOM into a multidimensional partial differential equation in phase space, we can define a new reaction coordinate, and the previous method can be generalized to the full quantum regime. The validity of the new method is demonstrated by using numerical examples, including the spin-Boson model, and the double well model for proton transfer reaction. The new method is found to resolve some key problems of the previous theory based on high temperature approximation, including possible numerical instability in long time simulation and wrong rate constant at low temperatures.     

Non-adiabatic molecular dynamics with quantum solvent effects

Three novel approaches extending quantum-classical non-adiabatic (NA) molecular dynamics (MD) to include quantum effects of solvent environments are described. In a standard NA-MD the solute subsystem is treated quantum mechanically, while the larger solvent part of a system is treated classically. The three novel approaches presented here are based on the Bohmian formulation of quantum mechanics, the stochastic Schrödinger equation for the evolution of open quantum systems and the quantized Hamilton dynamics generalization of classical mechanics. The approaches extend the standard NA-MD to incorporate the following quantum effects of the solvent. (1) Branching, i.e. the ability of solvent quantum wave packets to split and follow asymptotically diverging trajectories correlated with different quantum states of the solute. (2) Decoherence, i.e. loss of quantum interference within the solute subsystem induced by the diverging solvent trajectories. (3) Zero point energy that contributes to NA coupling and must be preserved during the energy exchange between solvent and solute degrees of freedom. The Bohmian quantum-classical mechanics, stochastic mean-field and quantized mean-field approximations incorporate the quantum solvent effects into the standard quantum-classical NA-MD in a straightforward and efficient way that can be easily applied to quantum dynamics of condensed phase chemical systems.     

Non-equilibrium spin-boson model: counting statistics and the heat exchange fluctuation theorem

We focus on the non-equilibrium two-bath spin-boson model, a toy model for examining quantum thermal transport in many-body open systems. Describing the dynamics within the noninteracting-blip approximation equations, applicable, e.g., in the strong system-bath coupling limit and/or at high temperatures, we derive expressions for the cumulant generating function in both the Markovian and non-Markovian limits by energy-resolving the quantum master equation of the subsystem. For a Markovian bath, we readily demonstrate the validity of a steady-state heat exchange fluctuation theorem. In the non-Markovian limit a "weaker" symmetry relation generally holds, a general outcome of microreversibility. We discuss the reduction of this symmetry relation to the universal steady-state fluctuation theorem. Using the cumulant generating function, an analytic expression for the heat current is obtained. Our results establish the validity of the steady-state heat exchange fluctuation theorem in quantum systems with strong system-bath interactions. From the practical point of view, this study provides tools for exploring transport characteristics of the two-bath spin-boson model, a prototype for a nonlinear thermal conductor.     

Communication: Decoherence in a nonequilibrium environment: an analytically solvable model

We describe an analytically solvable model of quantum decoherence in a nonequilibrium environment. The model considers the effect of a bath driven from equilibrium by, for example, an ultrafast excitation of a quantum chromophore. The nonequilibrium response of the environment is represented by a nonstationary random function corresponding to the fluctuating transition frequency between two quantum states coupled to the surroundings. The nonstationary random function is characterized by a Fourier series with the phase of each term starting initially with a definite value across the ensemble but undergoing random diffusion with time. The decay of the off-diagonal density matrix element is shown to depend significantly on the particular pattern of initial phases of the terms in the Fourier series, or equivalently, the initial phases of bath modes coupled to the quantum subsystem. This suggests the possibility of control of quantum decoherence by the detailed properties of an environment that is driven from thermal equilibrium.     

Quantum dynamical simulation of electron-transfer reactions in an anharmonic environment

The multilayer multiconfiguration time-dependent Hartree theory is applied to study the quantum dynamics of ultrafast electron-transfer reactions in a condensed-phase environment with anharmonic potential functions. Effects of anharmonicity for both the nuclear degrees of freedom of the environment and the intramolecular vibrational degrees of freedom are investigated. Whereas the former can in principle be mapped to a fictitious harmonic bath, the latter cannot be represented in this way and, thus, go beyond the commonly employed linear response approximation. Numerical examples are presented to illustrate these findings.     

Dynamically biased RRKM model of activated gas-surface reactivity: vibrational efficacy and rotation as a spectator in the dissociative chemisorption of CH4 on Pt(111)

The reactivity of CH(4) impinging on a Pt(111) surface was examined using a precursor-mediated microcanonical trapping model of dissociative chemisorption wherein the effects of rotational and vibrational energy could be explored. Dissociative sticking coefficients for a diverse range of non-equilibrium effusive beam, supersonic beam, and eigenstate-resolved experiments were simulated and an average relative discrepancy between theory and experiment of better than 50% was achieved by treating molecular rotations and translation parallel to the surface as spectator degrees of freedom, and introducing a dynamically-biased vibrational efficacy. The model parameters are {E(0) = 57.9 kJ mol(-1), s = 2, η(v) = 0.40} where E(0) is the apparent threshold energy for reaction, s is the number of surface oscillators participating in energy exchange within each gas-surface collision complex formed, and η(v) is the mean vibrational efficacy for reaction relative to normal translational energy which figures in the assembly of the active exchangeable energy which is available to surmount the activation barrier to dissociative chemisorption. GGA-DFT electronic structure calculations provided vibrational frequencies for the transition state for dissociative chemisorption. The asymmetry of the rotational state populations in supersonic and effusive molecular beam experiments allowed kinetic analysis to establish that taking rotation as a spectator degree of freedom is a good approximation. Surface phonons, rather than the incident molecules, are calculated to play the dominant role in supplying the energy required to overcome the activation barrier for dissociative chemisorption under the thermal equilibrium conditions relevant to high pressure catalysis. Over the temperature range 300 K ≤T≤ 1000 K, the thermal dissociative sticking coefficient is predicted to be well described by S(T) = S(0) exp(-E(a)/RT) where S(0) = 0.62 and E(a) = 62.6 kJ mol(-1).     

Natural coordinates for polyatomic reactions

The curvilinear natural coordinates (NRC) are considered to describe the polyatomic reaction dynamics. An invariant equation is found for the reaction path curve (RP), allowing its calculation in an arbitrary initial coordinate system. Apart from common properties, well known for triatomic reactions, the polyatomic NRC take explicit account of a smooth change of normal vibrational modes of a chemical system during its evolution along RP. This generates special terms in the kinetic energy operator that induce nonadiabatic vibrational transitions and do not vanish in a classical limit, no matter how small the curvature of the RP is. For the case when the amplitudes of vibrations perpendicular to the RP are small it is shown that, as a first approximation, the vibrational transitions can be treated independently of the rotational degrees of freedom. This conclusion is obtained as a generalization of Eckart theory stating the possibility of vibration—rotational separation for stable molecules.