New rules on potential surface topology and critical point search

The topology of potential energy surfaces provides a unified framework for the study of individual molecular properties, all conformational changes as well as chemical reactions. Molecular behavior, electronic and vibrational properties, conformational freedom, reactivity bond formation and bond breaking are all energy dependent, and the potential energy surface approach provides an elegant, conceptually convenient, although rather complicated representation of this energy dependence. Topology as a mathematical tool is exceptionally suitable for the extraction of the most essential features of complicated representations. By applying topological methods for potential surface analysis, a new, global perspective of many aspects of chemistry emerges. Some of these topological results also have important practical, computational significance. A family of new topological rules and symmetry relations will be adapted for applications in low dimensional relaxed cross-sections of configuration spaces, with a special emphasis on their role in the search for critical points, primarily energy minima and saddle points of transition structures of potential energy surfaces and hypersurfaces.     

Simulated annealing to locate various stationary points in semiempirical methods

In this work, an algorithm was developed to study the potential energy surfaces in the coordinate spaces of molecules by a nonlocal way, in contrast to classic energy minimizers as the BFGS or the DFP method. This algorithm, based on the specificities of semiempirical methods, mixes simulated annealing and local searches to reduce computation costs. By this technique, the global energy minimum can be localized. Moreover, local minima that are close in energy to the global minimum are also obtained. If the search is not only for minima but for all stationary points (minima, saddle points…), then the energy is replaced by the gradient norm, which reaches its minimum values at stationary points. The annealing process is stopped before having accurately reached the global minimum and generates a list of geometries whose energies (respectively, whose gradients) are optimized by local minimizers. This list of geometries is shortened from the nearly equivalent geometries by a dynamic single-clustering analysis. The energy/gradient local minimizers act on the clustered list to produce a set of minima/stationary points. A targeted search of these points and reduction of the costs are reached by the way of several penalty functions. They eliminate—without energy calculation—most of the points generated by random walks on the potential energy surface. These penalty functions (on the total moment of inertia or on interatomic distances) are specific to the class of problem studied. They account for the nonrupture of either specific chemical bonds or rings in cyclic molecules, they assure that molecular systems are kept bonded, and they avoid the collapsing of atoms. © 1992 John Wiley & Sons, Inc.     

Accelerated molecular dynamics: a promising and efficient simulation method for biomolecules

Many interesting dynamic properties of biological molecules cannot be simulated directly using molecular dynamics because of nanosecond time scale limitations. These systems are trapped in potential energy minima with high free energy barriers for large numbers of computational steps. The dynamic evolution of many molecular systems occurs through a series of rare events as the system moves from one potential energy basin to another. Therefore, we have proposed a robust bias potential function that can be used in an efficient accelerated molecular dynamics approach to simulate the transition of high energy barriers without any advance knowledge of the location of either the potential energy wells or saddle points. In this method, the potential energy landscape is altered by adding a bias potential to the true potential such that the escape rates from potential wells are enhanced, which accelerates and extends the time scale in molecular dynamics simulations. Our definition of the bias potential echoes the underlying shape of the potential energy landscape on the modified surface, thus allowing for the potential energy minima to be well defined, and hence properly sampled during the simulation. We have shown that our approach, which can be extended to biomolecules, samples the conformational space more efficiently than normal molecular dynamics simulations, and converges to the correct canonical distribution.     

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.     

Stereodynamics of triethylamine: Molecular mechanics calculations on the direct rotational racemization of the C3-symmetric conformers

The direct interconversion of the two C₃-symmetric enantiomeric conformations of triethylamine, via C? N bond rotation, has been studied by molecular mechanics (MM 2) calculations. The MM 2 calculations have been used to characterize the minima (equilibrium geometries) and first-order saddle points (transition states) for this process. For one interconversion, there are five saddle points and six minima. The highest energy saddle point results from the uncoupled rotation of one ethyl group to eclipse the lone pair. Two of the barriers result from coupled rotation of two ethyl groups in close passage.     

The intricacies of the stacking interaction in a pyrrole–pyrrole system

Extensive calculations of potential energy surfaces for parallel-displaced configurations of pyrrole–pyrrole systems have been carried out by the use of a dispersion-corrected density functional. System geometries associated with the energy minima have been found. The minimum interaction energy has been calculated as ?5.38 kcal/mol. However, bonding boundaries appeared to be relatively broad, and stacking interactions can be binding even for ring centroid distances larger than 6 Å. Though the contribution of the correlation energy to intermolecular interaction in pyrrole dimers appeared to be relatively small (around 1.6 smaller than it is in a benzene–benzene system), this system’s minimum interaction energy is lower than those calculated for benzene–benzene, benzene–pyridine and even pyridine–pyridine configurations. The calculation of the charges and energy decomposition analysis revealed that the specific charge distribution in a pyrrole molecule and its relatively high polarization are the significant source of the intermolecular interaction in pyrrole dimer systems.     

Following the streambed reaction on potential-energy surfaces: a new robust method

A simple procedure with low computational efforts is proposed to follow the reaction path of the potential-energy hypersurface (PES) starting from minima or saddle points. The method uses a modification of the so-called “following the reduced gradient” [Quapp W, Hirsch M, Imig O, Heidrich D (1998) J Comput Chem 19:1087]. The original method connects points where the gradient has a constant direction. In the present article the procedure is replaced by taking iterative varying directions of the gradient controlled by the last tangent of the searched curve. The resulting minimum energy path is that valley floor gradient extremal (GE) which belongs to the smallest (absolute) eigenvalue of the Hessian and, hence, that GE which usually leads along the streambed of a chemical reaction. The new method avoids third derivatives of the PES and obtains the GE of least ascent by second-order calculations only. Nevertheless, we are able to follow the streambed GE uphill or downhill. We can connect a minimum with its saddles if the streambed leads up to a saddle, or we find a turning point or a bifurcation point. The effectiveness and the characteristic properties of the new algorithm are demonstrated by using polynomial test surfaces, an ab initio PES of H₂O, and the analytic potentials of Lennard-Jones (LJ) clusters. By tracing the streambeds we located previously identified saddle points for LJ _N with N=3, 7, 8, and 55. Saddles for LJ _N with N=15, 20, and 30 as presented here are new results. Received: 8 March 2000 / Accepted: 17 July 2000 / Published online: 24 October 2000     

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     

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.     

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.     

Stationary points on the H2CO potential energy surface: dependence on theoretical level

A total of 36 stationary points have been located on the H₂CO potential energy surface by means of gradient extremal following. These 36 points are believed to represent all the important stationary points on this surface. There is no indication that the structure of the surface becomes less complicated as the size of the basis set is enlarged at the Hartree-Fock level of theory, but many of the second- and third-order saddle points disappear when electron correlation is introduced. Of the ten first-order saddle points (transition structures) located, the majority have reaction paths entering the associated minima in a side-on approach, i.e. these cannot be located by uphill walking from the minimum. Received: 5 February 1998 / Accepted: 21 May 1998 / Published online: 29 July 1998     

Lowest singlet and triplet potential energy surfaces of S2N2

Forty four stationary points have been located on the lowest singlet and triplet potential energy surfaces of S(2)N(2). Ten minima and ten saddle points on the lowest singlet surface and eleven minima and thirteen saddle points on the lowest triplet surface were found. All saddle points were connected to minima or lower-order saddle points by following the intrinsic reaction coordinate. Renner-Teller effects in the linear isomers were studied by examining their bending curves. The S(2)N(2) polymerization mechanism was investigated by first locating the transition state corresponding to ring opening and then considering all species connected to it that are close in energy. The commonly accepted mechanism is problematic due to the number of species that would lead to dissociation to SN + SN. Other possible isomers that are consistent with the experimental evidence but do not connect to SN radicals in the dissociation limit were examined. A mechanism of polymerization to (SN)(x)() is proposed that involves excitation of the square planar singlet molecule to the triplet surface. The triplet species then undergoes a puckering, and polymerization occurs in a direction approximately perpendicular to the S(2)N(2) plane. Consideration of the predicted vibrational frequencies suggests the structure of the second isomer of S(2)N(2). This isomer has a trans-NSSN structure with a long SS bond. The energetics of trans-NSSN are consistent with the observed temperature effects in the dimerization of SN. Analysis of the bending curves of linear NSSN and NSNS indicates that trans-NSSN is the only isomer which has a small yet significant barrier to that dimerization.     

Exploring the similarities between potential smoothing and simulated annealing

Simulated annealing and potential function smoothing are two widely used approaches for global energy optimization of molecular systems. Potential smoothing as implemented in the diffusion equation method has been applied to study partitioning of the potential energy surface (PES) for N‐Acetyl‐Ala‐Ala‐N‐Methylamide (CDAP) and the clustering of conformations on deformed surfaces. A deformable version of the united‐atom OPLS force field is described, and used to locate all local minima and conformational transition states on the CDAP surface. It is shown that the smoothing process clusters conformations in a manner consistent with the inherent structure of the undeformed PES. Smoothing deforms the original surface in three ways: structural shifting of individual minima, merging of adjacent minima, and energy crossings between unrelated minima. A master equation approach and explicit molecular dynamics trajectories are used to uncover similar features in the equilibrium probability distribution of CDAP minima as a function of temperature. Qualitative and quantitative correlations between the simulated annealing and potential smoothing approaches to enhanced conformational sampling are established. © 2000 John Wiley & Sons, Inc. J Comput Chem 21: 531–552, 2000     

Lower and upper bounds for the number of critical points on energy hypersurfaces

Lower and upper bounds are derived for the number of critical points of all indices λ(λ=0, ?n), for potential energy hypersurfaces defined over a subset S of the nuclear configuration space ⁿR. The number of minima (λ = 0) and saddle points of transition states (λ = 1) are important in analysis of conformational changes and chemical reactions.     

How good is a Broyden–Fletcher–Goldfarb–Shanno-like update Hessian formula to locate transition structures? Specific reformulation of Broyden–Fletcher–Goldfarb–Shanno for optimizing saddle points

Based on a study of the Broyden–Fletcher–Goldfarb–Shanno (BFGS) update Hessian formula to locate minima on a hypersurface potential energy, we present an updated Hessian formula to locate and optimize saddle points of any order that in some sense preserves the initial structure of the BFGS update formula. The performance and efficiency of this new formula is compared with the previous updated Hessian formulae such as the Powell and MSP ones. We conclude that the proposed update is quite competitive but no more efficient than the normal updates normally used in any optimization of saddle points. © 1998 John Wiley & Sons, Inc. J Comput Chem 19: 349–362, 1998     

Density functional theory studies on the reagent Ph3PBr2

We report density functional theory geometry optimizations at the B3LYP/6-311G(d,p) level of theory for the title reagent. Four stationary points on the molecular potential energy surface were located and characterized. Three of these stationary points are energy minima, one a saddle point. The minima correspond to the conventional Ph3PBr2 (three-fold Br-P-Br axis with twisted phenyl rings), the ion-pair [Ph3PBr]+Br- and a four-coordinated Ph3PBr2 spoke structure that can best be described as charge transfer on account of the substantial charge transfer from the Ph3P fragment to Br2 (as determined by a standard Mulliken population analysis and other considerations). The particular saddle point found corresponds to a three-fold Br-P-Br structure with coplanar phenyl rings. Single point B3LYP/6-311+g(3d,2p) calculations were done at the stationary point geometries in order to investigate possible deficiencies in the basis set. Solvent effects for the three solvents water, dichloroethane and cyclohexane were modelled using the self consistent reaction field Onsager method at the single point B3LYP/6-311+g(3d,2p) level of theory. In the gas phase, the charge transfer complex is the most stable of the four; in solution it is the least stable.     

Correlated ab initio study of the ground electronic state of the O2-HF complex

In this paper, we present the first correlated ab initio investigations on the ground electronic state of the O(2)-HF complex. Calculations were performed using the CCSD(T) method with the aug-cc-pVDZ and aug-cc-pVTZ basis sets. The results show that there are two equivalent minimum energy hydrogen-bonded structures of planar bent geometry, where the minima correspond to exchange of the oxygen atoms. For each minimum the length of the O-H hydrogen bond is 2.16 A. The best calculated value of D(e) of the equivalent minima is 271 cm(-1). The T-shaped geometry of the complex, with oxygen perpendicular to the axis connecting the center of masses of O(2) and the HF molecule, represents a barrier to tunneling between the equivalent minima. The best estimated value of that barrier height is 217 cm(-1). The linear O-O-HF geometry of the complex represents a saddle point. The calculated geometrical parameters of the minimum energy structure of the complex are in reasonable agreement with the previously reported spectroscopic results. However, results of the current calculations suggest that a full understanding of the fine structures of the observed infrared spectrum of the complex requires the development of an effective Hamiltonian that takes the effects of tunneling into account.     

Conformational analysis and flexibility of carbohydrates using the CICADA approach with MM3

The potential energy hypersurfaces (PES) of several carbohydrate molecules were studied with a new algorithm for conformational searches, CICADA (Channels in Conformational Space Analyzed by Driver Approach) interfaced with the molecular mechanics program MM3(92). The method requires (1) one or a few low-energy conformations as starting points; and (2) designation of the torsion angles important for understanding the conformational behavior of the molecule. The PES is explored by driving separately each selected torsion angle (in both directions) with a concomitant full-geometry optimization at each increment (except for the driven angle). When a minimum has been detected, the molecule is freely optimized, and the minima so detected are then stored if not encountered previously. The new minima serve as starting structures for further explorations. The results from CICADA permit prediction of relative and absolute flexibility and conformational softness for both the entire molecule as well as for individual group rotations and local minima. The carbohydrates analyzed were Me-α-D -glucopyranoside, β-D -GlcNAc(1-2)α-D -Man, and α-D -GalNAc(1-3)[α-L -Fuc(1-2)]Gal-O-Me. All the low-energy conformers along with the transition states and flexibilities features were characterized. CICADA found all minima and low-energy conversion pathways for the disaccharide that were found by a traditional grid search. In contrast to the grid search method, CICADA concentrates mostly on the exploration of the low-energy regions of the PES, thereby saving a significant amount of computational time. The performance of the method opens new routes for exploring conformational space of larger molecules, such as oligosaccharides. © 1995 by John Wiley & Sons, Inc.