Automatic identification of slip pathways in ductile inorganic materials by combining the active learning strategy and NEB method

Introduction

Recently discovered ductile inorganic non-metallic materials have attracted widespread attention due to their significant potential in flexible electronics applications^1,2,3,4,5. Despite the promise of these materials, the number of identified ductile inorganic semiconductors remains limited, with only a handful having been experimentally and theoretically characterized^6,7,8,9,10. The scarcity of ductile inorganic semiconductors has drawn significant attention to the need for systematic exploration and discovery of new inorganic materials with desirable mechanical properties, where high-throughput computational screening has emerged as a powerful strategy¹¹. By high-throughput screening with first-principles calculations, one recent attempt by Gao et al. has been successful in identifying several potential candidates with plastic deformability from binary 2D van der Waals crystals¹². In the high-throughput screening process, the most critical aspect, especially when it comes to understanding and predicting ductile behavior in inorganic materials, is the calculation of slip energy^6,13. The key to calculating slip energy lies in identifying the optimal slip pathway. In Gao’s work, the minimum energy slip pathway and the corresponding energy barrier were determined manually, which becomes inefficient for inorganic materials with complex slip pathways. Overall, to date, there is still a lack of an automated workflow for the determination of slip pathways or slip systems in complex ductile inorganic non-metallic materials.

Slip is the fundamental mechanism of plastic deformation in crystalline materials. When a material undergoes mechanical stress, atoms in the crystal lattice move relative to one another along specific crystallographic planes, known as slip planes. The combination of a slip direction and plane is referred to as a “slip system”^14,15. The energy required to drive this movement is termed the generalized stacking fault energy (GSFE)¹⁶, or simply, the slip energy. Calculating the slip energy provides insight into how easily a material deforms plastically and, by extension, its ductility. For brittle materials, the energy barriers are often too high to allow significant slip, leading to fracture instead of plastic deformation. Conversely, in ductile materials, the slip pathways offer low energy barriers, allowing atoms to slide past each other more readily under stress. Usually, the more slip systems a material has, the easier it undergoes slip, leading to better ductility¹⁷.

In metals, due to the relatively simple arrangement of atoms, the potential slip systems have been identified¹⁸. Digital Image Correlation (DIC) combined with Scanning Electron Microscope (SEM) technology has been successfully applied in the automated identification of activated slip systems during the process of plastic deformation^{19,20,21,22,23}. However, the complex elemental compositions and atomic arrangements complicate the slip pathways in inorganic materials, rendering this method ineffective. Compared to traditional experimental methods, energy-based approaches offer an alternative and effective means of determining slip systems²⁴. The grown density functional theory (DFT)^25,26 has greatly facilitated the construction of slip surfaces^27,28,29. Furthermore, the advent of machine learning (ML) has reduced the costs and accelerated the speed of the process^30,31,32,33. Studies have shown that slip systems typically represent minimum energy pathways (MEPs) upon the slip surface³⁴. Finding MEPs is a common and important problem in theoretical chemistry. Many different methods have been presented to find MEPs^35,36. Among these, the nudged elastic band (NEB) method, which has been extensively employed, is robust, efficient, and easily parallelizable^37,38. Hence, the combination of ML and NEB offers a potential solution to the automatic identification of the slip pathways in complex inorganic materials.

In this paper, for complex inorganic material systems, we propose a workflow based on interlayer slip potential energy surfaces to automatically and rapidly identify potential slip pathways. Active learning strategies and the climbing image nudged elastic band (CI-NEB) method are utilized in the workflow. Our model and workflow have been tested in several archetypical materials. The results confirm the reliability of our workflow. Our work solves the problem of identifying slip pathways in complex inorganic materials and offers valuable support for the further high-throughput screening of ductile inorganic materials.

Results

Workflow

Figure 1 shows the schematic of our workflow, which is based on the implementation of “interlayer slip potential energy” (ISPE). While similar concepts have been reported in other literature^39,40,41,42, our work uniquely incorporates ISPE. In a perfect crystal, the rigid slip of one group of atoms relative to another leads to variations in the system’s energy. We define ISPE as the energy computed from the post-slip structural configuration. Our methodology based on ISPE involves two modules. Module 1 implements an active learning strategy that effectively samples the slip plane and learns the ISPE surfaces. Module 2 leverages the trained ISPE surfaces to automatically detect potential slip pathways using the CI-NEB method, ultimately selecting the final pathway based on the energy profile. The following sections will provide a comprehensive explanation of our process.

Automatic identification of slip pathways in ductile inorganic materials by combining the active learning strategy and NEB method — **Fig. 1: The schematic of the workflow for the automated identification of slip pathways.**

In Module 1, the initial structure must be a layered configuration capable of slip. Considering the definition of the slip system, we constrain the slipping to occur within a single plane, disregarding any slipping perpendicular to the plane. To simplify subsequent calculations, we fix the slip plane within the ab plane of the unit cell, with the c-axis oriented perpendicular to the ab plane. Additionally, to prevent atoms from becoming too close during the slip process, which could lead to unrealistic structures, we perform an initial check. If any slip configurations result in interlayer atomic distances that are less than 0.5 Å, we will expand the interlayer distance by inserting a vacuum layer with a length of around 1.0 Å between the layers before performing the slip process. The rationalization for this operation is discussed in detail in the Supplementary Materials.

Once the initial slip structure is prepared, we uniformly sample the slip plane along the ab axes. To obtain a sufficiently fine grid of sampling points, we divide the slip plane into 100 × 100 sampling points, partitioning it into 99 segments along both the a and b axes. A subset of these sampling points is then extracted for DFT static calculations using a coarser grid, such as 5 × 5, to generate the initial training set. Based on this initial training set, we train a Gaussian process regression (GPR) model to produce a preliminary model.

Using this preliminary model, we predict the ISPE at the remaining sampling points, excluding those in the training set, and calculate the uncertainty at each point. The points with the highest uncertainties are ranked, and the 6 points with the greatest uncertainty are selected for further DFT static calculations to determine their energies. These points form the test set, and the model is evaluated using the r² score. If the score is below the predefined threshold (such as 0.999), the test points are added to the training set, the model is retrained, and new high-uncertainty points are selected for the next iteration. This process is repeated until the model converges. During this process, we actively learn the ISPE surfaces based on uncertainty, aiming to achieve a high-quality model using minimal computational resources.

In Module 2, based on the trained ISPE surface model, our first task is to locate the global minimum energy point on the ISPE surface. In most cases, the initial point of the slip process corresponds to a local energy minimum on the ISPE surface, but it is not necessarily the global minimum. This depends on the specific system. For example, in the case of InSe with the space group P6₃/mmc, we found that the global energy minimum on its ISPE surface does not coincide with the initial position. A recent study has reported similar findings⁴³. To address this issue, we have implemented a two-step procedure in our code to accurately locate the global energy minimum. This involves first using a trained model to roughly identify an initial point, followed by optimization from this point to ensure the global energy minimum is found.

The initial state and final state of the slip process is then fixed at the global minimum energy points of the neighboring unit cells. We then explore possible slip pathways in different directions. Specifically, we search within neighboring unit cells along both the a-axis and b-axis. By using linear interpolation, two initial pathways are generated. These two pathways are then refined using the CI-NEB method to identify the two MEPs.

The NEB is a chain-of-state method. The process of projecting the force components perpendicular to the path while allowing the spring force to act only along the path is known as “nudging”³⁷. The force on the ith image is

$$overrightarrow{{F}_{i}}=-vec{nabla }Vleft(overrightarrow{{R}_{i}}right)|_{perp}+overrightarrow{{F}_{i}^{s}}cdot widehat{{tau }_{{{||}}}}widehat{{tau }_{{{||}}}}$$

Where (overrightarrow{{R}_{i}}) denotes as the position vector of the ith image in the scalar field (V), (hat{{tau }_{{left||right.}}}) is the unit tangent vector to the path and (vec{nabla }Vleft(overrightarrow{{R}_{i}}right)|_{perp }=vec{nabla }Vleft(overrightarrow{{R}_{i}}right)-vec{nabla}Vleft(overrightarrow{{R}_{i}}right)cdot widehat{{tau }_{{||}}}widehat{{tau }_{{||}}}). This decoupling of the path relaxation from its discrete representation is critical for ensuring convergence to the MEP⁴⁴. The CI-NEB method constitutes a small modification to the NEB method⁴⁵. Specifically, after a few iterations with the regular NEB, the full force of the image with the highest energy is expressed as:

$$overrightarrow{{F}_{{i}_{max}}}=-vec{nabla }Vleft(overrightarrow{{R}_{{i}_{max }}}right)+2vec{nabla }V(overrightarrow{{R}_{{i}_{max }}})cdot widehat{{tau }_{{{||}}}}widehat{{tau }_{{{||}}}}$$

The subscript ({i}_{max }) refers to the index of the image along the path that corresponds to the highest energy. In this way, information about the shape of the MEP is retained, but a rigorous convergence to a saddle point is also obtained. This additional feature does not add any significant computational effort.

Since we are using the CI-NEB method, the choice of spring constants and the number of images is not particularly stringent. To further ensure the global minimum energy pathway, we compare the energy barriers along each of the two MEPs and select the path with the lower energy barrier as the final slip pathway. The structure corresponding to the energy barrier along this pathway is used to compute the maximum slip energy of this material.

Model training

Feature vector

For the active learn algorithms to learn the ISPE related to their corresponding slip structures, a numerical representation of the slip structures is required, as opposed to the character strings we normally use to identify different slip structures. Since our focus is on ISPE, we propose using the fractional slip vector (FSV) of the crystal as the preferred numerical representation of the slip structures. The concept of FSV has been widely employed in visualizing the stacking energy surface of van der Waals layers. It bears resemblance to the well-established concept of the “Burgers vector”^46,47,48,49. However, compared to “Burgers vector”, FSV focuses on the interlayer slip behaviors within a single unit cell, and is derived from the lattice vector, denoted as (vec{a}), (vec{b}) and (vec{c}). For example, the FSV of (f_a, f_b, f_c) means that a portion of the atoms in the unit cell slips a distance of ({f}_{a}cdot left|vec{a}right|) along the a-axis, ({f}_{b}cdot |vec{b}|) along the b-axis, and ({f}_{c}cdot left|vec{c}right|) along the c-axis, relative to another portion of atoms. Since the slip process is constrained to the ab plane, a two-dimensional FSV of (f_a, f_b) is employed as the feature vector. However, in specific scenarios where the energy range of the ISPE surfaces exceeds 40 eV, we have observed that the FSV becomes ineffective. Therefore, we have also developed two additional feature vectors to address these exceptional cases. See the discussion and Fig. S1 in Supplementary Materials for details.

ML model

After extensive testing and a thorough literature review^50,51, we find that the GPR model, with the advantages of uncertainty quantification, flexibility for small datasets, smoothness, and robustness^52,53, is the most suitable for the active learning of ISPE. However, the performance of the GPR model is highly correlated with the kernel functions. In this study, we have tested the performance of the model using the Rational Quadratic (RQ) kernel⁵⁴, Pairwise (PW) kernel, and Matérn (MT) kernel⁵⁵.

Figure 2 shows our results. DFT static calculations are initially performed in a dense 51 × 51 grid across the entire slip surface to obtain ISPEs for the test set. Subsequently, a uniform 11 × 11 subset from this grid is extracted to train the model. Model testing is conducted on the Ag₂S, InSe, and CrAuTe₄ systems. The energy span for the Ag₂S system is approximately 20 eV, about 0.5 eV for InSe, and around 5 eV for CuAuTe₄. The results indicate that, under the same kernel function, as the energy span of the system increases, the model’s performance deteriorates, with a gradual decrease in the r² score and an increase in the root mean squared error (RMSE).

**Fig. 2: Predicted ISPEs within the GPR model using different kernel functions.**

Furthermore, we evaluated three different kernel functions for each material system. The results demonstrate that, regardless of the system, the MT kernel consistently performs the best. The RQ kernel exhibits slightly lower performance compared to the MT kernel but remains competitive, while the PW kernel shows the lowest performance. Overall, the MT and RQ kernels perform at a comparable level. Notably, the parameter ν in the MT kernel is crucial; our detailed test reveals that a ν value of approximately 2.5 optimizes model performance. In summary, for various tested systems, the selected feature vectors and kernel functions consistently yield model scores exceeding 0.99, demonstrating the reliability and effectiveness of our model. We have discussed the influence of these three kernel functions on CI-NEB results in the supplementary materials.

Case studies

To validate our workflow, we first apply it to the InSe system with a relatively simple slip pathway. The algorithm for Module 1 is deployed on the AiiDA platform^56,57, where calculations are automatically submitted and managed through AiiDA for model training. Thanks to the integration of active learning algorithms, we are able to efficiently select sampling points across the ISPE surface, resulting in rapid convergence of the calculations. Subsequently, we test the model using the pre-calculated 51 × 51 grid points. The results show an r² of 0.999998 and an RMSE of 9.4 × 10⁻⁵ eV, indicating that the final trained model performed exceptionally well.

Building on the strategy from Module 2, we then conduct CI-NEB simulations based on the learned ISPE surface. Figure 3 presents our results. In this case, the initial pathways are generated by evenly interpolating a straight line between the two lowest energy points in adjacent cells along both the a-axis and b-axis directions. Figure 3a, c shows the initial pathways along the a-axis and b-axis. Figure 3b, d shows the identified pathways after CI-NEB running, which locate MEPs correctly. Figure 3e shows the energy of each image along the MEPs. The resulting energy profiles indicate that the energy barriers along both pathways are identical, confirming that the two pathways are completely equivalent. From these curves, we can then identify the location of the energy barrier during the slip process, which corresponds to the 6th image along the pathway of the a-axis, representing the FSV at (0.844, 0.156).

**Fig. 3: Identification of slip pathways for InSe.**

Furthermore, we have tested our workflow on the CrAuTe₄ system. Compared to InSe, CrAuTe₄ has a more complex interlayer structure, leading to a more intricate slip plane, and its ISPE surface operates on a larger energy span than that of InSe. However, thanks to the efficient active learning strategy, the calculations converged after approximately 20 iterations. We then tested the trained model using the previously computed 51 × 51 grid points, obtaining an r² of 0.9972 and an RMSE of 0.050 eV, indicating a reasonably high accuracy of the model.

Subsequently, CI-NEB simulations are performed utilizing the learned ISPE surface model, with search directions along both the a-axis and b-axis. Figure 4 displays our findings. Figure 4a, c illustrates the initial pathway configurations analogous to the InSe system. Figure 4b, d illustrates the slip pathways, or MEPs, obtained from the simulations, corresponding to the a-axis and b-axis, respectively. We can observe that, in contrast to the InSe system, the slip pathways in this system are notably more complex, no longer following a simple linear trajectory. Figure 4e presents the energy associated with each image along the slip pathways, highlighting a notable difference in energy barriers. The pathway along the a-axis exhibits a considerably higher energy barrier than that along the b-axis. Consequently, we identify the MEP along the b-axis as the optimal slip pathway, with the 13th image located at the position of maximum slip energy, corresponding to the FSV at (0.928, 0.496).

**Fig. 4: Identification of slip pathways for complex CrAuTe₄.**

It is worth noting that, due to the complexity of the slip plane, Fig. 4b indicates that a potential slip pathway may have been missed, located on the opposite side of the identified pathway. To resolve this, additional initial pathways should be set to ensure that all possible slip pathways are captured. See Fig. S2 in supplementary materials for our test results.

Discussion

In summary, we have developed a robust automated workflow for the identification of slip pathways in ductile inorganic materials, successfully integrating an active learning strategy with the CI-NEB method. Our workflow comprises two key modules: the first leverages uncertainty-based sampling to construct the ISPEs within the unit cell, while the second employs the ISPE surfaces to perform the slip pathway identification using the CI-NEB method. This innovative approach allows for the rapid and efficient identification of MEPs and facilitates the determination of slip structures associated with varying energy barriers.

The performance of our machine learning descriptors and models was rigorously evaluated across several material systems, including Ag₂S, InSe, and CrAuTe₄. Notably, the ML model demonstrated exceptional sensitivity to the energy span of the ISPE surfaces, maintaining a high degree of accuracy (r² > 0.99) for spans within 20 eV. Our kernel function analysis further reinforced the reliability of our approach, with both the RQ and MT kernels yielding satisfactory fitting results.

The comprehensive testing of our computational workflow has been conducted on InSe and CrAuTe₄ systems with simple and complicated slip pathways, respectively, confirming its efficacy and reliability in identifying slip pathways. The outcomes highlight the potential of our methodology in advancing the understanding of slip mechanisms in ductile inorganic materials. This work paves the way for future high-throughput screening research to explore new inorganic non-metal materials with plastic deformability.

Methods

ML model test

GPR is employed as our ML model. It is a non-parametric model that utilizes a Gaussian Process (GP) prior for regression analysis of data. A GP is characterized as a collection of random variables, where any finite subset exhibits a multivariate normal distribution. In GPR, we initially define a prior GP that encapsulates our initial comprehension of the shape of the function. Upon observing the data, we update this prior through Bayesian inference to obtain a posterior distribution, which reflects our comprehension of the function given the known data. Assuming we have a set of training data (X, y), where X represents the feature vector and y denotes the corresponding target variables, follows GP:

$$y,{{sim }},{mathcal{G}}{mathcal{P}}(mu left(Xright),sigma (X,,X^{prime} ))$$

Where (mu (X)) represents the mean function, while (sigma (X,{X}^{prime} )) denotes the kernel function (or covariance function). Once the distribution is established, to predict new input X^*, we can obtain the posterior distribution through Bayesian inference:

$${y}^{* }left|{X}^{* },X,y,{{sim }}{mathscr{N}}({mu }^{* },{Sigma }^{* })right.$$

Where ({mu }^{* }) represents the predicted mean value, and ({Sigma }^{* }) denotes the predicted variance value. The predicted ({mu }^{* }) provides an estimate for the new input data, while the ({Sigma }^{* }) quantifies the uncertainty associated with this prediction. A larger ({Sigma }^{* }) indicates lower confidence in the prediction, whereas a smaller variance signifies higher confidence. Our active learning process is precisely based on this uncertainty.

For the construction of the ML model, two-dimensional FSV was selected as X, and the ISPE was selected as y. GPR algorithm provided by scikit-learn package⁵⁸ was used. Before training, the built-in regularization method of the algorithm was applied to regularize the ISPEs.

To evaluate the kernels, we computed the ISPE surfaces within 51 × 51 grid points for the Ag₂S, InSe, and CrAuTe₄ systems, respectively. For Ag₂S, we selected a two-layer zigzag slab structure containing 24 atoms, the construction of which has been detailed in our previous work¹³. Given that, in certain slip configurations, the interlayer atomic distances could become extremely short, causing a sharp increase in system energy and model failure, we inserted a vacuum layer with a length of 1.0 Å between the two atomic layers and used this modified structure as our initial configuration. To obtain slip structures for the calculations of ISPE, atoms in the bottom layer were fixed, while atoms in the top layer displaced along the a-axis and b-axis in a certain displacement within the unit cell. We uniformly divided the displacements along the a-axis and b-axis of the unit cell into 50 equal segments, thereby constructing a 51 × 51 grid. Pymatgen package⁵⁹ was used to generate the slip structures for each grid point. For InSe and CrAuTe₄, a two-layer slab structure was constructed based on the structure from The Materials Project (mp-20485 and mp-12743)⁶⁰. The way to construct slip structures was the same as Ag₂S.

For these slip structures, DFT static calculations were performed to obtain sufficient data, which forms our dataset. Vienna Ab initio Simulation Package (VASP) with a plane-wave basis set and projector-augmented wave (PAW)^61,62 pseudopotentials were used. For the three systems, the exchange-correlation functional was approximated by the Perdew–Burke–Ernzerhof (PBE) formulation of the generalized gradient approximation (GGA)⁶³. In these calculations, van der Waals corrections were considered using the DFT-D3 method with Becke-Johnson damping function⁶⁴. For Ag₂S, InSe, and CrAuTe₄, the plane-wave cutoff energy used in the calculations for all was 400 eV, and a gamma-centered k-point with the mesh of 3 × 3 × 1, 11 × 11 × 1 and 9 × 5 × 1 was utilized, respectively.

The calculated 51 × 51 ISPEs formed our test set. We uniformly extracted 11 × 11 ISPEs from the test set to create the training set for GPR model training. We initially experimented with a variety of kernel functions; however, only three of them provided relatively reasonable results. Therefore, we focused our testing on these three models. They were the Rational Quadratic kernel, Pairwise kernel, and Matérn kernel.

Setups for the case studies

For InSe, we utilized the AiiDA platform^56,57 in combination with an active learning strategy to train the ISPE surface. The entire slip plane was divided into a 100 × 100 grid. A sparse uniform 5 × 5 grid of points was chosen for DFT calculations, forming the initial training set of the model. During the active learning process, the 6 grid points with the highest uncertainty were selected for each new round of DFT calculations. Throughout the iterations, we set a convergence threshold of r² = 0.999. For the CI-NEB simulation, the method built in ASE⁶⁵ was employed. The k of springs was set 1.0 eV/Å², and 21 images were used to identify MEP. The MDMin algorithm built in the software was used to optimize the pathways. In the optimization process, we set the maximum number of optimization steps to 3000. The optimization was terminated when the maximum force on these images fell below 0.01 eV/Å, at which point the results were output. For CrAuTe₄, the parameters remained the same as InSe except for the use of 25 images.