"In this tutorial , we introduce the computational framework ARISE (<ins>Ar</ins>tificial <ins>I</ins>ntelligence based <ins>S</ins>tructure <ins>E</ins>valuation) for crystal-structure recognition in single- and polycyrstalline systems[1]. ARISE can treat more than 100 crystal structures, including one-, two-dimensional, and bulk materials - in robust and threshold-independent fashion. The Bayesian-deep-learning model yields not only a classification but also uncertainty estimates, which are principled (i.e., they approximate the uncertainty estimates of a Gaussian process) [2,3]. These uncertainty estimates correlate with crystal order.\n",
"\n",
"\n",
"\n",
"For additional details, please refer to \n",
"\n",
"[1] A. Leitherer, A. Ziletti, and L. M. Ghiringhelli, arXiv:2103.09777, (2021)\n",
"\n",
"ARISE is part of the code framework *ai4materials* (https://github.com/angeloziletti/ai4materials).\n",
"\n",
"\n",
"The outline of this tutorial is as follows:\n",
"\n",
"* Quickstart (jump here if you want to directly use ARISE)\n",
"* Single-crystal classification\n",
"* Polycyrstal classification\n",
"* Unsupervised learning / explanatory analysis of the trained model"
"Please specify the geometry files that you want to analyze in the list 'geometry_files'.\n",
"\n",
"ARISE can be used in two modes: Local and global (controlled via the keyword 'mode'). If mode = 'local', then the strided pattern matching (SPM) framework introduced in [1] is employed, requiring the definition of stride and box size (standard setting: $4 \\overset{\\circ}{\\mathrm {A}}$ stride in all directions and $16 \\overset{\\circ}{\\mathrm {A}}$ box size)."
"Typically, one uses the global mode for single crystals an the local mode for large (polycrystalline) samples (to zoom into a given structure and detect structural defects such as grain boundaries). However, one can also mix this up and investigate, for instance, the global assignments for polycrystals (see for instance the analysis of STEM graphene images with grain boundaries in Fig. 3 of Ref. [1]).\n",
"\n",
"The following sections provide more details on the internal workings of ARISE and SPM. "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## The Bayesian-deep-learning model"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Given an unknown atomic structure, the goal of crystal-strucutre recognition is - in general terms - to find the most similar prototype (that is currently known to be found in nature, for example). In the figure below, this list includes fcc, bcc, diamond, and hcp symmetry but as we will see later, our framework allows to treat a much broader range of materials.\n",
"The initial representation of a given atomic structure (single- or polycrystalline) are atomic positions and chemical species symbols (as well as lattice vectors for periodic systems).\n",
"\n",
"The first step is to define how this input will be treated in the ARISE framework, i.e., how we arrive at a meaningful prediction."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Classification of mono-species single crystals"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Considering exemplarily the body-centered cubic system (see http://www.aflowlib.org/prototype-encyclopedia/A_cI2_229_a.html and also [4]), the prediction pipeline for single crystals is sketched as follows:\n",
""
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Each of these steps will be explained in the following.\n",
"\n",
"We first load the geometry file using ASE (https://wiki.fysik.dtu.dk/ase/index.html). We typically use the FHI-aims format for the geometry files (while via ASE's geometry file parser, compatibility with many other formats is provided). "
"To avoid dependence on the lattice parameters, we isotropically scale a given structure using the function 'get_nn_distance' from the *ai4materials* package. By default, we compute the radial distribution function (as approximated by the histogram of nearest neighbor distances) and then choose the center of the maximally populated bin (i.e., the mode) as the nearest neighbor distance:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"scale_factor = get_nn_distance(structure)\n",
"print('Scale factor for fcc structure = {}'.format(scale_factor))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Given an atomic structure, scaling of atomic positions and unit cell is summarized in the function 'scale_structure': "
"Using atomic positions $\\mathbf{r}_i$ labeled by chemical species $c_i$ ($i = 1, ..., \\text{N}_{\\text{atoms}}$) directly as input to a classification model introduces several issues. Most importantly, physically meaningful invariants that we know to be true are not respected (invariance with respect to translations, rotations and permutations of idential atoms). Also the input dimension would depend on the number of atoms.\n",
"\n",
"A well-known and -tested materials science descriptor that is by definition invariant to the above mentioned symmetries is the so-called smooth-overlap-of-atomic-positions (SOAP) descriptor [5-7]. \n",
"\n",
"In this tutorial (and in [1]) we employ the SOAP implementation that is being made available via the quippy package (https://github.com/libAtoms/QUIP) for non-commercial use.\n",
"\n",
"In short, given an atomic structure with N atoms, we obtain N SOAP vectors (respecting the mentioned invariants) that represent the local atomic environments of each atom. The local atomic environments of an atom is defined as a sphere (centered at that particular atom) that has a certain tunable cutoff radius $R_C$. Each atom within that sphere is represented by a Gaussian function (with tunable width $\\sigma$). The sum of these Gaussians defines the local density:\n",
"Applying this procedure to every atom, gives a collection of SOAP vectors, which one may average. The default behavior in quippy (that is employed in [1]) is to first average the coefficients and then perform the rotational average.\n",
"\n",
"The output of quippy is adapted using a typical averaging procedure (while additional details such as the treatment of \"cross-correlation terms\" is discussed in the next section and also the supplementary material of [1]).\n",
"\n",
"First, we have to define some parameters necessary for the calculation of SOAP, summarized in the string \n",
"descriptor_options that is used to create the SOAP descriptor object:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"cutoff = 4.0\n",
"central_weight = 0.0 # relative weight of central atom in atom density\n",
"The object quippy_SOAP_descriptor above has standard settings that allow to reproduce the results of [1], but in principle one may choose different values - in particular for the SOAP parameters:"
"In [1], ARISE employs a Bayesian neural network as a classification model, which it does not only yield classification probabilities, but also allows to quantify the predictive uncertainty. \n",
"\n",
"Specifically, we employ so-called Monte Carlo (MC) Dropout introduced in [2,3], which employs the stochastic regularization technique dropout. In dropout, individual neurons are randomly dropped in each layer, usually only during training to avoid over-specialization of individual units and this way control overfitting. One can show that employing dropout also at test time yields a probabilistic model whose uncertainty estimates approximate those of a Gaussian process. In particular, neural-network optimization with dropout is equivalent to performing variational inference (see [3] for more details, especially sections 1, 2, and 3.3). \n",
"\n",
"In practice, for a given input, the model output is computed for several iterations (100-1000 typically suffice), yielding a collection of differing probability vectors. In contrast, standard neural networks are deterministic and thus all predictions are identical. Sampling the output layer of the Bayesian neural network actually corresponds to sampling from (an approximated version of) the true probability distribution of outputs. Deterministic networks only provide a point estimate of this distribution, leading to overconfident predictions (see section 1.5 of [3]). Averaging the forward-passes gives an average classification probability vector, and the predicted class label can be inferred by selecting the most likely class (i.e., computing the argmax). Additional statistical information is contained in the classification probability samples. In [1], we employ mutual information (as implemented in *ai4materials*) to obtain a single uncertainty-quantifying number from the output samples. \n",
"\n",
"For the bcc structure from above, prediction and uncertainty quantification is performed in the following cell:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# reshape data array such that it fits model input shape\n",
"print('Maximal classifcation probability with value {:.4f} for prototype {}'.format(prediction[0][argmax_prediction],\n",
" predicted_label))\n",
"plt.plot(prediction[0])\n",
"plt.xlabel('Class (integer encoding)')\n",
"plt.ylabel('Probability')\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The dictionary \"uncertainty\" contains quantifiers of uncertainty, in particular the mutual information employed in [1]:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"print('Mutual information = {:.4f}'.format(uncertainty[\"mutual_information\"][0]))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"***In summary***: All of these steps - isotropic scaling, descriptor calculation, neural network predictions (classification probabilities + uncertainty) - are summarized in the following function, to provide a quick and easy usage. In particular, you may simply pass a list of geometry files, here we choose the fcc (Cu) and bcc (W) prototypes from AFLOW:"
"To analyze the predictions, in particular the ranking by classification probability and the uncertainty quantifiers, we can use the function 'analyze_predictions':"
"For the predictions in the fcc, non-zero probability is also assigned to a tetragonal prototype (identifier 'bctIn_In_A_tI2_139_a', where the AFLOW identifier is specified in this case as well, see [4] for more details). This assignment is justified since it is a slightly distorted (tetragonal) version of fcc (http://aflowlib.org/prototype-encyclopedia/A_tI2_139_a.In.html). Moreover, the uncertainty is non-zero, indicating that there is more than one candidate for the most similar prototype. Interestingly, also the low probability candidates are still meaningful. Further analysis of this classification-probability induced ranking is provided in [1], specifically Fig. 3 (for experimental scanning transmission electron microscopy data of graphene) and Supplementary note 2.3 (for out-of-sample single-crystal structures from AFLOW and NOMAD). This study becomes particularly important if one wants to investigate atomic structures that are not include in the training set. In the mentioned figures of [1], we provide guidelines on how one can still use ARISE in this situation and interpret its predictions. \n",
"\n",
"In the next section, we will first extend to multi-species systems and then provide the code for this out-of-sample studies, where you may also try different structures (upload the geometry files to './data/ARISE/calculations_path'). "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Classification of multi-species single crystals"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"One of the most important characteristics of ARISE is that it can treat much more classes than available methods (see [1] for a benchmarking study).\n",
"\n",
"So far we only considered mono-species systems. We consider the rock salt structure (NaCl) as an example for the extension of ARISE to multi-component systems.\n",
"\n",
"The usual way in treating multiple chemical species is to calculate [8]\n",
"where the coefficients $c_{b_1 lm}^{\\alpha}, c_{b_1 lm}^{\\beta}$ are obtained from separate densities for species $\\alpha,\\beta$. Then one may simply append all of these components, or average over certain combinations. We take the latter route, where it turns out that it suffices (for our purposes) to only select coefficients for $Z_1 = Z_2$ and ignore cross terms for $Z_1 \\neq Z_2$. Then, we consider all substructures, i.e., we compute SOAP according to \n",
"where we sit on all Na atoms and consider only Na atoms, sit on all Na atoms and consider only Ca atoms etc. This procedure yields 4 vectors of equal length, which are averaged. Generalization to structures with arbitrary number of chemical species is straightforward - and so we end up with a descriptor that is independent of 1. the number of atoms, 2. number of chemical species, and 3. incorporates all relevant symmetries (rotational invariance, translational invariance, and invariance with respect to permutations of identical atoms).\n",
"\n",
"Coming back to our implementation, an important keyword is 'scale_elment_sensitive' which is by default set to 'True' - this way, we isotropically scale each substructure separately (one may also not do this, but the model we provide her is trained for 'scale_element_sensitive'==True):"
"One can see that NaCl is correctly predicted with low mutual information."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"You may try to choose a different prototype and analyze its predictions in the following cells. You can choose one of the prototypes of [1] (see also Supplementary Table 4-6 for a complete list), or other examples from AFLOW (including those explored in Supplementary Figure 2).\n",
"\n",
"First we define some functions for loading (no need to understand this in detail, so you can skip for now):"
"* When executing the code below, don't be surprised if the predictions are not identical between different runs. Due to the stochastic nature of the model, this may very well happen, especially for out-of-sample predictions. Still, the final class prediction (via the argmax operation) remains constant for the parameters being used here. \n",
"* To arrive at more constant classification probabilties (and uncertainty quantifications), one may increase the stochastic forward passes (see 'n_iter' argument in the above cell). The increase of computational cost for larger values of n_iter is typicall mild (unless you go to extreme values with billions/trillions of forward passes). \n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Run the following cell to investigate the training prototypes used in [1]. Click 'Run Interact' to start a calculation:"
"You may also want to have a look at the top predicted prototypes. When executing the next cell, an interactive menu will open up where you can change via \"Ranking_idx\" which of the top predictions you want to have a look at and via \"Supercell\" which replicas you want to inspect - exception being nanotubes which have fixed length. Executing the cell after the following one visualizes the selected structure."
"Please run the following cell if you want to study out-of-sample structures (including the prototypes investigated in [1]). If you want to inspect your own, new structures, just upload the to './data/ARISE/calculations_path' and rerun the code cell before 'Some practical tips'). "
"To classify polycrystals, we introduce the strided pattern matching (SPM) framework, sketched in the following figure:\n",
"\n",
"\n",
"\n",
"***The first line*** depicts this procedure for slab-like systems and ***the second line*** for bulk materials. \n",
"\n",
"In both scenarios, a box of certain size is scanned across the whole crystal volume with a certain stride. For each stride, the classification model is applied, yielding both predictions and uncertainty values. These values are then rearranged in 2D or 3D maps. One can construct these for each of the classification probabilities and the uncertainty (here as quantified by mutual information), or also for the most similar prototypes (i.e., the argmax prediction). \n",
"\n",
"\n",
"This way one can discover the most prevalent symmetries (via classification maps) and detect defective regions. In particular, the uncertainty can indicate when a particular region is defective or far outside the known training examples, thus allowing for instance to discover grain boundaries, or, as it is depicted above, a cubic-shaped precipitate."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"As introduced in the 'Quickstart' section, the local analysis can be obtained by simply changing the mode to 'local', while one also has to provide stride and box size (standard settings stride=1.0, box_size=12.0). \n",
"\n",
"We investigate exemplarily the mono-species structure shown in the first line of the above plot (and discussed in more detail in Fig. 2 of [1]):"
"We choose a box size that is roughly equal to the slab thickness and perform the SPM analysis. We choose here a very coarse stride and the smaller we will choose it the better the resolution will get. We select a stride along the z-axis that is larger than the slab thickness, a setting for which no stride in this direction will be made.\n",
"\n",
"A folder (in './data/ARISE/calculations_path') in which all relevant information is saved is created automatically. The name of the folder is given by the structure file being passed and the options for SPM (box size and stride).\n",
"Note that if you run the cell below multiple times, a new folder (with \"_run_i\", i=1,... appended to the file name) will be created, so no data will be lost."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"box_size = [16.0]\n",
"# you may also pass a list for both stride and box size if you have more than one geometry file.\n",
"stride = [[15.0, 15.0, 20.0]] \n",
"# choosing the stride in z larger as the extension of the material (here as 20.0, exceeding \n",
"# the slab thickness by about 4 angstrom) will prohibit that a stride is made in z direction. \n",
"\n",
"previous_level = logger.getEffectiveLevel()\n",
"logger.setLevel(0) # switch off logging, you may change this.\n",
"Within the automatically created folder, the results are saved in the folder 'results_folder' (in this case, the probabilities as 'four_grains_elemental_solid.xyz_probabilities.npy' and the uncertainty as 'four_grains_elemental_solid.xyz_mutual_information.npy'), where each of the numpy arrays has in general the shape (#classes, z_coordinate, y_coordinate, x_coordinate), where here #classes=108 and the coordinates refer to the spatial positions of the box. In this case of a slab structure, the array shape is (#classes, y_coordinate, x_coordinate) while in general it is (#classes, z_coordinate, y_coordinate, x_coordinate)."
"The finer you make the stride, the better the resolution will get. In [1], we employ a stride of $1 \\overset{\\circ}{\\mathrm {A}}$, resulting in the following high-resolution picture, revealing further details at the grain boundary (see also supplementary figure 7 and 8 of [1]):\n",
"One important aspect of [1] is the application of unsupervised learning techniques to analyze the internal neural network representations. This allows to *explain* the trained model by showing that regularities in the neural-network representation space correspond to physical characteristics. Using this strategy, one can conduct exploratory analysis, with the goal to discover patterns in high-dimensional representation space and to relate them to the crystal structure - leading to the discovery of defective regions (e.g., grain boundaries or regions that are distorted in similar/distinct fashion, cf. Fig. 2 and 4 of [1]).\n",
"\n",
"Specifically, we apply clustering (HDBSCAN https://hdbscan.readthedocs.io/en/latest/how_hdbscan_works.html) to identify groups in the high-dimensional neural-network representations. Furthermore the manifold-learning method UMAP (https://umap-learn.readthedocs.io/en/latest/index.html) is employed to obtain a low-dimensional embedding (here 2D) that can be interpreted by human eye. \n",
"\n",
"The most important parameters for the unsupervised techniques employed in this work are the following:\n",
"\n",
"* HDBSCAN: min_cluster_size (int) - determines the minimum number of points a cluster must contain. In more detail, HDBSCAN employs (agglomerative) hierarchical clustering (https://www.displayr.com/what-is-hierarchical-clustering/) to construct a hierarchical cluster tree, where only those clusters are considered that contain as many data points as defined via min_cluster_size. Cluster assignments are extracted by computing a cluster stability value and the most stable clusters are used to construct a final, so-called flat clustering.\n",
"\n",
"\n",
"* UMAP: n_neighbors (int) - determines the trade-off between capturing local (small n_neighbors) vs global (large n_neighbors) relationships in the data. Specifically, UMAP constructs a topological representation of the data (in practice corresponding to a specific type of k-neighbor graph), while some assumptions are made (e.g., the data being uniformly distributed on a Riemannian manifold). Then a low-dimensional representation is found by matching its topological representation with the one of the original data (where the matching procedure amounts to the definition of a cross-entropy loss function and the application of stochastic gradient descent). The parameter n_neighbors influences the high-dimensional topological representation. The distances in the low-dimensional representation can also be tuned via 'min_dist', while this parameter is typically considered more as only tuning the visualization / the distances in the low-dimensional embedding - see also https://pair-code.github.io/understanding-umap/ for a discussion with hands-on examples."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In the cell below, we consider again the mono-species polycrystal with four crystalline grains (the slab structure depicted above and discussed in Fig. 2 of [1], where we chose a stride of 1 $ \\overset{\\circ}{\\mathrm {A}}$ in x,y and a box size of 16 $ \\overset{\\circ}{\\mathrm {A}}$). \n",
"\n",
"You can change the parameters to get a feeling for the importance of min_clsuter_size, n_neighbors. \n",
"\n",
"Also you may choose the color scale according to ARISE's predictions or uncertainty (mutual information), or the HDBSCAN cluster labels. \n",
"\n",
"\n",
"***Remark:*** For the HDBSCAN clustering, there are two options: Either you can choose the 'cluster_labels', which is the cluster assignment that arises from the standard (flat) clustering procedure, or the 'argmax_clustering' color scale, which is determined from the soft-clustering feature of HDBSCAN. Specifically, HDBSCAN allows to calculate a vector for each point, with component i capturing the \"probability\" that the point is a member of cluster i. Then, we can infer a cluster assignment for points that would normally be considered outliers by taking the cluster whose membership probability is maximal (while we consider the point an outlier if none of its cluster probabilities exceeds a certain threshold - here, we choose 10 %). This procedure allows to obtain reasonable clusterings also in the presence of higher levels of noise (as encountered in atomic electron tomography data[9-12], which is analyzed using ARISE in Fig. 4 of [1]). Here, it will not make a qualitative difference in performance, i.e., the basic structure of the polycrystal (four grains + grain boundaries) can be recovered by both clustering procedures."
" cbar.ax.set_yticklabels([int_to_target[_] for _ in range(len(unique_targets))]) \n",
"\n",
"\n",
" ax.axis('off')\n",
" if not color_scale == 'ARISE_predictions':\n",
" ax.set_title('Crystal map (real space)')\n",
"\n",
"\n",
" plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"A consistent picture with four main clusters corresponding to the four crystalline grains separated by grain boundaries arises for min_cluster_size > 250 (and also n_neighbors > 250, while also at lower values reasonable embeddings are obtained, i.e., points corresponding to different grains are separated). For more details on the interpretation of these figures, we refer to [1]."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Conclusion"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"You have reached the end of this tutorial. \n",
"\n",
"*Please let us know if you have any questions, wishes, or suggestions for improvement. Feel free to reach out to us, e.g., via mail (leitherer@fhi-berlin.mpg.de, andreas.leitherer@gmail.com).*"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# References"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"[1] ARISE: A. Leitherer, A. Ziletti, and L. M. Ghiringhelli, arXiv:2103.09777, (2021)\n",
"\n",
"[2] Gal, Y. & Ghahramani, Z. Dropout as a Bayesian approximation: Representing model uncertainty in deep learning. In international conference on machine learning, 1050-1059 (2016)\n",
"\n",
"[3] Gal, Y. Uncertainty in deep learning. Ph.D. thesis, University of Cambridge (2016)\n",
"\n",
"[4] Mehl, M. J. *et al.* The AFLOW library of crystallographic prototypes: part 1. Comput. Mater. Sci. 136, S1–S828 (2017).\n",
"\n",
"[5] A. P. Bartok *et al.* Physical Review Letters 104, 136403 (2010)\n",
"\n",
"[6] A. P. Bartok *et al.* Physical Review B 87, 184115 (2013)\n",
"\n",
"[7] The quippy software is available for non-commercial use from www.libatoms.org (or https://github.com/libAtoms/QUIP)\n",
"\n",
"[8] De, S., Bartók, A. P., Csányi, G. & Ceriotti, M. Comparing molecules and solids across structural and alchemical space. Phys. Chem. Chem. Phys. 18, 13754–13769 (2016)\n",
"\n",
"[9] Chen, C.-C. *et al.* Three-dimensional imaging of dislocations in a nanoparticle at atomic resolution. Nature 496, 74–77 (2013)\n",
"\n",
"[10] Miao, J., Ercius, P. & Billinge, S. J. Atomic electron tomography: 3D structures without crystals. Science 353, aaf2157–aaf2157 (2016).\n",
"\n",
"[11] Yang, Y. *et al.* Deciphering chemical order/disorder and material properties at the single-atom level. Nature 542, 75–79 (2017).\n",
"\n",
"[12] Zhou, J. *et al.* Observing crystal nucleation in four dimensions using atomic electron tomography. Nature 570, 500–503 (2019)."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Please refer to [1] for additional information and references."
