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.

The dimer method6 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,

$$\dot{x} = -\nabla V(x),$$
$ẋ=−∇V(x),$
(1)

where

$x \in \mathbb {R}^{n}$
$x∈Rn$ is the position, 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,

\begin{align}\dot{x} &= F^{\perp } - F^{\parallel }, \nonumber\\F^{\parallel } &= -\hat{\tau }^{\intercal } \nabla V(x) \hat{\tau }, \\F^{\perp } &= -\nabla V(x) - F^{\parallel }, \nonumber\end{align}
$ẋ=F⊥−F‖,F‖=−τ̂⊺∇V(x)τ̂,F⊥=−∇V(x)−F‖,$
(2)

where

$\hat{\tau }$
$τ̂$ is the unit eigenvector of the Hessian matrix H(V) associated with the smallest eigenvalue.6
$\hat{\tau }$
$τ̂$
is also called the min-mode direction, since it is aligned along the minimum mode of the Hessian. The name of the dimer method comes from the way
$\hat{\tau }$
$τ̂$
is calculated. Two nearby images (a dimer of configurations) are rotated to find the lowest mode. This finite-difference approximation avoids explicit construction and diagonalization of the Hessian matrix. The center of the dimer 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
$\hat{\tau }$
$τ̂$
, so that Eq. (2) converges. All other types of critical points of Eq. (1) are unstable under Eq. (2).

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

$\hat{\tau }$
$τ̂$⁠) is stable under Eq. (2), the dimer is unable to cross a local ridge. The heart of the disconnected-saddle problem is that while the local ridge follows the true ridge in the harmonic region of a saddle point, the two tend to diverge away from the saddle point. Even more problematic is that the local ridge is not nessecarily a closed surface; holes in the local ridge allow the dimer to escape from the initial state. Basins of attraction associated with disconnected saddles can penetrate the initial state through these holes and attract the dimer across the true ridge.

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, A1. 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

$\hat{\tau }$
$τ̂$ and the force −∇V(x) directions.

FIG. 1.

A two dimensional potential shown in the black contour lines has four minima M1-4 and four saddle points S1-4. The colors indicate the dot product between the minimum mode and force directions,

$|\hat{\tau }^{\intercal } \hat{\nabla V}|$
$|τ̂⊺∇V̂|$⁠. In the shadowed areas, the lowest eigenvalue of the Hessian matrix is positive, and in the magenta circles, both eigenvalues of the Hessian are negative.

FIG. 1.

A two dimensional potential shown in the black contour lines has four minima M1-4 and four saddle points S1-4. The colors indicate the dot product between the minimum mode and force directions,

$|\hat{\tau }^{\intercal } \hat{\nabla V}|$
$|τ̂⊺∇V̂|$⁠. In the shadowed areas, the lowest eigenvalue of the Hessian matrix is positive, and in the magenta circles, both eigenvalues of the Hessian are negative.

Close modal

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

$\hat{\tau }$
$τ̂$ abruptly switches from one eigenvector to the other.

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.

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,

$$\dot{x} = F^{\perp } - \gamma F^{\parallel }; (0 < \gamma \le 1).$$
$ẋ=F⊥−γF‖;(0<γ≤1).$
(3)

We note that the activation-relaxation technique (ART-neuveau)8 and hybrid eigenvector following7 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.

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.

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(x1, x2). At any point, there is a tangent direction of the contour,

$\hat{c}$
$ĉ$⁠, which can be used to define an angle, θ, between the contour and a fixed reference.
$\hat{c}$
$ĉ$
is a unit vector perpendicular to ∇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,

$$\kappa = \frac{d\theta }{ds} = -\frac{\nabla _{\hat{c}} (\nabla _{\hat{c}} V)}{\Vert \nabla V \Vert } = -\frac{\hat{c}^{\intercal } \mathbf {H} \hat{c}}{\Vert \nabla V \Vert }.$$
$κ=dθds=−∇ĉ(∇ĉV)‖∇V‖=−ĉ⊺Hĉ‖∇V‖.$
(4)

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

$V(x_1, x_2) = a x_1^2 - b x_2^2$
$V(x1,x2)=ax12−bx22$ with a > 0 and b > 0 near the saddle point (0, 0), then

$$\kappa = \frac{b x_2^2 - a x_1^2}{\sqrt{(x_1/a)^2 + (x_2/b)^2}} \Rightarrow \lim _{(x_1, x_2) \rightarrow (0, 0)} \kappa = 0.$$
$κ=bx22−ax12(x1/a)2+(x2/b)2⇒lim(x1,x2)→(0,0)κ=0.$
(5)
FIG. 2.

The attraction regions of the dimer method using (a) Eq. (2) and (b) Eq. (3) with γ = 0.2. The colored regions indicate the following: deep blue: both eigenvalues of the Hessian are positive; green: the attraction region of S1; light blue: S2; orange: S3; yellow: S4. The red lines are the local ridges and the white lines are the true ridges.

FIG. 2.

The attraction regions of the dimer method using (a) Eq. (2) and (b) Eq. (3) with γ = 0.2. The colored regions indicate the following: deep blue: both eigenvalues of the Hessian are positive; green: the attraction region of S1; light blue: S2; orange: S3; yellow: S4. The red lines are the local ridges and the white lines are the true ridges.

Close modal
FIG. 3.

Regions of positive (light blue) and negative (dark blue) isopotential curvature κ are separated by the κ = 0 restraining boundary for the κ-dimer method. The evolutions of the dimer (orange) and κ-dimer (yellow) methods are compared from the same initial position, x0.

FIG. 3.

Regions of positive (light blue) and negative (dark blue) isopotential curvature κ are separated by the κ = 0 restraining boundary for the κ-dimer method. The evolutions of the dimer (orange) and κ-dimer (yellow) methods are compared from the same initial position, x0.

Close modal

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

$\hat{n}$
$n̂$⁠:
$\nabla _{\hat{n}} V(x_1, x_2) |_{ (x_1, x_2) \in S } = 0$
$∇n̂V(x1,x2)|(x1,x2)∈S=0$
. The normal
$\hat{n}$
$n̂$
is also parallel to the contour because the contour is perpendicular to the force. The sign of κ is then determined by the sign of
$\hat{n}^{\intercal } \mathbf {H} \hat{n} \equiv \nu$
$n̂⊺Hn̂≡ν$
. When ν > 0, κ < 0 at the true ridge, and constraining the dimer in the κ < 0 region does not necessarily avoid the true ridge. Since ν is the second derivative along the normal direction of the true ridge, having both ν > 0 and a normal gradient of zero means V(x1, x2) 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

$\mathbb {R}^{n}$
$Rn$⁠, the isopotential contours form an (n − 1)-dimensional hyperplane. The direction perpendicular to the force is therefore not unique, and our choice of κ is the minimum directional curvature,

$$\kappa = -\min _{\hat{c} \perp \nabla V} \frac{ \hat{c}^{\intercal } \mathbf {H} \hat{c} }{\Vert \nabla V \Vert },$$
$κ=−minĉ⊥∇Vĉ⊺Hĉ‖∇V‖,$
(6)

where

$\hat{c}$
$ĉ$ is a unit vector and
$\hat{c} \perp \nabla V$
$ĉ⊥∇V$
constrains
$\hat{c}$
$ĉ$
on the tangent plane of the isopotential hyperplane. Denoting λ1 ⩽ λ2 ⩽ ⋅⋅⋅ as the eigenvalues of H and ν1 as
$\min _{\hat{c} \perp \nabla V} ( \hat{c}^{\intercal } \mathbf {H} \hat{c} )$
$minĉ⊥∇V(ĉ⊺Hĉ)$
, we still have λ1 ⩽ ν1 ⩽ λ2. On the isopotential hyperplane, there is always a vector
$\hat{c}_x$
$ĉx$
on the plane spanned by the eigenvectors associated with λ1 and λ2. Thus,
$\hat{c}_x^{\intercal } \mathbf {H} \hat{c}_x \le \lambda _2$
$ĉx⊺Hĉx≤λ2$
. Combined with
$\nu _1 \le \hat{c}_x^{\intercal } \mathbf {H} \hat{c}_x$
$ν1≤ĉx⊺Hĉx$
, we can prove that λ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.

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

\begin{align}\dot{x} &= \gamma _2 F^{\perp } - \gamma _1 F^{\parallel }, \nonumber\\\gamma _2 &= 1- \frac{1}{1 + e^{\beta \kappa }} , \\\gamma _1 &= \frac{2}{1 + e^{ \beta \kappa }} - 1, \nonumber\end{align}
$ẋ=γ2F⊥−γ1F‖,γ2=1−11+eβκ,γ1=21+eβκ−1,$
(7)

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.

1. When κ is sufficiently negative,

$\dot{x} = -F^{\parallel }$
$ẋ=−F‖$ 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.

2. Near κ = 0 where γ1 = γ2 > 0 we recover the regular dimer method of Eq. (2).

3. When κ is sufficiently positive,

$\dot{x} = F^{\perp } + F^{\parallel } = -\nabla V$
$ẋ=F⊥+F‖=−∇V$ and the dimer moves down the potential away from the region containing the true ridge.

FIG. 4.

Values of γ1 and γ2 as a function of κ, with β = 5.

FIG. 4.

Values of γ1 and γ2 as a function of κ, with β = 5.

Close modal

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.

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.

FIG. 5.

Basins of attraction of the κ-dimer method, as described by Eq. (7).

FIG. 5.

Basins of attraction of the κ-dimer method, as described by Eq. (7).

Close modal

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.

FIG. 6.

A counter example, where the disconnected problem still exists with the κ = 0 restraint applied. (a) Regions divided by the κ = 0 surface; the white lines are true ridges and the white points are saddles. (b) The basins of attraction for the κ-dimer method. The deep blue regions have two positive Hessian eigenvalues and the other colors are the basins of attraction of the contained saddle. The red circle is the boundary of the attractor to S1 of the regular dimer method.

FIG. 6.

A counter example, where the disconnected problem still exists with the κ = 0 restraint applied. (a) Regions divided by the κ = 0 surface; the white lines are true ridges and the white points are saddles. (b) The basins of attraction for the κ-dimer method. The deep blue regions have two positive Hessian eigenvalues and the other colors are the basins of attraction of the contained saddle. The red circle is the boundary of the attractor to S1 of the regular dimer method.

Close modal

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.

FIG. 7.

Comparison of the dimer and κ-dimer methods in terms of the distribution of barriers found and number of force calls required to find them. The inset figure in (b) shows the geometry of the reactant minimum containing a Pt7 island on Pt(111). The magnitude of the initial Gaussian displacements used to initiate the searches is (a) relatively conservative, 0.3 Å; (b) relatively aggressive, 0.5 Å.

FIG. 7.

Comparison of the dimer and κ-dimer methods in terms of the distribution of barriers found and number of force calls required to find them. The inset figure in (b) shows the geometry of the reactant minimum containing a Pt7 island on Pt(111). The magnitude of the initial Gaussian displacements used to initiate the searches is (a) relatively conservative, 0.3 Å; (b) relatively aggressive, 0.5 Å.

Close modal

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.

Table I.

Success ratio for finding connected saddles for the dimer and κ-dimer methods for two displacement amplitudes σ. The average number of force calls are measured in terms of all the successful searches

$\langle \rm {FC} \rangle _1$
$⟨ FC ⟩1$ and successful searches that found saddles below 1.4 eV
$\langle \rm {FC} \rangle _2$
$⟨ FC ⟩2$
.

Methodσ (Å)Success
$\langle \rm {FC} \rangle _1$
$⟨ FC ⟩1$
$\langle \rm {FC} \rangle _2$
$⟨ FC ⟩2$
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$
$⟨ FC ⟩1$
$\langle \rm {FC} \rangle _2$
$⟨ FC ⟩2$
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

$\langle \rm {FC} \rangle _1$
$⟨ FC ⟩1$⁠. In terms of the cost to find a saddle below 1.4 eV
$\langle \rm {FC} \rangle _2$
$⟨ FC ⟩2$
, the improvement is even more dramatic. These large-displacement results are particularly important because they systematically increase the probability of finding relative higher-energy saddles, as shown in Fig. 7. Two additional aspects of the κ-dimer method can be seen from this figure. First, there is a better chance of finding all saddles in the energy window of interest, such as the saddles around 1.2 eV. Second, while the κ-dimer costs twice as many force calls per rotation, the cost per saddle search is typically very similar to the regular dimer, indicating that the κ-dimer is taking a more direct path to the saddles as a result of the κ restraint.

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 Au55 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 Au55 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.

FIG. 9.

Deformation pathways of the Au55 cuboctahedron. (a) Mackay transformation to the icosahedron. The red lines indicate the front edges. (b) Four possible pathways to deform the triangle connected by red lines to the same product structure. The black arrows indicate the movement of atoms from the reactant to the saddle structure. The red dots emphasize the differences between the saddle points in the four mechanisms.

FIG. 9.

Deformation pathways of the Au55 cuboctahedron. (a) Mackay transformation to the icosahedron. The red lines indicate the front edges. (b) Four possible pathways to deform the triangle connected by red lines to the same product structure. The black arrows indicate the movement of atoms from the reactant to the saddle structure. The red dots emphasize the differences between the saddle points in the four mechanisms.

Close modal
Table II.

Success ratio for finding connected saddles for the dimer and κ-dimer methods. The average number of force calls are measured in terms of all the successful searches

$\langle \rm {FC} \rangle _1$
$⟨ FC ⟩1$ and successful searches that found saddles below 0.5 eV
$\langle \rm {FC} \rangle _2$
$⟨ FC ⟩2$
.

Methodσ (Å)Success
$\langle \rm {FC} \rangle _1$
$⟨ FC ⟩1$
$\langle \rm {FC} \rangle _2$
$⟨ FC ⟩2$
dimer 0.3 71.6% 5085 20 115
κ-dimer 0.3 97.6% 4965 7276
Methodσ (Å)Success
$\langle \rm {FC} \rangle _1$
$⟨ FC ⟩1$
$\langle \rm {FC} \rangle _2$
$⟨ FC ⟩2$
dimer 0.3 71.6% 5085 20 115
κ-dimer 0.3 97.6% 4965 7276
FIG. 8.

Comparison of the dimer and κ-dimer methods for a Au 55 atom cuboctahedral cluster. The magnitude of the initial Gaussian displacements used to initiate the searches is 0.3 Å.

FIG. 8.

Comparison of the dimer and κ-dimer methods for a Au 55 atom cuboctahedral cluster. The magnitude of the initial Gaussian displacements used to initiate the searches is 0.3 Å.

Close modal

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.

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.

1.
H.
Eyring
,
J. Chem. Phys.
3
,
107
(
1935
).
2.
G. H.
Vineyard
,
J. Phys. Chem. Solids
3
,
121
(
1957
).
3.
G.
Henkelman
and
H.
Jónsson
,
J. Chem. Phys.
115
,
9657
(
2001
).
4.
L.
Xu
and
G.
Henkelman
,
J. Chem. Phys.
129
,
114104
(
2008
).
5.
F. H.
Stillinger
and
T. A.
Weber
,
J. Phys. Chem.
87
,
2833
(
1983
).
6.
G.
Henkelman
and
H.
Jónsson
,
J. Chem. Phys.
111
,
7010
(
1999
).
7.
L. J.
Munro
and
D. J.
Wales
,
Phys. Rev. B
59
,
3969
(
1999
).
8.
R.
Malek
and
N.
Mousseau
,
Phys. Rev. E
62
,
7723
(
2000
).
9.
H.
Jónsson
,
G.
Mills
, and
K. W.
Jacobsen
, in
Classical and Quantum Dynamics in Condensed Phase Simulations
, edited by
B. J.
Berne
,
G.
Ciccotti
, and
D. F.
Coker
(
World Scientific
,
Singapore
,
1998
), pp.
385
404
.
10.
D.
Asenjo
,
J. D.
Stevenson
,
D. J.
Wales
, and
D.
Frenkel
,
J. Phys. Chem. B
117
,
12717
(
2013
).
11.
J. J.
Koenderink
and
A. J.
van Doorn
,
Opt. Acta
27
,
981
(
1980
).
12.
J. J.
Koenderink
and
W.
Richards
,
J. Opt. Soc. Am. A
5
,
1136
(
1988
).
13.
L. J.
van Vliet
,
Grey-scale Measurements in Multi-dimensional Digitized Images
(
Delft University Press
,
1993
).
14.
D.
Sheppard
and
G.
Henkelman
,
J. Comput. Chem.
32
,
1769
(
2011
).
15.
G.
Henkelman
,
G.
Jóhannesson
, and
H.
Jónsson
, in
Progress on Theoretical Chemistry and Physics
, edited by
S.
Schwartz
(
,
New York
,
2000
), pp.
269
299
.
16.
S. T.
Chill
,
M.
Welborn
,
R.
Terrell
,
L.
Zhang
,
J.-C.
Berthet
,
A.
Pedersen
,
H.
Jónsson
, and
G.
Henkelman
,
Model. Simul. Mater. Sci. Eng.
22
,
055002
(
2014
).
17.
A.
Pedersen
,
S. F.
Hafstein
, and
H.
Jónsson
,
SIAM J. Sci. Comput.
33
,
633
(
2011
).
18.
Y.
Mikami
,
A.
Dhakshinamoorthy
,
M.
Alvaro
, and
H.
García
,
Catal. Sci. Tech.
3
,
58
(
2013
).
19.
G.
Schmid
,
Chem. Soc. Rev.
37
,
1909
(
2008
).
20.
T.
Çağin
,
Y.
Kimura
,
Y.
Qi
,
H.
Li
,
H.
Ikeda
,
W. L.
Johnsonb
, and
W. A.
Goddard
,
MRS Proceedings
(
Cambridge University Press
,
1998
), Vol.
554
, p.
43
.
21.
A.
Mackay
,
Acta Crystallogr.
15
,
916
(
1962
).