Locating saddle points of any index on potential energy surfaces by the generalized gentlest ascent dynamics

The system of ordinary differential equations for the method of the gentlest ascent dynamics (GAD) is tested to determine the saddle points of the potential energy surface of some molecules. The method has been proposed earlier [E and Zhou in Nonlinearity 24:1831 (2011)]. We additionally use the metric of curvilinear internal coordinates. By a number of examples, we explain the possibilities of a GAD curve; it can find the transition state of interest by a gentlest ascent, directly or indirectly, or not. A GAD curve can be a model of a reaction path, if it does not contain a turning point for the energy. We further discuss generalized GAD formulas for the search of saddle points of a higher index. We calculate diverse examples.


INTRODUCTION
The potential energy surface (PES) [1,2,3] and the chemical reaction path [4] are the basic concepts for the theories of chemical dynamics. The PES is a continuous function with respect to the coordinates of the nuclei in an IR N , thus it is an N -dimensional hypersurface. It should have continuous derivatives up to a certain order. If Cartesian coordinates are used, then N = 3 n, however, we may also use N = 3 n − 6 nonredundant, internal coordinates for an natomic molecule. Both cases need some special care.
The PES can be seen as formally divided in catchments associated with local minimums [1]. The first order saddle points (SPs) or transition states (TSs) are located at the deepest points of the boundary of the basins. TS and minimum correspond to stationary points of the PES. Two minimums of the PES can be connected through a TS via a continuous curve in the Ndimensional coordinate space, describing the coordinates of the nuclei. The curve characterizes a usual reaction path. However, if one looks for dynamical trajectories then also SPs of a higher index can be touched. One can define many types of curves satisfying the requirement of a general reaction path. In this paper we test the 'gentlest ascent' pathways proposed by E and Zhou [5] and already used by Samanta and E [6], by Li et al. [7], and the authors of this paper in a preliminary, theoretical letter [8], and by Zeng et al. [9] in a comparison. The name 'gentlest ascent' was coined in a very early paper [10].
Of course, the location of a TS of index one is the second oldest problem for the investigation of a PES, besides the search of minimums. Many methods were developed in the last 80 years. For reviews see refs. [11,12]. The GAD method is a further ansatz. In this paper we test its applicability for molecular PES. A comparison with other search methods for the TS search or for the determination of general reaction pathways is a next task. TSs of higher index emerge often [13,14,15,16,17,18] on the PES of medium and larger molecules. They may be important for dynamical trajectories, as well as for the comparison with conical intersections [19]. To follow GAD-trajectories can be used to determine such TSs.
Let be V (q) an N -dimensional surface over IR N , let grad(q) be its gradient vector, ∇ q V (q), and let H(q) be the matrix of its second derivatives, the Hessian ∇ q grad T . The reason for introduction trajectories by a 'gentlest ascent' dynamics (GAD) is the purposeful search for SPs of index one, in a first step, or for SPs of a higher, fixed index, for example for 'summits' of index two. The index-1-GAD curves are trajectories along the following system of 2N coupled differential equations [5,8] with the initial conditions and v(t o ) = grad(t o ) .
The meaning of the two first equations is the following: If v is an eigenvector of H then the right hand side of eq.(2) is zero. Thus no change to v emerges. If additionally the eigenvector is the one belonging to the smallest eigenvalue near the SP, and the gradient is parallel to this eigenvector, then the ascent to the SP in eq.(1) is along this direction. This is the search along the SP-col. And at the SP the gradient disappears, so the dynamical system converges there. In ref. [5] is shown that the index one SPs are stable fixed points of the system (1,2). Further aspects of the action of GAD are explained in SubSect. 6.1. Usually, if one starts near a minimum in the corresponding minimum well, then the dynamical system turns itself the v vector to the smallest eigenvector; and we find convergence to a corresponding SP. SPs of index one are usually the attractors of the dynamical system (1,2). But we will also report other cases for the radicals CH 3 O and CH − 5 in Sect. 6. The use of the guide vector v in GAD is a generalization of the method of eigenvector following [20,21].
Recently, an application of the GAD to Physical Chemistry is demonstrated for the point defect activity on a surface [6,7]. We will apply GAD to single molecules in an empty space. The description of the PES of an n-atomic molecule is usually done in two different versions: (i) Cartesian coordinates are used. But here emerge the so-called zero eigenvalues of the Hessian. It means that the movement of the molecule as a whole does not change the force constants of the Hessian. If parts of the corresponding eigenvectors point into the direction of v there then the system of the GAD differential equations (1,2) can become stiff. The way out is a projector which is given in Sect. 2.
(ii) One only uses 3n − 6 internal coordinates. They have to be nonredundant, however, they are curvilinear. The inclusion of the curvilinear metric of nonredundant internal coordinates is given in Sect. 3. Formulas for the application of the GAD to SPs of higher index are given in Sect. 4. For every next index we need N additional equations of the GAD system of differential equations, and a further guide vector. In Sect. 5 we report the numerical methods used, and in Sect. 6 we present results of calculations. GAD finds the known test SPs of different index, for analytical surfaces and for different molecules: Lennard-Jones (LJ 7 ) cluster, and the radicals H 3 CO, NH 5 , and CH − 5 . Thus we demonstrate that the GAD method works; however we show also some drawbacks. In a last Section we discuss the results.

PROJECTION OUT OF THE ZERO EIGENVECTORS
Let z i , i = 1, ..., 6, be the six so-called zero-eigenvectors of the 3 n×3 n Hessian matrix. The z i are the fix directions under translations and over-all rotations of the molecule. We form the dyadic products with the 3 n × 3 n unit matrix I. The six single projectors can form the full projection matrix The ansatz P v does not contain any parts of the zero-eigenvectors; it projects out these parts out of vector v. Thus, if we use the 'projected' Hessian [22] in eq.(2) then the GAD system avoids roundabout ways into the overall translation or rotation directions. In the application below in Sect. 6.2 to Lennard-Jones clusters, we simply calculate the six zero-eigenvectors of the Hessian at an initial node, and we use these for the projection. It works, although the rotational eigenvectors depend on the current node.

METRIC
We assume N=3n-6 nonredundant internal coordinates for an n-atomic molecule. Now the corresponding metric relations are used [11,23], because normally, the PES is represented in curvilinear coordinates. The calculations are made with the modified ansatz (g kj ) is the covariant metric matrix, but (g ij ) is the contravariant one. In contrast to eqs.(1,2) we use the summation convention: over a double index, an upper and lower one, we have to sum up. Thus when an index variable appears twice in a single term it implies the summation of this term over all the possible values of the index. The summation over a set of indexed terms in a formula is used to achieve a notational brevity. Note that the upper indices are not exponents but are indices depicting a contravariant vector or tensor. We assume that q and v are contravariant vectors, however, grad is a covariant vector, and H is (an approximation) of a twofold covariant tensor (at least near the SP of the PES) [11,24]. Thus, we can form the scalar products v j grad j in eq. (1), as well as v k H kj v j in eq.(2), but for the norm in the denominator in eqs. (1,2) we have to include the metric matrix. For the first summand on the right hand side in eq.(8), the gradient component, we also have to 'change' the character to adapt it to the left hand side vector using the inverse metric matrix. This also happens formally for the vector H kj v j in eq.(2), correspondingly in eq.(9), which is covariant like the gradient.
Eq.(9) contains the (contravariant) eigenvector problem for the covariant tensor H: To make the characters equal, we have to treat the problem for the eigenvalue λ In (2) follows that for an eigenvector v the steplength becomes zero. Contrary in (1), if v becomes parallel to the contravariant gradient vector g ij grad j = kv i , or the covariant component of v becomes parallel: grad j = kg ij v i , (which is the usual behavior on the col at the SP) then the step becomes kv i along the corresponding eigenvector with steplength k =(norm of the gradient), and it converges along the valley to the SP. From the point of view of Algebra, eq.(1) is the 'mirror' transformation of the gradient to a plane which is orthogonal to the v-direction.
Note that the 'gentlest ascent' dynamics does not always ascend! This happens, e.g. in Fig.2 in ref. [8]. There v is furthermore the col-vector direction to the smallest eigenvalue, and near the eigenvector, but the dynamical trajectory comes from the ridge near the SP downhill. Then the eq.(2) brings again zero steps for the v-part of the system, but in eq.(1) the second summand disappears, and the first summand, the minus gradient, now steps downhill the ridge to the SP of index-1. Here, the examples of SubSect. 6.2 (Figs. 3 and 4) are also of this type. They show turning points of the energy profile.

SEARCH OF SPs OF A HIGHER INDEX
SPs of index-2 are interesting for Chemistry [13,14,15,16,17,18]. The field of SPs of index two is a vibrant and active one. Nowadays, Minyaev et al. [16] found that an SP of index two can be part of an IRC downhill from an usual SP of index one. In the Gamess-US program [21] one can use the commands nneg=2 and method=schlegel to search for an SP of index two, for example. It uses the quasi-Newton-Raphson optimizer of H.B. Schlegel [21]. E and Zhou [5] made a propose to search SPs of higher index. We use the SP-of-index-2 ansatz from ref. [8], again modified by the nonlinear metric tensors: always for i = 1, ..., N = 3n − 6. It is a system of 3N differential equations.
The initial values are where w o is not unique. The reasoning for the introduction of such a system of equations will come true with its success: The two guide directions v and w are now used analogously to search a summit with its two negative eigendirections of the H matrix. This GAD system searches purposefully an SP of index two. The idea is from the Lanczos diagonalization procedure which is a well established diagonalization algorithm, see ref. [25] and references therein with respect to the numerical problems. Note, that sometimes this GAD 2 does not converge to an index-2 SP if the vector w turns into direction v. The two eqs. (12) and (13) are uncoupled; sometimes they do not lead to correctly orthogonal vectors v and w. This is avoidable if one uses a Schmidt re-orthogonalization (SO) of w against v at every step. Results are discussed in Sect. 6. A more sophisticated ansatz is the following: We define the projection operator to vector v by a dyadic product forming an (N × N )-matrix An analogous definition should hold for P w . Be G −1 the inverse metric matrix (g ij ), G the usual metric matrix (g ij ), and I be the unit matrix, then we put the eqs. (11), (12), and (13) to and where again SO implies a Schmidt orthogonalization in the metric sense to previous guide vectors. Correspondingly the initial vectors (14) have to be orthogonal with and the value for the Rayleigh-Ritz relation for v has to be in relation to that of w by ev 1 ≤ λ v (q) ≤ ev 2 ≤ λ w (q) ≤ ev 3 for every current node q of the GAD 2 trajectory. In diverse tests, see Sect. 6.1, we sometimes found cases where these relations are not allways fulfilled; however, along the pathway they usually stabilize by themselves and converge with the correct relation of v and w.
If one searches SPs of a still higher index s, the system has to be extended correspondingly up to the (s + f irst) line for the guide vectors (21) A final remark on the above set of equations to the generalized GAD method: If we are interested in the location of a stationary point of index N where the number is the maximum of the dimension of the space then the generalized GAD equations are reduced to the steepest ascent equation. The proof is trivial, we take either eq. (9) or eq. (13), In eq. (22) the resolution of identity has been used. The generalized GAD method to locate a stationary point of index N is independent of the set of We summarize, the generalized GAD is an algorithm to locate a stationary point of any index, say s, on a PES, with 1 ≤ s ≤ N. Given a set of s orthog- Note that the contravariant version of the gradient vector has to be used, and the orthogonality relations are to be taken with the metric. The generalized GAD curve, If the GAD s method converges then it stops at an SP of indes s. If s = N the algorithm is just the steepest ascent method.

NUMERICAL METHODS
The GAD equations are of a form that is amenable to a direct numerical solution: The calculation of the GAD trajectories is done with a high order Runge-Kutta method like in ref. [8]. The performance of such methods is very good and they converge as fast as possible to a correct solution. The step length along the pathway is automatically selected. It usually costs some additional effort. However, on the other hand, in some regions the step length becomes large; and one avoids unnecessary steps. The method of choice here is the explicit Runge-Kutta method of order 8(5,3) [26], subroutine DOP853. Because this paper is a first test of the GAD method for molecules, we are interested in exact pathways. We use here this very exact procedure. That often causes a high number of steps. The same reason holds for not working with updates of the Hessian. In later, more 'practical' calculations, one should use update steps, as well as a coarser Runge-Kutta method.
For the PES of interesting molecules, the GamessUS program is used with functions of the B3YLP 6-311G(d,p) type [21].
6 N -DIMENSIONAL ANALYTICAL SURFACES AND CHEMICAL SYSTEMS

N -dimensional Surfaces
Before dealing with chemical systems, we test the generalized GAD method as described in Sect. 4 using analytic surfaces as a representation of an Ndimensional PES. For this aim the test will be done on the Rastrigin and Ackley functions [27]. Note that for the two examples the G matrix is the unit matrix in the eqs. (16), (17), and (18).

Rastrigin Energy Surface
The Rastrigin function is a non-convex function used as a performance test problem for optimization algorithms. It is an example of a non-linear function. The 2-dimensional case was proposed by Rastrigin [27] and was generalized by Mühlenbein et al. [28]. The search of stationary points of this function is a difficult problem due to the large search space and the large number of stationary points. In Fig.1 we show the 2-dimensional version of the Rastrigin surface. The Rastrigin function is where N is the dimension of the system, q T = (q 1 , . . . , q N ), 1 T = (1 1 , . . . , 1 N ), and the vectorial function c (q) T = (cos (2πq 1 ) , . . . , cos (2πq N )). The function is usually used in the hypercubic domain q i ∈ [−5.0, 5.0] for i = 1, . . . , N . It is easy to calculate gradient vector and Hessian matrix. The Hessian is diagonal at any point q of the domain. Each pair of the set of eigenpairs of the Hessian matrix, for i = 1, . . . , N . Note that the eigenvectors are the columns of the unit matrix. This result is due to the fact that the Rastrigin function does not have coupling terms between variables, like q i q j .

Fig. 1 2-dimensional Rastrigin surface
The global minimum is located in the point q gm = 0. In this point V (q gm ) = 0.0, and the Hessian matrix is a degenerated matrix of order N , because The properties of the Rastrigin function make it easy to analyze the behavior of the GAD algorithm to locate a stationary point of an index higher than one. If we are interested to locate a stationary point of index, say, s, where, 1 ≤ s ≤ N , then the set of the initial guide vectors, v 0 i s i=1 , are selected to be equal to the subset of eigenvectors of the Hessian matrix showing the highest contribution in the absolute value of the initial gradient vector, when the gradient is represented by the set of eigenvectors. In more detail, if at the initial point, q 0 , the absolute value of the gradient components in the eigenvector space representation satisfies the sequence then we use the subset of eigenvectors z 1 , . . . , z s for the initial vectors v 0 1 , . . . , v 0 s .
With the initial conditions for the set of guide vectors, the GAD equations take their simplest form. The initial guide vectors are eigenvectors of the Hessian, and they are constant in the domain of interest. Consequently, the equality dv i /dt = 0 holds for i = 1, . . . , s, and it is satisfied during the GAD search process. As noted in Sect. 4, the tangent vector, dq/dt, is a function of the guide vectors through the Householder transformation, eqs. (11) or (16), but in the present problem these guide vectors are s columns of the unit matrix at every point, and it is The search direction ascends in the first s coordinates because the gradient components have the + sign, and it descends in the remainder N − s components, where the − sign multiplies the corresponding components.
Using the numerical Runge-Kutta algorithm [8,26] described in Sect. 5 we will locate a stationary point of index 5 on the 100-dimensional Rastrigin surface starting from a node near the global minimum, q gm . The selected initial point is The located point by the Runge-Kutta algorithm is q i = 0.5025 for i = 1, . . . , 5 where the remainder of the q vector elements are zero, as expected. The gradient vector is zero. The eigenvalues of the Hessian at this point are h − = −392.734 and h + = 396.784 which are multiply degenerate in the orders 5 and 95 respectively. The result confirms that the located point is a stationary point of index 5. Note that for this special example, the GAD curve does not have a turning point as defined in ref.
between the initial and the converged point but not including these two points. In this derivation eq. (26) has been used. The result indicates that the GAD curve follows the steepest ascent curve and this is due to both the special structure of the Rastrigin function, eq. (23), and the selected initial conditions.

Ackley Surface
The Ackley function was proposed for the same purpose as the Rastrigin function, as a test of new optimizer algorithms [27]. It is defined as where N is the dimension of the problem and a, b and c are parameters. We use the parameters a = 20.0, b = 0.2 and c = 2π. The meaning of the q, and 1 vectors is the same as the one given for the Rastrigin function, see subsection 6.1.1 and it is c (q) T = (cos (cq 1 ) , . . . , cos (cq N )). The function is characterized by a nearly flat outer region, and a large hole at the center. The Ackley function has a useful search region in the hypercubic domain q i ∈ [−2 a, 2 a] for i = 1, . . . , N . It is easy to calculate gradient vector and Hessian matrix. Opposite to the Rastrigin surface the Hessian of the Ackley surface has coupling terms between the variables.
The global minimum is located at q gm = 0. There holds V (q gm ) = 0, however, the gradient vector becomes singular grad T (q gm ) = (∞, . . . , ∞), as well as all the elements of the Hessian are indeterminate. This implies that the Ackley surface, eq. (28), is a continuous function but not derivable in the point q gm = 0. Around the minimum point the function has the form of a cone with the vertex located at 0. In Fig.2 we show the 2-dimensional version of the Ackley surface. Regarding Fig. 2 we may see that the topography of this surface is close to that proposed as the PES model of a protein folding rearrangement, see e.g. [29]. If we take the Ackley surface as a model, the native conformation of the protein is located in the global minimum, q gm , and the minimums around the deep valley are related to the conformations associated to the denatured states of the protein. These minimums are connected through a reaction path to the global minimum. They have direct access to the native conformation. The remainder of denatured conformations is related to other minimums. They are called dead-end or kinetic-traps. To reach the global minimum related to the native conformation, any conformation beginning in the dead-end region of the PES must first pass through the conformations of denatured states that are located around the deep native valley. The features of the PES of a protein folding [29] are well represented by the model of the Ackley surface given in eq. (28).
The study of the Ackley surface using the GAD method reveals some limitations in the selection of the initial conditions. Using the same procedure in the selection of the initial conditions as that described in the previous subsection, we attempt to locate a stationary point of index 2 on the 4-dimension Ackley surface. We start at the initial node q and v 0 2 are z 0 1 and z 0 3 , respectively. However, the degeneracy of z 0 3 eigenvector precludes to use eqs. (17) and (18), because they force a partition of the eigenvector space which is not correct.
Note that the generalized GAD equations to find a stationary point of index say, s, imply a subspace spanned by the selected guide vectors, {v i } s i=1 , and the corresponding orthogonal subspace. If at the initial point the selected guide vectors are some eigenvectors of the Hessian matrix, the selection should take into account the possible degeneracy of the eigenpairs. The partition of the guide vectors should be in accord to the degeneracy of the Hessian matrix at the converged stationary point. If at the converged stationary point a guide vector, which is already an eigenvector of the Hessian matrix is degenerate with another eigenvector not included in the subspace of the guide vectors, then the GAD procedure has been converged to a stationary point of index s + 1. The reason is that at the initial node the Rayleigh-Ritz quotient is stationary with respect to the guide-vectors and during the GAD search procedure this stationarity is preserved as much as possible. The present example shows the last situation. The coordinates of the converged stationary point are, q 0 1 = q 0 2 = q 0 3 = 0.6635, and, q 0 4 = −1.1635. The eigenvalues of the Hessian matrix are, h 1 = −8.8893, h 2 = h 3 = −2.9014, h 4 = 1.8146. One easily concludes that the converged point is a stationary point of index 3 rather than 2 as it was required in the construction of the GAD search equations.
Finally, we present a case where we locate a stationary point of a desired index starting from a point such that the Hessian matrix is highly degenerate. The surface is the 100-dimension Ackley function. We are now interested in the location of a stationary point of index 5 starting from the node with coordinates

Lennard-Jones cluster
Before coming to molecules, we choose a cluster of n = 7 abstract LJ kernels (σ = 1, = 1) for the GAD experiment. It has an analytical PES. Energy, Gradient, and Hessian can be explicitely calculated, without approximation. Geometries of LJ 7 for the start and for comparisons are taken from the Doye homepage [30], for the formulas see an Appendix. We use 21 Cartesian coordinates in the full 3n-dimensional space. Thus, we have to use the projected Hessian matrix (7) for the 42-dimensional GAD system (2). A test starts near the minimum, MIN 1 , (in the Doye numbering) with small initial deviations for particle 4 (δx = 0.002, δy = −0.001, δz = 0.000 3). GAD finds the index-1 SP (number 3 in the Doye numbering), see Fig.3. Up to the cumulative steplength R = 1 one needs 80 steps and obtains an accuracy |grad| max = 0.16 of the largest component of the gradient. To get any accuracy, one has to carry on the calculation: for example, for |grad| max = 1.6 * 10 −5 , GAD needs 517 steps, and additionally 14 steps for step length corrections. Less than the tenth number of the steps leads to the neighborhood of the SP, the remainder is used for the refinement of the SP. So, if one stops GAD 'near' the SP of interest, and if one jumps there to a Newton-Raphson step, one can avoid the long convergence process of the DOP853 routine.
Note that bullets on the profile of Fig. 3 are drawn at fixed R steps. They do not reflect the steps of the DOP853 routine. Of course, there are also many start values near the minimum which lead to 'dissociative' geometries of the clusters: LJ 6 +LJ 1 , or LJ 5 +LJ 2 , for many combinations of the seven particles. Such structures are not really interesting here. We search for an 'isomerization' like in Fig. 3. A 'good' choice of the initial deviation from the minimum is a critical aspect of the GAD method. Shown is the energy profile over the reaction coordinate, R A next test for the GAD method is the search of SPs of a higher index, along the eqs. (11) to (13). The metric here is the trivial unit matrix; and the additional vector w is orthogonalized to v at every step of the procedure. (Or of vectors v i , 2 ≤ i ≤ s, if s is the index of the searched SP. We then use (s − 1) additional equations like eq.(13).) If we look for an index-2 SP, we need one more vector w, and the index-2-GAD system has 3 * 3n = 63 dimensions.
For small LJ clusters, the description of the numbers and characters of diverse SPs is given by Doye and Wales [20]. So, we can stop here our test after demonstrating that an index-s-GAD procedure can find an SP of index s. We start near a corresponding SP of a deeper index, here the SP i1 , after randomly perturbing the coordinates of the initial node. The initial choice of the second guide vector w can be an eigenvector being orthogonal to the direction of the saddle valley. We automatically perform so many iterated searches up to a success: here we needed 11 iterations for an SP i2 with E = −14.723 units. At R = 2.85 units we get a convergence of |grad| max = 1.4 * 10 −4 . The trajectory has 971 nodes and additionally 47 steps are done for step length calculations. Here, nearly the full trajectory, up to R = 2.25 units, is used to meander over the PES. There is |grad| max = 0.1.
A further index-2-GAD test is done for the deflation ansatz of E and Zhou [5], eqs. (16,17) there. Because the Hessian is symmetric, we can use here 3 * 3n = 63 dimensions only with the two guide vectors v, w. We get the same SP of index 2, in a third iteration, at R = 1.47 units with the accuracy result |grad| max = 0.1, and we need 121 steps. The ansatz works without an additional Schmidt orthogonalization of the test vectors v, w. And so on: an index-4-GAD can find an SP of index 4, see Fig. 4, at E = −13.708 units. Now the search space of the Runge-Kutta method has the dimension 5 * 3n = 105. The energy profile over the trajectory is shown. It enrolls the complicate search with some turning points of energy. Again the region of the SP i4 is already touched at R ≈ 0.6 units, but the remainder of the calculations, up to R = 6 units, is done for the exact convergence to the SP with a threshold of |grad| max < 2 * 10 −5 . The full trajectory makes 3560 steps; note that the bullets in Fig. 4 are drawn at fixed R-steps. They do not reflect the much more steps of the GAD trajectory.
A conclusion of this test is: the generalized Gad can find SPs of any index of the PES. It is a determined, purposeful method. To search for an SP of index s, one should start near an SP of the lower index (s − 1). However, any other initial nodes are possible.

H 3 CO, methoxy radical
The radical H 3 CO is an important, prototypical part of combustion chemistry [31]. A deep exploration of its PES was done by Knyazev 2002 [32], see also references therein. The radical consists of a methyl group bound to oxygen [33]. Its ground state has C s symmetry, thus, the umbrella symmetry C 3v of CH 3 is slightly changed. It is an electron-donating group and can cause an organic compound to become less acidic. Although this group exerts a negative inductive effect, its positive mesomeric effect is, however, greater and thus dominates, and so this group is an electron donating group. The radical has two radical electronic states. In the σ-bond the electron is localized in an orbital collinear with respect to the CO-bond, and in the π-electronic state the electron is localized in the orbital (SOMO) orthogonal to the CO-bond. Note that we do not treat the isomerization to the hydroxymethyl radical [34] here in this work The used method is UHF, and the basic set is B3LYP 6-311G(p,d). The minimum of the methoxy radical is given (in ANG,DEG internal coordinates) in The H 2 CO+H van der Waals minimum is given in Table 2. Note that the B3LYP basis set does not well describe such a van der Waals minimum, however, this is meaningless for the test of the GAD method. The SP number 1 (SP1) of index 1 is given in Table 3. It is the threshold for an escaping H-atom. The SP number 2 (SP2) of index 1 is given in Table  4. The eigenvector of this SP to the col direction is the lowering of the bond of the second escaping H-atom of the full radical. Thus, the SP2 forms an unsymmetric structure for the dissociation of an H 2 part.  Three index-1-GAD trajectories are shown in Fig. 5 where we used formulas (8) and (9) for the GAD method. The trajectories show the three H-atoms for three different GAD runs for the radical H 3 CO. One trajectory finds SP1, one finds SP2, and one goes to the valley of the van der Waals minimum at H 2 CO+H without crossing an SP. The right triple of curves (bullets) depicts the distances of the three H-atoms for a trajectory from the minimum to the SP1 of index one. In the representation of Fig. 5 the bullets now depict the real nodes of the calculated trajectory. One can clearly see the advance of the used Runge-Kutta method DOP853: the steplength is variable. At curved pieces the steps are denser. On a straight pathway are large steps, but near the SP the steps become very short; and the SP is exactly located. The method stops there after convergence to a given threshold.
The region of a valley-ridge inflection (VRI) is indicated here as well. It is a bifurcation (bif) of a valley between the two SPs which is not directly attracted by GAD trajectories. It is detected by Newton trajectories, see ref. [35]. A next search trajectory (thin lines), however, overcomes the barrier between minimum and dissociation to H+H 2 CO, but does not meet an SP. The left depicted minimum in the Figure is H 2 CO but the corresponding third, the single dissociating H escapes over the frame. Thus, this GAD trajectory misses its aim to converge at an SP. It connects two minimum regions without meeting an SP. The calculation is stopped after the trajectory leaves the frame of the Figure. The case opens a gap of GAD: if it converges to a node, then it is a saddle point, but it converges not in every case.
The middle trajectory (dashes + bullets) meets the SP2 of index one. It is the pathway to the (unsymmetric) dissociation to HCO+H 2 . The trajectory shows the advance of the Runge-Kutta method even more. Near the SP the steps become very short; and the trajectory can 'spirally' include the SP which again is exactly located. (In the case of a fixed step length of the Runge-Kutta method, however, the GAD trajectory would circumvent the SP2 in an infinite, false 'limit cycle'.) Note that also the SP2 is a result of a direct GAD search for SPs of index one, compare ref. [32] where this was not possible. The GAD trajectory goes across the 'ascending' ridge from SP1 over which the direct dissociation should go to HCO+H 2 . We find the SP2 for this reaction by a direct trajectory which 'misses' the SP1. (Compare ref. [36] for another 2dimensional picture of such an SP constellation.) There are now the SPs of index one for the reactions and with ref. [32] the rearrangement which is also found coming from the minimum H 3 CO. However, a direct symmetric dissociation should go over a putative SP of index two, corresponding to Fig. 5 of ref. [37]. It forms a 'summit' in the corresponding 2-dimensional section. It is a task for GAD to search this higher index SP. Diverse GAD 2 search tests, however, go to regions with an extremely negative eigenvalue of the Hessian, however, the norm of the gradient does not converge to zero: this is a hint for a conical intersection (CI) point nearby   [19]. Usually, a CI point has two negative eigenvalues, like an SP of index two. Thus, the search along a 'full two-dimensional direction' may find CI points, as well [19]. The difference to an SP is that at the CI the gradient is not zero and the GAD system will not converge. But a strong indicator for the corresponding region is, behind the large negative first eigenvalue, that the Runge-Kutta procedure extremely shortens the steplength. The reason is the large eigenvalue indicating a strong curvature of the PES. This causes a short steplength. For example, we find a CI point at a node with coordinates given in Table 5. The Hessian at the CI node has two negative eigenvalues: λ 1 = −5.42, and λ 2 = −0.748. The first one is larger than 'normal' by a factor of 20. Three index-2-GAD trajectories leading to CI points are shown in Fig. 6. None of the trajectories cross the region of the VRI bifurcation because only one H-atom moves in this region, however, the VRI concerns all three H-atoms. Usually, the CI points form a high-dimensional seam [19] on the intersecting PESs, so the CI points which we met may be only representants of this seam, or even of different seams.
A next test with the formulas of E and Zhou [5], eqs. (16), (17) there, was done for a search of SPs of index two: Here the vectors v and w stay orthogonal, however, we again do not find any stationary points. In contrast, again CI points are localized.
Note: using the GAD ansatz, Eqs. (11)-(13), we found that sometimes the vectors v and w do not stay strictly orthogonal, but the search can be successful, nevertheless, without a strict orthogonalization of both vectors. The Schmidt orthogonalization (SO) makes the GAD method safer; however, it may also disturb the steplength search of procedure DOP853. That can cause a bigger effort, besides the effort of the orthogonalization itself.
A last observation in this Subsection: if the GAD trajectory crosses a region of a VRI bifurcation then at least one curvature of the PES is nearly zero: To 'feel' the true pathway, in this case, the DOP853 procedure also uses very tiny steplengths. Thus, one needs a large effort to cross such a region. It is quite the same like in the contrary case of large curvatures at CI nodes. But here it seems to be a gap of the DOP853 procedure.
Conclusion: Apparently, an index-2-SP of the PES does not exist in the interesting region, compare also ref. [37]. On the contrary, many Conical Intersections (CI) emerge there. They can also be indicated by the GAD method. In contrast to the LJ clusters in SubSect. 6.2 where one single, unique PES exists, the relations in real molecules are quite more complicated.

NH 5
Is the NH 5 structure possible? NH + 4 is formed by the spare pair of electrons on the nitrogen atom bonding to an H + ion. This is a co-ordinate bond where both electrons come from the nitrogen. There are now no electrons left on the nitrogen to form either another regular covalent or a co-ordinate bond. To add another hydrogen, one had to form a co-ordinate bond from the hydrogen, which needed two electrons there, a hydride (H − ) ion. Even when one can get these to interact, this forms an ionic bond instead of a pure covalent bond. It is an 'ammonium hydride' salt, rather than a true NH 5 compound. Unfortunately this salt is stable only at a high pressure. When one decreases the pressure, it simply decomposes back into regular ammonia and hydrogen again.
Here we treat the NH 3 ...H 2 complex under a pyramidal inversion of the ammonia molecule, see Ref. [15]. The used method to calculate the PES is RHF, and the basic set is 6-311G(3d,2p). The van der Waals minimum NH 3 ...H 2 is given by its z-matrix (in ANG,DEG internal coordinates) in Table 6. A calculation along the GAD 2 ansatz eqs. (11)-(13) starts anywhere nearby to successfully search the SP of index two [15], see Fig. 7. Note: we circumvent possible SPs of index one by passing directly from the minimum region to the SP 2 . There is a recent paper on this radical, see ref. [16]. The used method is RHF, and the basic set is 6-311G(3d,2p). In Table 7 we report an SP (in ANG,DEG internal coordinates) with D 3h -symmetry [38].  In a first result we find this high lying SP 1 of index one if we start midway on the ridge after a valley-ridge inflection (VRI) uphill to the SP 1 , see ref. [16]. We start here above the VRI region to avoid endless tiny steps, see the former Subsect. 6.3 for an explanation.
A second result for CH − 5 is a so-called roundaboutway: We start anywhere in the plane of a kind of van der Waals minimum of the molecule. However, the trajectory comes back to a minimum structure in the 'plane' of a dissociating H from the molecule CH 4 . There is a turning point (TP) at the emerging mountains of the PES; and then the trajectory returns, like in the strange case of the MB test [8]. At the TP we find the first three eigenvalues of H (in Bohr,Radiant system) -0.009, 0.001, and 0.018 . Thus: the GAD trajectory 'finds' a region with one negative eigenvalue but turns back. GAD goes over the TP further and further, but then back to the minimum. GAD is a dynamical system method. It is programmed that the system trajectory locally holds a given direction. GAD may not necessarily converge globally [5].
A third result concerns the search of an SP of index two. We again start anywhere in the plane of a kind of a van der Waals minimum of the molecule. GAD 2 works successfully, it ends at the SP i2 given in Table 8. At the SP the first three eigenvalues of H are -0.004, -0.001, and 0.001 . The small first eigenvalues show that the PES of this radical is extremely flat, see [16]. The run needs the following effort: function calls 8478, steps 671, accepted steps 365, but many steplength searches: 306.

DISCUSSION and CONCLUSION
We have shown the applicability of the GAD method to small molecules. Two different kinds of coordinate systems are used: first Cartesian coordinates including the projection out of the zero eigenvectors of the Hessian in the GAD ansatz, or, second, internal coordinates where the curvilinear metric has to be included in the GAD formulas.
The GAD method is a good 'growing string' method for the exploration of a PES. We found all SPs of index 2 of interest, or the method indicates a CI case. The sticking point is the handicap of the initial node, and the initial guide vectors. Thus, one has to guess an abstract region of attractivity of the corresponding SP (which is, of course, usually unknown). To looking to an SP of interest we have to aim in the correct direction. In ref. [7] the authors worked with an 'available good initial guess': a coarse minimum energy path is determined before by a double-ended string method. If no SP of index s is found, then we cannot say that there is no SP of the searched index. Because GAD is 'only' a growing string method, we have no indication of how many SPs of a given index it overlooks. In the case of Newton trajectories [11] we know that they connect SPs of an index difference of one in the regular case. This does not fit here. We report a GAD 2 trajectory from the minimum region to an SP of index 2.
In our tests not in every case the searched kind of SP is really found; the GAD trajectory can have turning point(s) and can go 'wrong' on a way back. However, if one knows what to look for then the method can be useful. If the true region of the SP is found, then an exact localization of the SP is usually not very effective. One may jump automatically to another method, in such a case. The convergence of the dynamical system of GAD is slow, compare also ref. [9]. However, if GAD converges, then it converges to every given threshold.

Appendix: Hessian matrix of the Lennard-Jones cluster
The potential energy surface (PES) of an n-atomic Lennard-Jones (LJ) cluster is: V (x 1 , y 1 , z 1 , ..., x n , y n , z n ) = 4 n i=1 n j=i+1 σ 12 ( 1 r ij ) 12 − σ 6 ( 1 r ij ) 6 where the double sum can be shortened to i<j , and is the square of the distance between atoms i and j. , σ are parameters of the current cluster.
V (q) is a 3n-dimensional surface over IR 3n . Let grad(q) be its gradient vector, ∇ q V (q), and let H(q) be the matrix of its second derivatives, the Hessian ∇ q grad T . We could not find the calculation of the Hessian in the older literature; thus, we give it here in the appendix. We get for the distance between two atoms A single component of the gradient becomes If we put the abbreviation T 1(ij) =: −12σ 12 ( 1 r ij ) 14 + 6σ 6 ( 1 and if we set for simplification of the summation T 1(ii) = 0, then we can write We further define the abbreviation T 2(ij) =: 12 * 14σ 12 ( 1 r ij ) 16 − 6 * 8σ 6 ( 1 r ij ) 10 .
Again we set T 2(ii) = 0. Then single components of the Hessian are for the same atom and and analogously one gets H(x i , z i ) and H(y i , z i ) and their symmetric counterparts. If different atoms with numbers i and j are involved, then we get and and analogously one gets H(x i , z j ) and H(y i , z j ) and their symmetric counterparts.