On the automatic restricted-step rational-function-optimization method

The rational function optimization algorithm is one of the widely used methods to search stationary points on surfaces. However, one of the drawbacks of this method is the step reduction procedure to deal with the overstepping problem. We present and comment on a method such that the step obtained from the solution of the rational function equations possesses the desired correct length. The analysis and discussion of the method is mainly centered on the location and optimization of transition states. Received: 18 June 1998 / Accepted: 17 September 1998 / Published online: 23 November 1998     

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.     

Quadratic string method for determining the minimum-energy path based on multiobjective optimization

Based on a multiobjective optimization framework, we develop a new quadratic string method for finding the minimum-energy path. In the method, each point on the minimum-energy path is minimized by integration in the descent direction perpendicular to path. Each local integration is done on a quadratic surface approximated by a damped Broyden-Fletcher-Goldfarb-Shanno updated Hessian, allowing the algorithm to take many steps between energy and gradient calls. The integration is performed with an adaptive step-size solver, which is restricted in length to the trust radius of the approximate Hessian. The full algorithm is shown to be capable of practical superlinear convergence, in contrast to the linear convergence of other methods. The method also eliminates the need for predetermining such parameters as step size and spring constants, and is applicable to reactions with multiple barriers. The effectiveness of this method is demonstrated for the Muller-Brown potential, a seven-atom Lennard-Jones cluster, and the enolation of acetaldehyde to vinyl alcohol.     

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.     

A growing string method for determining transition states: comparison to the nudged elastic band and string methods

Interpolation methods such as the nudged elastic band and string methods are widely used for calculating minimum energy pathways and transition states for chemical reactions. Both methods require an initial guess for the reaction pathway. A poorly chosen initial guess can cause slow convergence, convergence to an incorrect pathway, or even failed electronic structure force calculations along the guessed pathway. This paper presents a growing string method that can find minimum energy pathways and transition states without the requirement of an initial guess for the pathway. The growing string begins as two string fragments, one associated with the reactants and the other with the products. Each string fragment is grown separately until the fragments converge. Once the two fragments join, the full string moves toward the minimum energy pathway according to the algorithm for the string method. This paper compares the growing string method to the string method and to the nudged elastic band method using the alanine dipeptide rearrangement as an example. In this example, for which the linearly interpolated guess is far from the minimum energy pathway, the growing string method finds the saddle point with significantly fewer electronic structure force calculations than the string method or the nudged elastic band method.     

Superlinearly converging dimer method for transition state search

Algorithmic improvements of the dimer method [G. Henkelman and H. Jonsson, J. Chem. Phys. 111, 7010 (1999)] are described in this paper. Using the limited memory Broyden-Fletcher-Goldfarb-Shanno (L-BFGS) optimizer for the dimer translation greatly improves the convergence compared to the previously used conjugate gradient algorithm. It also saves one energy and gradient calculation per dimer iteration. Extrapolation of the gradient during repeated dimer rotations reduces the computational cost to one gradient calculation per dimer rotation. The L-BFGS algorithm also improves convergence of the rotation. Thus, three to four energy and gradient evaluations are needed per iteration at the beginning of a transition state search, while only two are required close to convergence. Moreover, we apply the dimer method in internal coordinates to reduce coupling between the degrees of freedom. Weighting the coordinates can be used to apply chemical knowledge about the system and restrict the transition state search to only part of the system while minimizing the remainder. These improvements led to an efficient method for the location of transition states without the need to calculate the Hessian. Thus, it is especially useful in large systems with expensive gradient evaluations.     

Mode‐tracking based stationary‐point optimization

In this work, we present a transition‐state optimization protocol based on the Mode‐Tracking algorithm [Reiher and Neugebauer, J. Chem. Phys., 2003, 118, 1634]. By calculating only the eigenvector of interest instead of diagonalizing the full Hessian matrix and performing an eigenvector following search based on the selectively calculated vector, we can efficiently optimize transition‐state structures. The initial guess structures and eigenvectors are either chosen from a linear interpolation between the reactant and product structures, from a nudged‐elastic band search, from a constrained‐optimization scan, or from the minimum‐energy structures. Alternatively, initial guess vectors based on chemical intuition may be defined. We then iteratively refine the selected vectors by the Davidson subspace iteration technique. This procedure accelerates finding transition states for large molecules of a few hundred atoms. It is also beneficial in cases where the starting structure is very different from the transition‐state structure or where the desired vector to follow is not the one with lowest eigenvalue. Explorative studies of reaction pathways are feasible by following manually constructed molecular distortions. © 2015 Wiley Periodicals, 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.     

Efficiency considerations in the construction of interpolated potential energy surfaces for the calculation of quantum observables by diffusion Monte Carlo

A modified Shepard interpolation scheme is used to construct global potential energy surfaces (PES) in order to calculate quantum observables--vibrationally averaged internal coordinates, fully anharmonic zero-point energies and nuclear radial distribution functions--for a prototypical loosely bound molecular system, the water dimer. The efficiency of PES construction is examined with respect to (a) the method used to sample configurational space, (b) the method used to choose which points to add to the PES data set, and (c) the use of either a one- or two-part weight function. The most efficient method for constructing the PES is found to require a quantum sampling regime, a combination of both h-weight and rms methods for choosing data points and use of the two-part weight function in the interpolation. Using this regime, the quantum diffusion Monte Carlo zero-point energy converges to the exact result within addition of 50 data points. The vibrationally averaged O-O distance and O-O radial distribution function, however, converge more slowly and require addition of over 500 data points. The methods presented here are expected to be applicable to both other loosely bound complexes as well as tightly bound molecular species. When combined with high quality ab initio calculations, these methods should be able to accurately characterize the PES of such species.     

Excited state tracking during the relaxation of coordination compounds

The ability to locate minima on electronic excited states (ESs) potential energy surfaces both in the case of bright and dark states is crucial for a full understanding of photochemical reactions. This task has become a standard practice for small- to medium-sized organic chromophores thanks to the constant developments in the field of computational photochemistry. However, this remains a very challenging effort when it comes to the optimization of ESs of transition metal complexes (TMCs), not only due to the presence of several electronic ESs close in energy, but also due to the complex nature of the ESs involved. In this article, we present a simple yet powerful method to follow an ES of interest during a structural optimization in the case of TMCs, based on the use of a compact hole-particle representation of the electronic transition, namely the natural transition orbitals (NTOs). State tracking using NTOs is unambiguously accomplished by computing the mono-electronic wave function overlap between consecutive steps of the optimization. Here, we demonstrate that this simple but robust procedure works not only in the case of the cytosine but also in the case of the ES optimization of a ruthenium nitrosyl complex which is very problematic with standard approaches. © 2019 The Authors. Journal of Computational Chemistry published by Wiley Periodicals, Inc.     

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.     

How to find a saddle point

This paper explains a method for finding saddle points on a multidimensional surface and shows how it may be used to define saddle-point seeking curves that have properties similar to well-known orthogonal trajectories. It is shown that a gradient extremal is a special case of one of these curves, and its chemical significance as the path, defined by local criteria, which starts from a stable structure and leads to a transition state, is discussed briefly in relation to the intrinsic reaction coordinate. It is emphasized that this theory gives a natural method for locating points that have Hessians of similar structure to those of transition states.     

Integration of chemical stiff ODEs using exponential propagation method

In this paper, we study the numerical long time integration of large stiff systems of differential equations arising from chemical reactions by exponential propagation methods. These methods, which typically converge faster, use matrix-vector products with the exponential or other related function of the Jacobian that can be effectively approximated by Krylov sabspace methods. We equip these methods to an automatic stepsize control technique and apply the method of order 4 for numerical integration of some famous stiff chemical problems such as Belousov-Zhabotinskii reaction, the Chapman atmosphere, Hydrogen chemistry, chemical Akzo-Nobel problem and air pollution problem.     

Single‐ended transition state finding with the growing string method

Reaction path finding and transition state (TS) searching are important tasks in computational chemistry. Methods that seek to optimize an evenly distributed set of structures to represent a chemical reaction path are known as double‐ended string methods. Such methods can be highly reliable because the endpoints of the string are fixed, which effectively lowers the dimensionality of the reaction path search. String methods, however, require that the reactant and product structures are known beforehand, which limits their ability for systematic exploration of reactive steps. In this article, a single‐ended growing string method (GSM) is introduced which allows for reaction path searches starting from a single structure. The method works by sequentially adding nodes along coordinates that drive bonds, angles, and/or torsions to a desired reactive outcome. After the string is grown and an approximate reaction path through the TS is found, string optimization commences and the exact TS is located along with the reaction path. Fast convergence of the string is achieved through use of internal coordinates and eigenvector optimization schemes combined with Hessian estimates. Comparison to the double‐ended GSM shows that single‐ended method can be even more computationally efficient than the already rapid double‐ended method. Examples, including transition metal reactivity and a systematic, automated search for unknown reactivity, demonstrate the efficacy of the new method. This automated reaction search is able to find 165 reaction paths from 333 searches for the reaction of NH₃BH₃ and (LiH)₄, all without guidance from user intuition. © 2015 Wiley Periodicals, Inc.     

Reliable and efficient reaction path and transition state finding for surface reactions with the growing string method

The computational challenge of fast and reliable transition state and reaction path optimization requires new methodological strategies to maintain low cost, high accuracy, and systematic searching capabilities. The growing string method using internal coordinates has proven to be highly effective for the study of molecular, gas phase reactions, but difficulties in choosing a suitable coordinate system for periodic systems has prevented its use for surface chemistry. New developments are therefore needed, and presented herein, to handle surface reactions which include atoms with large coordination numbers that cannot be treated using standard internal coordinates. The double‐ended and single‐ended growing string methods are implemented using a hybrid coordinate system, then benchmarked for a test set of 43 elementary reactions occurring on surfaces. These results show that the growing string method is at least 45% faster than the widely used climbing image‐nudged elastic band method, which also fails to converge in several of the test cases. Additionally, the surface growing string method has a unique single‐ended search method which can move outward from an initial structure to find the intermediates, transition states, and reaction paths simultaneously. This powerful explorative feature of single ended‐growing string method is demonstrated to uncover, for the first time, the mechanism for atomic layer deposition of TiN on Cu(111) surface. This reaction is found to proceed through multiple hydrogen‐transfer and ligand‐exchange events, while formation of H‐bonds stabilizes intermediates of the reaction. Purging gaseous products out of the reaction environment is the driving force for these reactions. © 2017 Wiley Periodicals, Inc.     

Algorithmic tools in the study of semiempirical potential surfaces

The most efficient optimization methods implemented in the semiempirical package AMPAC are presented. They concern the minimization of the energy or of the gradient norm by either pseudo-Newton or quadratic procedures, the search for transition states, and the intrinsic reaction coordinate in conjunction with variational transition-state theories. Nonlocal methods such as simulated annealing are also introduced. © 1992 John Wiley & Sons, Inc.     

Spectral convergence of the quadrature discretization method in the solution of the Schrodinger and Fokker-Planck equations: comparison with sinc methods

Spectral methods based on nonclassical polynomials and Fourier basis functions or sinc interpolation techniques are compared for several eigenvalue problems for the Fokker-Planck and Schrodinger equations. A very rapid spectral convergence of the eigenvalues versus the number of quadrature points is obtained with the quadrature discretization method (QDM) and the appropriate choice of the weight function. The QDM is a pseudospectral method and the rate of convergence is compared with the sinc method reported by Wei [J. Chem. Phys., 110, 8930 (1999)]. In general, sinc methods based on Fourier basis functions with a uniform grid provide a much slower convergence. The paper considers Fokker-Planck equations (and analogous Schrodinger equations) for the thermalization of electrons in atomic moderators and for a quartic potential employed to model chemical reactions. The solution of the Schrodinger equation for the vibrational states of I2 with a Morse potential is also considered.     

Interpolated potential energy surfaces: How accurate do the second derivatives have to be?

A global potential energy surface for the water dimer is constructed using the modified Shepard interpolation scheme of Collins et al. According to this interpolation scheme, the energy at an arbitrary geometry is expressed as a weighted sum of Taylor series expansions from neighboring data points, where the energy and derivative data required are obtained from ab initio calculations. For some ab initio methods, errors are introduced into the second derivative matrix, either by numerical differencing of ab initio energies or numerical integration during the ab initio calculation. Therefore, we test the accuracy required of the second derivative data by truncation of the exact second derivatives to a series of approximate second derivatives, and assess the effect on the results of a quantum diffusion Monte Carlo (QDMC) simulation. Our results show that the calculated zero-point energy and wave function histograms converge to within the numerical uncertainty of the QDMC simulation by inclusion of either three significant figures or three decimal places in the second derivatives.     

Smooth solvation method for d-orbital semiempirical calculations of biological reactions. 1. Implementation

The present paper describes the extension of a recently developed smooth conductor-like screening model for solvation to a d-orbital semiempirical framework (MNDO/d-SCOSMO) with analytic gradients that can be used for geometry optimizations, transition state searches, and molecular dynamics simulations. The methodology is tested on the potential energy surfaces for separating ions and the dissociative phosphoryl transfer mechanism of methyl phosphate. The convergence behavior of the smooth COSMO method with respect to discretization level is examined and the numerical stability of the energy and gradient are compared to that from conventional COSMO calculations. The present method is further tested in applications to energy minimum and transition state geometry optimizations of neutral and charged metaphosphates, phosphates, and phosphoranes that are models for stationary points in transphosphorylation reaction pathways of enzymes and ribozymes. The results indicate that the smooth COSMO method greatly enhances the stability of quantum mechanical geometry optimization and transition state search calculations that would routinely fail with conventional solvation methods. The present MNDO/d-SCOSMO method has considerable computational advantages over hybrid quantum mechanical/molecular mechanical methods with explicit solvation, and represents a potentially useful tool in the arsenal of multi-scale quantum models used to study biochemical reactions.     

New optimization method for conformational energy calculations on polypeptides: Conformational space annealing

A new optimization method is presented to search for the global minimum-energy conformations of polypeptides. The method combines essential aspects of the build-up procedure and the genetic algorithm, and it introduces the important concept of “conformational space annealing.” Instead of considering a single conformation, attention is focused on a population of conformations while new conformations are obtained by modifying a “seed conformation.” The annealing is carried out by introducing a distance cutoff, D_cut, which is defined in the conformational space; D_cut effectively divides the whole conformational space of local minima into subdivisions. The value of D_cut is set to a large number at the beginning of the algorithm to cover the whole conformational space, and annealing is achieved by slowly reducing it. Many distinct local minima designed to be distributed as far apart as possible in conformational space are investigated simultaneously. Therefore, the new method finds not only the global minimum-energy conformation but also many other distinct local minima as by-products. The method is tested on Met-enkephalin, a 24-dihedral angle problem. For all 100 independent runs, the accepted global minimum-energy conformation was obtained after about 2600 minimizations on average. © 1997 John Wiley & Sons, Inc. J Comput Chem 18: 1222–1232