Mathematical foundation of a global strategy for searching reaction paths

For the determination of reaction paths and critical points on the potential energy hypersurface of chemical reactions, a rigorous mathematical background for the theory of a global searching procedure based on the catchment regions of the gradient field is given.     

A procedure for determining reaction paths and saddle points

To determine reaction paths, a method that does not require the previous location of extrema is presented and illustrated by an example. The procedure is based on a local symmetry property of the potential surface.     

A method for determining reaction paths in large molecules: Application to myoglobin

An algorithm is described for determining reaction paths between two known structures with many degrees of freedom. The method uses first-derivative techniques to optimize the entire path between the two end forms subject to certain constraints. The stability and convergence properties of the method are illustrated by applications to structural transitions in two test systems (cyclohexane and dialanine) and to a conformational change involving all degrees of freedom in the protein, myoglobin, with 1531 atoms.     

An algorithm for determining dynamically defined reaction paths (DDRP)

Summary A numerically stable and well-parallelizable curve variational algorithm is described for determining tangent curves of vector fields between two given stationary points. In particular, the method is suitable for finding reaction paths and saddle points on potential energy hypersurfaces (PHS). The stability of the procedure is illustrated by an artificial mathematical function, showing phases of following the reaction on the PHS.Dedicated to Professor Zoltán G. Szabó, the great teacher and scientist in reaction kinetics and in many other fields of physical chemistry, on his 84th birthday.     

Coordinate reduction for exploring chemical reaction paths

The potential energy surface for the reaction of a typical molecular system composed of N atoms is defined uniquely by 3N-6 coordinates. These coordinates can be defined by the Cartesian coordinates of the atomic centers (minus overall translation and rotation), or a set of internally defined coordinates such as bond stretches, angle bends, and torsions. By applying principal component analysis to the geometries along a reaction path, a reduced set of coordinates, d ≪ 3N-6, can be obtained. This reduced set of coordinates can reproduce the changes in geometry along the reaction path with chemical accuracy and may help improve the efficiency of reaction path optimization algorithms.     

Intrinsic reaction coordinate calculations for reaction paths possessing branching points

A 3-21+G energy surface corresponding to the proton transfer reaction in the hydroperoxyl anion solvated by one water molecule presents interesting topological features. In particular the intrinsic reaction coordinate that begins at the transition state does not lead to a minimum but to a saddle point of second order passing through two branching points. A new strategy to obtain the true reaction path in these cases is proposed.     

Survival paths for reaction dynamics in fluctuating environments

We study rate processes in general Gaussian fluctuating environments using a path integral formalism. We derive a variational equation for the dominant survival path when the fluctuations relax exponentially or according to a stretched exponential law. In the case of a slowly varying barrier, the equilibrium regression approximation which is used by Frauenfelder and coworkers emerges. In this approximation, the survival path follows the ordinary law of relaxation to equilibrium. If the rate coefficients vary rapidly with environmental variables, however, the dominant survival paths exhibit more complex behaviour. Many phenomena analogous to geometrical optics occur. These include reflection off of rapid variations in rate constant, as well as refraction, giving paths very different from the equilibrium relaxation properties. A model with a piece-wise linear rate exhibits the basic phenomena, and the survival path equation is exactly solved for the general stretched exponential relaxing environment.     

Tunneling reactions with two reaction paths

The WKB method and Taylor series expansions are used to derive an approximate expression relating the tunneling rates with two reaction paths to the tunneling rates of the individual paths. The two reaction paths originate in the same reactant well, travel over different potential barriers, and finish in the same product well. The validity of the approximate expression and the applicability of the WKB method to this problem are investigated.     

Exploring SCC-DFTB paths for mapping QM/MM reaction mechanisms

A new first-order procedure for locating transition structures (TS) that employs hybrid quantum mechanical/molecular mechanical (QM/MM) potentials has been developed. This new technique (RPATh+RESD) combines the replica path method (RPATh) and standard reaction coordinate driving (RCD) techniques in an approach that both efficiently determines reaction barriers and successfully eliminates two key weaknesses of RCD calculations (i.e., hysteresis/discontinuities in the path and the sequential nature of the RCD procedure). In addition, we have extended CHARMM's QM/MM reaction pathway methods, the RPATh and nudged elastic band (NEB) methods, to incorporate SCC-DFTB wave functions. This newly added functionality has been applied to the chorismate mutase-catalyzed interconversion of chorismate to prephenate, which is a key step in the shikimate pathway of bacteria, fungi, and other higher plants. The RPATh+RESD barrier height (DeltaE=5.7 kcal/mol) is in good agreement with previous results from full-energy surface mapping studies (Zhang, X.; Zhang, X.; Bruice, T. C. Biochemistry 2005, 44, 10443-10448). Full reaction paths were independently mapped with RPATh and NEB methods and showed good agreement with the final transition state from the RPATh+RESD "gold standard" and previous high-level QM/MM transition states (Woodcock, H. L.; Hodoscek, M.; Gilbert, T. B.; Gill, P. M. W.; Schaefer, H. F.; Brooks, B. R. J. Comput. Chem. 2007, 28, 1485-1502). The SCC-DFTB TS geometry most closely approximates the MP2/6-31+G(d) QM/MM result. However, the barrier height is underestimated and possibly points to an area for improvement in SCC-DFTB parametrization. In addition, the steepest descents (SD) minimizer for the NEB method was modified to uncouple the in-path and off-path degrees of freedom during the minimization, which significantly improved performance. The convergence behavior of the RPATh and NEB was examined for SCC-DFTB wave functions, and it was determined that, in general, both methods converge at about the same rate, although the techniques used for convergence may be different. For instance, RPATh can effectively use the adopted basis Newton-Raphson (ABNR) minimizer, where NEB seems to require a combination of SD and ABNR.     

Using swarm intelligence for finding transition states and reaction paths

We describe an algorithm that explores potential energy surfaces (PES) and finds approximate reaction paths and transition states. A few (≈6) evolving atomic configurations ("climbers") start near a local minimum M1 of the PES. The climbers seek a shallow ascent, low energy, path toward a saddle point S12, cross over to another valley of the PES, and climb down to a new minimum M2 that was not known beforehand. Climbers use both energy and energy derivatives to make individual decisions, and they use relative fitness to make team-based decisions. In sufficiently long runs, they keep exploring and may go through a sequence M1-S12-M2-S23-M3 ... of minima and saddle points without revisiting any of the critical points. We report results on eight small test systems that highlight advantages and disadvantages of the method. We also investigated the PES of Li(8), Al(7)(+), Ag(7), and Ag(2)NH(3) to illustrate potential applications of this new method.     

Thermochemistry for enthalpies and reaction paths of nitrous acid isomers

Recent studies show that nitrous acid, HONO, a significant precursor of the hydroxyl radical in the atmosphere, is formed during the photolysis of nitrogen dioxide in soils. The term nitrous acid is largely used interchangeably in the atmospheric literature, and the analytical methods employed do not often distinguish between the HONO structure (nitrous acid) and HNO₂ (nitryl hydride or isonitrous acid). The objective of this study is to determine the thermochemistry of the HNO₂ isomer, which has not been determined experimentally, and to evaluate its thermal and atmospheric stability relative to HONO. The thermochemistry of these isomers is also needed for reference and internal consistency in the calculation of larger nitrite and nitryl systems. We review, evaluate, and compare the thermochemical properties of several small nitric oxide and hydrogen nitrogen oxide molecules. The enthalpies of HONO and HNO₂ are calculated using computational chemistry with the following methods of analysis for the atomization, isomerization, and work reactions using closed‐ and open‐shell reference molecules. Three high‐level composite methods G3, CBS‐QB3, and CBS‐APNO are used for the computation of enthalpy. The enthalpy of formation, ΔH^o_f(298 K), for HONO is determined as ?18.90 ± 0.05 kcal mol^?1 (?79.08 ± 0.2 kJ mol^?1) and as ?10.90 ± 0.05 kcal mol^?1 (?45.61 ± 0.2 kJ mol^?1) for nitryl hydride (HNO₂), which is significantly higher than values used in recent NO_x combustion mechanisms. H‐NO₂ is the weakest bond in isonitrous acid; but HNO₂ will isomerize to HONO with a similar barrier to the HO? NO bond energy; thus, it also serves as a source of OH in atmospheric chemistry. Kinetics of the isomerization is determined; a potential energy diagram of H/N/O₂ system is presented, and an analysis of the triplet surface is initiated. © 2007 Wiley Periodicals, Inc. Int J Chem Kinet 39: 378–398, 2007     

Finding reaction paths using the potential energy as reaction coordinate

The intrinsic reaction coordinate curve (IRC), normally proposed as a representation of a reaction path, is parametrized as a function of the potential energy rather than the arc-length. This change in the parametrization of the curve implies that the values of the energy of the potential energy surface points, where the IRC curve is located, play the role of reaction coordinate. We use Caratheodory's relation to derive in a rigorous manner the proposed parametrization of the IRC path. Since this Caratheodory's relation is the basis of the theory of calculus of variations, then this fact permits to reformulate the IRC model from this mathematical theory. In this mathematical theory, the character of the variational solution (either maximum or minimum) is given through the Weierstrass E-function. As proposed by Crehuet and Bofill [J. Chem. Phys. 122, 234105 (2005)], we use the minimization of the Weierstrass E-function, as a function of the potential energy, to locate an IRC path between two minima from an arbitrary curve on the potential energy surface, and then join these two minima. We also prove, from the analysis of the Weierstrass E-function, the mathematical bases for the algorithms proposed to locate the IRC path. The proposed algorithm is applied to a set of examples. Finally, the algorithm is used to locate a discontinuous, or broken, IRC path, namely, when the path connects two first order saddle points through a valley-ridged inflection point.     

Stereochemistry of reaction paths at carbonyl centres

Results of recent experimental and theoretical studies of nucleophilic addition to carbonyl groups are described. The reaction paths found by different methods for different nucleophiles show some striking similarities that appear to be characteristic for the reaction type.     

Chemical reaction paths and calculus of variations

The reaction path is an important concept of theoretical chemistry. We analyze different forms of reaction pathways in the light of the abstract theory of calculus of variations such as steepest descent from saddle point, the intrinsic reaction coordinate (IRC), Newton trajectory, variationally optimized reaction paths and others. The paper is both a mathematical review and a pointer to future research. Besides the theoretical definitions, we shortly discuss hints at the numerical effect of the definitions.     

Local Hammond postulate and a local model for monitoring reaction paths

A local model for the following reaction path is suggested. Exploratory calculations on simple isomerization reactions are undertaken. A local version of Hammond's postulate is proposed and tested on some model systems.     

An algorithm for the location of branching points on reaction paths

A branching point is a point on a reaction path leading from reactants to products (via a transition state) at which it is energetically favorable for the system to break symmetry. Such a point can be defined in terms of normal modes along the reaction path and corresponds to zero curvature (a zero Hessian eigenvalue) along a symmetry-breaking mode. An effective method for the location of such points is presented and realized in an efficient, practical algorithm designed for use in the ab initio program package Gaussian 82.     

MINDO-Forces calculation of molecular geometries and reaction paths

A new method for the minimization of molecular energies is described, based on the Murtagh-Sargent procedure and the MINDO/2(3) semiempirical MO method. The derivative of energy is calculated according to the Pulay's Force method. This method was applied to the calculation of heats of formation and geometries of different organic molecules. The results agree well with the experimental and theoretical values known in the literature. The MINDO-Forces method was applied to the calculation of the internal rotation barrier of the benzyl carbonium ion. The calculated value 18.79 kcal/mol agrees well with the theoretical values known in the literature. The MINDO-Forces method was applied to the calculation of the internal rotation barrier of cyclopropyl carbinyl cation. The calculated rotation barrier, 21.28 kcal/mol is in agreement with the known theoretical values and with the expectation based on the NMR measurement of the barrier height in cyclopropyl dimethyl-carbinyl cation. Part of the Ph.D. thesis, S. M. Khalil, University of Baghdad 1976.     

Analysis of molecular shape changes along reaction paths

In this work, a detailed discussion and illustrative example are given for a new conceptual approach to the notion of similarity between nuclear configurations. We analyze the changes in molecular shape, as described by a topological characterization of molecular envelope surfaces, along reaction paths. This approach provides an independent definition of formal species along the path, defined in terms of molecular shape, rather than based on the curvature properties of the potential energy function. The two approaches are compared in this work, using the [1,3] hydrogen-shift in formic acid as an illustrative example. The reaction is studied at SCF and MP 2 ab initio levels. Although the actual reaction paths are different, the molecular shape changes encountered along them are similar, indicating that the essential shape changes provide chemical information on a more fundamental level than do detailed reaction paths.     

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.