next up previous contents
Next: A Revised Connection Algorithm Up: Results Previous: Doubly Nudged Elastic Bands   Contents

Comparison of the DNEB/L-BFGS and DNEB/SQVV Methods for Permutational Isomerisations of LJ$ _7$

It is sometimes hard to make a direct comparison of different double-ended methods for a particular rearrangement because the calculations may converge to different paths. Another problem concerns the choice of a consistent termination criterion: the RMS force usually converges to some finite system-dependent value, which in turn may depend on the number of images and other parameters. A low-energy chain of NEB images does not necessarily mean that a good pathway has been obtained, since it may arise because more images are associated with regions around local minima, rather than the higher energy transition state regions. Here we present the results of DNEB/L-BFGS, DNEB/SQVV and, where possible, NEB/SQVV calculations for all the distinct permutational rearrangements of the global minimum for the LJ$ _7$ cluster (see Figure 2.7 for the endpoints and nomenclature).

Figure: Structures of the most stable isomers for (a) LJ$ _7$, (b) LJ$ _{38}$ and (c) LJ$ _{75}$ clusters, which were used as endpoints in the NEB calculations. The first endpoint was the global minimum in each case. For LJ$ _{38}$ and LJ$ _{75}$ the second endpoint was chosen to be second-lowest minimum shown on the right in parts (b) and (c), respectively, while a permutational isomer of the global minimum was used as the second endpoint in all the LJ$ _7$ calculations. The notation 1-2 denotes an LJ$ _7$ rearrangement where the second endpoint is structure (a) with atoms 1 and 2 swapped. The structures and numbering employed for LJ$ _{38}$ and LJ$ _{75}$ are defined at http://www-wales.ch.cam.ac.uk/$ \sim $sat39/DNEBtests/. Picture of LJ$ _7$ was generated using XMakemol program written by Dr. Matthew Hodges [139]. Structures of LJ$ _{38}$ and LJ$ _{75}$ clusters were visualised using a Mathematica [140] notebook for making and manipulating triangulated polyhedra written by Dr. David Wales [141].
\begin{psfrags} \psfrag{11} [bc][bc]{$1$} \psfrag{22} [bc][bc]{$2$} \psfra... ...\includegraphics[width=.38\textheight]{dneb/ep_pics/ep_pics.eps}} \end{psfrags}

It is possible to draw a firm conclusion as to how well the NEB represents the pathway when the corresponding stationary points and steepest-descent paths are already known. We therefore base our criterion for the effectiveness of an NEB calculation on whether we obtain good estimates of all the transition states. By considering several systems of increasing complexity we hope to obtain comparisons that are not specific to a particular pathway.

Connections between two minima are defined by calculating an approximation to the two steepest-descent paths that lead downhill from each transition state, and two transition states are considered connected if they are linked to the same minimum via a steepest-descent path. We will say that minima are `connected' if there exists a path consisting of one or more transition states and intermediate minima linking them. Permutational isomers of the same minimum are distinguished in these calculations. We refer to the chain of images produced by the NEB calculation as `connected' if going downhill from each transition state using steepest-descent minimisation yields a set of minima that contains the endpoints linked together.

For NEB/SQVV calculations we used the NEB formulation defined in Reference henkelmanj00. DNEB is different from the above method because it includes an additional component in the NEB gradient, as described by Equation 2.39 and Equation 2.40. In addition, for the following DNEB calculations we did not remove overall rotation and translation (ORT), because we believe it is unnecessary when our gradient modification is used. To converge transition state candidates tightly we employed EF optimisation, limiting the maximum number of EF iterations to five with an RMS force convergence tolerance of $ 10^{-5}$. (Standard reduced units for the Lennard-Jones potential are used throughout this work.) Initial guesses for all the following calculations were obtained by linear interpolation between the endpoints. To prevent `atom-crashing' from causing overflow in the initial guess we simply perturbed such images slightly using random atomic displacements of order $ 10^{-2}$ reduced units.

In each case we first minimised the Euclidean distance between the endpoints with respect to overall rotation and translation using the method described in Reference rhee00.[1] Structure alignment methods that are based on finding a pseudorotation matrix that satisfies Eckart axis conditions might equally well be used for this purpose [147,145,143,144,146]. SQVV minimisation was performed with a time step of $ 0.01$ and a maximum step size per degree of freedom of $ 0.01$. This limit on the step size prevents the band from becoming `discontinuous' initially and plays an important role only during the first $ 100$ or so iterations. The limit was necessary because for the cases when the endpoints are permutational isomers linear interpolation usually yields bands with large gradients, and it is better to refrain from taking excessive steps at this stage. We did not try to select low energy initial guesses for each rearrangement individually, since one of our primary concerns was to automate this process. For the same reason, all the L-BFGS optimisations were started from guesses preoptimised using SQVV until the RMS force dropped below $ 2.0$.

Table 2.1 shows the minimum number of images and gradient calls required to produce a connected pathway using the DNEB/L-BFGS and DNEB/SQVV methods. These calculations were run assuming no prior knowledge of the path. Normally there is no initial information available on the integrated path length or the number of intermediate minima between the endpoints, and it takes some experimentation to select an appropriate number of images. Our strategy is therefore to gradually increase the number of images to make the problem as computationally inexpensive as possible. Hence we increment the number of images and maximum number of NEB iterations in each calculation until a connected path is produced, in the sense defined above. The permitted image range was $ 2 \leqslant N_i \leqslant 20$ and the maximum number of NEB iterations ranged from $ 1 \leqslant l \leqslant 3,000N_i$. We were unable to obtain connected pathways for any of the four LJ$ _7$ rearrangements using the NEB/SQVV approach.


Table: The minimal number of images and total number of gradient calls (in parentheses) are shown for degenerate rearrangements of LJ$ _7$. The image range was $ 2 \leqslant N_i \leqslant 20$ and the iteration range was $ 1 < \ell \leqslant 3,000N_i$. Each SQVV calculation was started from the guess produced using linear interpolation, while guesses for L-BFGS runs were preoptimised using DNEB/SQVV until the RMS force dropped below 2.0. Every iteration the images that satisfy $ V_i > V_{i\pm 1}$ were optimised further using eigenvector-following. The transition state candidates that converged to a true transition state within five iterations were used to generate the connected minima using energy minimisation. If this procedure yielded a connected pathway the calculation was terminated and the rest of the parameter range was not explored. Otherwise, the number of images was incremented and the procedure repeated. The number of gradient calls is a product of the number of images and the total number of iterations. For the L-BFGS calculations the number of iterations includes the SQVV preoptimisation steps (100 on average) and the actual number of L-BFGS steps. Dashes signify cases where we were unable to obtain a connected pathway.
Method 1-2 2-3 3-4 4-5
DNEB/L-BFGS 5(1720) 18(30276) 11(2486) 18(8010)
DNEB/SQVV 16(21648) - 10(14310) -

Table 2.2 presents the results of analogous calculations where we keep the number of images fixed to 50. Unlike the performance comparison where the number of images is kept to a minimum (Table 2.1), these results should provide insight into the performance of the DNEB approach when there are sufficient images to resolve all the transition states. All the optimisations for a particular rearrangement converged to the same or an enantiomeric pathway unless stated otherwise. The energy profiles that correspond to these rearrangements are shown in Figure 2.8.

Figure: The potential energy, $ V$, as a function of the integrated path length, $ s$,for four degenerate rearrangements of LJ$ _7$. These profiles were constructed using energy minimisation to characterise the paths connected to transition states obtained by EF refinement of candidate structures obtained from DNEB calculations [19].
\begin{psfrags} \psfrag{2} [bc][bc]{$2$} \psfrag{4} [bc][bc]{$4$} \psfrag{... ...terline{\includegraphics[width=.8\textheight]{dneb/profiles.eps}} \end{psfrags}


Table: The minimal number of iterations needed to produce connected pathways for four degenerate rearrangements of LJ$ _7$ using a 50-image NEB. The strategy of this calculation is identical to the one described in the caption to Table 2.1, except that the number of images was fixed.
Method 1-2 2-3 3-4 4-5
DNEB/L-BFGS 131$ ^a$ 493 171 326
DNEB/SQVV 1130 15178 2777 23405$ ^b$
NEB/SQVV 11088$ ^b$ - 30627 -

$ ^a$ The number of iterations is the sum of the SQVV preoptimisation steps (100 on average) and the actual number of iterations needed by L-BFGS minimiser. $ ^b$ This value is not directly comparable since DNEB converged to a different path that contains more intermediate minima. Dashes signify cases where we were unable to obtain a connected pathway.


From Table 2.1 and Table 2.2 we conclude that in all cases the DNEB/L-BFGS approach is more than an order of magnitude faster than DNEB/SQVV. It is also noteworthy that the DNEB/SQVV approach is faster than NEB/SQVV because overall rotation and translation are not removed. Allowing the images to rotate or translate freely can lead to numerical problems, namely a vanishing norm for the tangent vector, when the image density is very large or the spring force constant is too small. However, when overall rotation and translation are not allowed there is less scope for improving a bad initial guess, because the images are more constrained. This constraint usually means that more images are needed or a better initial guess is required. Our experience is that such constraints usually slow down convergence, depending on which degrees of freedom are frozen: if these are active degrees of freedom (see above) the whole cluster must move instead, which is usually a slow, concerted multi-step process.

next up previous contents
Next: A Revised Connection Algorithm Up: Results Previous: Doubly Nudged Elastic Bands   Contents
Semen A Trygubenko 2006-04-10