Separability of tight and roaming pathways to molecular decomposition

Recent studies have questioned the separability of the tight and roaming mechanisms to molecular decomposition. We explore this issue for a variety of reactions including MgH(2) → Mg + H(2), NCN → CNN, H(2)CO → H(2) + CO, CH(3)CHO → CH(4) + CO, and HNNOH → N(2) + H(2)O. Our analysis focuses on the role of second-order saddle points in defining global dividing surfaces that encompass both tight and roaming first-order saddle points. The second-order saddle points define an energetic criterion for separability of the two mechanisms. Furthermore, plots of the differential contribution to the reactive flux along paths connecting the first- and second-order saddle points provide a dynamic criterion for separability. The minimum in the differential reactive flux in the neighborhood of the second-order saddle point plays the role of a mechanism divider, with the presence of a strong minimum indicating that the roaming and tight mechanisms are dynamically distinct. We show that the mechanism divider is often, but not always, associated with a second-order saddle point. For the formaldehyde and acetaldehyde reactions, we find that the minimum energy geometry on a conical intersection is associated with the mechanism divider for the tight and roaming processes. For HNNOH, we again find that the roaming and tight processes are dynamically separable but we find no intrinsic feature of the potential energy surface associated with the mechanism divider. Overall, our calculations suggest that roaming and tight mechanisms are generally separable over broad ranges of energy covering most kinetically relevant regimes.     

Comparison of methods for finding saddle points without knowledge of the final states

Within the harmonic approximation to transition state theory, the biggest challenge involved in finding the mechanism or rate of transitions is the location of the relevant saddle points on the multidimensional potential energy surface. The saddle point search is particularly challenging when the final state of the transition is not specified. In this article we report on a comparison of several methods for locating saddle points under these conditions and compare, in particular, the well-established rational function optimization (RFO) methods using either exact or approximate Hessians with the more recently proposed minimum mode following methods where only the minimum eigenvalue mode is found, either by the dimer or the Lanczos method. A test problem involving transitions in a seven-atom Pt island on a Pt(111) surface using a simple Morse pairwise potential function is used and the number of degrees of freedom varied by varying the number of movable atoms. In the full system, 175 atoms can move so 525 degrees of freedom need to be optimized to find the saddle points. For testing purposes, we have also restricted the number of movable atoms to 7 and 1. Our results indicate that if attempting to make a map of all relevant saddle points for a large system (as would be necessary when simulating the long time scale evolution of a thermal system) the minimum mode following methods are preferred. The minimum mode following methods are also more efficient when searching for the lowest saddle points in a large system, and if the force can be obtained cheaply. However, if only the lowest saddle points are sought and the calculation of the force is expensive but a good approximation for the Hessian at the starting position of the search can be obtained at low cost, then the RFO approaches employing an approximate Hessian represent the preferred choice. For small and medium sized systems where the force is expensive to calculate, the RFO approaches employing an approximate Hessian is also the more efficient, but when the force and Hessian can be obtained cheaply and only the lowest saddle points are sought the RFO approach using an exact Hessian is the better choice. These conclusions have been reached based on a comparison of the total computational effort needed to find the saddle points and the number of saddle points found for each of the methods. The RFO methods do not perform very well with respect to the latter aspect, but starting the searches further away from the initial minimum or using the hybrid RFO version presented here improves this behavior considerably in most cases.     

The curvature of the conical intersection seam: an approximate second-order analysis

We present a method for analyzing the curvature (second derivatives) of the conical intersection hyperline at an optimized critical point. Our method uses the projected Hessians of the degenerate states after elimination of the two branching space coordinates, and is equivalent to a frequency calculation on a single Born-Oppenheimer potential-energy surface. Based on the projected Hessians, we develop an equation for the energy as a function of a set of curvilinear coordinates where the degeneracy is preserved to second order (i.e., the conical intersection hyperline). The curvature of the potential-energy surface in these coordinates is the curvature of the conical intersection hyperline itself, and thus determines whether one has a minimum or saddle point on the hyperline. The equation used to classify optimized conical intersection points depends in a simple way on the first- and second-order degeneracy splittings calculated at these points. As an example, for fulvene, we show that the two optimized conical intersection points of C2v symmetry are saddle points on the intersection hyperline. Accordingly, there are further intersection points of lower energy, and one of C2 symmetry--presented here for the first time--is found to be the global minimum in the intersection space.     

Optimization methods for finding minimum energy paths

A comparison of chain-of-states based methods for finding minimum energy pathways (MEPs) is presented. In each method, a set of images along an initial pathway between two local minima is relaxed to find a MEP. We compare the nudged elastic band (NEB), doubly nudged elastic band, string, and simplified string methods, each with a set of commonly used optimizers. Our results show that the NEB and string methods are essentially equivalent and the most efficient methods for finding MEPs when coupled with a suitable optimizer. The most efficient optimizer was found to be a form of the limited-memory Broyden-Fletcher-Goldfarb-Shanno method in which the approximate inverse Hessian is constructed globally for all images along the path. The use of a climbing-image allows for finding the saddle point while representing the MEP with as few images as possible. If a highly accurate MEP is desired, it is found to be more efficient to descend from the saddle to the minima than to use a chain-of-states method with many images. Our results are based on a pairwise Morse potential to model rearrangements of a heptamer island on Pt(111), and plane-wave based density functional theory to model a rollover diffusion mechanism of a Pd tetramer on MgO(100) and dissociative adsorption and diffusion of oxygen on Au(111).     

The Diels-Alder reaction of ethene and 1,3-butadiene: an extended multireference ab initio investigation.

The concerted and stepwise mechanisms of the Diels-Alder reaction between 1,3-butadiene and ethene have been investigated using highly correlated multireference methods (MRAQCC) and extended basis sets. Full MRAQCC geometry optimizations have been performed in all cases. The best estimate for the energy barrier of the Diels-Alder reaction is 22 kcalmol(-1). Anti- and gauche-out minima for the biradical structures and corresponding fragmentation saddle points have been determined. The biradical anti fragmentation saddle point is located 6.5 kcalmol(-1) above the concerted saddle point. The gauche-in structure does not correspond to a local minimum, but leads on geometry optimization directly to cyclohexene.     

Oxygen island formation on Pt(111) studied by dynamic Monte Carlo simulation

The formation of oxygen islands on the Pt(111) surface has been studied as a function of temperature by low energy electron diffraction (LEED) experiments and dynamic Monte Carlo (DMC) simulations. By raising the temperature, the (2 x 2) LEED spot intensity increases gradually and decays after a peak at around 255 K (T(p)) with full width of half maximum of 160 K. This behavior is interpreted by DMC simulations with the kinematical LEED analysis. In the DMC simulation, an oxygen atom hops to the neighboring site via the activation barrier of the saddle point. The potential energies at initial, saddle, and final points are changed at each hopping event depending on the surrounding oxygen atoms. By comparing the observed T(p) with the simulated one, the interaction energy E of oxygen atoms on Pt(111) was determined to be 25+/-3 meV at 2a(0). The DMC simulations visualize how the oxygen islands are formed and collapse on Pt(111) with increase of the temperature and well reproduce the surface configurations observed by scanning tunneling microscopy.     

Reaction Pathways and Convexity of the Potential Energy Surface: Application of Newton Trajectories

The reaction path is an important concept of theoretical chemistry. We employ the definitions of the intrinsic reaction coordinate (IRC), the gradient extremal (GE), and the Newton trajectory (NT). The usual imagination in chemistry is that a minimum energy path is in a convex region of the potential energy surface. We describe different schemes of convexity to handle the situation. It comes out that NTs are the best ansatz for the problem: NTs, which monotonically increase (or monotonically decrease), are automatically strictly pseudo-convex throughout, and they go throughout along a valley between minimum and saddle point.     

Cyclic C3 structures

Alternative cyclic C₃ structures are not competitive energetically with the linear ¹∑_g⁺ ground state, 1. The D₃h triplet, 4, is a local minimum on the potential energy surface. Cyclopropynylidene, a C_2v singlet (3) is a saddle point for the degenerate isomerization which permutes the order of carbon atoms in 1. The activation energy for this transformation is predicted to be 29 kcal/mole.     

The Structures and Stability of HNOS Isomers

Potential energy surface of HNOS system is investigated by means of MP2 method with 6-311 G(d,p) basis set.The energy for each minimum and saddle point on the potential energy surface is corrected at the QCISD(T)/6-311 G(3df,2p) level of theory with zero-point vibrational energy included.As a result ,eighteen isomers are theoretically predicted and cis-HNSO is found to be global minimum on the potential energy surface,Wherein,fourteen isomers are considered as kinetically stable species,and should be experimentally observed.Comparisons are made for HNOS system with its analogues,HNO2 and NHS2.The nature of bonding and isomers‘ stability of HNOS system are similar to HNS2.The obvious similarities and discrepancies among HNOS,HNO2 and HNS2 are attributed to the hypervalent capacity of sulfur,oxygen and nitrogen atoms.     

Theoretical study of potential energy surface and vibrational spectra of ArF2 system

An ab initio potential energy surface (PES) of ArF2 system has been obtained by using MP4 calculation with a large basis set including bond functions. There are two local minimums on the PES: one is T-shaped and the other is L-shaped. The L-shaped minimum is the global minimum with a well depth of -119.62 cm- 1 at R = 0.3883nm. The T-shaped minimum has a well depth of -85.93cm -1 at R = 0.3486 nm. A saddle point is found at R = 0.3486 and θ = 61° with the well depth of -61.53 cm-1. The vibrational energy levels have been calculated by using VSCF-CI method. The results show that this PES supports 27 vibrational bound states, and the ground states are two degenerate states assigned to the L-type vibration.     

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.     

Structures and vibrational frequencies of the A′and A″ states of the CH3O radical

Analytic gradient and second-derivative methods are used to obtain the equilibrium geometries and vibrational frequencies of the methoxy radical CH₃O. The analytic methods are found to be superior in both speed and accuracy to numerical techniques which for this molecule are qualitatively and quantitatively incorrect. The A′ state is found to be a minimum on the potential-energy surface, whereas the A″ state is a saddle point.     

Potential energy surface of the (H2)2 dimer: an MP2 study

Fifteen structures of the (H₂)₂ dimer have been investigated at the MP2/[4s3p] level. The SCF and MP2 (2nd order Møller-Plesser treatment) interaction energies have been corrected for the basis set superposition error. Only the T-shaped structure has been established as a minimum on the potential energy surface. Two equivalent T-shaped structures are connected by a saddle point with a rhomboid structure.Dedicated to Professor J. Koutecký on the occasion of his 65th birthday     

A method for finding the ridge between saddle points applied to rare event rate estimates

A method is presented for finding the ridge between first order saddle points on a multidimensional surface. For atomic scale systems, such saddle points on the energy surface correspond to atomic rearrangement mechanisms. Information about the ridge can be used to test the validity of the harmonic approximation to transition state theory, in particular to verify that second order saddle points--maxima along the ridge--are high enough compared to the first order saddle points. New minima along the ridge can also be identified during the path optimisation, thereby revealing additional transition mechanisms. The method is based on a string of discretisation points along a path between the first order saddle points and using an iterative optimisation which requires only the force acting on the atoms. At each iteration during the optimisation, the force is inverted along an unstable eigenmode perpendicular to the path. The method is applied to Al adatom diffusion on the Al(100) surface to find the ridge between 2-, 3- and 4-atom concerted displacements and hop mechanisms. A correction to the harmonic approximation of transition state theory was estimated by direct evaluation of the configuration integral along the ridge.     

Analysis of the concept of minimum energy path on the potential energy surface of chemically reacting systems

Some confusion regarding the properties of minimum energy paths is evident in the literature. We show that a way of steepest descent on a potential surface can be defined independently upon the choice of the coordinate systems. The result is applied to mass-weighted coordinates and their use is critically reviewed. Fukui's IRC appears to be a special case of the steepest descent path starting from a saddle point. The impossibility to define a general ascent path is illustrated and the relations of IRC to real trajectories are discussed.     

Effect of a complex formation on the calculated low-pressure rate constant of a bimolecular gas-phase reaction governed by tunneling

Many important bimolecular hydrogen-transfer processes that take place in the atmosphere proceed via a potential energy minimum (hydrogen-bonded complex) that precedes along the minimum energy path the unique saddle point of the reaction, the one corresponding to the hydrogen transfer. It is clear that the one-step low-pressure rate constant of such a reaction does not depend on the existence of any complex along the minimum energy path below the reactant if the reaction takes place by thermal activation over a transition state that lies quite above the reactants (for instance 10 kcal/mol). However, we have quantitatively shown in this article that the scenario notoriously changes if the reaction involves significant tunneling. In this work, we have theoretically calculated the rate constants and their temperature dependence for the reaction HO+HOH→HOH+OH by means of a canonical variational transition state theory and a canonical unified statistical theory (when necessary). Multidimensional tunneling effects have been included with a semiclassical transmission coefficient. Two kinds of modified potential energy surfaces (PESs), obtained from an original ab initio potential energy surface, previously calculated by us, have been used. The Eckart-modified PESs serve to model the hydrogen-abstraction profiles with no complexes along the path, while the Gaussian-modified PESs model the energy profiles with two complexes along the path symmetrically distributed at each side of the abstraction saddle point. Our results show that the existence of those complexes reduces the thickness of the classically forbidden region for energies below the adiabatic barrier, and then tunneling is promoted and the reaction is accelerated. The effect of the complex formation in several kinetic magnitudes, as the Arrhenius parameters and the kinetic isotope effect has also been analyzed. © 1999 John Wiley & Sons, Inc. J Comput Chem 20: 1685–1692, 1999     

Searching for saddle points of potential energy surfaces by following a reduced gradient

The old coordinate driving procedure to find transition structures in chemical systems is revisited. The well-known gradient criterion, ∇E( x )= 0 , which defines the stationary points of the potential energy surface (PES), is reduced by one equation corresponding to one search direction. In this manner, abstract curves can be defined connecting stationary points of the PES. Starting at a given minimum, one follows a well-selected coordinate to reach the saddle of interest. Usually, but not necessarily, this coordinate will be related to the reaction progress. The method, called reduced gradient following (RGF), locally has an explicit analytical definition. We present a predictor–corrector method for tracing such curves. RGF uses the gradient and the Hessian matrix or updates of the latter at every curve point. For the purpose of testing a whole surface, the six-dimensional PES of formaldehyde, H₂CO, was explored by RGF using the restricted Hartree–Fock (RHF) method and the STO-3G basis set. Forty-nine minima and saddle points of different indices were found. At least seven stationary points representing bonded structures were detected in addition to those located using another search algorithm on the same level of theory. Further examples are the localization of the saddle for the HCN⇌CNH isomerization (used for steplength tests) and for the ring closure of azidoazomethine to 1H-tetrazole. The results show that following the reduced gradient may represent a serious alternative to other methods used to locate saddle points in quantum chemistry. © 1998 John Wiley & Sons, Inc. J Comput Chem 19: 1087–1100, 1998     

A spline for your saddle

Pinpointing extrema on a multidimensional hypersurface is an important generic problem with a broad scope of application in statistical mechanics, biophysics, chemical reaction dynamics, and quantum chemistry. Local minima of the hypersurface correspond to metastable structures and are usually the most important points to look for. They are relatively easy to find using standard minimizing algorithms. A considerably more difficult task is the location of saddle points. The saddle points most sought for are those which form the lowest barriers between given minima and are usually required for determining rates of rare events. We formulate a path functional minimum principle for the saddle point. We then develop a cubic spline method for applying this principle and locating the saddle point(s) separating two local minima on a potential hypersurface. A quasi-Newton algorithm is used for minimization. The algorithm does not involve second derivatives of the hypersurface and the number of potential gradients evaluated is usually less than 10% of the number of potential evaluations. We demonstrate the performance of the method on several standard examples and on a concerted exchange mechanism for self-diffusion in diamond. Finally, we show that the method may be used for solving large constrained minimization problems which are relevant for self-consistent field iterations in large systems.     

Potential energy surface and molecular dynamics of MbNO: existence of an unsuspected FeON minimum

Ligands such as CO, O(2), or NO are involved in the biological function of myoglobin. Here we investigate the energetics and dynamics of NO interacting with the Fe(II) heme group in native myoglobin using ab initio and molecular dynamics simulations. At the global minimum of the ab initio potential energy surface (PES), the binding energy of 23.4 kcal/mol and the Fe-NO structure compare well with the experimental results. Interestingly, the PES is found to exhibit two minima: There exists a metastable, linear Fe-O-N minimum in addition to the known, bent Fe-N-O global minimum conformation. Moreover, the T-shaped configuration is found to be a saddle point, in contrast to the corresponding minimum for NO interacting with Fe(III). To use the ab initio results for finite temperature molecular dynamics simulations, an analytical function was fitted to represent the Fe-NO interaction. The simulations show that the secondary minimum is dynamically stable up to 250 K and has a lifetime of several hundred picoseconds at 300 K. The difference in the topology of the heme-NO PES from that assumed previously (one deep, single Fe-NO minimum) suggests that it is important to use the full PES for a quantitative understanding of this system. Why the metastable state has not been observed in the many spectroscopic studies of myoglobin interacting with NO is discussed, and possible approaches to finding it are outlined.     

Atomistic simulations of rare events using gentlest ascent dynamics

The dynamics of complex systems often involve thermally activated barrier crossing events that allow these systems to move from one basin of attraction on the high dimensional energy surface to another. Such events are ubiquitous, but challenging to simulate using conventional simulation tools, such as molecular dynamics. Recently, E and Zhou [Nonlinearity 24(6), 1831 (2011)] proposed a set of dynamic equations, the gentlest ascent dynamics (GAD), to describe the escape of a system from a basin of attraction and proved that solutions of GAD converge to index-1 saddle points of the underlying energy. In this paper, we extend GAD to enable finite temperature simulations in which the system hops between different saddle points on the energy surface. An effective strategy to use GAD to sample an ensemble of low barrier saddle points located in the vicinity of a locally stable configuration on the high dimensional energy surface is proposed. The utility of the method is demonstrated by studying the low barrier saddle points associated with point defect activity on a surface. This is done for two representative systems, namely, (a) a surface vacancy and ad-atom pair and (b) a heptamer island on the (111) surface of copper.