Sampling algorithms
The most primitive way to determine transition states would be to systematically sample configuration space point by point. However, by definition transition states are metastable points in configuration space. Therefore, they cannot be obtained by a minimization of the potential energy, and more advanced optimization techniques such as Newton’s algorithm have to be used. While the latter may stop on any stationary state, we are only interested in those with exactly one unstable direction so that many unnecessary computations have to be done.
A better approach is to take advantage of the vanishing of the forces at the transition state and define an auxiliary function (Stillinger and Weber, 1984b)
Ψ = (∇Φ)2=kFk2. (A.11)
Hence, any point in the potential energy landscape with vanishing force is a minimum of Ψ, thus replacing the task of finding the transition state with the much easier problem of finding local minima. However, there are two serious drawbacks: First, the minimization will end on any minimum of Ψ, which may either be higher order equilibria or even inflection points of the potential, i.e. local minima in kFk. Typically, transition states represent only a small subset thereof, and many excess computations have to be done. To identify the desired states, in a second step the spectrum of the minimum configuration has to be inspected. Another issue is that during the minimization second derivatives have to be evaluated, which may be computationally expensive.
154
A.4 Algorithms to locate transition states
An approach inspired by the separating nature of the transition state was proposed by Berry et al. (1988). The idea is, in a first step to generate a path connecting the minima, e.g. by MD-simulations. Thereafter, initial conditions are sampled randomly in the vicinity of the configuration with maximum potential energy, which presumably should be close to the transition state. For each initial condition, a steepest descent quench is performed.
Close to a transition state, forces become very small and the minimization algorithm stays in its vicinity for several iterations. Thus, by stopping the algorithm at moderate forces kFk, quenches can be interrupted close to the transition state before they relax towards a minimum. By sampling several initial points, chances to hit the transition state during the minimization are increased. In contrast to the previous method, the algorithm will not stop at local force minima. However, the co-dimension of the point still has to be determined since the algorithm will terminate at any force equilibrium. Furthermore, due to its unsystematic approach, the algorithm cannot guarantee to find the transition state corresponding to the transition in question.
Ridge algorithm
An algorithm which is conceptually close to the edge tracking algorithm was proposed by Ionova and Carter (1993). Here, we start with two sets of initial conditions, one on either side of the transition state. Hence, the potential energy has a maximum on the line connecting the two configurations. Instead of computing the complete trajectories, each configuration is shifted in the downhill direction, such that the energy maximum on the line connecting the two resulting configurations is lower than the previous maximum. In the next step, two new initial conditions for the next iteration are determined, each one situated on the connecting line and in a distance pfrom the maximum. Since they are on different sides of a maximum, they evolve towards two different minima. Thereafter, the algorithm repeats itself, gradually approaching the transition state, which by definition has the lowest maximum energy. The algorithm does not check the minima corresponding to the intermediate initial conditions. Thus, when investigating a transition between two minima consisting of a chain of several transitions and intermediate minima, the algorithm will end up on any one of the transition states. Therefore, the algorithm is not suitable for a systematic search for transition paths.
Step and slide algorithm
Similar to the ridge algorithm, the step and slide method (Miron and Fichthorn, 2001) considers two replicas of the system as well. The algorithm incorporates two steps. It starts with a configuration close to each minimum, set to a potential energy Φtrial. In the slide phase, both configurations are displaced along their equipotential lines such that the distance between the ensembles is minimized. Thereafter, in the step phase, both configurations are shifted towards a higher potential energy level, and the algorithm repeats itself. In case that the new level was set too high and both ensembles collide, an intermediate value of Φtrial is tested. Unlike other algorithms, it thereby brackets the potential energy of the transition state both from below and above. However, it is unsuitable to determine a chain of several subsequent transitions connecting two minima:
Appendix
The potential energy levels will be increased until only the energetically highest transition state is obtained, ignoring all intermediate transitions.
Nudged elastic band method
In the Nudged Elastic Band (NEB) method (J´onsson et al., 1998), not only one or two but many copies of the system are considered, spaced between e.g. two minimum con-figurations. Each for itself will relax towards the closest minimum, giving no additional information. To prevent this, the copies are connected with each other, forming a band of configurations in the potential energy landscape. Along the connecting lines, spring forces are invoked so that neighboring copies will stay close together. The forces constituted by the potential are then projected on the direction perpendicular to the band, whereas the spring forces are projected on the tangent. In consequence, the copies will form an elastic string, connecting minima via transition states and approximating the minimum energy path. A slight modification is made in theClimbing Image Nudged Elastic Band (CI-NEB) algorithm (Henkelman et al., 2000). In the last step, in order to pin the transition state the copy with the largest potential energy is allowed to move uphill. To this end, forces along the line connecting it to its neighbors are inverted. Another modification is the Doubly Nudged Elastic Band (DNEB) method, which allows for spring forces perpendicular to the band, as long as they are perpendicular to the potential forces as well (Trygubenko and Wales, 2004). It speeds up convergence, but also leads to deviations from the minimum energy path as it cuts corners in the potential energy landscape. A broad overview over the different realizations of the NEB-algorithm can be found in Sheppard et al. (2008).
An algorithm which is conceptually closely related to the NEB-method is the string method (E et al., 2002). The difference is that the spacing between copies of the system is set explicitly at each iteration, and not determined implicitly by spring forces. More re-cently, this method has been improved by the application of splines to interpolate between images (E et al., 2007).
All those methods are very efficient. Since the system has to be duplicated several times, memory demand increases likewise, and the implementation is less straightforward as for example of edge tracking. A downside of the NEB algorithms is that with only few images taken into account, the methods may lead to deviations from the optimum path or, more severely, fail to converge.
Activation-relaxation technique & eigenvector following
Another variety of algorithms use a bottom-up approach: They start at a minimum and use different types of hill climbing procedures to find the transition state. The second minimum is obtained by a subsequent minimization.
The activation-relaxation technique (ART) (Barkema and Mousseau, 1996; Mousseau and Barkema, 1998; Barkema and Mousseau, 2001) implements the hill climbing by in-verting the force component parallel to the displacement from the starting minimum
∆x=x−x0,
g=f −(1 +α)hF,e∆xie∆x. (A.12)
156
A.5 Analytical solution for the singular force in an infinite 2d-solid
Here, g is the adjusted force, α is a parameter which determines the strength of the inverted force, and e∆x is the unit vector of ∆x. Thus, the system will be minimized in all but one direction, at the same time driving it away from the minimum towards a saddle point. At the transition state, the overlap hF,e∆xi will reverse its sign, thus providing a stopping criterion.
It is crucial to define ∆x properly. If all particles are considered in ∆x, the method is likely to fail, as the displacement quickly spreads across all particles such that the few particles to exchange their positions fall back upon their initial coordinates. However, when restricting ∆x to only a few particles, for example to those which should exchange their positions, this method gives a good approximation of the transition state, i.e. the force is non-zero, and the associated energy and unstable eigenvalue are slightly off.
Obviously, this method is only able to detect a single transition at a time. Thereafter, the next minimum has to be determined, and the process can be restarted. However, in a chain of several transitions, it is a non-trivial task to identify the next search direction.
Moreover, it becomes unreliable the more particles are considered in ∆x.
A more involved method is the eigenvector following (Wales, 1994; Munro and Wales, 1999). As the name suggests, the direction of the least stable mode is evaluated, and the system is then displaced uphill along this direction and optimized in the perpendicular space. Just as ART, it is a single ended method, less suitable for finding prescribed connections between minima.
Dimer method
A method which can be seen as closely related to both the band methods and eigenvector following is the dimer-method (Henkelman and J´onsson, 1999). It uses two copies of the system, with the distance between them kept fixed. This dimer is then rotated such that its orientation is parallel to the direction of lowest curvature. In the picture of an energy landscape, it is oriented along a valley. This is equivalent to a minimization of the combined potential energy of the constituting particles, while keeping the center of mass fixed.
Thereby, potentially expensive evaluations of second derivatives are avoided. Afterwards, the potential force is applied, where the force component parallel to the orientation of the dimer is inverted, subsequently pushing the dimer uphill towards a transition state.
This last step is conceptually the same in the eigenvector following, where the direction is directly derived from the spectrum. This method is especially interesting for systems where the evaluation of the spectrum is computationally expensive. The crucial point is an effective way to implement the rotation of the dimer.