A preconditioning scheme for Minimum Energy Path finding methods

Popular methods for identifying transition paths between energy minima, such as the nudged elastic band and string methods, typically do not incorporate potential energy curvature information, leading to slow relaxation to the minimum energy path for typical potential energy surfaces encountered in molecular simulation. We propose a preconditioning scheme which, combined with a new adaptive timestep selection algorithm, substantially reduces the computational cost of transition path finding algorithms. We demonstrate the improved performance of our approach in a range of examples including vacancy and dislocation migration modelled with both interatomic potentials and density functional theory.


I. INTRODUCTION
In computational chemistry, structural biology, materials science and engineering, the time taken for processes is often dominated by transitions between energy minima in a potential energy landscape. The computational evaluation of the Minimum Energy Path (MEP) of the transition is a familiar technique used to find the energy barrier ∆E of such a transition 1 . The objective is the evaluation of the transition rate to leading order which is given by ν ∼ ν 0 exp −∆E / kBT 2,3 , where the attempt rate ν 0 may be estimated using Eyring's heuristic derivation 2 , or approximated with Harmonic Transition State Theory 4 , k B is the Boltzmann constant and T is the temperature of the system. Knowing the transition rate enables the simulation of the transition on the mesoscale using, for example, the kinetic Monte Carlo method 5 .
We restrict our focus to 'double ended' cases where both energy minima are known. The most notable techniques in this case are the string method 6-8 and the Nudged Elastic Band (NEB) method 9,10 . Both methods find the MEP by iteratively relaxing a discretised path, of N images, until convergence to an approximate MEP is achieved. Typically, the path is evolved in the energy landscape via a steepest descent-like optimisation technique, which may converge slowly when the geometry of the potential is ill-conditioned, that is, the Hessian matrix of the potential along the path has a large condition number 11 . Such a situation arises, for example, in large computational domains or if bonds with significant stiffness variations are present. Preconditioning is commonly used in linear algebra and optimisation to effectively reduce the condition number and thus improve the rate of convergence of an iterative scheme 11 .
It has been shown for example in Refs. 12-14 how to construct and invert effective preconditioners for the potential energy landscape of materials and molecules at a cost comparable to the evaluation of an interatomic potential and much lower than the cost of evaluating a DFT model. When used correctly, preconditioning leads to a substantial reduction in the number of force calls and thus is expected to significantly improve computing times 12,15 .
In this paper we introduce a simple yet effective way to precondition the standard NEB and string methods to obtain efficient and robust algorithms for computing MEPs in ill-conditioned geometries. Our scheme is further enhanced by a novel adaptive step length selection method for optimising the MEP. We demonstrate the effectiveness of this combination on a range of material modelling examples.

II. THE NEB AND STRING METHODS
Let x ∈ R M , M ∈ N, be a state, or configuration, of the dynamical system in question. We denote by V (x) the potential energy of x and assume that V is twice differentiable and that it has at least two local minima, which we denote by x A and x B , separated by a single saddle point x S of Morse index 1 (to ensure that there is a unique direction of steepest descent at x S 16 ). An MEP of the transition from x A to x B is defined as the intrinsically parametrised path x * (s), s ∈ [0, 1], satisfying with end points at the local minima where x = dx ds . (We note that, strictly speaking ∇ ⊥ V depends on x as well as x but for the sake of simplicity of notation we will only write ∇ ⊥ V (x).) The NEB and string methods discretise a path x(s) by interpolating N discrete points {x n } N n=1 . In the present work we will employ cubic spline interpolation 17 , imposing the "not-a-knot" boundary condition, but the methods we discuss can be readily extended to other interpolation schemes as well.
To evolve the discrete path to equilibrium we introduce a pseudo-temporal coordinate τ and writeẋ = dx dτ . The evolution of x n (τ ) is then described by the system of ODEsẋ where η = 0 leads to the string method, while the NEB method introduces elastic interactions between adjacent arXiv:1810.02705v1 [physics.comp-ph] 5 Oct 2018 images along the path by adding the term The system (2) can be solved with any ordinary differential equation (ODE) numerical integrator. Most commonly, Euler's method 7 is used, which yields an update step of the form where η k n = η((x k n ) , (x k n ) ) and α k is the timestep at iteration k.
While for NEB the presence of the elastic interaction η enforces an approximate equidistribution of the nodes along the path, the string method reparametrises the path after each iteration to ensure that the images remain equidistant with respect to a suitable metric. In the continuous limit, as N → ∞ a converged discretised path tends to the correct MEP, independently of the choice of the reparametrisation metric 8 . We initially use the standard 2 -norm defined by x 2 = x · x, but we will introduce a different notion of distance later on.
To summarise, the updating relations are given by: where, for the string method only, there is an additional redistribution of the images after the update step. We follow precisely the approach described in Eq. 12 in Ref. 7, but for simplicity of presentation do not make this step explicit.
The updating steps Eq. (4) for the string and NEB methods as well as the subsequent analysis were defined in terms of total derivatives of the path variable x (i.e. in terms of x and x ), as they are motivated from the respective laws of classical dynamics. This information is available at each iteration at no extra cost as we use cubic spline interpolation to find an expression for x(s) 6,9 .

III. PRECONDITIONING
The NEB and string methods have slow convergence rates when they are subjected to ill-conditioned energy landscapes V . However, a suitable preconditioner P ∈ R M × R M that is cheap to compute can be used to reduce the condition number of the Hessian ∇∇V along the path. In steepest descent optimisation, preconditioning has related but distinct interpretations: (a) as an approximation of the hessian, P ≈ ∇∇V , in analogy to Newton's scheme or (b) as a coordinate transformation in the state space, x → P 1/2 x, that captures information of the local curvature of the potential landscape (mapping hyperellipsoids to balls) 11 .
We will now describe a preconditioning technique for NEB and string methods. The same preconditioners used in geometry optimisation of interatomic potentials 12 are expected to be valid for the purposes of preconditioning each image separately. We first present our construction of the preconditioned string method which has a simpler updating step.

A. Preconditioned String Method
Let us first consider the simple case where P is constant in x. Starting from the coordinate transformation ij . The string method in the transformed space has updating stepx k+1 .
Reversing the coordinate transformation we obtain an equivalent formulation in the original coordinates with updating step where care needs to be taken to normalise the tangents x with respect to the P-norm, y P = (y · Py) 1/2 , instead of the usual 2 -norm, y = (y · y) 1/2 . Expressing the reparametrisation step in terms of coordinates in the configuration space is trivial, as it suffices to replace the usual 2 -norm with the P-norm, due to linearity of the d ds operator. The systems of interest, however, are described by preconditioners that are not constant in the configuration space 12 , which leads to a Riemannian metric framework and in particular Eq. (6) involves the evaluation of ∇P 1/2 (x k n ) which is computationally expensive. We circumvent these issues entirely by dropping these terms. Preliminary tests (which we do not discuss here) showed that this does not lead to any loss of performance. Thus, we obtain the preconditioned string method where we defined the quantity in terms of the P k n = P(x k n ). We are left to specify how to re-parametrise the path. Recall that in the continuous limit, we are free to use any parametrisation for the path. In our setting, the premise is that · P is a more natural notion of distance than the standard 2 -norm · , hence we will use the following notion of distance along the path: We note that d P is not a metric in the technical sense, as it does not satisfy the triangle inequality. However, it is an approximation (discretisation) of the geodesic distance on the Riemannian manifold induced by the preconditioner P, hence it is reasonable to expect that it can be used for the reparametrisation of the path. In practise, we have not encountered any difficulties related to this issue. The details of the preconditioned reparametrisation algorithm are given in Appendix A.

B. Preconditioned NEB method
An entirely analogous argument yields the preconditioned NEB method, where Notice that this class of preconditioning schemes disregards the interactions between images and therefore, the preconditioner aids the convergence of the path only in the transverse direction. This is justified when the main source of ill-conditioning is due to the potential energy landscape, which is the case when only few images are used as is often done in practise. To summarise, the preconditioned updating relations are given by where, in analogy to our earlier notation, (η P ) k n = 0 for the string method and (η P ) k n = (η neb,P ) k n for NEB.

C. ODE solvers and steepest descent
The optimisation steps Eq. (4) were derived by applying Euler's method to the first order differential equation (2), but any ODE solver can be used instead. Here, we use an adaptive ODE solver based on Ref. 18 to allow for some adaptivity in the step selection mechanism.
The user supplies an absolute and a relative tolerance atol and rtol, which control the accuracy of the solution. Choosing these two parameters is more intuitive and more robust than choosing the step length of the static method. We will demonstrate that the ODE solver improves the convergence of the string and NEB methods significantly.
We modify an adaptive ODE solver, ode12 18 . To begin we compute a trial step x k+1 n using Eq. (11) with a given step-length α k . Next, we use x k+1 n to compute a secondorder solution to the underlying ODE system, viã where f k n = −∇ ⊥ V P (x k n ) + (η P ) k n is the driving force on image n at timestep k. We can then use the differencẽ x k+1 n − x k+1 n , or equivalently the difference f k n − f k+1 n as an error indicator.
Taking this as a starting point and following, for example, Ref. 19 to implement an adaptive time-stepping algorithm we obtain an algorithm that underestimates the local error in the neighbourhood of equilibria and in particular will not converge as k → ∞. To overcome this, we add a second step-length selection mechanism based on minimising the residual. In essence, the adaptive ODE step selection should be used in the preasymptotic regime while minimising the residual is a suitable mechamism in the asymptotic regime.
This leads to the following step-length selection algorithm, which we label ode12r : we define the re-scaled residual error and local error where the index j denotes vector components. We then accept the proposed x k+1 n if the scaled residual error satisfies either one of the two following conditions: 2) R k+1 ≤ R k c 2 AND E k+1 ≤ rtol, for contraction and growth parameters c 1 and c 2 ∈ R.
Whether the step is accepted or rejected, we now compute two step-length candidates using (1) the adaptive solver and (2) a simple line-search procedure.
If the step x k+1 is rejected, then the new step-length candidate starting from x k is α k = max 1 10 α k , min 1 4 α k , α k+1 ls , α k+1 ode12 . Figure 1 demonstrates how ode12 effectively selects appropriate step lengths in the pre-asymptotic regime, but stagnates in the asymptotic regime for the case of vacancy migration in tungsten modelled with the EAM4 class of the Embedded Atom Model (EAM) interatomic potential proposed by Marinica et al. 20 . The convergence rate of the modified ode12r agrees with the results of ode12 in the pre-asymptotic regime but successfully converges upon reaching the asymptotic regime.

IV. RESULTS
We tested our preconditioning scheme for a variety of examples. In the following tables we compare the number of force evaluations needed to converge to 'coarse' and 'fine' target accuracies (maximum force less than 10 −1 eV/Å and 10 −3 eV/Å, respectively) using unpreconditioned and preconditioned schemes with either static or adaptive ode12r step selection. The criterion for convergence is the magnitude of the residual error R k+1 as defined in Eq. (12).

A. Vacancy Migration
First we consider the diffusion of a vacancy in a two dimensional 60-atom triangular lattice and governed by a Lennard-Jones potential V (r) = 4 [(σ/r) 12 − 2(σ/r) 6 ] with parameters = 1.0, σ = 2 − 1 6 . The vacancy is located at the centre of the cell initially and migrates in the y direction by one lattice spacing. Periodic boundary conditions are imposed in the x and y directions. Table I shows the computational cost for convergence. The exponential preconditioner introduced in Packwood et al. 12 was used to treat the ill-conditioning of the system allowing convergence beyond the 10 −3 tolerance, which the unpreconditioned case could not achieve within a reasonable number of iterations. The latter came as a suprise to us, as on the contrary to the real vacancy migration systems that we study next, this artificial set up exhibits more severe ill-conditioning.

2D Vacancy
Step selection static ode12r solver Next, we considered a three dimensional system containing a vacancy, specifically a 107 atom Cu fcc supercell in a fixed cell with periodic boundary conditions. Interactions were modeled with a Morse potential with parameters A = 4.0, = 1.0 and nearest neighbour distance r 0 = 2.55266Å for Cu and interactions between atoms expressed by V (r) = (e −2A(r/r0−1) − 2e −A(r/r0−1) ). The exponential preconditioner introduced in Packwood et al. 12 was used. Table II shows the number of force evaluations per image needed for convergence to two preset tolerance limits. This example demonstrates how the ode12r solver can aid the performance of the string and NEB methods if a static step is not suitable. Preconditioning gave a 2-fold speedup for the higher accuracy results, but only marginal improvements for the lower acuracy.

Vacancy in Cu supercell
Step selection static ode12r solver Tol / eV/Å 10 −1 10 −3 10 −1 10 −3 A 53-atom W bcc supercell modelled with the EAM4 potential described in Ref. 20 was examined as well. Periodic boundary conditions were imposed. A force field (FF) preconditioner was constructed from the Hessian of the EAM force field as proposed in Mones et al. 13 giving up to 16 times faster convergence for higher accuracies as shown in Table III.
We studied the same 53-atom W vacancy system with density functional theory (DFT), as implemented in the Castep 21 software. The exchange correlation functional
was approximated by the Perdew, Burke and Ernzerhof (PBE) generalised gradient approximation (GGA) 22 , with a planewave energy cut-off of 500 eV and a 2 × 2 × 2 Monkhorst-Pack grid to sample the Brillouin zone.
Step selection was achieved with the ode12r step selection scheme. The same FF preconditioner was used, based on the EAM Hessian. Traversing the path in subsequent iterations of the NEB and string methods was performed in an alternating order, allowing efficient reuse of previous electronic structure data to start the next optimisation step. Unlike the EAM case above, the preconditioner we used for the DFT model does not describe the potential energy surface of the DFT model exactly, but nevertheless gives a speed-up of a factor of two for an accuracy of ∼ 10 −2 eV/Å and furthermore allows accuracies of the order of ∼ 10 −3 eV/Å to be achieved, unlike the unpreconditioned case, as shown in Fig. 2. The results of Table III suggest that constructing a better preconditioner would improve these results further. Notice further that the number of force evaluations needed for convergence and the time needed for convergence are in agreement (by comparison of the upper and lower panes of Fig 2), confirming that the computational cost of constructing the preconditioner model is negligible compared to the cost of computing DFT forces, justifying our earlier assumptions. We note that the gain of preconditioning would be expected to further increase with system size 12 .

B. Screw Dislocation
In the final example we study a 1 2 111 screw dislocation in a 562-atom W bcc structure confined in a cylinder of radius equal to 20Å and surrounded by an 11Å cylindical shell of clamped atoms, with periodic boundary conditions along the dislocation line (z) direction. The system is simulated with the same EAM4 potential. The dislocation advances by one glide step. Table IV shows the computational costs for converging the MEP with the NEB and string methods, using either static or ode12r step length selection. A force field preconditioner built from the same EAM potential was used for geome-  Upon preconditioning, we observed a 5-fold speed up for the static case for low accuracies but no significant speed up for the ode12r case. For a higher accuracy, a speed up of a factor of 6 was observed with static steps and a speed up of a factor of 3 for the ode12r step selection. Furthermore, there was an overall speed up of a factor of at least 2 from using the ode12r step selection over the static step selection. We investigated this system further, focussing on the NEB implementation to allow comparison with the widely used Limited memory Broyden -Fletcher -Goldfarb -Shanno (LBFGS) 23 optimisation algorithm, which can be used with the NEB implementation 10 in the Atomic Simulation Environment (ASE) 24 . This required fixing the endpoints of the path at the minima as is done in the ASE code. The comparison was carried out on systems of two sizes. A force field preconditioner was used as before for the preconditioned cases. Figure 3 shows the convergence rate of the various NEB schemes for a radius of 20Å in the upper panel (a) and for a radius of 40Å in the lower panel (b). Although LBFGS gave good convergence for the smaller system, it lacks robustness as exemplified by the large spikes along the curve in the larger system. LBFGS's behaviour for the bigger system can be explained by the fact that the force field of the NEB algorithm is not conservative, violating one of LBFGS's assumptions. LBFGS constructs a Hessian matrix corresponding to a scalar field, failing to capture the effects of the transport terms of the NEB force field. Moreover, the lack of the energy function prevents the use of line search required to ensure the method's stability; in the ASE LBFGS implementation a heuristic is instead used to impose a maximum step length of 0.04Å. Convergence of NEB variants for a screw dislocation in a 562-atom W cylindrical structure (a) and a 1489atom W cylindrical structure (b) modeled with the EAM4 Marinica potential 20 . Atoms outside outer radii of R = 20Å and R = 40Å respectively were clamped, with periodic boundary conditions along the dislocation line. The path was discretised with 7 images (excluding the minima at each end, which were held fixed).

V. CONCLUSIONS
We have demonstrated that MEP finding techniques such as the NEB and the string method can exhibit slow convergence rates due to poor search direction and steplength selection during the optimisation procedure. We have introduced a new optimisation technique combining an adaptive time-stepping scheme with preconditioning to address ill-conditioning of the energy landscape in directions transverse to the path and to allow faster convergence to the minimum energy path.
We observed that our new scheme gives a significant speed up and improved robustness over currently used approaches for a range of systems using both force fields and DFT. Moreover, it allows higher accuracies to be reached than existing methods.
However, our preconditioning scheme targets transverse ill-conditioning only. The longitudinal terms, (e.g. the NEB spring interactions) are unaffected by the preconditioner, suggesting that our scheme provides a baseline for further improvements.
An open source implementation of our technique is available at https://github.com/cortner/ SaddleSearch.jl.

Define
This algorithm does not ensure that images will be equidistributed according to d P . However it does ensure that images remain bounded away from one another, which is the key property required for the string method.