Rate constants from the reaction path Hamiltonian. I. Reactive flux simulations for dynamically correct rates

As ab initio electronic structure calculations become more accurate, inherent sources of error in classical transition state theory such as barrier recrossing and tunneling may become major sources of error in calculating rate constants. This paper introduces a general method for diabatically constructing the transverse eigensystem of a reaction path Hamiltonian in systems with many degenerate transverse frequencies. The diabatically constructed reaction path Hamiltonian yields smoothly varying coupling constants that, in turn, facilitate reactive flux calculations. As an example we compute the dynamically corrected rate constant for the chair to boat interconversion of cyclohexane, a system with 48 degrees of freedom and a number of degenerate frequencies. The transmission coefficients obtained from the reactive flux simulations agree with previous results that have been calculated using an empirical potential. Furthermore, the calculated rate constants agree with experimental values. Comparison to variational transition state theory shows that, despite finding the true bottleneck along the reaction pathway, variational transition state theory only accounts for half of the rate constant reduction due to recrossing trajectories.     

A hierarchical transition state search algorithm

A hierarchical transition state search algorithm is developed and its implementation in the density functional theory program deMon2k is described. This search algorithm combines the double ended saddle interpolation method with local uphill trust region optimization. A new formalism for the incorporation of the distance constrain in the saddle interpolation method is derived. The similarities between the constrained optimizations in the local trust region method and the saddle interpolation are highlighted. The saddle interpolation and local uphill trust region optimizations are validated on a test set of 28 representative reactions. The hierarchical transition state search algorithm is applied to an intramolecular Diels-Alder reaction with several internal rotors, which makes automatic transition state search rather challenging. The obtained reaction mechanism is discussed in the context of the experimentally observed product distribution.     

Efficient exploration of reaction paths via a freezing string method

The ability to efficiently locate transition states is critically important to the widespread adoption of theoretical chemistry techniques for their ability to accurately predict kinetic constants. Existing surface walking techniques to locate such transition states typically require an extremely good initial guess that is often beyond human intuition to estimate. To alleviate this problem, automated techniques to locate transition state guesses have been created that take the known reactant and product endpoint structures as inputs. In this work, we present a simple method to build an approximate reaction path through a combination of interpolation and optimization. Starting from the known reactant and product structures, new nodes are interpolated inwards towards the transition state, partially optimized orthogonally to the reaction path, and then frozen before a new pair of nodes is added. The algorithm is stopped once the string ends connect. For the practical user, this method provides a quick and convenient way to generate transition state structure guesses. Tests on three reactions (cyclization of cis,cis-2,4-hexadiene, alanine dipeptide conformation transition, and ethylene dimerization in a Ni-exchanged zeolite) show that this "freezing string" method is an efficient way to identify complex transition states with significant cost savings over existing methods, particularly when high quality linear synchronous transit interpolation is employed.     

Rate constants from the reaction path Hamiltonian. II. Nonseparable semiclassical transition state theory

For proton transfer reactions, the tunneling contributions to the rates are often much larger than thermally activated rates at temperatures of interest. A number of separable tunneling corrections have been proposed that capture the dependence of tunneling rates on barrier height and imaginary frequency size. However, the effects of reaction pathway curvature and barrier anharmonicity are more difficult to quantify. The nonseparable semiclassical transition state theory (TST) of Hernandez and Miller [Chem. Phys. Lett. 214, 129 (1993)] accounts for curvature and barrier anharmonicity, but it requires prohibitively expensive cubic and quartic derivatives of the potential energy surface at the transition state. This paper shows how the reaction path Hamiltonian can be used to approximate the cubic and quartic derivatives used in nonseparable semiclassical transition state theory. This enables tunneling corrections that include curvature and barrier anharmonicity effects with just three frequency calculations as required by a conventional harmonic transition state theory calculation. The tunneling correction developed here is nonseparable, but can be expressed as a thermal average to enable efficient Monte Carlo calculations. For the proton exchange reaction NH2 + CH4 <==> NH3 + CH3, the nonseparable rates are very accurate at temperatures from 300 K up to about 1000 K where the TST rate itself begins to diverge from the experimental results.     

Proton transfer in a polar solvent from ring polymer reaction rate theory

We have used the ring polymer molecular dynamics method to study the Azzouz-Borgis model for proton transfer between phenol (AH) and trimethylamine (B) in liquid methyl chloride. When the A-H distance is used as the reaction coordinate, the ring polymer trajectories are found to exhibit multiple recrossings of the transition state dividing surface and to give a rate coefficient that is smaller than the quantum transition state theory value by an order of magnitude. This is to be expected on kinematic grounds for a heavy-light-heavy reaction when the light atom transfer coordinate is used as the reaction coordinate, and it clearly precludes the use of transition state theory with this reaction coordinate. As has been shown previously for this problem, a solvent polarization coordinate defined in terms of the expectation value of the proton transfer distance in the ground adiabatic quantum state provides a better reaction coordinate with less recrossing. These results are discussed in light of the wide body of earlier theoretical work on the Azzouz-Borgis model and the considerable range of previously reported values for its proton and deuteron transfer rate coefficients.     

A hybrid elastic band string algorithm for studies of enzymatic reactions

A common challenge in theoretical biophysics is the identification of a minimum energy path (MEP) for the rearrangement of a group of atoms from one stable configuration to another. The structure with maximum energy along the MEP approximates the transition state for the process and the energy profile itself permits estimation of the transition rates. In this work we describe a computationally efficient algorithm for the identification of minimum energy paths in complicated biosystems. The algorithm is a hybrid of the nudged elastic band (NEB) and string methods. It has been implemented in the pDynamo simulation program and tested by examining elementary steps in the reaction mechanisms of three enzymes: citrate synthase, RasGAP, and lactate dehydrogenase. Good agreement is found for the energies and geometries of the species along the reaction profiles calculated using the new algorithm and previous versions of the NEB and string techniques, and also those obtained by the common method of adiabatic exploration of the potential energy surface as a function of predefined reaction coordinates. Precisely refined structures of the saddle points along the paths may be subsequently obtained with the climbing image variant of the NEB algorithm. Directions in which the utility of the methods that we have implemented can be further improved are discussed.     

TSNet: predicting transition state structures with tensor field networks and transfer learning

Transition states are among the most important molecular structures in chemistry, critical to a variety of fields such as reaction kinetics, catalyst design, and the study of protein function. However, transition states are very unstable, typically only existing on the order of femtoseconds. The transient nature of these structures makes them incredibly difficult to study, thus chemists often turn to simulation. Unfortunately, computer simulation of transition states is also challenging, as they are first-order saddle points on highly dimensional mathematical surfaces. Locating these points is resource intensive and unreliable, resulting in methods which can take very long to converge. Machine learning, a relatively novel class of algorithm, has led to radical changes in several fields of computation, including computer vision and natural language processing due to its aptitude for highly accurate function approximation. While machine learning has been widely adopted throughout computational chemistry as a lightweight alternative to costly quantum mechanical calculations, little research has been pursued which utilizes machine learning for transition state structure optimization. In this paper TSNet is presented, a new end-to-end Siamese message-passing neural network based on tensor field networks shown to be capable of predicting transition state geometries. Also presented is a small dataset of S_N2 reactions which includes transition state structures – the first of its kind built specifically for machine learning. Finally, transfer learning, a low data remedial technique, is explored to understand the viability of pretraining TSNet on widely available chemical data may provide better starting points during training, faster convergence, and lower loss values. Aspects of the new dataset and model shall be discussed in detail, along with motivations and general outlook on the future of machine learning-based transition state prediction.

Transition states are among the most important molecular structures in chemistry, critical to a variety of fields such as reaction kinetics, catalyst design, and the study of protein function.     

Thermal decomposition of decalin: an ab initio study

Density functional theory calculations (B3LYP and BH&HLYP functionals) of the potential energy surface have been performed to investigate the mechanisms of decalin breakdown, and the Rice-Ramsperger-Kassel-Marcus and transition state theory methods have been used to compute the high-pressure limit thermal rate constants for the new reaction pathways. The new pathways connect decalin to five primary monoaromatic species: benzene, toluene, styrene, ethylbenzene, and xylene. The reactions used for the new routes are carbon-carbon bond cleavage reaction, dissociation reaction, and hydrogen abstraction and addition reactions. A kinetic analysis was performed for pyrolytic conditions, and benzene, toluene, and xylene were identified as major products.     

Model for a three-center decomposition reaction. I. Perfluorodiazirine

A semiempirical approach has been used to evaluate rate parameters for a three-center decomposition reaction from the point of view of transition state theory, with perfluorodiazirine serving as the prototype molecule. Several activated complex models are considered in which the reaction coordinate is chosen as the ? NCN bending mode. The constraints imposed include the principle of concerted bond-order conservation in passing from the initial to the final state, and use is made of empirical bond order–bond length and bond order–force constant relationships. The geometric configuration of the transition state sought is one which conforms with the lowest energy path and is also consistent with the observed entropy of activation. The potential energy of activation is taken as the optimum difference in binding energies (based on the INDO method) between the transition and initial states, and the critical energy is obtained by applying a correction for the zero-point energy difference, derived from normal coordinate analysis. Satisfactory agreement is found in the case of the activated complex model for which the total bond order is conserved and bonds undergoing rupture are assigned a fractional bond order (FBO) of 2/3, derived from the postulate (FBO) = α/β whe re α(=2) is the number of bonds breaking, and β(=3) is the number of bonds undergoing change in the ring opening.     

Dependence of the rate constants on the treatment of internal rotation modes: the reaction OH + CH3SH --> CH3S + H2O as an example

The title reaction has been used as an example to test the importance of using a hindered rotor treatment instead of a harmonic oscillator model for calculating vibrational partition functions corresponding to low-frequency internal rotation modes. First, a normal-mode analysis according to the Ayala and Schlegel's algorithm has been used to identify the internal rotation modes of methanethiol and the transition state structure. Then, after calculation of the energy barrier for each internal rotation, the corresponding hindered rotor partition functions have been calculated following the CW scheme of Chuang and Truhlar. The results show that the anharmonic treatment produces a rather modest improvement of the rate constants at room temperature or below.     

Geometric parameters of the transition state in radical reactions of antioxidants

An algorithm for calculating interatomic distances in a transition state of radical abstraction reactions through the enthalpy of a reaction is developed. The algorithm is based on a combination of quantum-chemical calculations with the calculation by the intersecting-parabolas method. Using this method and experimental data (enthalpies and activation energies of reactions), interatomic distances in a reaction center are calculated for the reactions of antioxidants with a symmetrical transition state of the N...H...N and O...H...O types and a nonsymmetrical transition state of the N...H...C, O...H...C, and O...H...N types. These distances are compared with other characteristics of the reactions and reactants. The influence of such factors as the enthalpy of the reaction, triplet repulsion in the transition state, electronegativity of atoms in the reaction center, and steric repulsion of reactants on the geometric parameters of the transition state is examined. The results obtained are used to calculate the increments characterizing the influence of various factors on interatomic distances of radical reactions considered.Translated from Kinetika i Kataliz, Vol. 45, No. 6, 2004, pp. 877–896.Original Russian Text Copyright © 2004 by Denisova, Denisov.     

Acetylene cyclotrimerization by early second-row transition metals in the gas phase. A theoretical study

The acetylene cyclotrimerization reaction mediated by the left-hand-side bare transition metal atoms Y, Zr, Nb, and Mo has been studied theoretically, employing DFT in its B3LYP formulation. The complete reaction mechanism has been analyzed, identifying intermediates and transition states. Both the ground spin state and at least one low-lying excited state have been considered to establish whether possible spin crossings between surfaces of different multiplicity can occur. Our results show that the overall reaction is highly favorable from a thermodynamic point of view and ground state transition states lie always below the energy limit represented by ground state reactants. After the activation of two acetylene molecules and formation of a bis-ligated complex, the reaction proceeds to give a metallacycle intermediate, as the alternative formation of a cyclobutadiene complex is energetically disfavored. All the examined reaction paths involve formation of a metallacycloheptatriene intermediate that in turn generates a metal-benzene adduct from which finally benzene is released. Similarities and differences in the behaviors of the considered four metal atoms have been examined.     

An analytic hindered rotor model for calculating microcanonical variational unimolecular rate constants from reaction path properties

A model is proposed for performing microcanonical variational transition state theory calculations which incorporates ideas from vibrator and flexible variational transition state models. Vibrational frequencies, moments of inertia, and potential energy for the variational transition state are found by reaction path following as for the vibrator model. However, the transitional modes are treated as hindered rotors using an analytic potential and an analytic density of states, which are fit to barriers for hindered rotation determined from reaction path following. The model proposed here differs from the flexible transition state model in that the density of states for the transitional modes is analytic and transitional modes and external rotational angular momenta are uncoupled. For the H + CH₃ ? CH₄ system, rate constants calculated with this new model are only 6–23% smaller on average from those of the flexible transition state model for values of total angular momentum which correspond to average rotational temperatures of 0–2000 K. Harmonic frequencies calculated for the transitional modes from the hindered rotor Hamiltonian are in good agreement with the exact values found by a reaction path analysis. © 1994 John Wiley & Sons, Inc.     

A new perspective on transition states: χ1 separatrix

We present a new definition of the transition state for chemical reactions, named the χ(1) separatrix. In contrast to previous transition state definitions which depend on the choice of reaction coordinates, the χ(1) separatrix is defined by choosing a time scale for observation and is connected to exact rate constants in the high friction limit. We demonstrate that this separatrix appears in the isomerization of alanine dipeptide as a stationary population in quasi-equilibrium, without assuming a particular coordinate system or reactant and product surfaces.     

Effect of absolute laser phase on reaction paths in laser-induced chemical reactions

Potential surfaces, dipole moments, and polarizabilities are calculated by ab initio methods [unrestricted MP2(full)/6-311++G(2d,2p)] along the reaction paths of the F+CH4 and Cl+CH4 reaction systems. It is found that in general dipole moments and polarizabilities exhibit peaks near the transition state. In the case of X=F these peaks are on the products side and in the case of X=Cl they are on the reactants side indicating an early transition state in the case of fluorine and a late transition state in the case of chlorine. An analysis of the geometric changes along the reaction paths reveals a one-to-one correspondence between the peaks in the electric properties and peaks in the rate of change of certain internal geometric coordinates along the reaction path. Interaction with short infrared intense laser fields pulses leads to the possibility of interferences between the dipole and polarizability laser-molecule interactions as a function of laser phase. The larger dipole moment in the Cl+CH4 reaction can lead to the creation of deep wells (instead of energy barriers) and new strongly bound states in the transition state region. This suggests possible coherent control of the reaction path as a function of the absolute phase of the incident field, by significant modification of the potential surfaces along the reaction path and, in particular, in the transition state region.     

The role of the putative catalytic base in the phosphoryl transfer reaction in a protein kinase: first-principles calculations

Protein kinases are important enzymes controlling the majority of cellular signaling events via a transfer of the gamma-phosphate of ATP to a target protein. Even after many years of study, the mechanism of this reaction is still poorly understood. Among many factors that may be responsible for the 1011-fold rate enhancement due to this enzyme, the role of the conserved aspartate (Asp166) has been given special consideration. While the essential presence of Asp166 has been established by mutational studies, its function is still debated. The general base catalyst role assigned to Asp166 on the basis of its position in the active site has been brought into question by the pH dependence of the reaction rate, isotope measurements, and pre-steady-state kinetics. Recent semiempirical calculations have added to the controversy surrounding the role of Asp166 in the catalytic mechanism. No major role for Asp166 has been found in these calculations, which have predicted the reaction process consisting of an early transfer of a substrate proton onto the phosphate group. These conclusions were inconsistent with experimental observations. To address these differences between experimental results and theory with a more reliable computational approach and to provide a theoretical platform for understanding catalysis in this important enzyme family, we have carried out first-principles structural and dynamical calculations of the reaction process in cAPK kinase. To preserve the essential features of the reaction, representations of all of the key conserved residues (82 atoms) were included in the calculation. The structural calculations were performed using the local basis density functional (DFT) approach with both hybrid B3LYP and PBE96 generalized gradient approximations. This kind of calculation has been shown to yield highly accurate structural information for a large number of systems. The optimized reactant state structure is in good agreement with X-ray data. In contrast to semiempirical methods, the lowest energy product state places the substrate proton on Asp166. First-principles molecular dynamics simulations provide additional support for the stability of this product state. The latter also demonstrate that the proton transfer to Asp166 occurs at a point in the reaction where bond cleavage at the PO bridging position is already advanced. This mechanism is further supported by the calculated structure of the transition state in which the substrate hydroxyl group is largely intact. A metaphoshate-like structure is present in the transition state, which is consistent with the X-ray structures of transition state mimics. On the basis of the calculated structure of the transition state, it is estimated to be 85% dissociative. Our analysis also indicates an increase in the hydrogen bond strength between Asp166 and substrate hydroxyl and a small decrease in the bond strength of the latter in the transition state. In summary, our calculations demonstrate the importance of Asp166 in the enzymatic mechanism as a proton acceptor. However, the proton abstraction from the substrate occurs late in the reaction process. Thus, in the catalytic mechanism of cAPK protein kinase, Asp166 plays a role of a "proton trap" that locks the transferred phosphoryl group to the substrate. These results resolve prior inconsistencies between theory and experiment and bring new understanding of the role of Asp166 in the protein kinase catalytic mechanism.     

Experimental and theoretical studies on the isomerization of allyl thiocyanate to allyl isothiocyanate

The mechanism of isomerization of allyl thiocyanate to allyl isothiocyanate has been investigated both experimentally and theoretically. The kinetic study indicates that the reaction is unimolecular and is not ionic. The entropy of activation suggests strongly that the mechanism involves a cyclic transition state. The rate of reaction was retarded to a small extent in polar solvents relative to that in nonpolar solvents. Ab initio MO calculations indicate, in agreement with the experimental results, that the reaction proceeds through a cyclic transition state, one in which the SCN moiety is almost linear. Thus, this is a [3,3] sigmatropic rearrangement. The charge separation in the transition state was substantial. The retardation of the reaction in polar solvents was attributed to the difference in solvation in the original state and in the transition state. © 1997 John Wiley & Sons, Inc.     

Force reversed method for locating transition states

To identify the transition state accurately and efficiently on a high-dimensional potential energy surface is one of the most important topics in kinetic studies on chemical reactions. We present here an algorithm to search the transition state by so-called force reversed method, which only requires a rough reaction direction instead of knowing the initial state and final state. Compared to the nudged elastic band method and the dimer method that require multiple images, the present algorithm with only single image required saves significantly the computational cost. The algorithm was implemented in the first-principle periodic total energy calculation package and applied successfully to several prototype surface processes such as the adsorbate diffusion and dissociation on metal surfaces. The results indicate that the force reversed method is efficient, robust to identify the transition state of various surface processes.