An Algorithm to Locate Optimal Bond Breaking Points on a Potential Energy Surface for Applications in Mechanochemistry and Catalysis

Josep Maria Bofill, Jordi Ribas-Ariño, Sergio Pablo García, Wolfgang Quapp Departament de Química Inorgànica i Orgànica and IQTCUB, Universitat de Barcelona, Martí i Franquès 1, 08028 Barcelona, Spain. E-mail: jmbofill@ub.edu Departament de Ciència de Materials i Química Física and IQTCUB, Universitat de Barcelona, Martí i Franquès 1, 08028 Barcelona, Spain. E-mail: j.ribas@ub.edu Mathematisches Institut, Universität Leizpig, Augustus-Platz PF 100920, D-04009 Leipzig, Germany. E-mail: quapp@uni-leipzig.de


I. INTRODUCTION
][3][4][5][6][7][8][9][10][11][12][13][14][15] Several experimental methods allow us to apply mechanical forces to molecules.On the one hand, it is now well-established that ball milling and grinding techniques can be used not only for crushing solids but also for the initiation of chemical reactions. 16,17On the other hand, recent developments in single molecule force spectroscopy (especially in atomic force microscopy [18][19][20][21] ), sonochemical techniques, [22][23][24] and molecular force probes 25 have enabled the application of tensile forces to single molecules.7][28][29][30][31][32][33][34][35][36][37][38][39] From a conceptual point of view, the phenomenon of mechanical activation can be understood on the basis of the fact that the potential energy surface (PES) of a given reactive system changes when this system is subjected to tensile stress.As a result of these force-induced changes of the PES, the barriers between the minimums and the saddle points change. 38For this reason, when a force is applied to a molecular system, the reactivity is either enhanced or suppressed. 39The extent to which a barrier can be modified depends on the direction and on the magnitude of the external a) Electronic mail: jmbofill@ub.edub) Electronic mail: j.ribas@ub.educ) Electronic mail: quapp@uni-leipzig.deforce.If we consider the case in which a constant stretching force is applied to a molecular system, the resulting modified PES (which can be called force-transformed PES 31 or force-modified PES 33 ), V f (x), is obtained via 1,7,13,31-33,40-42 where the superscript T means transposition of a vector or a matrix.According to Eq. (1), the effective potential results from a subtraction of a linear potential with the force vector, f, and the displacement, (x x 0 ), from the original potential energy, V (x).x 0 is any anchor point.Equation (1) is the basis of some mechanochemical methods for the calculation of deformed molecular structures due to an external mechanical force.The dimension of the vectors, f, x, and x 0 , is N = 3n 6 internal coordinates.A proof of the invariance of Eq. (1) under a coordinate transformation is given in Appendix A. An alternative proof is given in Ref. 43.
On the effective potential, V f (x), the stationary points are located at different points with respect to the unperturbed potential, V(x), where it holds ∇ x V (x) = g(x) = 0.The stationary points on the effective potential have to satisfy the analogous condition, ∇ x V f (x) = 0. Since V f (x) is that given in Eq. (1), it follows that its stationary points (minima and SPs) should satisfy One searches a point where the gradient of the original PES, g(x), has to be equal to the constant mechanochemical force, f, being the force that induces the chemical process.If the mechanical stress in a defined direction is f = F l with a fixed unit vector, l, and with magnitude, F, then we get from Eq. (1) at the stationary points: l (||g(x)|| − F) = 0 since g(x) = l||g(x)|| and F = ||g(x)||.The magnitude of F coincides with the gradient norm, F = ||g(x)||, at the stationary points of the effective PES.From a geometrical point of view, Eq. ( 2) means that the tangential hyperplane of the original PES, characterized by the gradient, g(x), is equal to the hyperplane of the pulling force, f T (x − x 0 ), of Eq. (1).In the absence of any tensile force (i.e., in the situation that F = 0), Eq. (1) reduces to the pure thermal limit.
Let us consider the case that the force f changes its magnitude F in a continuous way but its direction l is constant.This would give rise to a set of effective PESs, V f (x), each of them with their corresponding stationary points.This set of stationary points defines a continuous curve on the original PES, V (x).Under our construction, at each point of the curve, the equation Fl = l ||g(x)|| = g(x) holds, see Eq. (1).Thus, in each point of this path, the gradient has the same direction.A curve with this feature is known as Newton trajectory (NT) [44][45][46][47][48][49][50][51] or more recently as the curve of the force displaced stationary points (FDSPs). 38,39It is known that every NT satisfies the Branin equation, 38,52,53 where t is the parameter that characterizes the curve and the adjoint matrix A(x) is computed as a product of the determinant of the Hessian matrix, H(x), of the original PES, V (x), with its inverse, A(x) = Det(H(x))(H(x)) 1 .Equation ( 3) is also a way to generate the FDSP curve.The slope of the effective PES, V f (x), along the FDSP curve is given by In the derivation of Eq. ( 4), the directional derivative concept and Eqs.(1) and ( 3) have been used.Equation ( 4) tells us that one has to increase the magnitude of the force F in order to move along the FDSP path.If the force increases, a sub-arc of this path leads uphill from the minimum and another sub-arc leads downhill from the saddle point.Eventually, the two arcs will meet at the point where the gradient norm achieves its maximal value.
If we go along an NT and if we differentiate the square of the gradient norm of the PES, we get It is zero, either at the stationary points or at the points of the manifold Det(H(x)) = 0 of the original PES, V (x).At the points with Det(H(x)) = 0 emerges the bond breaking point (BBP). 38,54At the latter type of points, the gradient norm reaches a turning point, and the effective V f (x) PES presents a shoulder on the FDSP path, 38,39 see Eq. ( 4).As a result, the difference in energy between the minimum and saddle point configurations (i.e., the energy barrier that must be overcome for the reaction to occur) vanishes on the effective PES, V f (x).This is so because, in the evolution of the FDSP path from the reactant minimum, the point where Det(H(x)) = 0 is located before the saddle point of the original PES.The gradient of the original PES at a given BBP gives the magnitude of the force that should be applied to a molecule in order to mechanically induce a chemical transformation along the pulling direction associated with this gradient.5][56][57] For this theory, a BBP represents a catastrophe of the changing PES function, V (x), being unfolded by a force through the additional perturbation term, f T • (x − x 0 ), of Eq. (1).This gives us the general structure of all possible effective PESs, V f (x).They have a Hessian matrix at the BBP with zero eigenvalue.This simple mathematical model can be used to predict the effect of a pulling force acting over a molecular system, at least as a first approximation to describe the phenomenon.In fact, the catastrophe theory has already been used to compute BBPs and the corresponding rupture forces for benzocyclobutene derivatives in a simple pulling scenario in which a pair of opposing forces act on two terminal atoms. 58he BBPs define a manifold or a curve in the twodimensional case.A way to generate a curve with condition Det(H(x)) = 0 is a predictor-corrector method with predictor steps along the tangent, t, by the solution of It is a homogeneous system of N linear equations.A more elaborated treatment of the derivation in Eq. ( 6) is given in Appendix B. Within the manifold of BBPs of a given PES, there is an optimal BBP, 38,39 see Figs. 2-4 and 6 of Ref. 38.This BBP defines the lowest force in magnitude and the corresponding pulling direction that should be applied to a molecular system in order to mechanically promote a given chemical transformation.The optimal BBP satisfies the equation, 38 where g(x) 0. Thus at the optimal BBP, the gradient is an eigenvector of the Hessian matrix with null eigenvalue.The optimal BBP point coincides with a point of the gradient extremal (GE) [59][60][61] exactly at the intersection point with the Det(H(x)) = 0-line. 38,39Equation ( 7) is an eigenvalue equation where the eigenvalue of the eigenvector g(x) is zero.The location of optimal BBPs is extremely important in the context of mechanochemistry because these points reveal which is the most efficient way to trigger a reaction by means of a mechanical force.It is thus of paramount importance to devise algorithms for the location of these special critical points.In this article, we present an algorithm to locate BBPs on the PES of a molecular system, which is based on the Gauss-Newton method.
This paper is organized as follows.First, we demonstrate that locating an optimal BBP is equivalent to minimizing or finding the zeros of a positively defined function and that such a minimization can be performed using a Gauss-Newton method (Sec.II).Then, the algorithm based on the Gauss-Newton method and the rational function optimization technique will be worked out in detail (Sec.III).After this, we show that the devised algorithm works efficiently in well-known test 2D functions (Sec.IV).Subsequently, we show that the algorithm is also able to efficiently locate the optimal BBP on the multidimensional PES of a real chemical example (Sec.V).Finally, Sec.VI of the article is the conclusions.

II. THE σ-FUNCTION OF THE OPTIMAL BARRIER BREAKDOWN POINTS
In a theory of mechanochemistry, the importance of the optimal BBP resides in the fact that it gives the optimal direction and magnitude of the pulling force.This necessitates an algorithm to find this type of points.0][61] However, this is very expensive, and the algorithms to locate this type of paths are often not numerically stable, except for low dimensional systems.Sometimes the GE has an avoided crossing. 60At this point, we can take advantage of the fact that we are not interested in the whole GE path but only in a specific point of such a path, namely, the point where the Hessian has a zero eigenvalue.Locating this point is equivalent to finding the zero of the σ-function, where the vector s(x) is The function σ(x) is a sum of squares of nonlinear functions, s(x); thus, it is a non-negative function.It is clear that a minimum x min of the function σ(x) for which σ(x min ) = 0 is a desired solution since this can only arise if x min satisfies s(x min ) = 0, which is the optimal BBP condition given in Eq. (7).To find σ(x min ) = 0 of Eq. ( 8) falls into the class of the so-called nonlinear least squares problems. 62For this type of problems, the derivatives of σ(x) with respect to x are given, and where s i (x) is the ith-element of the s(x) vector and is the N × N Jacobian matrix, an element of this matrix has the form, J ij (x) = ∂s j (x)/∂x i , for i, j = 1, . . ., N. Thus, by differentiation of s(x) given in Eq. ( 9), with respect to x, the Jacobi matrix, J(x), in the present problem has the form, Here I is the identity matrix of dimension N × N, g(x) is the normalized gradient vector, and T(x)g(x) is the matrix of the contraction of the gradient vector and the tensor of the third energy derivatives with respect to the coordinates, x.Thus it is again an N × N matrix.One can find the minimum using the formulas of Eqs. ( 10) and ( 11) in conjunction with the Newton method.Nevertheless, it is necessary to evaluate the set of Hessian matrices, , one for each element of the s(x)-vector.In addition, these matrices involve fourth-order energy derivatives with respect to x.However, because we are interested in the zero of the s(x) vector in the least square sense, the set of elements s i (x) N i=1 are small.This suggests that the second term of Eq. ( 11) can be neglected.We get the approximate expression for the second derivatives of the σ-function, which is equivalent to a linear approximation with respect to x to each element of the s(x) vector.In this way, the first and approximate second derivatives of the σ-function can be determined by using only s(x) and J(x), Eqs. ( 10) and ( 14), respectively.When the second derivative is approximated as given by Eq. ( 14), the basic Newton method becomes the Gauss-Newton method or the generalized least squares method. 62Note that the Gauss-Newton method can fail or can converge very slowly.To improve this, we use the restricted step algorithm.In the ith-iteration, the modified Gauss-Newton is where , and α (i) is a parameter.

III. AN ALGORITHM TO LOCATE OPTIMAL BBPs
Many different algorithms based on Eq. ( 15) have been suggested.The algorithm used in the present context is based on the rational function optimization technique, [63][64][65][66] where ν (i) and ∆ x (i) are obtained by the solution of the eigenvalue equation, 0 taking the eigenvector of the lowest eigenvalue.The last N 1 rows of Eq. ( 16) correspond to Eq. (15a) ordered in a different way, and they permit to compute ∆ x (i) , whereas the first row permits to obtain ν (i) .In this way, the decrease goes with the second order of ∆ σ(x To obtain a better efficiency of the Gauss-Newton method, a line search is performed once ∆ x (i) has been obtained.The algorithm determines x (i+1) = x (i) + α (i) ∆ x (i) , by using a line search strategy to find the optimal α (i) .The matrix J (i) is updated following the Broyden formula. 67The line search determines α (i) opt to minimize σ(x (i+1) ) = σ(x (i) + α (i) opt ∆ x (i) ) = s (i+1)T s (i+1) and σ(x (i) ) = s (i)T s (i) for α (i) = 0.
Defining the difference vector y (i) = s (i+1) s (i) , the matrix J (i+1) can be calculated by the update formula, This update formula satisfies the Newton condition only for the current i-iteration, J (i+1)T ∆ x (i) α (i) opt = y (i) . 62The Broyden update formula is a special case of the Barnes update formula. 68The algorithm flow is illustrated in the following steps: Step 1. Set i = 0 and set x (1) .Calculate the J (1) matrix with where the e i vector is the i-th column of the unit matrix.Calculate σ(x (1) ) = σ (1) .
Step 2. Set i = i + 1 and form the vector q (i) = J (i) s (i) and the matrix J (i) J (i)T .
Step 3. Form the matrix 0 q (i)T q (i) J (i) J (i)T and diagonalize it.
Step 4. Let t (i) be the eigenvector of the smallest eigenvalue of the matrix of Step 3. Set t (i) = t (i) /t (i) 1 , where t (i) 1 is the first component of the t (i) vector.Set ∆ x (i) = t (i) N , where t (i) N is the vector formed by the last N components of the t (i) vector.
Step 5. Perform a line search to determine α (i) opt to minimize σ(x Step 6.If max 1≤j ≤N |s (i+1) j | ≤ , exit or otherwise go to Step 7.
Step 7. Compute the vectors y (i) and v (i) and the matrix J (i+1) by Step 8.If i + 1 >, maximal number of iterations then exit, otherwise go to Step 2. The algorithm only needs the parameters, h, , l , and the maximal number of iterations.
As the start point for the above algorithm, we take the point with the highest ||g(x)|| value of the IRC, the steepest-descent path from the transition state to the reactant minimum.In this way we ensure that the optimal BBP found by the algorithm is located in the desired reaction valley or near the reaction path.
Finally the located point is characterized by computing the Hessian and the gradient in this point.The Hessian matrix is diagonalized and the null eigenvector should overlap with the normalized gradient vector within a given criterium.
A main problem of the present algorithm is the computational cost by the computation of the Hessian matrix at each iteration.We did some work to reduce this computational effort by introducing some type of update such that the stability of the algorithm remains using this formula.Results should be reported in a subsequent paper.

A. The Rosenbrock surface
A well-known test function for optimization methods is the Rosenbrock function. 69,70Over a two-dimension plane, it is FIG. 1. Contours of Rosenbrock's function 69 given in Eq. ( 18).Inner contours are in steps of 1 energy unit, but the outer contours are in steps of 5 energy units.The green curve is the Det(H(x, y)) = 0 line and the gray curve is the GE. 59The black point is the minimum at (1,1).
The contours of the function are shown in Fig. 1.The global minimum at point (1, 1) is inside a long, narrow, parabolically shaped flat valley.Normally it is a difficult task to converge to this global minimum.The line of points where Det(H(x)) = Det(H(x, y)) = 0 is characterized by the line The line is described by the green parabola in Figs.1-3.At each point of this line, the gradient is parallel to the unit vector, Thus the gradients evaluated at the points belonging to the line Det(H(x, y)) = 0 point to the same direction, which is independent of the position.The optimal BBP, the point that belongs to this line and also satisfies Eq. ( 7), is x T = (x, y) FIG. 2. Contours of the σ(x, y) function, Eq. ( 8), of the Rosenbrock function, Eq. ( 18).The green curve is the Det(H(x, y)) = 0 line.The optimal BBP is located at the point (0.250, 0.0675) (orange).Here is the zero minimum of the σ function.Thus the symmetric shape of the level lines is a delusion.FIG. 3. Behavior of the algorithm to locate the optimal BBP on the Rosenbrock function, Eq. ( 18).The set of points is a selected subset of all points of the optimization process, focused on the domain −0.26 ≤ x ≤ −0.2 and 0.04 ≤ y ≤ 0.1.The initial point, (0.250, 0.0675), is out of this region.The selected points of the optimization process are located in the deep valley of the σ-function, Eq. ( 8).The green curve is the Det(H(x, y)) = 0 line.The optimal BBP is at the point (0.250, 0.0675) corresponding to the label x 11 .
= (0.250, 0.0675).Note that the green line and the valley branch of the GE are not to tell apart in Fig. 1.Equidistant contours of the function σ(x, y), Eq. ( 8), for the Rosenbrock function are shown in Fig. 2. The function seems to be planar in the region −0.3 ≤ x ≤ 0.3; however, it is caused by very steep "mountains" at the both sides.The line Det(H(x, y)) = 0 is characterized by the parabola of Eq. ( 19).It is located in a deep valley of the function σ(x, y).Starting at the point (0.250, 0.0675), the above algorithm reaches the optimal BBP with a behavior that is shown in Fig. 3.The set of points generated by the algorithm is also located in the deep valley of the σ-function, Eq. ( 8).The points are near or on the Det(H(x, y)) = 0-line.
FIG. 4. The Müller-Brown surface. 71The IRC paths are the red curves.The green curves are the set of points Det(H(x, y)) = 0.The pink points are the optimal BBPs located in reaction path valleys.The black points are the minima and transition states.

B. The M üller-Brown surface
The Müller-Brown surface 71 has been used several times to test new algorithms to find both stationary points and reaction paths.This surface can be used to describe a reaction mechanism that involves a rate-determining step followed by an equilibrium step or vice-versa, that is to say, a pre-equilibrium step followed by the rate-determining step.We use this 2-dimensional surface to test the algorithm, as well as to give a procedure for the selection of a start point.The Müller-Brown PES is displayed in Fig. 4, together with the intrinsic reaction coordinate (IRC) 72  minima through the corresponding saddle points of index one.The optimal BBPs located in the reaction valley are also shown.
The valley-ridge inflection (VRI) points of the surface define a special sort of NTs, the singular NTs (Fig. 5).A regular NT directly leads from a minimum to an SP of index one (index theorem 73 ).A singular NT bifurcates at the respective VRI point.Two bifurcating branches can connect two neighboring SPs of index one, or they can connect two minimums, or they can connect a minimum with an SP of index two (which here does not exist).The singular NTs form the borders of families FIG. 7. The σ(x, y)-function of the Müller-Brown surface. 71The orange points are the optimal BBPs located in the reaction valleys.The black points are the stationary points associated to the minima and saddle points of index one on the MB surface.The green curves are the set of points with Det(H(x, y)) = 0. Apart from the marked optimal BBPs, the surface has more points of this type that are not located in the reaction valleys. of regular NTs in between.Such families are named chemical corridors for a pulling. 38,39,74he main direct chemical corridor for an enhanced pulling from minimum R to SP 1 is shown in Fig. 6, compare Ref. 73.Later the corridor follows the full reaction valley from reactant, R, to intermediate, I, and product, P. The borders of the corridor are the singular NTs through the VRI points 1 and 4. The borders are drawn in blue and orange color.An optimal BBP is a red point.VRIs are black points.
Particularly important for the present discussion is the form of the σ(x, y)-function of the Müller-Brown surface.
FIG. 8. Evolution of the optimization process for the optimal BBPs of the Müller-Brown surface. 71The depicted level lines are the σ (x, y)-function of the Müller-Brown potential energy near the optimal BBPs distributed along the reaction valley.The blue points are the iterated points with the corresponding iteration number, x i .The green line is the set of the points where the equation, Det (H (x, y)) = 0, is satisfied.The thin gray curves are some iso-contours of the σ-function.The evolution of the optimization process to locate BBP 1 , BBP 2 , BBP 3 , and BBP 4 is displayed in the panels (a)-(d), respectively.This function is reported in Fig. 7 where the position of the set of optimal BBPs located in valley regions is marked in orange color.Apart from these optimal BBPs, there exist other extremal BBPs not located in these regions.These points are not shown in Fig. 7.They have to be located at the intersection of the green line, the line of points with Det (H (x, y)) = 0, with the deeper regions of the σ(x, y)-function since in these points this non-negative function takes the value zero.The optimal BBPs on the reaction valley are the points (−0.946, 1.040), (−0.575, 0.507), (0.163, 0.393), and (0.269, 0.094), labeled as BBP 1 , BBP 2 , BBP 3 , and BBP 4 , respectively, in Fig. 4.These optimal BBPs were located using the above algorithm starting at the BBP of the corresponding IRC path.
The behavior of the algorithm during the location process is depicted in Fig. 8.In the location process of the optimal BBP 1 in Fig. 8(a), all iterated points are in the catchment of a deep valley of the σ-function where this point is found.
Starting the algorithm at a point of the IRC curve near the SP 1 , the algorithm can converge to other optimal BBPs that are not in the reaction valley or far from the IRC, namely, (−1.170, 0.742), (−1.046, 0.339), and (−1.098, 0.648).This is because the valley of the σ-function where the optimal BBP 1 is located is a very narrow valley with some further zero minima.
The location process of the optimal BBP 2 and BBP 3 does not strongly depend on the start point of the IRC curve.The elliptic form of the iso-contour curves of the σ-function around these points is the reason 62 why the algorithm converges very fast to the desired optimal BBP.
The last optimal BBP, namely, BBP 4 , is the most difficult to locate due to the nonlinear form of the σ-function around this point, see Figs. 7 and 8(d).Starting the localization of the BBP 4 at any point near the IRC arc in between SP 2 and P can converge to the point (−0.258,0.081) being an extremal BBP far away from the IRC curve.But usually for the initial guess of a point of the IRC, the algorithm converges to the BBP 4 .In this case the convergence goes faster just when the location process enters the region where the set of iso-contours curves has an elliptic structure, see Fig. 8(d).

V. A CHEMICAL EXAMPLE: THE 1,2-SIGMATROPIC H-SHIFT REARRANGEMENT OF CYCLOPENTADIENE
After showing that the algorithm of the minimization of the σ-function is able to efficiently locate optimal BBPs in the case of well-known 2D test functions, we now turn our attention to the performance of the algorithm for a location of the optimal BBPs of a multidimensional PES associated with a real chemical transformation.We have chosen the 1,2sigmatropic H-shift rearrangement of cyclopentadiene as a model system to gauge the performance of our algorithm.This reaction was computationally studied by Jiao and Schleyer. 75ecently the same reaction has been used as a model for the simplest corridor of type one in the classification of forces acting in mechanochemistry and catalysis. 74However, in this former analysis of acting forces, the quasi optimal BBP was taken as the point of the highest value of the gradient norm along the IRC curve.Note that this is a simplification and this particular point does not need to coincide with the true optimal BBP, 38 compare also BBP 4 of the MB surface in Fig. 4. In order to locate the optimal BBP associated with the reactant valley of this reaction, our algorithm was interfaced with the TURBOMOLE code. 76The electronic structure calculations carried out to locate the optimal BBP were done using the B3LYP functional 77 and the Triple Zeta Valence plus Polarization (TZVP) basis set. 78rior to the minimization of the σ-function associated with the PES of the sigmatropic rearrangement, we ran an IRC calculation from the transition state of the reaction to the reactant.Then, we took the molecular configuration with the largest norm of the gradient along the IRC path (i.e., the BBP located on the IRC path) and used this particular geometry as a start point for our minimization algorithm.
As shown in Table I, the minimization converged in 19 steps.Figure 9 reveals that the geometry of the optimal BBP is quite similar to the geometry of the BBP on the IRC path.Thus it follows that taking the latter type of BBP as a start point for the location of the optimal BBP is usually a good strategy.The optimal pulling direction to mechanically induce the sigmatropic rearrangement is displayed in Fig. 10.

VI. CONCLUSIONS
We present an algorithm based on the Gauss-Newton method to locate optimal bond-breaking points on the PES of a molecular system.The algorithm minimizes a positively defined function (the so-called sigma-function) that features zeros at the points on the PES where the Hessian matrix has a zero eigenvalue and the corresponding zero eigenvector coincides with the gradient.We have demonstrated that our algorithm works efficiently for 2D test functions and for multidimensional chemical systems.Optimal BBPs reveal how a force should be applied to a molecular system in order to mechanically induce a reaction using a force with the minimal magnitude.Given the relevance of the concept of the optimal bond-breaking point, we hope that our algorithm will assist in the design of more efficient ways of harnessing mechanical forces in the activation of chemical reactions.e jk is the component k of the vector e j , and F(x) j is the matrix of the contraction of the third derivative tensor with the eigenvector e j .It is used to build the vector w j = F(x) j e j .Thus in our normal case of one zero eigenvalue of the Hessian, t is the orthogonal vector to w j .If we have the w j vector, then we can compute the vector t by the projector, t = c I − w j w T j w T j w j z, (B6) with an arbitrary vector z 0 and a normalization factor c. The matrix I is the unit matrix of dimension N × N. The calculation is not so expansive as the former iteration of the Gauss-Newton algorithm to find optimal BBPs.In Fig. 11 we show a tangent to a green line with Det(H(x(t))) = 0 in the main valley on the MB surface.It is calculated by Eq. (B6).Shown are also two eigenvectors.e 1 is the eigenvector to the zero eigenvalue.
Corollary.If during the integration of the curve of Det(H(x)) = 0, the equation e T j g(x)/||g(x)|| = 1 is satisfied, then this point is an optimal BBP.If in the contrary case, e T j g(x) = 0, then the point is a VRI point.In the last case, the analysis of the diagonalized matrix F(x) j gives the topography of the VRI point.
FIG. 5. Singular NTs through the four VRIs of the MB surface.The VRI points are marked by black bullets.

FIG. 9 .
FIG. 9. Molecular configuration of the optimal BBP of the 1,2-sigmatropic H-shift rearrangement of cyclopentadiene superimposed to the molecular configuration of the BBP along the IRC of the reaction.The configuration depicted with a darker gray scale corresponds to the optimal BBP.

FIG. 10 .
FIG. 10.Molecular configuration of the optimal BBP of the 1,2-sigmatropic H-shift rearrangement of cyclopentadiene.The arrows correspond to the components of the gradient at this point.

FIG. 11 .∂ 3 V
FIG.11.Tangent vector (red) to the line of the solution curve to Det(H(x(t))) = 0. e 1 and e 2 are the eigenvectors of the Hessian at the curve point of interest.e 1 belongs to the zero eigenvalue.

TABLE I .
Behavior of the algorithm described in Sec.III to locate the optimal BBP of the 1,2-sigmatropic H-shift rearrangement of cyclopentadiene.
a The units are au 2 • bohr 4 .b ∆x T ∆x .