Mode-Tracking Based Stationary-Point Optimization
MMode-Tracking Based Stationary-Point Optimization
Maike Bergeler † , Carmen Herrmann ‡∗ and Markus Reiher †∗ November 6, 2018
Abstract
In this work we present a transition-state optimization protocol based on the Mode-Tracking algorithm [
J. Chem. Phys. , , , 1634]. By calculating only the eigen-vector of interest instead of diagonalizing the full Hessian matrix and performing aneigenvector following search based on the selectively calculated vector, we can efficientlyoptimize transition-state structures. The initial guess structures and eigenvectors areeither chosen from a linear interpolation between the reactant and product structures,from a nudged-elastic band search, from a constrained-optimization scan, or from theminimum-energy structures. Alternatively, initial guess vectors based on chemical in-tuition may be defined. We then iteratively refine the selected vectors by the Davidsonsubspace iteration technique. This procedure accelerates finding transition states forlarge molecules of a few hundred atoms. It is also beneficial in cases where the start-ing structure is very different from the transition-state structure or where the desiredvector to follow is not the one with lowest eigenvalue. Explorative studies of reactionpathways are feasible by following manually constructed molecular distortions. Keywords:
Transition-state search, Mode-Tracking, Davidson subspace iteration,Eigenvector following, Reaction kinetics. † ETH Z¨urich, Laboratorium f¨ur Physikalische Chemie, Vladimir-Prelog-Weg 2, 8093 Z¨urich, Switzerland ‡ University of Hamburg, Institute of Inorganic and Applied Chemistry, Martin-Luther-King-Platz 6,20146 Hamburg, Germany ∗ corresponding authors: [email protected], [email protected] a r X i v : . [ phy s i c s . c o m p - ph ] D ec NTRODUCTION
The optimization of transition-state structures (TSs) is key to the understanding of mech-anisms and kinetics of chemical reactions on a computational basis. Transition states aredefined as first-order saddle-point structures located on the minimum (reaction) energy pathbetween reactants and products. First-order saddle points are characterized by one negativeeigenvalue of the matrix of second partial derivatives of the electronic energy with respectto the Cartesian nuclear coordinates, i.e. of the Hessian. Reactants and products are localminima on the Born–Oppenheimer potential energy surface (PES). The energy differencesbetween a TS and two minima of an elementary reaction are the activation energy barriers.They should in principle be evaluated from the Gibbs free energy, but are approximatedhere, as in most quantum-chemical studies, by the electronic energy at zero Kelvin (neglect-ing temperature and entropy contributions).Numerous methods have been developed to efficiently find TSs. Examples are interpola-tion methods , eigenvector following (EVF) , string methods , and the scaled hyperspheresearch method . The existing TS search methods can be divided into those that start fromone structure (often called single-ended methods) or those that require at least two startingstructures, usually reactant and product structures (double-ended methods). Double-endedTS search algorithms are often based on interpolation methods such as linear (LST ) orquadratic synchronous transit (QST ), string methods or nudged elastic band (NEB )algorithms. Since the double-ended methods usually show slow convergence near a TS ,they are mainly employed to find a guess structure close to the TS, which is then refinedby a more efficient single-ended method, such as EVF. Hence, it is beneficial to combinesingle-ended and double-ended methods for TS searches.In most of the EVF-based methods, the full Hessian of the transition-state guess structureis calculated to obtain the exact vibrational mode to follow. For large molecules, the com-plete Hessian calculation is computationally demanding as the calculation of the elementsof the Hessian matrix is very time consuming within a first-principles electronic-structuredescription. Therefore, several algorithms have been developed to circumvent the calculationof the full Hessian in structure-optimization algorithms. A quasi-Newton–Raphson method2as been introduced by Broyden. In this method, an approximate Hessian is built fromgradients only and then updated (according to Bofill and Powell ) by the gradients ofintermediate points obtained during the optimization. These methods reduce the computa-tional effort significantly, but for large molecules a further reduction of the computationalcost is desirable.Recently, Sharada et al. introduced an approximate-Hessian approach based on thetangent of the transition-state guess structure determined by an interpolation between reac-tant and product structures and by local curvature information. This approximate Hessianapproach combined with the growing string method turned out to be computationally lessexpensive than previous Hessian approximations. Since the efficiency of a TS search depends strongly on the initial Hessian, a main goalis to set up an approximate Hessian matrix that resembles the exact one as closely as pos-sible. In 2002, we proposed an algorithm based on Davidson subspace diagonalization forthe selective calculation of eigenvectors of the mass-weighted Hessian based on predefinedmolecular distortions. This so-called Mode-Tracking scheme turned out to be very efficientin vibrational spectroscopy . Because of the straightforward and flexible implementa-tion, Mode-Tracking was implemented in a semi-numerical fashion. At the same time,Deglmann and Furche presented an implementation of a fully analytical Davidson sub-space diagonalization of the Hessian for the optimization of its lowest eigenvalue requiredfor the identification of stationary points.Very recently, Sharada et al. described a semi-numerical Davidson subspace iterationmethod to obtain selected information of the Hessian spectrum, which is identical to Mode-Tracking . For transition-state optimizations, Sharada et al. extract the guess modefrom the coordinates along the pathway obtained from the freezing-string method (FSM).In contrast to the Hessian approach presented by Sharada et al. , we here develop a Mode-Tracking-based TS and minimum localization algorithm that can iteratively refine a specificeigenvector of interest, which does not have to be the one with lowest eigenvalue. Ouralgorithm can be executed in an explorative fashion as we can circumvent the NEB orFSM calculation by starting from only one minimum-energy structure and by followingseveral eigenvectors in one optimization in parallel. We will demonstrate these capabilities3t the example of the isomerization reactions of formaldehyde, which has been studied asa benchmark system for automated transition-state search algorithms . Subsequently,we investigate an internal proton-transfer reaction in a hydrazine complex, which is anintermediate in the Schrock N -fixation catalytic cycle. At this example, we examinethe applicability of our algorithm for finding TSs for large molecules (the Schrock catalystcontains 284 atoms). Although smaller model complexes can be generated, the smallest ones,which resemble the structure of the original catalyst, still comprise 41 atoms.We choose these examples to highlight the capabilities of the Mode-Tracking-based ap-proach to TS searches, which improves on existing methods rather than proposing a newTS search algorithm. Hence, validating the performance of our Mode-Tracking version ofexisting TS search algorithms at standard TS test sets is neither needed nor necessary.This paper is organized as follows: In the next section, the Mode-Tracking algorithmand the theoretical background of transition-state optimizations are described. After thesubsequent Computational Methodology section, results are reported for our benchmarkreactions.
THEORY
The main idea of the algorithm to be described is to find transition-state structures by follow-ing only certain eigenvectors of the Hessian matrix selectively calculated by Mode-Tracking.For stationary structures, the harmonic vibrational frequencies and the corresponding eigen-vectors of a system can be obtained by solving the following eigensystem, HQ k = λ k Q k , (1)where H is the mass-weighted Cartesian Hessian, λ k are the eigenvalues, and the eigenvectors Q k are the mass-weighted vibrational normal modes. For non-stationary structures, forwhich the length of the geometry gradient is nonzero, Eq. (1) cannot be related to thevibrational properties of a molecule, but the eigenpairs (eigenvalues and eigenvectors) of H still characterize the PES.In Mode-Tracking, the eigenpairs of interest are obtained through a Davidson-type sub-4pace iteration method , in which the original Hessian matrix H is transformed to the(reduced-dimensional) Davidson matrix ˜H i , ˜H i = ( B i ) T HB i ≡ ( B i ) T Σ i , (2)where i denotes the i -th iteration step. B i is a matrix whose columns contain collectivedisplacement vectors b l ( l = 1,..., i ) along the 3 M (mass-weighted) nuclear Cartesian basisvectors ( M is the number of atoms). In our semi-numerical implementation, Σ i contains allvectors σ l , which collect the (numerical) derivatives of the (analytical) Cartesian gradient g of the total electronic energy with respect to the corresponding collective displacementvector b l , σ l = Hb l = (cid:88) n H ,n b ln (cid:88) n H ,n b ln ... (cid:88) n H M,n b ln = ∂∂ b l ∂E el ∂R ∂∂ b l ∂E el ∂R ...∂∂ b l ∂E el ∂R M = ∂∂ b l g . (3)By solving ˜H i c ik = λ ik c ik , (4)for the eigenvectors c ik and eigenvalues λ ik . In the i -th iteration step, Mode-Tracking calcu-lates the approximate k -th normal mode Q ik as Q ik = i (cid:88) l =1 c ik,l b l . (5)New basis vectors b i +1 are generated from the residuum vector, r ik = [ ˜H i − λ ik ] Q ik , (6)after applying a preconditioner X i to it , b i +1 = X i r ik . (7)The initial guess mode b can be obtained from the LST, which linearly interpolatesbetween the reactant and product structures, or from other path methods such as NEB.5et R (nmw) j be the non-mass-weighted ’(nmw)’ Cartesian coordinates of a structure j on thePES, then the initial normalized, non-mass-weighted mode is constructed from the coordinatedifferences between this point and each of its neighboring points, R (nmw) j +1 = { R (nmw) k,j +1 } and R (nmw) j − = { R (nmw) k,j − } , b (nmw) , j = 12 (cid:34) ( R (nmw) j +1 − R (nmw) j ) | R (nmw) j +1 − R (nmw) j | + ( R (nmw) j − R (nmw) j − ) | R (nmw) j − R (nmw) j − | (cid:35) = (cid:110) b (nmw) k,j (cid:111) , k = 1 , . . . , M. (8)This mode is then mass-weighted, b j = b (nmw) k,j √ m k (cid:80) Mk =1 (cid:16) b (nmw) k,j √ m k (cid:17) , k = 1 , . . . , M, (9)where m k is the mass of the k -th atomic nucleus.In general, Mode-Tracking can either optimize the mode with largest overlap with theinitial guess vector or the one with largest overlap with the approximate eigenvector chosenin the last iteration (root-homing). If the initial guess vector differs strongly from the normalmode of the transition-state structure, the second option might be more suited to find a TS.The eigenvector following algorithm is then employed to steer the optimization into thedirection of the TS and to finally locate it. Newton–Raphson steps along the converged Mode-Tracking eigenvector, which is referred to as transition vector, are carried out to maximizethe energy in this direction, while in all directions orthogonal to the transition direction thestructure is relaxed . For this, we project out the gradient along the transition vector, g TS ,from the total molecular gradient, g = { g k } = (cid:26) ∂E∂R k (cid:27) , k = 1 , ..., M. (10)To obtain the components of the molecular gradient that are orthogonal to the eigenvec-tor, g ort , we subtract the gradient part along the transition vector from the original moleculargradient and obtain g (nmw)ort = g (nmw) − Q (nmw)TS Q (nmw) ,T TS g (nmw) (cid:124) (cid:123)(cid:122) (cid:125) g (nmw)TS = (1 − Q (nmw)TS Q (nmw) ,T TS ) g (nmw) , (11)6here Q (nmw)TS is the selected eigenvector calculated with Mode-Tracking, which approxi-mately points into the direction of the TS. This is done in no-mass-weighted coordinates.The corresponding eigenvalue is λ TS .Let R be the coordinates of the targeted stationary point, for which g ≡ { ( ∂E/∂R k ) R k = R ,k } vanishes component-wise, and H ≡ { ( ∂ E/∂R k ∂R l ) R k = R ,k ,R l = R ,l } . From a truncated Tay-lor series expansion of the potential energy around E = E ( R ) on the PES, E ( R ) = E + g T ∆R + 12 ∆R T H ∆R + O ( ∆R ) . (12)the coordinate displacement ∆R ≡ R − R that leads to a stationary point ( dE ( R )/ d R = ), ∆R = − g H , (13)can be derived. R is the position of a stationary structure, g its gradient and H itsHessian. ∆R can be split into a direction parallel to the transition vector, ∆R TS , and intoall other directions. The step in the direction of the transition vector reads ∆R TS = − g TS λ TS . (14)The energy in direction of the selected mode is maximized if λ TS is negative. If we do notstart the EVF procedure from a structure close to the TS, but, for instance, from a minimumstructure, we must ensure that the transition vector is still followed uphill. This can eitherbe accomplished by employing the absolute value of λ TS ∆R TS = g TS | λ TS | , (15)or by employing Eq. (20) described below.To improve on the convergence of the EVF optimization, Wales defined a Lagrangianwith Lagrangian multipliers κ k for each degree of freedom l : L = − E − g T ∆R − ∆R T H ∆R + 12 M (cid:88) l =1 κ l (∆ R l − c l ) . (16)Wales’ method employs the rational function by Banerjee , in which the Lagrangianmultipliers are defined by the eigenvalues λ k and the gradient components g k along theeigenvectors, κ l = 12 (cid:18) λ l ± (cid:113) λ l + 4 g l (cid:19) . (17)7t appears to be more efficient to modify the equation to the following one: κ k = λ l ± | λ l | (cid:32) (cid:115) g l λ l (cid:33) (18)where ’+’ is for maximization and ’ − ’ for minimization.Wales arrived at the following equation that describes the steps to be made along alldegrees of freedom l , ∆R l = ± g l | λ l | (cid:18) (cid:114) g l λ l (cid:19) , (19)where ’+’ leads to an uphill and ’ − ’ to a downhill energy step. For a TS search, an uphillstep along the desired mode (i.e., the approximate transition vector) is required, ∆R TS = +2 g TS | λ TS | (cid:18) (cid:114) g λ (cid:19) . (20) Computational Methodology
The
MTsearch program
We implemented the theory presented in the previous section in a computer program called
MTsearch . The program is based on the original Mode-Tracking program , whichis currently available in its latest release as part of the
MoViPac package . MTsearch is a parallelized meta-program that accesses standard quantum-chemical programs for thecalculation of gradients and electronic energies. The computational methodology for thegeneration of these raw data is described in detail in the next subsection. The algorithmicstructure of
MTsearch is sketched in Figure 1.[Figure 1 about here.]A set-up tool, called tsdefine , creates the necessary input files for an
MTsearch calculation. With tsdefine we read in initial guess structures and, if available, initialmodes. The initial guess modes and structures can either be created within
MTsearch ,from a LST or an NEB path, or by an external program, which provides guess structures and8odes, e.g., based on a constrained optimization scan. The LST or NEB path consists of sixto twelve nodes (that is, molecular structures on an (approximate) reaction path, includingreactants and products), which we found to be a reasonable number. The spring forces inan NEB calculation are set to 0.02 a.u., and such a calculation is considered converged whenthe difference between the gradient norm of the actual iteration and the previous one dropsbelow 1 × − a.u.The first step of the TS optimization procedure is the Mode-Tracking optimization of theinitial guess mode to produce the corresponding minimal mode. Mode-Tracking is assumedto have converged when the maximum element of the residuum vector drops below 5 × − a.u. and the change in the length of the residuum vector drops below 5 × − a.u. One mayalso choose the convergence criteria corresponding to the last-added basis vector contributionto the selected eigenvector or to the change in the eigenvalue between the last iterations.After the calculation of a specific mode with Mode-Tracking, an EVF step is performedbased on this converged mode. For the Newton–Raphson step along the transition vector, ∆R TS , we define a maximum step size of 0.2 ˚A/ √ amu, which is decreased to 0.1 ˚A/ √ amuwhen the norm of g TS drops below 3 × − a.u., and to 0.05 ˚A/ √ amu when the norm of g TS drops below 1 × − a.u. For TS searches starting from minimum-energy structures, thefirst four Newton–Raphson steps are set to a maximum length of 1.0 ˚A/ √ amu, wheneverthe Hessian eigenvalue is positive or close to zero, i.e., no imaginary frequency with a largeabsolute magnitude is obtained.After each Newton–Raphson step a predefined number of optimization steps orthogonalto the eigenvector is performed. As default, a maximum of three iterations is chosen, if nototherwise mentioned.If the norm of the total gradient is still above the threshold (default is 1 × − a.u.) afterthe predefined number of orthogonal optimization steps, another Mode-Tracking calculationis launched, for which the last converged Mode-Tracking eigenvector is chosen as defaultguess vector. By default, a root-homing scheme selects the eigenvectors during the Mode-Tracking calculation with respect to the largest overlap with the initial one. For comparison,we also employed a root-homing scheme in which the eigenvector is always compared to theprevious one. 9t would also be possible to reuse the same eigenvector for a predefined number of EVFsteps, and/or to perform more than one EVF step between orthogonal optimizations. Thishas not been explored in this study. It is also possible to supply more information about thetransition path direction to MTsearch than only the first eigenvector (e.g., a sequence ofstructures which can for example easily be generated by a haptic device ). The guess vec-tors for the first few Mode-Tracking calculations are then chosen according to the predefinedsequence of structures. It has to be specified how many times the initial guess structure pathshall be taken as reference for creating a guess mode, which is then refined by Mode-Tracking.The structures and normal modes were visualized with
Pymol and Jmol , respec-tively. Raw data generation
All energies and gradients which are read as raw data by
MTsearch were calculated withdensity functional theory employing the
Turbomole program package (version 6.3.1) with Ahlrichs’ def2-SV(P), def2-SVP and def2-TZVP basis sets . MTsearch launchesthese calculations by system calls. Restricted and unrestricted BP86 all-electron Kohn–Sham calculations in combination with the resolution-of-the-identity technique were carriedout. Self-consistent-field single-point calculations are considered to be converged when thetotal electronic energy difference between two iteration steps is less then 10 − Hartree, if nototherwise indicated. Molecular structure minimizations are considered converged when thenorm of the geometry gradient is below 10 − a.u. For the optimization of transition-statestructures a geometry-gradient threshold of 10 − a.u. is chosen. Reference calculations
For comparison, we performed
Turbomole (version 6.3.1) EVF calculations for compar-ison with the
MTsearch results. Starting structures were chosen from the LST, NEB,or constrained optimization paths. We carry out a single-point calculation on the startingstructure and continue with a calculation of all vibrational modes with
Turbomole . Then,we employ the trust-radius imaging method (the maximum radius and the trust radius are10hosen between 0.1 ˚A/ √ amu and 0.2 ˚A/ √ amu) to follow the lowest eigenvalue. We refer tothis procedure in the following as “standard EVF method”. In these Turbomole calcula-tions, the BP86 density functional is chosen with Ahlrichs’ def2-SV(P), def2-SVP anddef2-TZVP basis sets .Furthermore, we performed constrained optimizations by employing the Gaussian program (version 09, Revision C.1) to obtain transition-state guess structures. Essentially,one internal coordinate was kept fixed at defined values and for all other degrees of freedoma structure optimization was carried out. Furthermore, intrinsic reaction coordinates werecalculated with Gaussian . In these calculations we employed BP86 with the def2-SVP basisset.
We have chosen the default convergence criteria (scfconv=tight, which means thatthe energy difference between two SCF iterations was less than 10 − Hartree, and that thestructure optimizations were considered converged when the root-mean-square force actingon all atoms was below 3 × − a.u.).We should note that we provide data for the eigenvalues of the Hessian as ’frequencies’(reported in units of wave numbers). I.e., we take the square root of the eigenvalues, whichcorresponds to a harmonic vibrational frequency for a stationary structure, even for non -stationary structures and denote it a ’frequency’ for the sake of convenience (eventually,these data become harmonic frequencies upon convergence of the stationary-structure op-timization). Moreover, to highlight imaginary frequencies, we add a minus sign in front ofthem (this is possible as the square of such a frequency still yields the correct eigenvalue ofthe Hessian matrix). Note also that we use the term ’mode’ to denote an eigenvector of theHessian matrix. RESULTS
To study the capabilities of
MTsearch , we have chosen four intramolecular reactions in-volving molecules of different sizes (4 atoms, 8 atoms, 41 atoms, and 284 atoms; shown inFigure 2). [Figure 2 about here.]11e start with the investigation of the rotational barrier in the C H molecule, becausethe transition-state structure is well defined and the system is small, which allows us to inves-tigate the suitable settings and thresholds for MTsearch . Next, we analyze the possibilityof
MTsearch to optimize several transition-state structures starting from one minimum-energy structure using hydroxymethylene as an example.The last two reactions considered are possible side reactions of the Chatt–Schrock cycleof N fixation at a molybdenum containing catalyst , in which N is reduced to ammoniaunder acidic and reductive conditions. Under these reaction conditions, it is possible thatseveral unwanted intermediates are formed. Exemplarily, we have chosen one possible side-reaction pathway, where one proton of N H coordinated to molybdenum shifts to one of theamido nitrogens.The Schrock catalyst is ligated by a tetradentate hiptN N ligand (hipt = hexa- iso -propylterphenyl). It has been intensively studied both experimentally and theoretically .Because of the relatively large system size of the hiptN N ligated hydrazine molybdenumcomplex (278 atoms), several smaller generic model system of the catalyst, in which thearyl substituents have been substituted, e.g., by H atoms or CH groups, have been stud-ied. In the following, the small and large Schrock catalyst refer to the MeNCH CH Nor hipt ligated molybdenum catalyst, respectively, with a hydrazine ligand as shown in Fig-ure 2. Our focus is the optimization of transition-state structures for the proton-transferreaction in the small and large Schrock catalysts.
Benchmark Example: C H rotation We first calculated the transition-state structure of ethane rotation from staggered confor-mation to eclipsed conformation and back to staggered conformation. To obtain startingstructures and initial guesses for the Mode-Tracking scheme, we performed a linear syn-chronous transit in internal coordinates with six nodes on the path including reactant andproduct structures (both staggered). Since only one dihedral angle is changed during thetransition from one minimum structure to the other, the choice of internal coordinates isvery useful in this example. Due to the symmetry in the LST path, only three of the sixstructures are different. The first (minimum), second, and third structures of the LST path-12ay were chosen as starting structures for EVF procedures performed with
MTsearch and,for comparison,
Turbomole .Although the second structure of the LST pathway does not feature negative eigenvaluesof the Hessian, the EVF algorithm optimization started from these structures convergedtowards the transition-state structure with both
Turbomole and
MTsearch (see Sup-porting Information for details).If we start from the energy minimum structure, EVF relying on one negative eigenvalueis not able to find the TS, because the structure is far away from the quadratic regionaround the TS. Therefore, one would usually start from a guess structure closer to the TS.By contrast,
MTsearch locates the TS starting from the energy-minimum structure witha mode specified by the LST and by a manually chosen mode corresponding to the rotationof one CH group around the C − C axis. The initial-guess structures and converged TSs canbe found in the Supporting Information.We analyzed the effect of various parameters on the convergence of
MTsearch . First, weinvestigated the optimal length of the first Newton–Raphson step. If the starting structureis close to the energy minimum structure, the optimizer has to accomplish a larger step outof the minimum. We observed that a first Newton–Raphson step size of 1.0–1.5 ˚A/ √ amu isappropriate (see Supporting Information, Section 2, for details).Next, we adapted the number of orthogonal optimization steps performed until the nextmode is optimized by Mode-Tracking to values between 2 and 10. To generalize the algo-rithm, we defined a protocol which stops the orthogonal optimization if the norm of thegradient for the optimization orthogonal to the transition path drops below 1 × − a.u.,which means that the maximum number of orthogonal-optimization steps needs not to bereached. Then, the next Mode-Tracking calculation and Newton–Raphson step in the direc-tion of the converged transition vector is performed. Explorative Example: Isomerization of H CO In this section, we study the possibility to find several TSs with
MTsearch starting from oneminimum-energy structure only. We have chosen the isomerization reaction of formaldehydeto hydroxymethylene and a subsequent trans-/cis-isomerization of hydroxymethylene as an13xample (see Figure 3). The transition-state structures are well known . Since tworeaction pathways are possible from trans-hydroxymethylene, a selective way of choosing theeigenvector of interest is important.[Figure 3 about here.]The trans-hydroxymethylene structure, from which we start the explorative TS search,can either undergo an internal hydrogen transfer from the oxygen atom to the carbon atom(over
TS-1 ) that leads to formaldehyde or a rotation around the C-O axis that leads to cis-hydroxymethylene (over
TS-2 ). MTsearch is able to locate both transition-state structuresby following the modes which are shown in Figure 3 next to the arrows indicating the reactiondirection. The three lowest modes of the starting structure (obtained by a full Hessiancalculation) have the following frequencies: 1100 cm − , 1188 cm − , and 1317 cm − . Thefirst one leads to TS-2 and the second and third ones to
TS-1 . MTsearch can find
TS-1 and
TS-2 also by starting from guess modes which are based on chemical intuition (seeSupporting Information, Section 1).The standard EVF optimization from hydroxymethylene following the lowest vibrationalfrequency mode does not converge to a TS, but falls back to the minimum-energy structure.Already a minor distortion of the minimum-energy structure towards the TS can alreadylead to a successful location of
TS-2 (see Supporting Information, Section 4.2, for details).Besides the TS optimization starting from trans-hydroxymethylene, we have also carriedout TS localizations from formaldehyde and cis-hydroxymethylene. For cis-hydroxymethylene,
TS-2 was found by following the LST guess mode. For formaldehyde, neither a LST guessmode nor a guess mode based on chemical intuition led to convergence to
TS-1 . Intramolecular proton-transfer reaction in a hydrazine Mo complex
In this section, we study the hydrazine intermediate of Schrock’s nitrogen-reducing cata-lyst and a generic model complex with aryl groups substituted by methyl groups. For theSchrock hydrazine complex, the lowest-energy spin state is a doublet. All spin states withhigher multiplicity are at least 60 kJ/mol higher in energy. We investigated the transition-14tate structure for a proton shift reaction from the nitrogen atom of N H that ligates tomolybdenum to one of the amido nitrogen atoms.For the generic model complex, we carried out a constrained optimization scan along theN amido -H distance that changes from reactant to product in 12 steps including reactants andproducts of 0.2 ˚A step size to obtain a guess structure close to the TS. From this constrainedoptimization scan, we selected the highest-energy structure ( Scan12 , see Figure 4) and aninitial guess mode based on the structures
Scan11 and the product. This mode has beensimplified by retaining only those entries that refer to the transferring hydrogen atom. Thisrestriction to the “moving” part in the system improves convergence as other motions ofparts of the system are discarded. Moreover, it produces a guess mode that is transferablebetween homologous species (see the large complex below).[Figure 4 about here.]The frequency analysis of structures
Scan11 and
Scan9 revealed that the lowest eigen-value modes do not correspond to the desired transition vector. We observed that thestandard EVF algorithm often fails to find the TS in this situation (see Supporting Infor-mation for details). Only for structure
Scan12 , which is already very close to the TS, thestandard EVF optimization converges towards the TS. By contrast,
MTsearch was ableto find the TS also from
Scan11 and
Scan9 (see Supporting Information, Section 4.3).The root-mean-square deviation between the
MTsearch -optimized transition-state struc-tures and the one calculated with
Turbomole ’s EVF is only 0.04 ˚A, which means that thetwo algorithms converged to the same structure. The vibrational analysis revealed exactlyone imaginary frequency of -i1244 cm − , and the intrinsic reaction coordinates (IRC) con-nect the reactant and product structures, which confirms that we found the desired TS. Thestationary points calculated by MTsearch are shown in Figure 5.It is noteworthy that the initial guess modes cannot only be obtained from a constrainedscan, but also from a LST or NEB pathway or based on chemical intuition. For isomerizationreactions, in which only one atom re-positions, as in our example, it is straightforward tomanually choose an approximate transition pathway (cf. Supporting Information). However,the manual set-up of a proper molecular distortion that is likely to resemble a reaction15athway is possible and potentially useful also for other types of reactions.In a Mode-Tracking-based TS search, one should confirm whether the very first modeconverged with Mode-Tracking corresponds to the desired reaction pathway, since all follow-ing optimization steps are based on the direction of this initial mode. However, this can bedone automatically by calculating the overlap of the initial guess mode and the convergedone. If the initial mode is not close to the transition vector, the optimization may leadto a different TS than the desired one (in our case, the TS for a rotation of the terminalNH moiety of the N H ligand was often found when the initial guess mode was not clearlydominated by the shifting proton).[Figure 5 about here.]For the optimization of the analogous transition-state structure in the significantly largerhiptN N-ligated Schrock Mo catalyst, we rotated the coordinating N H ligand in the minimum-energy structure such that an initial guess structure comparable to the TS of the genericmodel complex is obtained. Then, we performed a constrained optimization (fixed atomsare: molybdenum, the proton that moves and the two nitrogen atoms to which the protonbinds in the reactant structures). We choose as an initial guess mode the converged modeof the TS in the generic model complex (after alignment of the large and small homolo-gous complexes, and choosing only those entries for atoms that occur in both complexes;all other entries are set to zero). Due to the significantly larger system size, the maximumnumber of orthogonal optimization steps performed within MTsearch is increased to 10.The stationary points with the TSs calculated by
MTsearch are displayed in Figure 5.With Mode-Tracking we obtain one imaginary frequency of -i1323 cm − for the TS, whichis similar to the one of the TS in the small model system. The mode is located on the protonthat shifts. To study the performance of our algorithm, we also calculated the full Hessianand obtained one negative frequency mode of -i1338 cm − .Since the complete Hessian calculation and diagonalization for this molecule consisting of284 atoms takes significantly longer than the Mode-Tracking calculation (about a week vs.2 hours on 12 cores on a blade system featuring two six-core AMD Opteron 2435 processors(i.e., a total of 12 cores)), the computational time needed for the TS search is considerably16educed. CONCLUSIONS
The search for multiple reaction pathways starting from one minimum structure is still a mainobstacle of current transition-state optimization programs. In general, chemical intuition isneeded to choose a suitable starting mode, which connects the reactants with the products.In this paper, we presented an algorithm that efficiently combines the calculation ofselected normal modes by the Mode-Tracking scheme and the eigenvector-following pro-cedure to locate and optimize transition-state structures. Since Mode-Tracking avoids thetime-consuming calculation of the complete Hessian matrix and instead only optimizes themodes of interest, MTsearch is particularly suitable for optimizing transition-state struc-tures of large reactive molecular systems. The search for several transition-state structuresis feasible and the starting structures for a search may lie outside the quadratic region of atransition-state structure.We investigated our algorithm at the example of four intramolecular reaction pathways:the rotational barrier of C H , the isomerization of H CO, a proton shift in the hydrazine-bound intermediate of the [Mo(hiptN N)] catalyst by Schrock as well as in a model systemwith methyl substituents instead of the hipt substituents. Initial guess modes for the Mode-Tracking procedure can be extracted either from an LST or NEB pathway or based on chem-ical intuition. A TS optimization can be started either from two or from only one minimum-energy structure. The potential energy surface can be explored in a customized way alongthe desired directions. Even for a large molecule such as the hydrazine-coordinating Mocomplex of Schrock and coworkers with more than 200 atoms, we were able to efficientlylocate a transition-state structure.By choosing different initial guess modes and/or branching off at certain structures dur-ing the optimization, one may scan the potential energy surfaces along different directionssimultaneously. 17
CKNOWLEDGMENTS
This work has been financially supported by ETH Zurich and the Schweizer Nationalfonds(SNF project 200021L 156598). MB gratefully acknowledges support by a fellowship of theFonds der Chemischen Industrie (FCI). CH acknowledges funding by the FCI.18 eferences
1. T. A. Halgren and W. N. Lipscomb, Chem. Phys. Lett. , 225 (1977).2. C. J. Cerjan and W. H. Miller, J. Chem. Phys. , 2800 (1981).3. J. Simons, P. Jorgensen, H. Taylor, and J. Ozment, J. Phys. Chem. , 2745 (1983).4. D. J. Wales, J. Chem. Soc., Faraday Trans. , 653 (1992).5. D. J. Wales, J. Chem. Soc., Faraday Trans. , 1305 (1993).6. F. Jensen, J. Chem. Phys. , 6706 (1995).7. G. Henkelman, B. P. Uberuaga, and H. J´onsson, J. Chem. Phys. , 9901 (2000).8. K. Ohno and S. Maeda, Chem. Phys. Lett. , 277 (2004).9. C. Peng and H. P. Schlegel, Isr. J. Chem. , 449 (1993).10. G. Henkelman and H. J´onsson, J. Chem. Phys. , 9978 (2000).11. F. Jensen, Introduction to Computational Chemistry (2 nd ed. Wiley & Sons, Chichester,1989).12. C. G. Broyden, Math. Comput. , 368 (1967).13. R. Fletcher, Practical Methods of Optimization. Unconstrained Optimization (Band 1.Wiley, New York, 1980).14. J. M. Bofill, J. Comput. Chem. , 1 (1994).15. M. J. D. Powell, Math. Prog. pp. 26–57 (1971).16. S. M. Sharada, P. M. Zimmermann, A. T. Bell, and M. Head-Gordon, J. Chem. TheoryComput. , 5166 (2012).17. P. M. Zimmerman, J. Chem. Phys. , 184102 (2013).18. P. M. Zimmerman, J. Comput. Chem. , 1385 (2013).199. M. Reiher and J. Neugebauer, J. Chem. Phys. , 1634 (2003).20. J. Neugebauer and M. Reiher, J. Phys. Chem. A , 2053 (2004).21. J. Neugebauer and M. Reiher, J. Comput. Chem. , 587 (2004).22. M. Reiher and J. Neugebauer, Phys. Chem. Chem. Phys. , 4621 (2004).23. M. Reiher, B. Le Guennic, and B. Kirchner, Inorg. Chem. , 9640 (2005).24. T. B. Adler, N. Bohro, M. Reiher, and M. A. Suhm, Angew. Chem. , 3518 (2006).25. C. Herrmann and M. Reiher, Surf. Science , 1891 (2006).26. C. Herrmann, J. Neugebauer, and M. Reiher, New. J. Chem. , 818 (2007).27. C. Herrmann, J. Neugebauer, and M. Reiher, J. Comput. Chem. , 2460 (2008).28. P. Deglmann and F. Furche, J. Chem. Phys. , 9535 (2002).29. S. M. Sharada, A. T. Bell, and M. Head-Gordon, J. Chem. Phys. , 164115 (2014).30. S. M. Sharada, A. T. Bell, and M. Head-Gordon, J. Chem. Phys. , 229902 (2014).31. H. L. Davis, D. J. Wales, and R. S. Berry, J. Chem. Phys. , 4308 (1990).32. K. Bondensg˚ard and F. Jensen, J. Chem. Phys. , 8025 (1996).33. W. Quapp, M. Hirsch, and D. Heidrich, Theor. Chem. Acc. , 285 (1998).34. S. Maeda, K. Ohno, and K. Morokuma, Phys. Chem. Chem. Phys. , 3683 (2013).35. D. V. Yandulov and R. R. Schrock, J. Am. Chem. Soc. , 6252 (2002).36. D. V. Yandulov and R. R. Schrock, Science , 76 (2003).37. R. R. Schrock, Acc. Chem. Res. , 955 (2005).38. J. Baker and F. Chan, J. Comput. Chem. , 888 (1996).39. Y. Zhao, N. Gonz´alez-Garc´ıa, and D. G. Truhlar, J. Phys. Chem. A , 2012 (2005).200. V. Guner, K. S. Khuong, A. G. Leach, P. S. Lee, M. D. Bartberger, and K. N. Houk, J.Phys. Chem. A , 11445 (2003).41. P. M. Zimmerman, J. Chem. Theory Comput. , 3043 (2013).42. D. J. Wales, J. Chem. Phys. , 3750 (1994).43. D. O. Neal, H. Taylor, and J. Simons, J. Phys. Chem. , 1510 (1984).44. A. Banerjee, N. Adams, J. Simons, and R. Shepard, J. Phys. Chem. , 52 (1985).45. D. J. Wales, Mol. Phys. , 1 (1991).46. J. Neugebauer, C. Herrmann, and M. Reiher, (2009).47. T. Weymuth, M. P. Haag, K. Kiewisch, S. Luber, S. Schenk, C. R. Jacob, C. Herrmann,J. Neugebauer, and M. Reiher, J. Comput. Chem. , 2186 (2012).48. K. H. Marti and M. Reiher, J. Comput. Chem. , 2010 (2009).49. M. P. Haag and M. Reiher, Faraday Discuss. , 89 (2014).50. M. P. Haag, A. C. Vaucher, M. Bosson, S. Redon, and M. Reiher, ChemPhysChem , 165(1989).54. F. Weigend and R. Ahlrichs, Phys. Chem. Chem. Phys. , 3297 (2005).55. A. D. Becke, Phys. Rev. A , 3098 (1988).56. J. P. Perdew, Phys. Rev. B , 8822 (1986).217. M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman,G. Scalmani, V. Barone, B. Mennucci, G. A. Petersson, et al., Gaussian 09 (Gaussian,Inc., Wallingford CT, 2009).58. A. Sch¨afer, H. Horn, and R. Ahlrichs, J. Chem. Phys. , 2571 (1992).59. A. Sch¨afer, C. Huber, and R. Ahlrichs, J. Chem. Phys. , 5829 (1994).60. B. LeGuennic, B. Kirchner, and M. Reiher, Chem. Eur. J. , 7448 (2005).61. F. Studt and F. Tuczek, Angew. Chem. Int. Ed. , 5639 (2005).62. S. Schenk, B. L. Guennic, B. Kirchner, and M. Reiher, Inorg. Chem. , 3634 (2008).63. S. Schenk and M. Reiher, Inorg. Chem. , 1638 (2009).64. S. Schenk, B. Kirchner, and M. Reiher, Chem. Eur. J. , 5073 (2009).65. D. V. Khoroshun, D. G. Musaev, and K. Morokuma, Mol. Phys. , 523 (2002).66. Z. Cao, Z. Zhou, H. Wan, and Q. Zhang, Int. J. Quantum Chem. , 344 (2005).67. F. Studt and F. Tuczek, Angew. Chem. , 5783 (2005).68. A. Magistrato, A. Robertazzi, and P. Carloni, J. Chem. Theory Comput. , 1708 (2007).69. L. Deng, T. Ziegler, and L. Fan, J. Chem. Phys. , 3823 (1993).70. J. Neugebauer, M. Reiher, C. Kind, and B. A. Hess, J. Comput. Chem. , 895 (2002).22 ) efficient calculation of transition states of large molecules B) explorative studie s MTsearch
The access of specific normal modes out of all vibrational modes opens up new possibilitiesfor the search of stationary points. Since the calculation of the complete Hessian matrixis avoided, transition-state structures of large molecules can be obtained. Because of thepossibility to start the transition-state search from one initial guess structure, an explorativeinvestigation of the potential energy surface is feasible.23 aw data generation (geometry gradient, electronic energy):standard quan-tum chemistry programs
MTsearch
Mode-TrackingEVF:
1) Newton-Raphson along TS mode 2) orthogonal opt. until max. number of ort. opt. steps
TS confirmation:
1) Mode-Tracking calc. of lowest frequency2) Intrinsic Reaction Coordinates tsdefine: set-up tool I nput preparation: TS-guess mode b and TS-guess structure R a) generation of b and R by LST or NEB from reactants b) read in external R and b not conv.? conv.? TS optimization: for i < max: TSTSTS
Figure 1: Overview of the
MTsearch meta-program structure.24 atoms 41 atoms 284 atoms
Figure 2: Molecular models considered in this work: H CO (left), ethane (second from left), asmall (third from left) and the full Schrock hydrazine tris(amido)amine Mo complex (right).Element color code: green, C; red, O; blue, N; cyan, Mo; white, H.25 nergy inkcal/mol
Reaction coordinate + TS-1 TS-2
Figure 3: BP86/RI/def2-SV(P) reaction path from formaldehyde to trans-hydroxymethyleneand cis-hydroxymethylene with transition-state structures optimized with
MTsearch . Themodes taken from a full vibrational analysis that lead to the TSs,
TS-1 and
TS-2 , are alsodepicted. A maximum number of three orthogonal optimization steps has been chosen in
MTsearch . Element color code: gray, C; red, O; white, H.26 can9 Scan11 Scan1268.9°68.6°76.2°1.9 Å 1.5 Å 1.3 Å TS69.5°1.3 Å
Figure 4: Initial guess structures chosen from a constrained optimization scan along oneH(N H )-N amido distance. Scan9 , Scan11 and
Scan12 are the 9th, 11th and 12th structurefrom a 12-step constrained optimization scan with the program
Gaussian (0.2 ˚A increase ofthe N amido -H distance in each step) along the N amido -H distance starting from the hydrazinebound complex. 27 .0 21.139.3Energy inkcal/mol0102030405060 Reaction coordinate0.0 13.127.9Energy inkcal/mol051015202530 Reaction coordinate
Figure 5: Internal proton shift reaction pathway of the generic Schrock model system (top)and the full Schrock hydrazine Mo complex (bottom) with BP86/RI/def2-SV(P) transition-state structures optimized with