An efficient algorithm for finding the minimum energy path for cation migration in ionic materials

The Nudged Elastic Band (NEB) is an established method for finding minimum-energy paths and energy barriers of ion migration in materials, but has been hampered in its general application by its significant computational expense when coupled with density functional theory (DFT) calculations. Typically, an NEB calculation is initialized from a linear interpolation of successive intermediate structures (also known as images) between known initial and final states. However, the linear interpolation introduces two problems: (1) slow convergence of the calculation, particularly in cases where the final path exhibits notable curvature; (2) divergence of the NEB calculations if any intermediate image comes too close to a non-diffusing species, causing instabilities in the ensuing calculation. In this work, we propose a new scheme to accelerate NEB calculations through an improved path initialization and associated energy estimation workflow. We demonstrate that for cation migration in an ionic framework, initializing the diffusion path as the minimum energy path through a static potential built upon the DFT charge density reproduces the true NEB path within a 0.2 Å deviation and yields up to a 25% improvement in typical NEB runtimes. Furthermore, we find that the locally relaxed energy barrier derived from this initialization yields a good approximation of the NEB barrier, with errors within 20 meV of the true NEB value, while reducing computational expense by up to a factor of 5. Finally, and of critical importance for the automation of migration path calculations in high-throughput studies, we find that the new approach significantly enhances the stability of the calculation by avoiding unphysical image initialization. Our algorithm promises to enable efficient calculations of diffusion pathways, resolving a long-standing obstacle to the computational screening of intercalation compounds for Li-ion and multivalent batteries.


INTRODUCTION
The nudged elastic band (NEB) method is an established technique for finding the minimum energy path (MEP) between the given initial and final states of a transition. 1,2 This method has been used in conjunction with Density Functional Theory (DFT) 3-6 and empirical potentials [7][8][9] for studying ion and molecule diffusion in a variety of systems such as semiconductors, metals, and organic molecules. NEB is also widely used for estimating transition states within the harmonic transition state theory (hTST) approximation. 10 In chemistry, the NEB method has been used to characterize transition paths and energetic profiles of reactions occurring on surfaces, 11 in enzymes 12 and solutions, 13 among others. NEB is the method of choice to study vacancy and defect diffusion in alloys and metals. [14][15][16] In the materials science a) Z. Rong and D. Kitchaev contributed equally to this work. b) Author to whom correspondence should be addressed. Electronic mail: gceder@berkeley.edu community and specifically in battery research, NEB has been applied successfully to address Li [17][18][19][20][21][22] and multivalent ion [23][24][25][26]66,67 diffusion in a multitude of cathode materials and ionic conductors. 27,28 In a NEB calculation, a group of images (replicas) of the system is interpolated between the initial and final states. A spring interaction between adjacent images is added to ensure continuity of the path, thus mimicking an elastic band. An optimization of the band, involving the minimization of the force acting on the images, relaxes the band to the MEP. More details of the NEB algorithm can be found in Ref. 50. Over the past two decades, a number of algorithmic improvements have been introduced to increase stability and accuracy. Henkelman et al. proposed the climbing image method 29 and the improved tangent estimate, 11 30 presented the adaptive nudged elastic band method, where NEBs are iteratively calculated to move the initial and final states closer to the saddle point. More recently, Sheppard et al. 31 generalized the NEB method to address solid-solid phase transitions. Crehuet and Field 32 expanded the NEB formalism to account for finite temperature effects. These efforts are accompanied by other work focusing on computational details to accelerate the optimization methods of finding the MEP. 33,34 Recent progress in high-throughput computational infrastructure has opened the door to efficient material innovation through the characterization of material properties by first-principles calculations, [35][36][37] where high-throughput computation is used to search for promising materials for specific applications. Currently, high-throughput firstprinciples calculations are being used to study an increasing range of problems, such as Li-ion battery optimization, [38][39][40][41] nano-porous material design, 42 the search for new catalysts, 43 crystal structure prediction, 44 and the study for surface phenomena. 45 However, despite the success of the NEB method in characterizing the dynamics of materials, the significant computational expense of NEB relative to standard DFT calculations, e.g., geometry relaxations and static energy calculations, has hampered its application in highthroughput work, where properties of thousands of structures are computed and analyzed. In order for the NEB method to be useful in screening materials in a high-throughput fashion, significant improvements are needed in both its stability and efficiency. In this work, we propose a workflow that resolves both issues, thereby enabling a high-throughput automation of NEB migration calculations.
The standard workflow for the NEB calculation consists of 3 main steps: • Step 1. Relax the initial and final state structures.
• Step 2. Linearly interpolate a number of images between the initial and final states. • Step 3. Apply the NEB algorithm to compute the MEP.
We find that the linear interpolation in Step 2 is the primary source of inefficiency and instability in the calculation procedure, especially if the final MEP displays substantial curvature from the initial linear interpolation. Furthermore, during the preparation of the NEB calculation in some systems, the linear interpolation can place atoms (of one image) at unreasonably close distances to one another, causing instability during the NEB relaxation (see, for example, the CaMoO 3 structure in the section titled "Discussion").
Here we present a new method to initialize the NEB interpolation close to the final relaxed band that we call PathFinder Algorithm. In the section titled "Methods," we discuss the idea behind the PathFinder algorithm and give details about its implementation. In the section titled "Results and Discussions," we test the PathFinder algorithm on a set of six materials, demonstrating its predictive capabilities and the computational runtime reduction it brings.
Along with the PathFinder algorithm, we have developed a new approximate method for characterizing the energetic profile of the MEP, hereafter referred as ApproxNEB. ApproxNEB is discussed in detail in the section titled "Methods." We show that ApproxNEB is able to predict migration barriers within an error of ∼20 meV from those obtained with traditional NEB calculations, while reducing the central processing unit (CPU) time by up to a factor of 5.

METHODS
First, we note that while our algorithm is general and independent of any given DFT implementation, in this paper we focus the discussion and implementations to the Vienna Ab Initio Simulation Package (VASP). 46 Nonetheless, we expect our analysis to be directly transferrable to the high-throughput calculation of cation migration barriers within other codes capable of outputting electrostatic potentials.

Path initialization
In the NEB algorithm, each image along the band is relaxed by two forces: the true force from the potential and the spring force (from the virtual springs) connecting adjacent images. Both forces are decomposed into components perpendicular and parallel to the path, and only the perpendicular component of the true force and parallel component of spring force are relaxed in the NEB procedure. The force projection is referred to as "nudging" and leads the chain of images to the MEP. To predict the MEP with fewer computational resources, we would like to imitate this relaxation process starting from a static potential. As the spring forces are very easy to simulate, the difficulty lies in finding a potential that is able to reproduce the true force from first-principles calculations.
The key idea behind PathFinder is that when an atom migrates inside a host structure, it moves to avoid atoms or bonds, as atomic charge density overlap with other species would correspond to reactions, or at least large changes in energy. Consequently, non-reactive migration paths should avoid concentrations of electronic charge density. Thus, we propose using the electronic charge density available from DFT as the potential landscape within which to estimate the migration MEP. In general, this potential will push migrating atoms to regions of diminishing charge density, corresponding to areas void of atoms or bonds, matching the intuition regarding the migration path geometry.
Based on this construction, each of the migrating images relaxes according to the sum of two forces, where ⇀ F1 is the true force acting on the migrating species at its current coordinates, which we approximate using the gradient of the static charge density "potential," and ⇀ F2 is a virtual spring force, which acts as a path arclength penalty, and thereby couples the images of a single migrating ion through the transition state path. In relation to conventional NEB, ⇀ F1 aims to approximate the true force that would otherwise be obtained from a DFT-derived energy gradient, while ⇀ F2 acts analogously to the virtual spring force of the NEB method, ensuring that the images interpolate a continuous path from the initial to the final state. In order to define the two forces ⇀ F1 and ⇀ F2, we rely on ρ, the DFT-derived scalar charge density normalized to the maximum charge density found in the cell, ⇀ r n , which is the position of image n in real space, and k PF , which is the spring force constant for the pathfinder, where all quantities are non-dimensionalized. The non-dimensional spring constant k PF = 0.17 is a constant fit to best reproduce paths from a full VASP NEB calculation with a default NEB spring constant of 5.0 eV/Å. 2 Finally, the positions of all non-migrating atoms can be interpolated linearly for the intermediate images, as their positions are nearly static and thus reasonably represented by the linear path.
An important detail in our method is that the forces used to relax the path are not projected along the path tangent and normal for the spring and real forces, respectively, as is done in traditional NEB. Instead, the force definition and implementation follow that of the simplified string method (also known as the zero-temperature string (ZTS) method 47,48 ) where the homogeneous image spacing is maintained by reparametrizing the path with a new set of evenly spaced images at every iteration of the relaxation algorithm. We choose this approach specifically because it demonstrates superior performance to NEB in cases when a large number of images are accesible, 33 as well as simplicity of implementation.
The PathFinder requires three inputs (see Fig. 1): 1. The initial state structure (structure of the atom-vacancy pair pre-migration jump). 2. The final state structure (structure of the atom-vacancy pair post-migration jump). 3. The charge density of the host structure with vacancies at both the initial and final locations of the migrating ion.
Using these data, we initialize the potential defined in Equations (1)-(3) and optimize the transition path using the steepest descent (gradient) method, until the average displacement of the images used to define the string falls below a predefined convergence threshold, in analogy to the traditional implementation of the ZTS method. 48 Good convergence is generally reached for average displacement values below 5 × 10 −6 Å together with a step size scaling factor on the displacement of ∼0.1. For illustration, Fig. 1 depicts the three inputs to compute Li diffusion paths in LiFePO 4 along the b axis, 49 where in the initial and final states Li-ions sit in the stable sites. The PathFinder algorithm relaxes intermediate Li images along the migration path to positions on the MEP. To initialize the PathFinder algorithm, we compute the charge density of the host structure using a static calculation with Γ−point sampling of reciprocal space, as we have found that the paths thus obtained are sufficiently converged for all test cases. The output of the PathFinder algorithm is the positions of the intermediate images which can then be used to initialize a NEB calculation. As the computational cost of the PathFinder itself is negligible compared to the full NEB calculation, we find that it is effective to use a large number of interpolated images in the PathFinder to ensure the optimal convergence of the string method, followed by a selection of a smaller subset of evenly spaced images to initiate the full NEB calculation.
The complete code set and an example for using the algorithm are available on the github code repository, 50 or as part of the MAST package, 68 and the code implementation depends on the Python Materials Genomics (pymatgen) library. 42

Static barrier estimation
While the PathFinder algorithm can provide a good approximation of the geometry of the MEP, it does not yield energetics along the path. For this reason, we have developed the complementary package ApproxNEB that allows one to estimate the energies of each image by decoupling the band into individual image calculations. The right panel depicts the path relaxation in the PathFinder algorithm, where the path is iteratively relaxed through the virtual potential derived from the DFT electronic charge density and the spring force. The potential is shown color-coded by magnitude with equipotential contours depicted by dashed lines, and with white arrows indicating the direction of relaxation. Finally, the lower left panel depicts the final relaxed Li MEP path produced by the PathFinder algorithm, which is in close agreement to the MEP obtained from a full NEB calculation (see Fig. 4(b)).

FIG. 2. A comparison between traditional NEB and ApproxNEB schemes.
Here it is assumed that 7 images are interpolated between the initial and final states (image 00 and image 08 are the initial and final states). While in both cases the relaxation of each image is iterative (taking steps 01, 01 ′ , 01 ′′ , etc.), in ApproxNEB the relaxation iterations of separate images are decoupled, decreasing the computational burden of the mean energy path (MEP) sampling.
The key idea behind ApproxNEB is that, if we fix the moving cation along the approximate MEP obtained from the PathFinder algorithm, and perform a single point relaxation image by image, we can access the missing energetics, thereby fully characterizing the MEP. The difference between ApproxNEB and NEB algorithms is depicted in Fig. 2. In general, the execution of the NEB algorithm, in first-principles or classical potential codes, requires communication between images, as they are connected by virtual springs. At the end of each ionic relaxation step, images communicate with each other to update spring forces, and a new step in the constrained potential energy is taken-this procedure is repeated iteratively until the NEB force and energy criteria are satisfied. The ApproxNEB method removes the spring force and estimates the migration barrier by fixing the positions of the moving ion and relaxing other atoms in each image. In order to constrain the translational degrees of freedom of the system, this procedure requires that the position of a reference atom that is farthest away from the moving ion in the unit cell to be fixed. This constraint prevents the whole cell from shifting uniformly to translate into the initial or final state. Because the framework ions are already very close to the local-minimum positions as they come from fully relaxed host structures, the quasi-Newton RMM-DIIS algorithm is sufficient to achieve fast convergence during ion relaxation. 63 Under these constraints, the energy of the independently relaxed images provides an approximate MEP trajectory.
As discussed earlier, in NEB calculations, the migrating ions are relaxed by a combination of virtual spring forces and true forces, while non-migrating atoms are relaxed only by the true forces. The spring forces serve to push the migrating ions to higher energy positions on the MEP. However, by knowing a priori the geometry of the MEP from the PathFinder method, the spring forces can be removed by fixing the moving cation on the MEP. From this perspective, ApproxNEB and NEB provide equivalent constraints on the system during relaxation.

RESULTS AND DISCUSSIONS
To assess the capabilities of PathFinder and ApproxNEB, we apply these methods to the migration of cations in a set of materials that are of practical interest in the field of batteries. As we expect the PathFinder and ApproxNEB methods to yield improved performance relative to standard NEB in cases where the migration paths deviate substantially from the straight-line paths, we report the curvature of the MEP as obtained by the NEB calculations.
The geometry error for the cation migration path for each benchmark material is given in Fig. 4. Specifically, for each material, we compare the error of the PathFinder path and the standard linear interpolation, with respect to the NEBconverged MEP, in order to understand which interpolation scheme can serve as a superior initialization. Note that while in the NEB calculations of the benchmark materials, seven images are used to interpolate the migration path, in the PathFinder algorithm, we use 21 images to ensure good performance of the string method. However, for consistency, in Fig. 4(b), we only show seven equally spaced images for visualization. As can be seen from Fig. 4(a), if the fully relaxed NEB path possesses a large degree of curvature, the PathFinder algorithm systematically provides a better initialization than the traditional linear interpolation. The migration path derived from the PathFinder algorithm falls within 0.2 Å of the NEB-derived MEP in all the test structures, which is a very small error in absolute terms, and is ≈5 times smaller than the typical error obtained from linear interpolation. Fig. 4(b) shows this agreement visually for LiFePO 4 and MgV 2 O 5 . In both structures, the PathFinder algorithm reliably yields migration path geometries very close to the true MEP structures, capturing the effect of nearby oxygens on the cation migration trajectories. In the cases where the MEP is linear, the linearly interpolated initial band usually shows a slightly smaller error than the PathFinder-derived path, as a linear interpolation is by circumstance already the optimal configuration. Nonetheless, the error of the PathFinder-derived path remains within the 0.2 Å bound observed earlier. This error is a sufficiently small absolute error that we can expect its effect on the NEB calculation speed, accuracy, and stability to be negligible, as compared to the traditional linear interpolation scheme. Thus, the PathFinder algorithm offers a robust estimate of cation migration MEPs, yielding a migration path within a small error of the true MEP for both linear and curved geometries, offering both an efficient estimate of MEP geometry and a reliable initialization for subsequent NEB calculations.
To characterize the computational efficiency gains through the PathFinder initialization, we compare the runtime of NEB calculations initialized using the PathFinder scheme versus the traditional linear interpolation. The computational resources are measured by the total CPU hours used on a Cray XC30 machine with a parallelization of 24 cores per image. To ensure a fair comparison, all computational parameters are kept the same for the two-initialization schemes. The results of our test are given in Fig. 5. As could be expected from our analysis of MEP geometry, initialization using the PathFinder algorithm does not significantly affect performance for structures with a linear MEP for migration, but does lead to consistent performance gains in cases where the MEP deviates substantially from a linear path.
Having established the PathFinder approach as a reliable method to efficiently estimate migration geometries, we turn to the ApproxNEB approach of characterizing the energetics of the MEP. To assess the validity of this approach, we compare the overall energy profile of the MEP and the migration barrier obtained from the ApproxNEB algorithm to those obtained from a traditional NEB scheme. As can be seen in Fig. 6, the two methods yield energy profiles and migration barriers within 20 meV of each other, suggesting that ApproxNEB is able to reproduce the results of NEB to good agreement across a variety of systems and migration geometries. As shown in Fig. 6(b), the barriers obtained from the ApproxNEB method are close but systematically higher than those obtained from NEB. This trend is to be expected in the ApproxNEB scheme, because the moving cation is fixed on the path provided by the PathFinder algorithm. By constraining the position of the diffusing species in each image, we reduce the number of degrees of freedom available during relaxation as compared to traditional NEB, such that any error in the MEP geometry obtained from the PathFinder translates to an increase of the migration barrier. However, just as the absolute error in the estimated MEP geometry remains within 0.2 Å across all tested systems, the error in the migration barrier remains within 20 meV, and represents a sufficiently small error margin for most high-throughput screening applications.
It is worth noting that the computational resources necessary for ApproxNEB are substantially lower than for traditional NEB, further justifying its use in high-throughput screening applications. As can be seen in Fig. 7, for both linear and curved paths, ApproxNEB is systematically faster than the NEB method. Notably, we find that in the case of structures with curved MEP geometries, ApproxNEB yields a speedup by a factor of up to 5 with respect to linearly initialized NEB, offering a significant improvement over even PathFinder-initialized NEB as discussed earlier. The reason for FIG. 7. CPU hours consumed by ApproxNEB and NEB methods. For the NEB method, the band is initialized from the linear interpolation. However, the performance gains of ApproxNEB are significantly higher than even PathFinder-initialized NEB shown in Fig. 5. While the ApproxNEB method requires an initial PathFinder calculation, we do not include the computational time required for PathFinder step as it is in all cases negligible. this improvement lies in the decoupling of image calculations from one another. Decoupled images experience a much simpler potential field that remains quasi-static throughout the relaxation, enabling efficient minima-searching during ionic relaxation.
Another issue is that of parallelization-in traditional NEB, because the position of the moving species must be communicated among images to update spring forces, every image must be at the same ionic relaxation step (see Fig. 2). This constraint limits the progress of the calculation since converged images have to wait until all other images reach convergence before the next NEB step is taken. Finally, error handling becomes much easier in ApproxNEB. In the traditional NEB scheme, if an image calculation fails due to a convergence issue, the whole calculation must be restarted. Given that in ApproxNEB each image is independent, only the failed images need to be recomputed. The improvements in both computation runtime and error handling make PathFinder and ApproxNEB suitable for scaling up to screen material properties in a high-throughput fashion.
The final advantage of the PathFinder initialization and ApproxNEB barrier estimation scheme is reflected in the improved calculation stability. One of the common issues in NEB calculations is that linear interpolation can yield highly unphysical initializations with image structures that are difficult to relax due to exceptionally high forces and instabilities occurring during the electronic minimization. PathFinder avoids this problem by biasing the migrating ion away from concentrations of electronic charge density, escaping unintended reactions in the intermediate images.
For example, when calculating the MEP of Ca inner-layer migration along the a axis in CaMoO 3 (see Fig. 8), we find that typical NEB with linear interpolation is unstable due to excessive forces in some images. The reason for this instability can be seen in Fig. 8-the initialization of the NEB calculation by the linear interpolation places one oxygen atom (colored in yellow) very close to a Ca ion in some images, an issue which is avoided by the PathFinder. The unphysically small Ca-O distance results in large inter-atomic forces, destabilizing the calculation. Conventionally, such instabilities are mitigated by the careful tuning of convergence and relaxation parameters or by "chemical intuition," resulting in a significant increase in runtime and human labor,

CONCLUSION
In this report, we have proposed a new scheme for estimating migration minimum-energy path (MEP) geometry and energetics. By testing our methodology against standard NEB calculations and literature values, we find that the PathFinder algorithm can reliably predict the geometry of cation migration MEP within 0.2 Å at negligible computational costs. Furthermore, we find that the ApproxNEB calculation scheme yields activation barriers for the migration within an error bound of 20 meV while using significantly fewer computational resources than traditional NEB schemes. We envision that our methods can be used to accelerate NEB calculations, as well as to provide a robust estimation criterion for migration barriers in ionic materials for high-throughput computational screening of materials.

SUPPLEMENTARY MATERIAL
See supplementary material for addressing the strategy adopted to implement the ApproxNEB method. In Figure S1 of the supplementary material we show the two implementation workflows to combine PathFinder and ApproxNEB methods.