Saddle Point Search Research Articles

Exploring potential energy surfaces to reach saddle points above convex regions.

Saddle points on high-dimensional potential energy surfaces (PES) play a determining role in the activated dynamics of molecules and materials. Building on approaches dating back more than 50 years, many open-ended transition-state search methods have been developed to follow the direction of negative curvature from a local minimum to an adjacent first-order saddle point. Despite the mathematical justification, these methods can display a high failure rate: using small deformation steps, up to 80% of the explorations can end up in a convex region of the PES, where all directions of negative curvature vanish, while if the deformation is aggressive, a similar fraction of attempts lead to saddle points that are not directly connected to the initial minimum. In high-dimension PES, these reproducible failures were thought to only increase the overall computational cost, without having any effect on the methods' capacity to find all saddle points surrounding a minimum. Using activation-relaxation technique nouveau (ARTn), we characterize the nature of the PES around minima, considerably expanding on previous knowledge. We show that convex regions can lie on activation pathways and that not exploring beyond them can introduce significant bias in the saddle-point search. We introduce an efficient approach for traversing the convex regions, almost eliminating exploration failures, while multiplying by almost 10 the number of identified unique and connected saddle points as compared to the standard ARTn, thus underlining the importance of correctly handling convex regions for completeness of saddle point explorations.

The Importance of Reaction Energy in Predicting Chemical Reaction Barriers with Machine Learning Models.

Improving our fundamental understanding of complex heterocatalytic processes increasingly relies on electronic structure simulations and microkinetic models based on calculated energy differences. In particular, calculation of activation barriers, usually achieved through compute-intensive saddle point search routines, remains a serious bottleneck in understanding trends in catalytic activity for highly branched reaction networks. Although the well-known Brønsted-Evans-Polyani (BEP) scaling - a one-feature linear regression model - has been widely applied in such microkinetic models, they still rely on calculated reaction energies and may not generalize beyond a single facet on a single class of materials, e. g., a terrace sites on transition metals. For highly branched and energetically shallow reaction networks, such as electrochemical CO2 reduction or wastewater remediation, calculating even reaction energies on many surfaces can become computationally intractable due to the combinatorial explosion of states that must be considered. Here, we investigate the feasibility of activation barrier prediction without knowledge of the reaction energy using linear and nonlinear machine learning (ML) models trained on a new database of over 500 dehydrogenation activation barriers. We also find that inclusion of the reaction energy significantly improves both classes of ML models, but complex nonlinear models can achieve performance similar to the simplest BEP scaling when predicting activation barriers on new systems. Additionally, inclusion of the reaction energy significantly improves generalizability to new systems beyond the training set. Our results suggest that the reaction energy is a critical feature to consider when building models to predict activation barriers, indicating that efforts to reliably predict reaction energies through, e. g., the Open Catalyst Project and others, will be an important route to effective model development for more complex systems.

COMPASS: Double-ended saddle point search as a constrained optimization problem.

We present an algorithm to find first order saddle points on the potential energy surface (PES). The algorithm is formulated as a constrained optimization problem that involves two sets of atomic coordinates (images), a time-varying distance constraint and a constraint on the energy difference. Both images start in different valleys of the PES and are pulled toward each other by gradually reducing the distance. The search space is restricted to the pairs of configurations that share the same potential energy. By minimizing the energy while the distance shrinks, a minimum of the constrained search space is tracked. In simple cases, the two images are confined to their respective sides of the barrier until they finally converge near the saddle point. If one image accidentally crosses the barrier, the path is split at suitable locations and the algorithm is repeated recursively. The optimization is implemented as a combination of a quasi-Newton optimization and a linear constraint. The method was tested on a set of Lennard-Jones-38 cluster transitions and a set of 121 molecular reactions using density functional theory calculations. The efficiency in terms of energy and force evaluation is better than with competing methods as long as they do not switch to single-ended methods. The construction of a continuous search path with small steps and the ability to focus on arbitrary subsegments of the path provide an additional value in terms of robustness and flexibility.

Scalar auxiliary variable approach in iterative minimization formulation for saddle point search

Scalar auxiliary variable approach in iterative minimization formulation for saddle point search

PARTn: A plugin implementation of the Activation Relaxation Technique nouveau that takes over the FIRE minimisation algorithm

Nowadays, the interoperability and interfacing of codes and libraries have become crucial aspects of software development and engineering, and the basis for enabling and simplifying the sharing of methods and tools, both within and among communities. One of the most important bottlenecks that arises when developing and maintaining an interface of a library with an already existing software, is to keep it aligned with the development route of the latter. This might include significant changes, such as changes in the data structures used by the library, which are communicated through the interface.In this paper, an approach for inserting a new algorithm into existing software is presented, through a minimally invasive interface, that takes over an already present algorithm, and thus changes its original purpose. The approach is applied to the well-established Activation-Relaxation Technique nouveau (ARTn) algorithm, that is revisited and re-engineered to bias and take over the FIRE minimization algorithm, as presently implemented in two community codes for atomistic simulations, namely Quantum ESPRESSO (PWscf) and LAMMPS. ARTn is a well established single-ended saddle-point search algorithm that allows for the exploration of potential energy surfaces. The resulting algorithm acts as a plugin, and is distributed in the form of an external library (pARTn). Program summaryProgram Title: plugin Activation Relaxation Technique nouveau (pARTn)CPC Library link to program files:https://doi.org/10.17632/xpwr56d2yk.1Developer's repository link:https://gitlab.com/mammasmias/artn-pluginLicensing provisions: Apache-2.0/GPLv3Programming language: Modern Fortran and C++Nature of problem: Original ARTn implementation was difficult to interface, port and maintain.Solution method: Full refactoring and re-engineering into a library, which utilizes and biases a minimization algorithm already implemented in the engine, to make it behave as ARTn. The biasing is done by reverse-engineering specific external conditions that drive the displacement of atoms.

Meso-timescale atomistic simulations on coalescence process of He bubbles in Fe by SEAKMC method

He bubbles are characteristic microstructures under fusion reactor conditions. They approach and coalesce through their own migration, which significantly impacts the microstructure and material properties. However, these processes, which involve multiple migrations of metal atoms, cannot be treated by molecular dynamics (MD) due to its timescale limitation. In this study, self-evolving atomistic kinetic Monte Carlo (SEAKMC) was used to expand the timescale and reproduce bubble coalescences in Fe. To enhance selections of events that led to the process by avoiding trivial events with an extremely low activation energy such as tiny vibrations of a He atom or short-range displacements of the Fe atom, we introduced two algorithms into SEAKMC, a two-step saddle point search for the former measure and setting a threshold for a displacement distance of the Fe atom for the latter. Furthermore, by adding another algorithm to set an upper bound for the activation energy to prevent selections of events with an impractically high activation energy, we succeeded to reproduce the change in the configuration from dumbbell to elliptical up to a simulated time of 10−1 s, 8 orders longer than MD timescales. The developed method is effective for analyzing microstructures of metallic materials containing light elements and is the only method that can reach timescales comparable to those of experiments.

Saddle point search with dynamic active volume

Sampling potential energy surface (PES) is critical for many problems in materials science, chemistry, physics, and biology and requires highly efficient saddle point searches (SPS). In the study, we introduce the concept of dynamic active volume (DAV) in addition to the active volume in self-evolving atomistic kinetic Monte Carlo (SEAKMC). The DAV method has further reduced the dimensionality of the PES at the elevation stage of a SPS. At the subsequent converging stage, the dynamic boundary in DAV is lifted to allow the system to converge to the right location of a saddle point. Coupled with the dimer method, the DAV method not only significantly reduces the time cost for a given search attempt, but also dramatically increases the probability of finding relevant saddle points for PES sampling. A Python software package within the framework SEAKMC (SEAKMC_py) with the DAV method has been developed.

Atomistic learning in the electronically grand-canonical ensemble

A strategy is presented for the machine-learning emulation of electronic structure calculations carried out in the electronically grand-canonical ensemble. The approach relies upon a dual-learning scheme, where both the system charge and the system energy are predicted for each image. The scheme is shown to be capable of emulating basic electrochemical reactions at a range of potentials, and coupling it with a bootstrap-ensemble approach gives reasonable estimates of the prediction uncertainty. The method is also demonstrated to accelerate saddle-point searches, and to extrapolate to systems with one to five water layers. We anticipate that this method will allow for larger length- and time-scale simulations necessary for electrochemical simulations.

Restricted multicanonical sampling for machine learning potential construction

The restricted multicanonical (rMUCA) ensemble method is developed and combined with the on-the-fly machine learning potential (MLP) generation scheme. The rMUCA simulation performs a random walk in the potential-energy subspace restricted by the selected collective variables and allows us to sample physically relevant configurations without being trapped in local energy minima. No preliminary simulation runs are required to construct a bias potential. The sample structures for training are collected dynamically from the simulations using the MLP itself where the simultaneous error estimation is utilized to judge whether an updated structure should be added to the sample data set or not. The rMUCA formula can be also used for the saddle-point search with minor modification. The method is applied to the oxidation of carbon monoxide on platinum surfaces. The results show that the rMUCA simulation provides an efficient and accurate way to sample rare events for the MLP construction.

Long-timescale transformations of self-interstitial atom clusters of Cu using the SEAKMC method: The effect of setting an activation energy threshold for saddle point searches

On-the-fly kinetic Monte Carlo (kMC), a computational technique for atomistic simulations, has attracted attention because it increases the simulation timescale beyond that of molecular dynamics (MD) simulations while maintaining atomistic fidelity. However, for most kMC methods, when events with high and low activation energies coexist in the event list, trivial events with extremely low activation energies that do not essentially affect the phenomena of interest, so-called flicker events, are frequently selected, making it challenging to observe the key dynamics. In this study, we use Self-Evolving Atomistic kMC (SEAKMC), one of the on-the-fly kMC methods, to model the unstable-to-stable transformations of irregular three-dimensional self-interstitial-atom (SIA) clusters in Cu generated through collision cascade. By setting an activation energy threshold once every five steps, transformations into stable configurations are enhanced. The algorithm renders the simulation timescales one or two orders of magnitude longer than those possible with MD simulations. Further, the probability of transformations into stable configurations is increased by 40 times compared to that of the original SEAKMC method. In addition, we find that the stable configurations obtained by the transformation of the SIA clusters are mostly Frank loops. In summary, this new algorithm for the SEAKMC method helps to resolve the inefficiency of kMC methods resulting from the selection of flicker events and will aid the study of meso-timescale atomistic dynamics.

Discovery of multimechanisms of screw dislocation interaction in bcc iron from open-ended saddle point searches

Dislocation motion and interactions determine mechanical properties in body-centered cubic (bcc) metallic materials. However, studying mechanisms for the screw dislocation interaction is fundamentally challenging since many underlying processes involve mesotimescales and atomistic resolution, currently inaccessible by either experimental techniques or continuum theoretical methods. In this paper, we develop a computational capability based on self-evolving atomistic kinetic Monte Carlo (SEAKMC) to sample the critical events and saddle point energies related to screw dislocations and their junctions. The method is first validated by calculating the stress dependence of Peierls barriers and formation energies of kink pairs and cross-slip kink pairs on a single screw dislocation in bcc iron. Then the method is applied to a binary junction of a pair of intersecting screw dislocations, the structure of which is crucial for low-temperature plastic deformation. We identify three important mechanisms: coplanar cross-slipping, jog-pinning, and a previously unknown unzipping mechanism during the evolution of the binary junction. The mechanisms are then further validated using classical molecular dynamics simulations. The computational capability developed in this paper provides an effective tool to evaluate screw dislocation related thermally activated events in complex stress conditions. The mechanisms discovered in this paper provide critical insights into temperature dependence of the anomalous slip, a breakdown of the Schmidt law, during the plastic deformation in bcc iron and can be generalized to other bcc metals.

Accelerating off-lattice kinetic Monte Carlo simulations to predict hydrogen vacancy-cluster interactions in [formula omitted]-Fe

We present an enhanced off-lattice kinetic Monte Carlo (OLKMC) model, based on a new method for tolerant classification of atomistic local-environments that is invariant under Euclidean-transformations and permutations of atoms. Our method ensures that environments within a norm-based tolerance are classified as equivalent. During OLKMC simulations, our method guarantees to elide the maximum number of redundant saddle-point searches in symmetrically equivalent local-environments. Hence, we are able to study the trapping/detrapping of hydrogen from up to five-vacancy clusters and simultaneously the effect hydrogen has on the diffusivity of these clusters. These processes occur at vastly different timescales at room temperature in body-centred cubic iron. We predict the diffusion pathways of clusters/complexes without a priori assumptions of their mechanisms, not only reproducing previously reported mechanisms but also discovering new ones for larger complexes. We detail the hydrogen-induced changes in the clusters’ diffusion mechanisms and find evidence that, in contrast to mono-vacancies, the introduction of hydrogen to larger clusters can increase their diffusivity. We compare the effective hydrogen diffusivity to Oriani’s classical theory of trapping, finding general agreement and some evidence that hydrogen may not always be in equilibrium with traps, when the traps are mobile. Finally, we compute the trapping atmosphere of meta-stable states surrounding non-point traps, opening new avenues to better understand and predict hydrogen embrittlement in complex alloys.

Saddle point sampling using scaled normal coordinates

The sampling of saddle points (SPs) on a potential energy surface (PES) is critical for describing the dynamics and transport properties of materials. Here, we propose a saddle point search (SPS) method that enables efficient sampling of practically meaningful SPs around a local minimum by converting the PES to the scaled normal coordinates (SNCs), together with the dimer method. We find that the pathway leading to a SP using SNCs is highly efficient and is independent of the system size. This results in the significantly increased SPS efficiency and the enhanced possibility of a complete catalog of the SPs. In addition, we perform SPSs using the SNCs for the diffusion of a vacancy and a dumbbell in body-centered cubic Fe and demonstrate (1) the required number of system force evaluations for a SPS decreases by at least an order of magnitude, (2) the effects of the system size on the number of force evaluations are greatly decreased, and (3) the effects of the number of atoms involved in the SPS on the probability of finding practically meaningful SPs is mostly eliminated. We also discuss the computational cost of introducing the SNCs. It should be highlighted that the SNCs could also be applied to other minimum-mode following methods, demonstrating the general versatility of the proposed method.

Projection Method for Saddle Points of Energy Functional in $$H^{-1}$$ Metric

Saddle points play important roles as the transition states of activated process in gradient systems driven by energy functional. However, for the same energy functional, the saddle points, as well as other stationary points, are different in different metrics such as the \(L^2\) metric and the \(H^{-1}\) metric. The saddle point calculation in \(H^{-1}\) metric is more challenging with much higher computational cost since it involves higher order derivative in space and the inner product calculation needs to solve another Possion equation to get the \(\Delta ^{-1}\) operator. In this paper, we introduce the projection idea to the existing saddle point search methods, gentlest ascent dynamics (GAD) and iterative minimization formulation (IMF), to overcome this numerical challenge due to \(H^{-1}\) metric. Our new method in the \(L^2\) metric can locate the saddle point in \(H^{-1}\) metric only by carefully incorporating a simple linear projection step. We show that our projection method maintains the same convergence speed of the original GAD and IMF, but the new algorithm is much faster than the direct method for \(H^{-1}\) problems. The numerical results of saddle points in the one dimensional Ginzburg-Landau free energy and the two dimensional Landau-Brazovskii free energy in \(H^{-1}\) metric are presented to demonstrate the efficiency of this new method.

Nudged Elastic Band Method for Molecular Reactions Using Energy-Weighted Springs Combined with Eigenvector Following.

The climbing image nudged elastic band method (CI-NEB) is used to identify reaction coordinates and to find saddle points representing transition states of reactions. It can make efficient use of parallel computing as the calculations of the discretization points, the so-called images, can be carried out simultaneously. In typical implementations, the images are distributed evenly along the path by connecting adjacent images with equally stiff springs. However, for systems with a high degree of flexibility, this can lead to poor resolution near the saddle point. By making the spring constants increase with energy, the resolution near the saddle point is improved. To assess the performance of this energy-weighted CI-NEB method, calculations are carried out for a benchmark set of 121 molecular reactions. The performance of the method is analyzed with respect to the input parameters. Energy-weighted springs are found to greatly improve performance and result in successful location of the saddle points in less than a thousand energy and force evaluations on average (about a hundred per image) using the same set of parameter values for all of the reactions. Even better performance is obtained by stopping the calculation before full convergence and complete the saddle point search using an eigenvector following method starting from the location of the climbing image. This combination of methods, referred to as NEB-TS, turns out to be robust and highly efficient as it reduces the average number of energy and force evaluations down to a third, to 305. An efficient and flexible implementation of these methods has been made available in the ORCA software.

Deep Constraint-Based Propagation in Graph Neural Networks.

The popularity of deep learning techniques renewed the interest in neural architectures able to process complex structures that can be represented using graphs, inspired by Graph Neural Networks (GNNs). We focus our attention on the originally proposed GNN model of Scarselli et al. 2009, which encodes the state of the nodes of the graph by means of an iterative diffusion procedure that, during the learning stage, must be computed at every epoch, until the fixed point of a learnable state transition function is reached, propagating the information among the neighbouring nodes. We propose a novel approach to learning in GNNs, based on constrained optimization in the Lagrangian framework. Learning both the transition function and the node states is the outcome of a joint process, in which the state convergence procedure is implicitly expressed by a constraint satisfaction mechanism, avoiding iterative epoch-wise procedures and the network unfolding. Our computational structure searches for saddle points of the Lagrangian in the adjoint space composed of weights, nodes state variables and Lagrange multipliers. This process is further enhanced by multiple layers of constraints that accelerate the diffusion process. An experimental analysis shows that the proposed approach compares favourably with popular models on several benchmarks.

Atomistic modeling of meso-timescale processes with SEAKMC: A perspective and recent developments

On-the-fly kinetic Monte Carlo (kMC) methods have recently garnered significant attentions after successful applications to various atomic-scale problems using a timescale outside the reach of classical molecular dynamics. These methods play a critical role in modeling atomistic meso-timescale processes, and it is therefore essential to further improve their capabilities. Herein, we review one of the on-the-fly kMC methods, Self-Evolving Atomistic kinetic Monte Carlo (SEAKMC) and propose two schemes that considerably enhance the efficiency of saddle point searches (SPSs) during the simulations. The performance of these schemes is tested using the diffusion of point defects in bcc Fe. In addition, we discuss approaches to significantly mitigate limitations of these schemes, which further improves their efficiencies. Importantly, these schemes improve the SPS efficiency not only for SEAKMC but also for other on-the-fly kMC methods, broadening the applications of on-the-fly kMC simulations to complex meso-timescale problems.

A Dual-Dimer method for training physics-constrained neural networks with minimax architecture

Data sparsity is a common issue to train machine learning tools such as neural networks for engineering and scientific applications, where experiments and simulations are expensive. Recently physics-constrained neural networks (PCNNs) were developed to reduce the required amount of training data. However, the weights of different losses from data and physical constraints are adjusted empirically in PCNNs. In this paper, a new physics-constrained neural network with the minimax architecture (PCNN-MM) is proposed so that the weights of different losses can be adjusted systematically. The training of the PCNN-MM is searching the high-order saddle points of the objective function. A novel saddle point search algorithm called Dual-Dimer method is developed. It is demonstrated that the Dual-Dimer method is computationally more efficient than the gradient descent ascent method for nonconvex–nonconcave functions and provides additional eigenvalue information to verify search results. A heat transfer example also shows that the convergence of PCNN-MMs is faster than that of traditional PCNNs.

Examining Saddle Point Searches in the Context of Off-Lattice Kinetic Monte Carlo

Examining Saddle Point Searches in the Context of Off-Lattice Kinetic Monte Carlo

The energy landscape governs ductility in disordered materials.

Based on their structure, non-crystalline phases can fail in a brittle or ductile fashion. However, the nature of the link between structure and propensity for ductility in disordered materials has remained elusive. Here, based on molecular dynamics simulations of colloidal gels and silica glasses, we investigate how the degree of structural disorder affects the fracture of disordered materials. As expected, we observe that structural disorder results in an increase in ductility. By applying the activation-relaxation technique (an open-ended saddle point search algorithm), we demonstrate that the propensity for ductility is controlled by the topography of the energy landscape. Interestingly, we observe a power-law relationship between the particle non-affine displacement upon fracture and the average local energy barrier. This reveals that the dynamics of the particles upon fracture is encoded in the static energy landscape, i.e., before any load is applied. This relationship is shown to apply to several classes of non-crystalline materials (oxide and metallic glasses, amorphous solid, and colloidal gels), which suggests that it may be a generic feature of disordered materials.