A combined explicit-implicit method for high accuracy reaction path integration

We present the use of an optimal combined explicit-implicit method for following the reaction path to high accuracy. This is in contrast to most purely implicit reaction path integration algorithms, which are only efficient on stiff ordinary differential equations. The defining equation for the reaction path is considered to be stiff, however, we show here that the reaction path is not uniformly stiff and instead is only stiff near stationary points. The optimal algorithm developed in this work is a combination of explicit and implicit methods with a simple criterion to switch between the two. Using three different chemical reactions, we combine and compare three different integration methods: the implicit trapezoidal method, an explicit stabilized third order algorithm implemented in the code DUMKA3 and the traditional explicit fourth order Runge-Kutta method written in the code RKSUITE. The results for high accuracy show that when the implicit trapezoidal method is combined with either explicit method the number of energy and gradient calculations can potentially be reduced by almost a half compared with integrating either method alone. Finally, to explain the improvements of the combined method we expand on the concepts of stability and stiffness and relate them to the efficiency of integration methods.     

Determination of complex reaction mechanisms. Analysis of chemical, biological and genetic networks

We present several methods of determining, not guessing, complex chemical reaction mechanisms and their functions. One method is based on the theory of correlation functions of measured time series of concentrations of chemical species; another is on measurements of temporal responses of concentrations to various perturbations of arbitrary magnitude; a third deals with the analysis of oscillatory systems; a fourth is on the use of genetic algorithms to determine functions of chemical reaction networks. All methods are applicable to chemical, biochemical, and biological reaction systems and to genetic networks and systems biology. The methods depend on the design of appropriate experiments on the whole system and corresponding theories for interpretation that lead to information on the causal chemical connectivity of species, on reaction pathways, on reaction mechanisms, on control centers in the system, and on functions of the system. The first three methods require no assumption of a model or hypothesis, nor extensive calculations, unlike the interpretation of measurements made on a gene network at only one time.     

Reaction path synthesis methodology for waste minimization

Being the core of whole process, a chemical process is a key step for reducing waste generation. Therefore, to organize a reasonable reaction path is very important for reducing waste emission. Because of the large amount of the feasible reaction path and other factors in a chemical reaction system, the opti-mizing for reaction path is difficult. Since the reaction path synthesis problem has presented, lots of re-searchers focused on it. At the beginning, only chem-ists took the research work,…     

Adapting the nudged elastic band method for determining minimum-energy paths of chemical reactions in enzymes

Optimization of reaction paths for enzymatic systems is a challenging problem because such systems have a very large number of degrees of freedom and many of these degrees are flexible. To meet this challenge, an efficient, robust and general approach is presented based on the well-known nudged elastic band reaction path optimization method with the following extensions: (1) soft spectator degrees of freedom are excluded from path definitions by using only inter-atomic distances corresponding to forming/breaking bonds in a reaction; (2) a general transformation of the distances is defined to treat multistep reactions without knowing the partitioning of steps in advance; (3) a multistage strategy, in which path optimizations are carried out for reference systems with gradually decreasing rigidity, is developed to maximize the opportunity of obtaining continuously changing environments along the path. We demonstrate the applicability of the approach using the acylation reaction of type A beta-lactamase as an example. The reaction mechanism investigated involves four elementary reaction steps, eight forming/breaking bonds. We obtained a continuous minimum energy path without any assumption on reaction coordinates, or on the possible sequence or the concertedness of chemical events. We expect our approach to have general applicability in the modeling of enzymatic reactions with quantum mechanical/molecular mechanical models.     

Examining the robustness of first-principles calculations for metal hydride reaction thermodynamics by detection of metastable reaction pathways

First principles calculations have played a useful role in screening mixtures of complex metal hydrides to find systems suitable for H(2) storage applications. Standard methods for this task efficiently identify the lowest energy reaction mechanisms among all possible reactions involving collections of materials for which DFT calculations have been performed. The resulting mechanism can potentially differ from physical reality due to inaccuracies in the DFT functionals used, or due to other approximations made in estimating reaction free energies. We introduce an efficient method to probe the robustness of DFT-based predictions that relies on identifying reactions that are metastable relative to the lowest energy reaction path predicted with DFT. An important conclusion of our calculations is that in many examples DFT cannot unambiguously predict a single reaction mechanism for a well defined metal hydride mixture because two or more mechanisms have reaction energies that differ by a small amount. Our approach is illustrated by analyzing a series of single step reactions identified in our recent work that examined reactions with a large database of solids [Kim et al., Phys. Chem. Chem. Phys. 2011, 13, 7218].     

Steepest descent reaction path integration using a first-order predictor-corrector method

The theoretical treatment of chemical reactions inevitably includes the integration of reaction pathways. After reactant, transition structure, and product stationary points on the potential energy surface are located, steepest descent reaction path following provides a means for verifying reaction mechanisms. Accurately integrated paths are also needed when evaluating reaction rates using variational transition state theory or reaction path Hamiltonian models. In this work an Euler-based predictor-corrector integrator is presented and tested using one analytic model surface and five chemical reactions. The use of Hessian updating, as a means for reducing the overall computational cost of the reaction path calculation, is also discussed.     

An efficient molecular dynamics simulation method for calculating the diffusion-influenced reaction rates

We present a molecular dynamics (MD) simulation method for calculating the diffusion-influenced reaction rates in the limit of low reactant concentrations. To calculate the reaction rate coefficient, we use MD trajectories of a nonreactive equilibrium system that are initiated with a pair of reactant molecules in reactive configuration. Hence reaction systems involving complicated reactant molecules with geometrically restricted reactivities can be treated with comparable efficiency as the simple hard-sphere reaction system. Compared to the similar MD method proposed by Van Beijeren, Dong, and Bocquet [J. Chem. Phys. 114, 6265 (2001)], the present method has a couple of advantages. First, reactions involving more general sink functions can be treated. Second, more accurate results can be obtained when the reaction probability upon collision is less than unity. As an application, we investigate the effects of nondiffusive dynamics and hydrodynamic interaction of reactants on the reaction rate.     

Multistructural microiteration technique for geometry optimization and reaction path calculation in large systems

We propose a multistructural microiteration (MSM) method for geometry optimization and reaction path calculation in large systems. MSM is a simple extension of the geometrical microiteration technique. In conventional microiteration, the structure of the non‐reaction‐center (surrounding) part is optimized by fixing atoms in the reaction‐center part before displacements of the reaction‐center atoms. In this method, the surrounding part is described as the weighted sum of multiple surrounding structures that are independently optimized. Then, geometric displacements of the reaction‐center atoms are performed in the mean field generated by the weighted sum of the surrounding parts. MSM was combined with the QM/MM‐ONIOM method and applied to chemical reactions in aqueous solution or enzyme. In all three cases, MSM gave lower reaction energy profiles than the QM/MM‐ONIOM‐microiteration method over the entire reaction paths with comparable computational costs. © 2017 Wiley Periodicals, Inc.     

Determination of initial reaction rates using wilkinson's relation

Uncertainty sometimes exists in determining initial reaction rates from experimental data. A method, originally proposed by Wilkinson [9] for estimating orders and rate constants for simple batch nth order reactions, has been generalized to complex kinetic systems. This method yields very accurate initial rates for all systems and extends the conversion range of experimental investigation of initial rates well beyond the “zero-order” region. Accurate initial rates are required in analytical methods used for screening alternate reaction mechanisms.     

Dynamics of the reaction glucose-catalase-glucose oxidase-hydrogen peroxide

Glucose-catalase-glucose oxidase-hydrogen peroxide reaction is one of the few known enzymatic systems studied in vitro in the field of nonlinear chemical dynamics. This reaction belongs to the family of oscillatory enzymatic reactions, which form a natural basis of oscillations in biological systems. A parametric study of dependence on mixing, temperature and initial concentrations of components in a batch stirred reactor was carried out. A newly proposed mathematical model of the reaction conforms to the obtained experimental data. Results of our experiments and simulations hint at further directions of research of non-linear dynamics in this reaction.     

Theoretical investigation of the gas-phase reaction of NiO+ with ethane

To explore the mechanisms for Ni-based oxide-catalyzed oxidative dehydrogenation (ODH) reactions, we investigate the reactions of C₂H₆ with NiO⁺ using density functional calculations. Two possible reaction pathways are identified, which lead to the formation of ethanol (path 1), ethylene and water (path 2). The proportion of products is discussed by Curtin-Hammett principle, and the result shows that path 2 is the main reaction channel and the water and ethylene are the main products. In order to get a deeper understanding of the titled reaction, numerous means of analysis methods including the atoms in molecules (AIM), electron localization function (ELF), natural bond orbital (NBO), and density of states (DOS) are used to study the properties of the chemical bonding evolution along the reaction pathways.

     

Efficiency analysis of reaction rate calculation methods using analytical models I: The two-dimensional sharp barrier

We analyze the efficiency of different methods for the calculation of reaction rates in the case of a simple two-dimensional analytical benchmark system. Two classes of methods are considered: the first is based on the free energy calculation along a reaction coordinate and the calculation of the transmission coefficient, the second on the sampling of dynamical pathways. We give scaling rules for how this efficiency depends on barrier height and width, and we hand out simple optimization rules for the method-specific parameters. We show that the path sampling methods, using the transition interface sampling technique, become exceedingly more efficient than the others when the reaction coordinate is not the optimal one.     

Exploring reaction pathways with transition path and umbrella sampling: application to methyl maltoside

The transition path sampling (TPS) method is a powerful approach to study chemical reactions or transitional properties on complex potential energy landscapes. One of the main advantages of the method over potential of mean force methods is that reaction rates can be directly accessed without knowledge of the exact reaction coordinate. We have investigated the complementary nature of these two differing approaches, comparing transition path sampling with the weighted histogram analysis method to study a conformational change in a small model system. In this case study, the transition paths for a transition between two rotational conformers of a model disaccharide molecule, methyl beta-D-maltoside, were compared with a free energy surface constrained by the two commonly used glycosidic (phi,psi) torsional angles. The TPS method revealed a reaction channel that was not apparent from the potential of mean force method, and the suitability of phi and psi as reaction coordinates to describe the isomerization in vacuo was confirmed by examination of the transition path ensemble. Using both transition state theory and transition path sampling methods, the transition rate was estimated. We have estimated a characteristic time between transitions of approximately 160 ns for this rare isomerization event between the two conformations of the carbohydrate. We conclude that transition path sampling can extract subtle information about the dynamics not apparent from the potential of mean force method. However, in calculating the reaction rate, the transition path sampling method required 27.5 times the computational effort than was needed by the potential of mean force method.     

Semiempirical QM/MM potential with simple valence bond (SVB) for enzyme reactions. Application to the nucleophilic addition reaction in haloalkane dehalogenase

 We present a method for the correction of errors in combined QM/MM calculations using a semiempirical Hamiltonian for enzyme reactions. Since semiempirical models can provide a reasonable representation of the general shape of the potential energy surface for chemical reactions, we introduce a simple valence bond-like (SVB) term to correct the energies at critical points on the potential energy surface. The present SVB term is not a stand-alone potential energy function, but it is used purely for introducing small energy corrections to the semiempirical Hamiltonian to achieve the accuracy needed for modeling enzymatic reactions. We show that the present coupled QM-SVB/MM approach can be parameterized to reproduce experimental and ab initio results for model reactions, and have applied the PM3-SVB/MM potential to the nucleophilic addition reaction in haloalkane dehalogenase. In a preliminary energy minimization study, the PM3-SVB/MM results are reasonable, suggesting that it may be used in free energy simulations to assess enzymatic reaction mechanism. Received: 1 November 2001 / Accepted: 6 September 2002 / Published online: 19 February 2003 Contribution to the Proceedings of the Symposium on Combined QM/MM Methods at the 222nd National Meeting of the American Chemical Society, 2001 Correspondence to: Lakshmi S. Devi-Kesavan e-mail: kesavan@chem.umn.edu Acknowledgments. The work is partially supported by the NIH and the NSF.     

Cope reaction families: to be or not to be a biradical

The cyclization reactions of polyunsaturated systems such as enediynes, enyne-allenes, as well as many others, can be structurally related to the Cope reaction and grouped into two branches: 6pi (e.g., Cope, allenyl-Cope) and [2sigma + 4pi] systems (e.g., Bergman, Myers-Saito). After extensive computional (BCCD(T)/cc-pVDZ//UBLYP/6-31G) examination of a large set of Cope-related reactions, we derive the following simple rule: biradical intermediates are involved along the reaction path if they are stabilized by allylic resonance or aromaticity.     

Automatic determination of reaction mappings and reaction center information. 1. The imaginary transition state energy approach

Chemical reactions transform the reactant molecules by deleting existing and forming new bonds. The identification of these so-called reacting bonds is important for studying the reaction mechanism and for applications in metabolomics, e.g. for interpreting substrate labeling experiments. Here, we introduce an approach which suggests the simplest possible reaction center at the heavy atom level, with high accuracy. In contrast to current methods the approach is motivated by a simple theoretical model based on a crude approximation of the reaction energetics, and takes the complete reacting system into account. Finally, it recovers all optimal solutions to the problem while removing all symmetry-related, redundant solutions. We apply the method on the complete KEGG database of biochemical reactions, and compare our approach with previous methods. The resulting reaction centers are represented as imaginary transition states, which are molecule-like representations of reaction mechanisms. We provide the statistics of the calculations on the KEGG database and discuss some examples for the different types of alternative solutions found.     

Parallel iterative reaction path optimization in ab initio quantum mechanical/molecular mechanical modeling of enzyme reactions

The determination of reaction paths for enzyme systems remains a great challenge for current computational methods. In this paper we present an efficient method for the determination of minimum energy reaction paths with the ab initio quantum mechanical/molecular mechanical approach. Our method is based on an adaptation of the path optimization procedure by Ayala and Schlegel for small molecules in gas phase, the iterative quantum mechanical/molecular mechanical (QM/MM) optimization method developed earlier in our laboratory and the introduction of a new metric defining the distance between different structures in the configuration space. In this method we represent the reaction path by a discrete set of structures. For each structure we partition the atoms into a core set that usually includes the QM subsystem and an environment set that usually includes the MM subsystem. These two sets are optimized iteratively: the core set is optimized to approximate the reaction path while the environment set is optimized to the corresponding energy minimum. In the optimization of the core set of atoms for the reaction path, we introduce a new metric to define the distances between the points on the reaction path, which excludes the soft degrees of freedom from the environment set and includes extra weights on coordinates describing chemical changes. Because the reaction path is represented by discrete structures and the optimization for each can be performed individually with very limited coupling, our method can be executed in a natural and efficient parallelization, with each processor handling one of the structures. We demonstrate the applicability and efficiency of our method by testing it on two systems previously studied by our group, triosephosphate isomerase and 4-oxalocrotonate tautomerase. In both cases the minimum energy paths for both enzymes agree with the previously reported paths.     

A partial-propensity formulation of the stochastic simulation algorithm for chemical reaction networks with delays

Several real-world systems, such as gene expression networks in biological cells, contain coupled chemical reactions with a time delay between reaction initiation and completion. The non-Markovian kinetics of such reaction networks can be exactly simulated using the delay stochastic simulation algorithm (dSSA). The computational cost of dSSA scales with the total number of reactions in the network. We reduce this cost to scale at most with the smaller number of species by using the concept of partial reaction propensities. The resulting delay partial-propensity direct method (dPDM) is an exact dSSA formulation for well-stirred systems of coupled chemical reactions with delays. We detail dPDM and present a theoretical analysis of its computational cost. Furthermore, we demonstrate the implications of the theoretical cost analysis in two prototypical benchmark applications. The dPDM formulation is shown to be particularly efficient for strongly coupled reaction networks, where the number of reactions is much larger than the number of species.     

RCARM: Reaction classification using automated reaction mapping

Detailed chemical kinetic modeling of gas‐phase reactions can result in automatically generated mechanisms that contain thousands of reactions. In this paper, we describe the development of a rule‐based expert system tool that organizes these reactions into classes such as hydrogen abstraction and beta scission. We have developed 29 simple classification rules, 20 complex (well‐skipping) classification rules, and four second‐stage classification rules. This greatly simplifies the task of the chemical kineticist who wishes to verify, analyze, and gain insights into the reactions comprising the mechanism. This system, which is based on the automated identification of the bonds that break and form in a chemical reaction (the reaction mapping problem), is used to classify reactions in three different mechanisms. © 2012 Wiley Periodicals, Inc. Int J Chem Kinet 45: 125–139, 2013     

Automatic integration of the reaction path using diagonally implicit Runge-Kutta methods

The diagonally implicit Runge-Kutta framework is shown to be a general form for constructing stable, efficient steepest descent reaction path integrators, of any order. With this framework tolerance driven, adaptive step-size methods can be constructed by embedding methods to obtain error estimates of each step without additional computational cost. There are many embedded and nonembedded, diagonally implicit Runge-Kutta methods available from the numerical analysis literature and these are reviewed for orders two, three, and four. New embedded methods are also developed which are tailored to the application of reaction path following. All integrators are summarized and compared for three systems: the Muller-Brown [Theor. Chem. Acta 53, 75 (1979)] potential and two gas phase chemical reactions. The results show that many of the methods are capable of integrating efficiently while reliably keeping the error bound within the desired tolerance. This allows the reaction path to be determined through automatic integration by only specifying the desired accuracy and transition state.