Predicting orientation-dependent plastic susceptibility from static structure in amorphous solids via deep learning

It has been a long-standing materials science challenge to establish structure-property relations in amorphous solids. Here we introduce a rotationally non-invariant local structure representation that enables different predictions for different loading orientations, which is found essential for high-fidelity prediction of the propensity for stress-driven shear transformations. This novel structure representation, when combined with convolutional neural network (CNN), a powerful deep learning algorithm, leads to unprecedented accuracy for identifying atoms with high propensity for shear transformations (i.e., plastic susceptibility), solely from the static structure in both two- and three-dimensional model glasses. The data-driven models trained on samples at one composition and a given processing history are found transferrable to glass samples with different processing histories or at different compositions in the same alloy system. Our analysis of the new structure representation also provides valuable insight into key atomic packing features that influence the local mechanical response and its anisotropy in glasses.


Supplementary Note 1: The optimization of rc and ∆ for the input features to CNN
During the optimization of CNN models, for 2D glasses we also tried rc = 10.0 and ∆ = 0.2 , or rc = 5.0 and ∆ = 0.1 with CNN containing 30 convolutional layers. We found that the accuracy can increase from 93.16% to 94.28% when using rc = 10.0 , ∆ = 0.2 , 30 convolutional layers compared to rc = 5.0 , ∆ = 0.2 , 25 convolutional layers, see Supplementary Fig. 1a. And it is very likely that the accuracy would be higher if there were more training data because the power of deeper CNN models has to be leveraged with larger training data and, as seen from Supplementary Fig. 23a, validation accuracy did not flatten until the training dataset size increases to ~19.6 million when training a CNN model containing 25 convolutional layers. (Note that for 2D glasses, we did not try 30 convolutional layers for fthres = 0.5%, because the accuracy achieved with 25 convolutional layer is already close to the corresponding ceiling and our training dataset size was only 1.96 million when fthres = 0.5%, which will not leverage the power of deeper CNN models). For 3D glasses, we tried rc =12.5 Å and ∆=0.5 Å with CNN containing 25 convolutional layers (see Supplementary Fig. 1b), and the possible reasons why it did not improve accuracy significantly are that the current dataset is not big enough and/or a very small mini-batch size of 32 was used due to limited computer memory. These suggest that the accuracy, especially for 3D glasses, can be further improved if more training data and memory are available. However, doubling our current datasets demands huge computation resources. We project that with ever increasing training dataset size the accuracy may approach the upper bound. An example along this line was ResNet 1 , which contains more than 100 convolutional layers.

Supplementary Note 2: The importance of including medium range structural information
We also tried to reduce rc from its optimized value of 5.0 (9.6 Å) for the 2D (3D) glass, while keeping the same CNN architecture and the same ratio of 2rc/∆. As shown in Supplementary Fig.   1, the validation accuracy decreases gradually. This demonstrates the importance of including surrounding particles beyond the first nearest neighbor shell in predicting plastic susceptibility, and in particular the environment in the medium range on the order of 1 nm.

Supplementary Note 3: The choice of fthres
We chose fthres = 0.5%, because a previous study 2 found that a higher predictive accuracy can be achieved with a higher D 2 min threshold. We confirmed this in Supplementary Fig. 14. The trend for CNN and GNN cannot reflect fully their dependency of accuracy on fthres, because large training datasets (see Supplementary Fig. 23) are needed. Here smaller fthres also means smaller training datasets for the same number of glass samples. If the training dataset size is the same for different fthres, the accuracy at small fthres is expected to be higher.

Supplementary Note 4: Training all species together versus training each species separately
We found that training each species separately could not improve significantly the accuracy achieved for our CNN model,see Supplementary Fig. 15. Although training single species seems more effective, it also means smaller dataset size. In addition, training all species together in a single model is also easier to implement and use. There are some infertile particles mislabeled as having top 0.5% D 2 min (green circles on the maps).
For the 2D glass, the GNN prediction is similar to CNN. In stark contrast, the SVM gives identical prediction for different loading protocols; as a result, too many particles are falsely predicted as belonging to the top 0.5% D 2 min group (green circles on the maps). There are also some fertile sites that are predicted incorrectly (blue circles).

Supplementary Note 6: Rotating and replicating 3D configurations for shearing along desired loading directions
In the LAMMPS package, the configurations with periodic boundary conditions in all three directions have to be cuboids, of which all the edges are parallel to the axes of the coordinate system. Thus, for the loading orientations of #7-24 listed in Supplementary Table 1, after appropriate rotation such that the shear direction and the normal direction of the shear plane are parallel respectively to X and Y axes of the coordinate system, we have to conduct appropriate replication to perform subsequent deformation simulations. Here we give an example to illustrate atoms, we have confirmed that the D 2 min of the replicated atoms is practically the same as that of the original atoms and we select atoms solely from the original atoms to construct training/validating/test datasets.

Supplementary Figures
Supplementary Fig. 1. Effect of cutoff (rc) and grid size (∆) on validation accuracy. (a) 2D glasses at shear strain of 3.0% with fthres of 5.0%. When rc =10.0 or ∆=0.1 , the input images will have 100×100 pixels and thus the corresponding CNN model contains 30 convolutional layers. (b) 3D glasses at shear strain of 5.0% with fthres of 0.5%. When 2rc/∆= 50, the input images will have 50×50×50 pixels and thus the corresponding CNN model contains 25 convolutional layers.

Supplementary Fig. 2. Deformation protocols for 2D model glasses.
Upper panel: four different loading protocols in simulations; for each case the corresponding effective shear orientation is shown in the lower panel, with an appropriate coordinate system set up to compute the SDM. White dash line marks the macroscopic elongation direction. Supplementary Fig. 3. The CNN architecture for 3D glasses. The convolutional neural network (CNN) model for the 3D model glasses contains 20 3D convolutional (conv.) layers. A 3D max-pooling layer (green box) is inserted after every 5 convolutional layers. The last convolutional layer is followed by the output layer which is a sigmoid neuron.
Supplementary Fig. 4. The difference in spatial distribution of surrounding L particles between S particles with the highest and those with the lowest plastic susceptibility. Each map in (a-h) shows the average SDM of the L particles in the neighborhood up to medium range for all S particles that have been identified to have large or small D 2 min when sheared to 3.0% strain, in 5,000 2D samples, for one of the four different loading protocols (from left to right: positive simple shear, negative simple shear, positive pure shear and negative pure shear). The SDM was calculated for these particles in the initial undeformed configurations. (a-d) are for central particles with top 0.5% of D 2 min while (e-h) for those with the lowest 0.5% of D 2 min. (i-l) show the difference in SDMs, which is calculated using the first row (a-d) minus the second row (e-h). The white dashed line in (i-l) represents the macroscopic elongation direction, . (m) shows orientational partial (S-centered) pair correlation function ( , ) for S-L correlation. Each ( , ) curve is an average, over all samples and all loading orientations. For the 3D Cu50Zr50 glass, we show two orientational pair correlation functions for atoms that showed extreme D 2 min upon straining to 5%. (n) and (o) show Cu-centric ( , ) and ( , ) curves, respectively, for Cu-Zr correlation. Supplementary Fig. 5. The difference in spatial distribution of surrounding S particles between L particles with the highest and those with the lowest plastic susceptibility. Each map in (a-h) shows the average SDM of the S particles in the neighborhood up to medium range for all L particles that have been identified to have large or small D 2 min when sheared to 3.0% strain, in 5,000 2D samples, for one of the four different loading protocols (from left to right: positive simple shear, negative simple shear, positive pure shear and negative pure shear). The SDM was calculated for these particles in the initial undeformed configurations. (a-d) are for central particles with top 0.5% of D 2 min while (e-h) for those with the lowest 0.5% of D 2 min. (i-l) show the difference in SDMs, which is calculated using the first row (a-d) minus the second row (e-h). The white dashed line in (i-l) represents the macroscopic elongation direction, . (m) shows orientational partial (L-centered) pair correlation function ( , ) for L-S correlation. Each ( , ) curve is an average, over all samples and all loading orientations. For the 3D Cu50Zr50 glass, we show two orientational pair correlation functions for atoms that showed extreme D 2 min upon straining to 5%. (n) and (o) show Zr-centric ( , ) and ( , ) curves, respectively, for Zr-Cu correlation. Supplementary Fig. 6. The difference in spatial distribution of surrounding L particles between L particles with the highest and those with the lowest plastic susceptibility. Each map in (a-h) shows the average SDM of the L particles in the neighborhood up to medium range for all L particles that have been identified to have large or small D 2 min when sheared to 3.0% strain, in 5,000 2D samples, for one of the four different loading protocols (from left to right: positive simple shear, negative simple shear, positive pure shear and negative pure shear). The SDM was calculated for these particles in the initial undeformed configurations. (a-d) are for central particles with top 0.5% of D 2 min while (e-h) for those with the lowest 0.5% of D 2 min. (i-l) show the difference in SDMs, which is calculated using the first row (a-d) minus the second row (e-h). The white dashed line in (i-l) represents the macroscopic elongation direction, . (m) shows orientational partial (L-centered) pair correlation function ( , ) for L-L correlation. Each ( , ) curve is an average, over all samples and all loading orientations. The distribution of L neighbors changes with mainly in the peak intensity, also indicating the asymmetric packing. For the 3D Cu50Zr50 glass, we show two orientational pair correlation functions for atoms that showed extreme D 2 min upon straining to 5%. (n) and (o) show Zr-centric ( , ) and ( , ) curves, respectively, for Zr-Zr correlation. Supplementary Fig. 7. Schematic illustrating the critical structural differences that influence the local mechanical response and its anisotropy. This insight is first revealed from the difference in SDMs between particles with the highest versus those with the lowest plastic susceptibility in 2D glasses and then further confirmed for both 2D and 3D glasses through contrasting orientational pair correlation functions, see Fig. 2 and Supplementary Figs. 4-6. Specifically, for particles with low plastic susceptibility, their neighbors tend to form circularly or spherically symmetric shells with almost uniform packing intensity, as schematically idealized using the blue circles in the schematic above. In contrast, packing tends to be more asymmetric surrounding fertile sites: the neighboring shells are more like ellipses, shown as the red elliptical circles, each with non-uniform particle distribution. For a fertile site to be activated upon a particular mechanical loading, the longest and shortest axes of neighboring shells/ellipses should be parallel to the elongation (the red dotted line) and contraction (the black dotted line) directions imposed by the loading, respectively. In other words, a fertile site is the most prone to be activated when its surrounding packing along the elongation direction is the loosest and the packing along the contraction direction is the densest.

Elongation direction Contraction direction
Supplementary Fig. 8. The orientational partial pair correlation function ( , ) of S particles that showed the top 0.5% D 2 min after being shear strained to 3.0% (colored), compared with that for the lowest 0.5% D 2 min (grey). The upper row is for the neighboring S particles, whereas the lower row is for L particles in the neighborhood. Each curve is the average over all surrounding S (or L) particles from all 5,000 samples. From left to right the columns correspond to positive simple shear, negative simple shear, positive pure shear and negative pure shear, respectively. The (a-d) are similar to each other and are averaged into Fig.  2m in the main text. And (e-h) are similar to each other and are averaged into Supplementary Fig. 4m.  Supplementary Fig. 9. The orientational partial pair correlation function ( , ) of L particles that showed the top 0.5% D 2 min after being shear strained to 3.0% (colored), compared with that for the lowest 0.5% D 2 min (grey). The upper row is for the neighboring S particles, whereas the lower row is for L particles in the neighborhood. Each curve is the average over all surrounding S (or L) particles from all 5,000 samples.  (c) are mirrored to each other along macroscopic contraction direction (green dash line); (a) and (d) will be the same after two mirror operations (along both elongation and contraction lines), which is equivalent to a rotation by 180° around the center of images. The 4 images can generate another 4 mirrored ones, across the plane of × for 3D configurations. These mirror operations will not change the packing density of local configurations along elongation or contraction orientations and thus those mirrored configurations are expected to exhibit similar mechanical response, which is confirmed in the Supplementary Fig. 11.  Supplementary Fig. 11. Mirrored samples show almost identical mechanical response. (a) and (b) are for a 2D glass strained to 2.5% in simple and pure shear, respectively. (c) is for 3D glasses strained to 5.0% in simple shear. Each symbol in the plots represents D 2 min of a single particle in one the of 3 (or 7) mirrored samples (see more details about mirror-symmetry in the Supplementary Fig. 10) against that in the original 2D (or 3D) glass. Almost all symbols locate on the diagonal dash lines, indicating that particles in mirrored samples show almost identical mechanical response as those in the original samples. Supplementary Fig. 12. (a) The shear stress-strain curves for 2D glasses, averaged over 5,000 post-loading samples for each of the four different protocols. For simple shear, the shear strain = 4 "# 4, the shear stress = 4 "# 4, and for pure shear, = 4 "" − ## 4, = $% !! &% "" $ ' . (b) The shear stress-strain curve for Cu50Zr50 3D glasses, averaged over 20 samples sheared to large (up to 20%) strain. Each sample was simple-sheared in 24 different orientations and thus the curve is the average of 480 curves. Supplementary Fig. 13. The cumulative distribution function (CDF) of D 2 min at different shear strains for the 2D glasses (a) and the 3D Cu50Zr50 glasses (b). Supplementary Fig. 15. Effect on validation accuracy, comparing training each species separately versus training all species together. (a) 2D glasses at shear strain of 3.0% with fthres of 5.0%. (b) 3D Cu50Zr50 glasses at shear strain of 5.0% with fthres of 0.5%. Supplementary Fig. 16. Spatial distribution map of prediction results for a given 2D glass sample. The particles with top 0.5% D 2 min are regarded as fertile sites after the model was loaded to shear strain of 3.0%. The results on the upper, middle and lower panels are from CNN, GNN, and SVM models, respectively. From left to right the columns represent four different loading protocols, i.e., positive simple shear, negative simple shear, positive pure shear and negative pure shear, all performed on the same glass. On the maps, fertile sites predicted correctly are marked with red color; non-fertile-sites predicted correctly are marked silver while those mislabeled as fertile are marked light green; fertile sites mislabeled as nonfertile are marked blue.

Supplementary Fig. 17. Comparing fully connected neural networks (NN) with CNN, GNN and SVM.
This plot shows the validation accuracy (symbols) reached, on 2D glasses at strain 3.0% with fthres = 5.0%, via the fully connected NNs with different architectures (each NN contains 2 or 4 hidden layers with different number of neurons and the number of neurons in different hidden layers is identical within a NN). For comparison, the accuracy achieved by using other machine learning methods is represented by lines with different colors (see legend). The input to these simple NNs is the same as the images fed to CNN, but flattened into vectors. Supplementary Fig. 19. Shear modulus (G) versus the CNN-predicted fraction of particles with class probability > 0.5. This fraction is denoted fpy=1 and used as the x axis, for (a) 2D L-J and (b) 3D Cu50Zr50 glass samples with different processing histories. Each data point is averaged over 50 (10) samples across 4 (24) loading protocols for the 2D (3D) glasses, sheared to a strain of 1% (2%). The dashed line serves as a guide to the eye. Supplementary Fig. 21. Accuracy achieved via the GNN model, when the hyperparameters are varied, for (a) 2D glasses at strain of 3.0% with fthres of 5.0% and (b) 3D Cu50Zr50 glasses at strain of 5.0% with fthres of 0.5%. The edge threshold is varied at a fixed number of unrolls of the network (nrec), and the nrec is varied at a fixed edge threshold.