Within the harmonic approximation to transition state theory, the rate of escape from a reactant is calculated from local information at saddle points on the boundary of the state. The dimer minimum-mode following method can be used to find such saddle points. But as we show, dimer searches that are initiated from a reactant state of interest can converge to saddles that are not on the boundary of the reactant state. These disconnected saddles are not directly useful for calculating the escape rate. Additionally, the ratio of disconnected saddles can be large, especially when the dimer searches are initiated far from the reactant minimum. The reason that the method finds disconnected saddles is a result of the fact that the dimer method tracks local ridges, defined as the set of points where the force is perpendicular to the negative curvature mode, and not the true ridge, defined as the boundary of the set of points which minimize to the reactant. The local ridges tend to deviate from the true ridge away from saddle points. Furthermore, the local ridge can be discontinuous and have holes which allow the dimer to cross the true ridge and escape the initial state. To solve this problem, we employ an alternative definition of a local ridge based upon the minimum directional curvature of the isopotential hyperplane, κ, which provides additional local information to tune the dimer dynamics. We find that hyperplanes of κ = 0 pass through all saddle points but rarely intersect with the true ridge elsewhere. By restraining the dimer within the κ < 0 region, the probability of converging to disconnected saddles is significantly reduced and the efficiency of finding connected saddles is increased.

## I. INTRODUCTION

The determination of saddle points is routinely used to find reaction mechanisms and calculate reaction rates in solid state systems and surface chemistry. The underlying approximation is that the bottleneck region of a reaction is localized to a transition state (TS) region in the neighborhood of a first-order saddle point on the potential energy surface. A harmonic expansion of the potential around the saddle point and reactant minimum can be used to quantify the reaction rate within the harmonic approximation to transition state theory (HTST).^{1,2}

In the context of exploring rare-event dynamics, quantification of the low-energy saddle points that lead from a reactant state gives the escape rate and branching ratio to possible product states. With this information, the adaptive kinetic Monte Carlo algorithm can be used to model the state-to-state dynamics over long time scales.^{3,4} Here, the definition of a state is the set of points which minimize to a single minimum, as followed by a steepest descent trajectory. In other words, the minimum geometry is the inherent structure of the state.^{5}

The focus of this paper is to improve the computational methods designed to determine saddle points that connect, in the sense of a steepest descent path, to a specified initial state.

## II. DIMER METHOD

The dimer method^{6} is one such minimum-mode following algorithm (see also Refs. 7,8) designed to find saddle points on a potential energy surface. In terms of a dynamical system, this is to find first-order saddle points of the gradient flow,

where

*V*(

*x*) is the potential, and −∇

*V*(

*x*) is the force. In the region of the potential with at least one negative curvature mode, which is the focus of this paper, the dimer method employs a modified dynamical system,

where

**H**(

*V*) associated with the smallest eigenvalue.

^{6}

*x*follows Eq. (2) to climb up the min-mode direction and descend in all perpendicular directions. If the min-mode direction is unstable (the associated eigenvalue is negative) under Eq. (1), then Eq. (2) stabilizes it. At a first-order saddle there is only one unstable direction, which is the min-mode direction

### A. Dimer escape from a local minimum

In the context of calculating the HTST rate from an initial state, the objective is to find low energy saddle points on the boundary of the state. Equation (2) guarantees dimer convergence to a first-order saddle, but it does not guarantee that the saddle is connected to the initial state, even if the dimer search is started well-inside in the initial state basin. Disconnected saddles are not directly useful for calculating the escape rate and are thus considered failures of the method when it is used to calculate an HTST escape rate.

To understand the nature of the disconnected-saddle problem, we need to define two types of ridges: the true ridge and the local ridge as seen by the dimer method. The true ridge is the boundary of a state, or alternatively the separatrix in the dynamical system of Eq. (1). If the dimer method did not cross the true ridge, it could not find disconnected saddles. The dimer evolution, however, is guided by local ridges, defined by the set of points with *F*^{‖} = 0 and at least one negative eigenvalue of **H**. Since the *F*^{‖} direction (along

The two-dimensional potential shown in Fig. 1 is used to illustrate the disconnected-saddle problem. The analytic form of the potential is defined in Ref. 6, modified so that the central Gaussian hill is turned into a basin by changing the sign of its amplitude, *A*_{1}. The contours of the potential are shown as black lines. There are four minima (M1-4) and four saddles (S1-4). The dark shaded regions around the minima have no negative curvature mode; in these regions Eq. (2) is not applied; instead, the dimer is moved directly up the lowest curvature mode to rapidly escape the convex region. The colors indicate the dot product between the min-mode

*V*(

*x*) directions.

Outside the dark shaded regions of the minima, the deep blue paths are the local ridges where the force and min-mode directions are perpendicular. The local ridges pass through each saddle point. It is also clear, however, that they are separated into segments and do not enclose the minima. Points at which the local ridges terminate are contained within the magenta contours, defined as regions where the two eigenvalues of the Hessian are both negative. At the termination points, the min-mode direction

The basins of attraction of the dimer, along with the local and true ridges, are plotted in Fig. 2(a). A first order quickmin (QM) optimizer,^{9,10} with a small time step of 1 fs (assuming a mass of 1 amu and potential units of eV and Å), and a maximum step size of 0.05 Å is used for the numerical implementation of Eq. (2), which gives similar basins of attraction as a steepest descent algorithm but converges more rapidly. Saddles S1, S3, and S4 line on the boundary of M1 (white lines) while S2 is disconnected. The basin of attraction of S2, however (light blue region), penetrates through the boundary of M1 and even lies close to the minimum. A similar problem exists for state M2. M2 has only one connected saddle, S2, but inside the boundary of M2 there are attractors to S4 (yellow) and S1 (green) and the area of these disconnected attractors is comparable to that of the connected saddle S2 (light blue). The boundaries between these basins of attraction intersect the termination of the local ridges (red lines), demonstrating that the penetration of attractors to disconnected saddles occur at breaks in the local ridge.

### B. Scaling the parallel and perpendicular forces

Since breaks in the local ridge allow the dimer to cross the true ridge, one strategy to prevent the dimer from finding disconnected saddle is to keep it away from the local ridge. This can be done by modifying the dimer evolution of Eq. (2) so that the dimer climbs up the potential along the min-mode direction more slowly than it descends the potential in perpendicular directions. Then, the dimer will be more likely to approach the saddle along the min-mode direction, where the true and local ridges are coincident, and less likely to find holes in the local ridge away from the connected saddles. This can be implemented with a scaling parameter γ between the perpendicular and parallel forces,

We note that the activation-relaxation technique (ART-neuveau)^{8} and hybrid eigenvector following^{7} methods make use of this idea by separating variables and relaxing in the perpendicular direction before moving up the potential along the parallel direction. Figure 2(b) shows the basins of attraction setting γ = 0.2 in Eq. (3) so that perpendicular relaxation is five times faster than activation parallel to the min-mode direction. Clearly the penetration of disconnected attractors decreases, but the problem is not solved because the holes in the local ridge remain, as well as clear channels of escape. Reducing γ below 0.2 was not found to improve the basins of attraction; it was also found to significantly increase the number of iterations required for the method to converge to a saddle.

## III. κ-DIMER METHOD

Constraining a dimer search to the basin of the reactant state is equivalent to keeping it away from the true ridge, except in the neighborhood of a saddle points. The true ridge is a nonlocal quantity and cannot be determined solely upon local information of the potential energy surface. It is possible, however, to define a region that contains the true ridge, with only local information. If we can define such a ridge-neighborhood, we can keep the dimer search from entering it and also prevent it from crossing the true ridge. Additionally, we will show that keeping the dimer search away from the true ridge can reduce the search space and increase the efficiency of finding saddle points. This will only work, however, if the boundary of the ridge region also intersects the true ridge at all first-order saddle points so that it does not prevent the dimer method from finding them. In the following, we present our choice of repulsive region based upon the local minimum directional curvature of the isopotential hyperplane κ, and show how to modify the dimer equations of motion to avoid this region and increase the efficiency of finding connected saddles.

### A. Isopotential curvature κ

The idea behind using the isopotential curvature is that regions around minima can be identified by contour lines curving around the minimum, whereas at ridges, the contours curve away. A similar concept to the isopotential curvature is widely used in image processing, where it is called the isophote curvature.^{11–13} For images, the quantity of interest is the luminance on a two-dimensional surface, while ours is the potential energy in a high dimensional space. The two-dimensional case will be discussed first so that results from image processing can be borrowed directly.

The black lines in Fig. 2 are the contours of the potential energy *V*(*x*_{1}, *x*_{2}). At any point, there is a tangent direction of the contour,

*V*. The distance along the contour is denoted as

*s*. The isopotential curvature κ is a description of how fast the contour tangent direction changes as moving along the isopotential contour. κ can be calculated from the Hessian matrix and the force,

For the two-dimensional potential discussed previously, regions indicating the sign of κ are plotted in Fig. 3. We find that the κ < 0 region (deep blue) never touches the true ridges (white lines) except at the saddle points. If κ > 0 is set as the repulsive region for dimer searches in this example, the two requirements for a good restraint are satisfied: (i) the κ < 0 region avoids the true ridge; (ii) the boundary κ = 0 intersects the true ridge at all saddle points. The second requirement can be proven in general. Near a saddle the potential energy can be approximated by a harmonic function and κ can be found explicitly. Taking

*a*> 0 and

*b*> 0 near the saddle point (0, 0), then

Requirement (i), however, is not always satisfied. On the true ridge *S*, the force is parallel to the tangent of *S*, by definition, and thus is zero along the normal direction

*V*(

*x*

_{1},

*x*

_{2}) reaches the minimum in the normal direction on the true ridge. Thus the true ridge must be of valley shape. This is a rare situation on typical potential energy surfaces, but it can exist in model potentials and may cause problems in some systems. An example of this scenario will be given below.

Although the κ < 0 constraint does not always satisfy requirement (i), it does strictly exclude regions with more than one negative eigenvalues of **H**, which reduces the possibility of leaving the initial basin and finding a disconnected saddle. This can be shown by taking the two eigenvalues of **H** as λ_{1} and λ_{2} where λ_{1} ⩽ λ_{2}. Then λ_{1} ⩽ ν ⩽ λ_{2} and in the permitted regions where κ < 0, we have λ_{2} ⩾ ν > 0.

In the extension of the isopotential curvature to an *n*-dimensional space

*n*− 1)-dimensional hyperplane. The direction perpendicular to the force is therefore not unique, and our choice of κ is the minimum directional curvature,

where

_{1}⩽ λ

_{2}⩽ ⋅⋅⋅ as the eigenvalues of

**H**and ν

_{1}as

_{1}⩽ ν

_{1}⩽ λ

_{2}. On the isopotential hyperplane, there is always a vector

_{1}and λ

_{2}. Thus,

_{1}⩽ ν

_{1}⩽ λ

_{2}. In high dimension, therefore, restraining the dimer to κ < 0 (ν

_{1}> 0) avoids regions where

**H**has more than one negative eigenvalue, as in the two-dimensional case.

Calculation of ν_{1} and κ can be done with the same dimer rotation algorithm used to find the minimum-mode direction τ, with one modification. The dimer is constrained to the isopotential hyperplane by projecting out the components along the force direction. The extra work to calculate κ approximately doubles the amount of computational work of the κ-dimer method as compared the regular dimer. Thus, the κ-dimer method must be twice as efficient at finding connected saddles, or reduce the number of iterations to find a saddle as a result of the κ restraint, to outperform the regular dimer method.

### B. Equation of motion of the κ-dimer method

To apply a restraint based on the sign of κ, two coefficients are added to scale the perpendicular and parallel forces of the dimer method,

where β is a constant determining how quickly the coefficients switch values around κ = 0. Figure 4 illustrates how the values of γ_{1} and γ_{2} change as a function κ from which three regions can be identified.

When κ is sufficiently negative,

$\dot{x} = -F^{\parallel }$$x\u0307=\u2212F\Vert $ and the dimer moves up potential along the lowest curvature mode. This is equivalent to the evolution of the regular dimer method in convex regions of the potential.Near κ = 0 where γ

_{1}= γ_{2}> 0 we recover the regular dimer method of Eq. (2).When κ is sufficiently positive,

$\dot{x} = F^{\perp } + F^{\parallel } = -\nabla V$$x\u0307=F\u22a5+F\Vert =\u2212\u2207V$ and the dimer moves down the potential away from the region containing the true ridge.

The logic behind these choices is as follows. For κ < 0 it is unlikely to find or cross the true ridge and we allow the dimer to aggressively climb up the potential. This is equivalent to rapidly escaping convex regions of the potential around minima to find a region with at least one negative curvature mode where saddles are present. For κ = 0 the regular dimer method is employed, but for κ > 0 there is a risk of crossing the true ridge and the κ-dimer is relaxed back towards κ = 0. In summary, the aim of the κ-dimer method is to restrain the dimer to the κ = 0 surface as it approaches a saddle. This improves upon the original dimer method which follows local ridges, and can escape to disconnected saddles. Our choices of setting β = 5 and flipping the sign of the parallel force at κ = 0 are arbitrary, yet sensible choices. Future calculations will demonstrate how universal these choices are or if system dependent tuning is required.

In Fig. 3, a trajectory of the κ-dimer is compared to that of the regular dimer initiated from the same configuration point. The κ-dimer converges to the connected saddle S1 following the κ = 0 surface, whereas the regular dimer escapes to the disconnected saddle S2.

### C. Two dimensional examples

First, we revisit the two-dimensional potential of Fig. 2. The new basins of attraction, following the κ-dimer motion of Eq. (7), are shown in Fig. 5. The disconnected problem is solved in this case; at the lower and upper center of the map, the boundaries of the green and light blue regions align with the true ridge. The yellow basin also increases in volume due to fact that the κ-dimer more rapidly climbs out of the minimum towards the κ = 0 surface.

Unfortunately, there are cases in which the κ = 0 surface intersects the true ridge and the κ-dimer is able to find disconnected saddles. Such a counter example is shown in Fig. 6. The potential, which is modified from Ref. 14, has four identical Gaussian basins in each quadrant and two Gaussian hills in the middle aligned along the *y*-axis. The lower hill is slightly larger than the upper one to break the mirror symmetry along the *x*-axis. Considering the initial state containing minimum M1, there are two saddles (S2 and S3) that lie on its boundary. The two hills make the horizontal true ridge points between them become minima along the normal directions of the ridges (valley shape), and thus the κ < 0 regions overlap the true ridges as the red arrows indicate in Fig. 6(a). As a result, Eq. (7) does not prevent S1 attracting some followers crossing the white lines from state M1. But compared with the red circle which is the attraction region of S1 from the regular dimer method, the constraint still shrinks the area of S1 within the boundary of state M1.

### D. Pt heptamer island

Next, we consider the performance of the κ-dimer method in a higher dimensional system, containing an island of seven adsorbed Pt atoms on a Pt(111) surface, as shown in the inset of Fig. 7. The interatomic interactions are described by a Morse potential, with parameters as described in Ref. 15. The heptamer island, as well as three substrate layers, are free to move, giving a total of 525 degrees of freedom. Three additional layers of substrate are frozen at the bottom of the slab to represent the bulk. Dimer searches are started from two sets of initial conditions. In both cases, 1000 initial configurations are generated by displacing one of the least-coordinated atoms on the corner of the island and any neighbors within 3.3 Å (nearest neighbors) of it. In the first set, relatively small Gaussian displacements with a standard deviation of 0.3 Å moved these atoms modestly away from the minimum; in the second set, displacements of 0.5 Å moved the atoms closer to the true ridge, but still within the reactant state.^{16,17} All saddle searches used the conservative QM optimizer with a time step of 1 fs and a maximum step size of 0.1 Å to avoid any discontinuities in the path which can be caused by more aggressive optimizers.^{10} Convergence was claimed when the magnitude of the total force dropped below 0.001 eV/Å. Both dimer methods were started from the same initial points to reduce statistical errors in the comparison.

The success ratio for finding connected saddles and the average number of force calls for the dimer and κ-dimer methods are listed in Table I. For both choices of initial displacements σ, the κ-dimer significantly increases the success ratio. The fact that the success ratio is close to unity indicates that the scenario in which the κ = 0 surface crosses the true ridge, as presented in the previous counter-example, appears to be unlikely in this high dimensional system.

Method . | σ (Å) . | Success . | $\langle \rm {FC} \rangle _1$ $\u27e8 FC \u27e91$
. | $\langle \rm {FC} \rangle _2$ $\u27e8 FC \u27e92$
. |
---|---|---|---|---|

Dimer | 0.3 | 88.1% | 2462 | 4985 |

κ-dimer | 0.3 | 99.5% | 3658 | 5449 |

Dimer | 0.5 | 65.5% | 4030 | 23 452 |

κ-dimer | 0.5 | 99.8% | 3508 | 6425 |

Method . | σ (Å) . | Success . | $\langle \rm {FC} \rangle _1$ $\u27e8 FC \u27e91$
. | $\langle \rm {FC} \rangle _2$ $\u27e8 FC \u27e92$
. |
---|---|---|---|---|

Dimer | 0.3 | 88.1% | 2462 | 4985 |

κ-dimer | 0.3 | 99.5% | 3658 | 5449 |

Dimer | 0.5 | 65.5% | 4030 | 23 452 |

κ-dimer | 0.5 | 99.8% | 3508 | 6425 |

For the conservative initial points (σ = 0.3 Å) the disconnected problem accounts for only 12% of the regular dimer searches and the additional cost of the κ-dimer method does not outweigh the fact that the κ-dimer has only 0.5% disconnected saddles. For larger displacements (σ = 0.5 Å), the disconnected saddle rate drops from 34.5% in the regular dimer method to 0.2% with the κ-dimer method, so that the later becomes more efficient in terms of the average cost to find a connected saddle

### E. Au_{55} cuboctahedral cluster

Our final example examines the deformation mechanisms of a 55-atom Au cluster, initially in a metastable cuboctahedral shape. Au nanoparticles have attracted interest for their catalytic activity in various chemical reactions.^{18,19} The determination of the morphology of small nanoparticles is tremendously challenging, experimentally. For a theoretical understanding of catalytic activity, however, it is crucial to know the structure and the corresponding stability of model particles. Here we show that a gas-phase Au_{55} cuboctahedral cluster easily deforms to other structures by overcoming small barriers. In this example, the κ-dimer again improves the success ratio of saddle searches. We also show a remarkable case where four low-energy saddle points connect the same reactant and product states via different reaction pathways.

The interatomic interaction between Au atoms is described by the quantum Sutton-Chen potential with a cutoff radius of 15 A to include all the atoms in the cluster.^{20} For the κ-dimer calculations, the forms of γ_{1} and γ_{2} are taken as in Fig. 4. To accelerate convergence near saddle points, when Δ*V* is less than 0.1 eV/Å, γ_{1} and γ_{2} are set to unity to recover the convergence properties of the regular dimer method. The six translational and rotational modes are projected out for both the dimer and κ-dimer calculations. 1000 searches are conducted from initial configurations generated by displacing all 55 atoms by a Gaussian distribution with a standard deviation of 0.3 Å.

The relative performance of the two methods is shown in Table II, with barrier distributions in Fig. 8. The results are similar to the Pt heptamer example: the κ-dimer decreases the ratio of disconnected saddles to near-zero and increases the probability of finding low energy saddles. Different from the Pt heptamer, the Au_{55} cluster is more complicated in terms of the number of escape pathways from the reactant, due to its high ratio of under coordinated surface atoms. The lowest barrier process, which involves all 55 atoms (162 degrees of freedom), is the Mackay cuboctahedron to icosahedron transformation.^{21} This concerted mechanism, shown in Fig. 9(a), has a remarkably low barrier in Au (ca. 0.3 eV). The second lowest barrier, shown in the first row of Fig. 9(b), involves more than eight atoms. Interestingly, there exist another three reaction pathways connecting the same reactant and product state, with slightly higher barriers, as shown in the following three rows of Fig. 9(b). This high dimensional system with complicated reaction mechanisms demonstrates the robustness of the κ-dimer method for finding connected saddles.

Method . | σ (Å) . | Success . | $\langle \rm {FC} \rangle _1$ $\u27e8 FC \u27e91$
. | $\langle \rm {FC} \rangle _2$ $\u27e8 FC \u27e92$
. |
---|---|---|---|---|

dimer | 0.3 | 71.6% | 5085 | 20 115 |

κ-dimer | 0.3 | 97.6% | 4965 | 7276 |

Method . | σ (Å) . | Success . | $\langle \rm {FC} \rangle _1$ $\u27e8 FC \u27e91$
. | $\langle \rm {FC} \rangle _2$ $\u27e8 FC \u27e92$
. |
---|---|---|---|---|

dimer | 0.3 | 71.6% | 5085 | 20 115 |

κ-dimer | 0.3 | 97.6% | 4965 | 7276 |

## IV. CONCLUSION

The original dimer method is demonstrated to converge to saddles that are disconnected from the initial state where the search is initiated. This problem is shown to be intrinsic to the method and related to how the dimer method follows local ridges across the true ridge bordering the initial state and into adjacent states. The concept of the isopotential curvature κ is introduced into the method because the κ = 0 surface tends to cross the true ridge only at saddle points. Restraining the dimer method to the κ = 0 surface is the basis of the κ-dimer method, which is shown to nearly eliminate the disconnected-saddle problem, and in cases where saddle searches are started away from the local minimum, significantly increases the efficiency of finding connected saddles.

## ACKNOWLEDGMENTS

This work was supported by the National Science Foundation under Grant No. CHE-1152342 and the Welch Foundation under Grant No. F-1841. We gratefully acknowledge the Institute for Pure and Applied Mathematics at the University of California in Los Angeles for hosting our participation in the long program “Materials for a Sustainable Energy Future” from 9 September 2013 to 13 December 2013. Calculations were done at the Texas Advanced Computing Center. We thank Keith Promislow for helpful discussions.