Scaled and Dynamic Optimizations of Nudged Elastic Bands

TL;DR

This paper introduces dyNEB, a dynamic optimization method for nudged elastic bands that reduces force calls by over 50% while accurately locating saddle points in solid state reactions.

Contribution

dyNEB is a novel, adaptive approach that selectively optimizes band regions based on force criteria, improving efficiency without sacrificing resolution.

Findings

Reduces force calls by more than 50%.

Effective for solid state reaction barriers.

Maintains band continuity and saddle point resolution.

Abstract

We present a modified nudged elastic band routine that can reduce the number of force calls by more than 50% for bands with non-uniform convergence. The method, which we call "dyNEB", dynamically and selectively optimizes states based on the perpendicular forces and parallel spring forces acting on that region of the band. The convergence criteria are scaled to focus on the region of interest, i.e., the saddle point, while maintaining continuity of the band and avoiding truncation. We show that this method works well for solid state reaction barriers---non-electrochemical in general and electrochemical in particular---and that the number of force calls can be significantly reduced without loss of resolution at the saddle point.

Tables

Table 1. Table 1: Force calls and normalized CPU time for a saddle point calculation of the Tafel reaction (Fig.  3 ) with eight interior states.

Method		$\Delta E^{‡}$ (eV)	Force calls	Normalized
				CPU time
Default parallel	no climbing	—	1,960	1.0
	climbing	0.564	224	1.0
${\alpha }_s=1$	no climbing	—	991	0.56
	climbing	0.563	93	0.49
${\alpha }_s=2$	no climbing	—	857	0.49
	climbing	0.566	70	0.39

Table 2. Table 2: Force calls and normalized CPU time for a semi-grand canonical saddle point calculation of the Heyrovsky reaction (Fig.  4 ) with eight interior states.

Method		$\Delta E^{‡}$ (eV)	Force calls	Normalized
				CPU time
Default parallel	no climbing	—	2,416	1.0
	climbing	0.601	304	1.0
${\alpha }_s=2$	no climbing	—	1,829	0.77
	climbing	0.607	46	0.14
${\alpha }_s=3$	no climbing	—	1,249	0.47
	climbing	0.602	90	0.22

Equations

$$\mathbf{F}_{i}=\mathbf{F}_{i}^{s}\big{|}_{\parallel}-\nabla V(\mathbf{R}_{i})\big{|}_{\perp}.$$

$$f_{\max }^i=f_{\max }^{i_{\max }}\cdot \left(1+\sqrt{\sum _{j=0}^N{\left({\mathbf{R}}_{i_{\max },j}-{\mathbf{R}}_{i,j}\right)}^2}\cdot {\alpha }_{\mathrm{s}}\right),$$

Full text

Reaction barrier estimates are pivotal in computational studies of chemical systems; reaction rates in homogeneous and heterogeneous catalysis, alloy stability Vej-Hansen et al. (2016) and current densities in electrochemistry Skúlason et al. (2007, 2010); Goodpaster et al. (2016); He et al. (2017); Singh et al. (2017), to name a few, all depend upon the height of the transition state separating two thermodynamically stable states. Kinetic information from ab initio calculations can be used to calculate rate constants via transition state theory Eyring (1935) and connect microscopic quantities to macroscopic observables. This information is contained on many-dimensional potential energy surfaces (PES)—the central theme in computational chemistry. Here, thermodynamically stable states are located in local minima, and the transition from one state to another depends upon the height of the transition state region separating the two basins. Transition states are first-order saddle points—regions on the PES with convex and concave curvature—and represent the lowest energy region connecting two minima. Saddle point calculations are typically categorized based on the information needed to initialize the calculation Jensen (2017). Some methods, notably the Dimer method Henkelman and Jónsson (1999), require only local information, while chain-of-states methods require two fixed endpoints. In the former, saddle points are located based on the lowest curvature out of the local minimum. While computationally efficient, these methods require statistical sampling to ensure that the relevant saddle point—and not just any saddle point—is found Henkelman and Jónsson (1999).

The nudged elastic band (NEB) Jónsson et al. (1998); Henkelman and Jónsson (2000); Henkelman et al. (2000) method is one of the most widely used algorithms for (first-order) saddle point calculations. This chain-of-states method connects two thermodynamically stable states by intermediate states and Hookean springs. These springs allow the intermediate states to traverse high-energy regions of the PES. The nudging force projection, after which the method is named, only includes spring forces parallel to the path and PES–derived forces perpendicular to the reaction pathway (Eq. 1) Henkelman and Jónsson (2000),

[TABLE]

The chain-of-states is conventionally initialized by linear interpolation between the initial and final states or by more sophisticated interpolation schemes Smidstrup et al. (2014), and the band is subsequently optimized to locate the minimum energy pathway (MEP) connecting the two basins. Although rigorous and widely used throughout computational chemistry, the NEB method is often hampered by its significant computational cost. This is due in part to size (all states must be simultaneously optimized) and the dimensionality of the optimization problem. Many methods have been proposed to alleviate the computational cost, ranging from adaptive bands Maragakis et al. (2002) and automated methods Kolsbjerg et al. (2016) to reflection symmetry operations Mathiesen et al. (2019) and machine-learning algorithms Peterson (2016); Koistinen et al. (2017); Garrido Torres et al. (2019).

Recent efforts in the electrocatalysis community have focused on developing constant-potential density functional theory (DFT) Hohenberg and Kohn (1964); Kohn and Sham (1965) methods. These electronically grand canonical methods Lozovoi et al. (2001); Otani and Sugino (2006); Taylor et al. (2006); Jinnouchi and Anderson (2008); Letchworth-Weaver and Arias (2012); Goodpaster et al. (2016); Sundararaman et al. (2017); Kastlunger et al. (2018); Bouzid and Pasquarello (2018); Melander et al. (2019) allow the number of electrons in the system to fluctuate along the reaction trajectory in order to keep the applied potential (work function) constant. Such schemes are required when calculating electrochemical reaction barriers, since finite-sized unit cells introduce an artificial potential bias along the reaction trajectory in constant-charge DFT Rossmeisl et al. (2008); Skúlason et al. (2010). The computational cost of NEB calculations can grow manifold when coupled with these semi-grand canonical methods, since the geometry of the band and the work function of each state along the reaction trajectory must be simultaneously optimized. In this emerging field, there is need for a computationally efficient approach to reaction barrier calculations.

Part of the allure of the NEB method is that it is embarrassingly parallel; computational resources can be efficiently partitioned to calculate all interior states simultaneously. This is certainly true for bands with uniform convergence. However, as we will show, convergence of bands is often highly non-uniform, and significant computational resources are often spent calculating states that are well below the convergence criterion. This issue is particularly severe for semi-grand canonical saddle point searches. Moreover, the efficiency of parallel schemes not only hinges on uniform convergence of interior states, but also on homogeneous computer architectures.

In this Letter, we introduce a methodology to improve saddle point convergence for non-uniform MEPs. Moreover, we demonstrate that an electronically grand canonical treatment makes convergence of the MEP inherently non-uniform, since charge equilibration is particularly important at or around the saddle point Kastlunger et al. (2018). This makes parallel NEB schemes unsuitable, since a large fraction of computational resources will idle for significant parts of the optimization while a small fraction is used to optimize the region of interest, i.e., the saddle point. We present a serial NEB implementation that dynamically and selectively optimizes the reaction pathway in NEB calculations. While this scheme is applicable to any NEB calculation, we note that it is particularly well-suited for semi-grand canonical saddle point searches.

Convergence of chains-of-states is often highly non-uniform; states in low-energy regions of the PES (in closer vicinity to local minima) generally converge faster than those close to the saddle point. Moreover, in all but a few applications, the important result of an NEB calculation is the saddle point. The interior states before and after the saddle point contribute to mapping out a reasonable minimum energy pathway and provide tangent estimates for the nudging force projection (Eq. 1), but are typically not used in subsequent analyses. Thus, the computational effort can be significantly reduced when optimizing states to different convergence criteria. We do this by scaling the convergence criteria along the reaction trajectory and dynamically optimizing the states. A state is not recalculated if the forces acting on that state are below the convergence criterion; if the forces rise above the convergence criterion—due to spring adjustments between the state and its neighbors—the state is recalculated, hence dynamic. This convergence check is performed after each optimization step. Fig. 1 shows an example of dynamic optimizations for three states in the associative desorption reaction of hydrogen on Au(111) when the climbing image Henkelman et al. (2000) implementation of NEB is invoked. For simplicity, we consider a simple convergence scaling based on state indices, $f_{\mathrm{max}}^{i}=f_{\mathrm{max}}^{i_{\mathrm{max}}}\cdot\left(1+|i_{\mathrm{max}}-i|\cdot\alpha_{\mathrm{s}}\right)$, where $i_{\mathrm{max}}$ is the state with the highest potential energy and $\alpha_{\mathrm{s}}$ is the scaling factor. In Fig. 1, the $i^{\mathrm{th}}$ state is the climbing image, and its geometry changes drastically when the state climbs up the gradient of the PES. This changes the parallel spring force, $\mathbf{F}_{i}^{s}|_{\parallel}$, between states $i$ and $i-1$ (although the climbing image does not feel the spring forces), and subsequently pulls state $i-1$ out of convergence. State $i-1$ is then recalculated, and the change in the parallel spring force, $\mathbf{F}_{i-1}^{s}|_{\parallel}$, between states $i-1$ and $i-2$ causes state $i-2$ to be recalculated. This implementation is therefore truly dynamic, since all states dynamically respond to any perturbations in the band. Convergence plots for all states of this example are shown in Figs. 2–3 of the supporting information.

The choice of convergence scaling is in principle arbitrary, but a very high scaling factor could disrupt the pathway (see Fig. 4 of the supporting information). This highlights the importance of dynamic optimizations, since “frozen” states can be perturbed away from the MEP if the geometries of their neighboring states change. This is particularly important when using the climbing-image Henkelman et al. (2000) implementation of NEB, since the geometry (and number of electrons in electronically grand canonical methods) changes significantly in this region of the PES. We find that high scaling factors only introduce modest changes in the number of force calls compared to a more conservative convergence scaling (Fig. 2 and table Scaled and Dynamic Optimizations of Nudged Elastic Bands), and could disrupt the pathway and delay convergence. Thus, a low scaling factor is recommended. Here, we choose to scale the convergence criteria based on a displacement metric,

[TABLE]

where $i$ loops over all interior states, $i_{\mathrm{max}}$ is the index of the interior state with the highest potential energy, $f_{\mathrm{max}}^{i_{\mathrm{max}}}$ is the convergence criterion given to the optimizer, $\mathbf{R}_{i}$ is the position matrix of the $i^{\mathrm{th}}$ state, $j$ is the atom index, $N$ is the number of atoms and $\alpha_{\mathrm{s}}$ is the scaling factor. The advantage of this implementation is that it is independent of sampling density along the band. We emphasize that this convergence scaling is dynamic; $\alpha_{\mathrm{s}}$ is assigned a priori, but $i_{\mathrm{max}}$ and $\mathbf{R}_{i}$ are updated after each optimization step. We have implemented this scaled and dynamic optimization method of nudged elastic bands in the Atomic Simulation Environment (ASE) Larsen et al. (2017).

As a simple example, we consider oxygen diffusion on Pt(111). Fig. 2 shows the MEP for the default NEB implementation (black) and “dyNEB” with different convergence scaling factors, as calculated with Effective Medium Theory Jacobsen et al. (1996). Here, the pathway leading to the saddle point differs between the two methods; the tails of the bands converge quickly, and are thus unnecessarily recalculated in the default implementation. However, the saddle point geometry and energy are identical in the two methods, and the reduction in force calls is significant; the highest scaling ($\alpha_{\mathrm{s}}=6$) of the convergence criteria reduces the number of force calls by 75%. We note that a dynamic relaxation without convergence scaling ($\alpha_{\mathrm{s}}=0$) reduces the number of force calls by 59% for this particular reaction barrier.

We use the method outlined above the calculate (i) non-electrochemical diffusion barriers and (ii) electrochemical reaction barriers at constant applied potential with the Solvated Jellium (SJ) method Kastlunger et al. (2018), and show that the number of force calls in NEB calculations can be significantly reduced without loss of resolution at the saddle point. However, the method is equally applicable to molecular saddle point searches. Each benchmark calculation is started from linear interpolation, and the climbing-image method is only invoked after the band has converged. All computational details are listed in the supporting information.

The non-electrochemical associative desorption of hydrogen (Tafel) of the hydrogen evolution reaction (HER) can be calculated without electronically grand canonical methods, since the charge transfer between the initial and final states (and any states along the reaction trajectory) is negligible. Fig. 3 shows the converged pathway on Au(111) for “dyNEB” with $\alpha_{\mathrm{s}}=2$ in black and the default implementation in gray. The implementations differ only in the convergence of states away from the saddle point; states far away from the saddle point converge faster and to a higher convergence criterion, and thus require less force calls. Note that the convergence scaling in Eq. 2 captures the asymmetry of the reaction pathway.

However, the geometry and energy of the saddle point are identical in both implementations. As shown in table Scaled and Dynamic Optimizations of Nudged Elastic Bands, this dynamic and selective optimization method rapidly locates the saddle point; “dyNEB” with $\alpha_{\mathrm{s}}=2$ reduces the total number of force calls and CPU time by 58% and 52%, respectively, compared to a default parallel implementation. Thus, the slightly higher CPU time–per–force call of “dyNEB”, which stems from the serial implementation, is offset by a significant reduction in the number of force calls. Moreover, the reduction in force calls is more pronounced when the climbing-image method is invoked; $\alpha_{\mathrm{s}}=2$ decreases the number of force calls by 69% compared to the default methodology.