Robust crystal structure identification at extreme conditions using a density-independent spectral descriptor and supervised learning

The effects of temperature must be accounted for when developing a robust crystal structure classifier suitable for materials at extreme conditions. In order to sample finite-temperature configurations for the database, we compute MD trajectories in the NPT ensemble for each material representing the different crystal structures. The simulation cell dimensions correspond to the equilibrium density of each material at 300 K and ambient pressure, and contain 864 atoms. Trajectories are integrated with a timestep of 1 fs, and the coupling parameters for the thermostat and barostat are set to 0.1 and 1.0 ps, respectively. NPT simulations are performed at ambient pressure while ramping up the temperature from 0 to 2/3 of each material’s melting temperature over a 500 ps time window. Configurations are extracted every 5 ps along each trajectory, leading to an ensemble of 100 configurations per crystal structure.

A macroscopic deformation gradient tensor ${\mathbf{F}}$ is applied to the entire system while remapping the $3N$ atoms coordinates $\mathbf{q}={\mathbf{r}}_{1}\oplus\ldots\oplus{\mathbf{r}}_{N}\in\mathbb{R}^{3N}$ (where ${\mathbf{r}}_{i}\in\mathbb{R}^{3}$ are the cartesian coordinates of the $i^{th}$ atom) into the deformed simulation cell. For every configuration extracted from the NPT trajectories:

where $\mathbf{r}_{i,\textrm{deformed}}\in\mathbb{R}^{3}$ stands for the cartesian coordinates of the $i^{th}$ atom subjected to a homogeneous deformation governed by ${\mathbf{F}}$. This deformation gradient tensor reads:

leading to $6$ independent variables describing all homogeneous deformation possibilities in any material. In the present work, we focus on the measure of deviatoric deformation $\epsilon_{\mathrm{eq}}=\sqrt{\frac{3}{2}\mathrm{dev}({\mathbf{E}}):\mathrm{dev}({\mathbf{E}})}$ as a threshold criterion. Here, $\mathrm{dev}{\mathbf{(E)}}={\mathbf{E}}-\frac{1}{3}\mathrm{tr}({\mathbf{E}}){\mathbf{I}}$ with ${\mathbf{E}}$ the Green-Lagrange strain tensor constructed from ${\mathbf{F}}$ with ${\mathbf{E}}=1/2({\mathbf{F}}^{T}{\mathbf{F}}-{\mathbf{I}})$ and ${\mathbf{I}}$ the identity. Only small strains, i.e. within the elastic deformation regime, are considered by imposing a threshold value on $\epsilon_{eq}$, i.e. $\epsilon_{eq}<0.05$. This way, subsets of atomistic configurations subjected to large local deformation such as dislocation-mediated plasticity or amorphous shear banding should emerge as outliers during the classification process.

In this work, we use the bispectrum SO4 [43] to map the local atomic density function into invariant representations. This descriptor has several advantages compared to using the Cartesian coordinates of the surrounding atoms, e.g., it has constant dimension and is invariant to atomic permutation, rotation, and translation. It has been widely used in the context of machine learning interatomic potentials - MLIP, including spectral neighbor analysis potentials (SNAP) [48] and other similar forms [49, 50, 51]. Bispectrum SO4 also has the capacity to provide an accurate description of the atomic neighborhood, suitable for advanced structural analysis [43, 41]. Below, we briefly recall the key concepts and the algebraic formalism used to compute this descriptor. The fully detailed mathematical definition is given in [48]. For all monoatomic systems employed in the present study, the atomic neighbor density around atom $i$ at location ${\mathbf{r}}_{i}$ reads:

where ${\mathbf{r}}_{ii^{\prime}}={\mathbf{r}}_{i}-{\mathbf{r}}_{i^{\prime}}$ is the distance between the central atom $i$ and the neighbor atom $i^{\prime}$, and the cutoff function $f_{c}$ ensures that the contribution of neighboring atoms smoothly decreases to zero at $r_{\mathrm{cut}}$. By mapping radial neighbor coordinates $r$ to an angular component $\theta_{0}=\theta_{0}^{\mathrm{max}}r/r_{\mathrm{cut}}$, the atomic neighbor density can be expanded in the basis functions of the unit 3-sphere, the 4D hyper-spherical harmonics $U^{j}_{m,m^{\prime}}(\theta_{0},\theta,\phi)$:

where the expansion coefficients $u^{j}_{m,m^{\prime}}$ are a sum over discrete values of the corresponding basis function evaluated at each neighbor position,

Finally, using the scalar triple products of these expansion coefficients, the real-value bispectrum components can be expressed as:

with $*$ the complex conjugation operator and where the constants $H\begin{subarray}{c}jmm^{\prime}\\ j_{1}m_{1}m^{\prime}_{1}\\ j_{2}m_{2}m^{\prime}_{2}\end{subarray}$ are the Clebsch-Gordan coefficients for the hyper-spherical harmonics. The final coefficient is invariant to rotation and permutation. The order of the expansion $J_{max}$ determines the accuracy of the geometrical representation of the atomic neighborhood, although bispectrum coefficients are not listed in order of importance. However, increasing the value of $J_{max}$ leads to better accuracy but also to a higher computational cost. In the following, we choose a value of the expansion parameter $J_{max}=4$ that represents a good compromise between the accuracy of geometrical description and computational cost [52, 41]. This leads to a bispectrum $\mathbf{B}$ with $55$ real components, i.e. $\in\mathbb{R}^{55}$, which corresponds to the dimensionality of feature space.

Concerning the cutoff parameter $r_{\mathrm{cut}}$ used to define the neighborhood of a central atom for which the bispectrum is computed, two strategies are proposed. Firstly we consider a fixed cutoff radius, as for the calculation of the potential which is the common procedure in the context of MLIP. Each neighbor in the cutoff sphere is included in the bispectrum calculation, weighted by the cut-off function which smoothly switches from 1 for distances lower than $r_{\mathrm{cut}}$ to exactly zero above $r_{\mathrm{cut}}$. The number of neighbors $N_{\mathrm{neigh}}$ of a central atom is determined by the magnitude of $r_{\mathrm{cut}}$, resulting in a larger number of neighbors for denser materials.

The second alternative is to compute the bispectrum of an atomic environment containing a fixed number of neighbors $N_{\mathrm{neigh}}$. Its main advantage is to remove the dependence of the crystal structure analysis on the density of the material. Hence different materials at varying densities (even locally) could be mapped to an equivalent descriptor representation.

A simple way to achieve the selection of the $N_{\mathrm{neigh}}$ neighbors would be to sort them by their distance to the central atom and consider only the $N_{\mathrm{neigh}}$ first ones while choosing the cutoff radius as the distance between the central atom and its farthest neighbor. More formally, using a basic dichotomy algorithm, one can compute the optimal cutoff radius $r_{\mathrm{cut}}$ that satisfies these requirements by defining the following equality:

where $W_{i}(r_{\mathrm{cut}})$ is the total weight factor associated with the central atom $i$, $H$ is the Heavyside function, $N_{\mathrm{neigh}}$ is the target number of neighbors and $M$ is the actual number of neighbors, related to $r_{\mathrm{cut}}$, required to satisfy $W_{i}(r_{\mathrm{cut}})=N_{\mathrm{neigh}}$. Using the Heavyside weight function systematically leads to the solution with $M=N_{\mathrm{neigh}}$ and $r_{\mathrm{cut}}$ equal to the distance to the $N^{\mathrm{th}}$ neighbor. However, at finite temperature, the BSO4 descriptor strongly depends on local thermal fluctuations since the position (and the identity) of the last neighbor may vary from one step to another. A way to regularize the neighborhood construction by limiting thermal temperature effects is to use a smooth weight function such as tanh:

The $\delta$ parameter ensures the smooth transition for the weights from 1 to 0 and is set to $0.3\text{\,}\mathrm{\SIUnitSymbolAngstrom}$ in the present work, close to thermal fluctuations of atomic positions. In Algorithm 2 we present an algorithm able to compute a general $\mathbf{B}$ descriptor for a constant numbers of neighbours.

In Figure 1, we highlight the distinctions between the standard bispectrum with a fixed cutoff radius, denoted as $\mathbf{B}_{\textrm{cut}}$, and the proposed approach that utilizes a fixed number of neighbors, denoted as $\mathbf{B}^{\mathrm{Heavyside}}_{nn}$ or $\mathbf{B}^{\mathrm{tanh}}_{nn}$ (depending on the weight function employed). We set up NPT simulations of BCC Fe at 100 K, during which the pressure was gradually increased to around 40 GPa (corresponding to a volumetric compression of approximately 20%). The variations of the $\mathbf{B}_{cut}$ and $\mathbf{B}_{nn}$ components with respect to the volumetric ratio $V/V_{0}$ are illustrated in Figure 1-a, b, c, respectively.

The two ways of computing the bispectrum lead to very different results. Indeed, $\mathbf{B}_{nn}$ is almost insensitive to density while the $\mathbf{B}_{cut}$ components evolve significantly with it. Thus, the $\mathbf{B}_{nn}$ should be better at identifying crystal structures even during MD simulations involving a large change in material density. Besides, one can notice a difference between $\mathbf{B}^{\mathrm{Heavyside}}_{nn}$ and $\mathbf{B}^{\mathrm{tanh}}_{nn}$. Indeed, even if it does not depend on local density, $\mathbf{B}^{\mathrm{Heavyside}}_{nn}$ displays some slight evolution that is caused by thermal fluctuations around the sharp step of the Heavyside function. On the contrary, $\mathbf{B}^{\mathrm{tanh}}_{nn}$ appears satisfyingly stable with density, exhibiting minor sensitivity. In the following, we consider BSO4 descriptors computed using either a fixed cutoff radius or a fixed target number of neighbors using the $\tanh$ regularization. The two descriptors will be respectively labelled $\mathbf{B}^{N}_{nn}$ and $\mathbf{B}^{R}_{\mathrm{cut}}$ with $N$ and $R$ the values of the corresponding target number of neighbors or cutoff radius.

As described above, configurations included in the database are selected at different temperatures and after distinct instantaneous deformations. The combination of these two effects makes the local environment of each atom of the supercell unique. Hence, there is no need to apply a peculiar sparsification procedure to avoid redundancy in the database.

For each of the four different crystal structures considered, 100 snapshots containing 864 atoms each are extracted from the NPT trajectory up to 2/3 $T_{m}$. Then, the 6 non-zero components of the deformation gradient tensor ${\mathbf{F}}$ are sampled using the Latin Hypercube Sampling with Multi-Dimensional Uniformity (LHSMDU) [53], in order to obtain 100 draws to be applied to the different snapshots, while ensuring the condition $\epsilon_{eq}<0.05$. Following this procedure, our total database contains $N_{\mathrm{atoms}}$ x $N_{\mathrm{snapshots}}$ = 864 x 100 = 86400 bispectrum vectors $\mathbf{B}\in\mathbb{R}^{55}$, for each crystal structures, giving a total of $M=345600$ $\mathbf{B}$ vectors.

Here we perform an analysis of NPT trajectories in BCC Iron, FCC Aluminium, HCP Zirconium, and c-DIA Silicon, where pressure was maintained at 0 GPa as the temperature is increased up to each material’s melting point $T_{m}$. These simulations are distinct from those of the training database.

We define the accuracy score of a crystal structure classifier algorithm as the number of atoms identified as crystalline over the total number of atoms in the simulation, assuming that the analyzed materials are fully crystalline at $\frac{2}{3}T_{m}$. Figure 4 reports the evolution of the accuracy score as a function of the average temperature in the simulation cell for the four different crystal structures. The SL-CSA classifier clearly outperforms a-CNA, and PTM for BCC Fe, FCC Al, and HCP Zr. The conventional methods shift from 100 % accuracy before reaching $\frac{2}{3}T_{m}$, while the SL-CSA retains more than 98 % accuracy along the whole NPT trajectory. The IDS tool used for c-DIA Si appears not very sensitive to thermal noise and provides comparable results with SL-CSA with an accuracy above 99 % for this structure along the entire trajectory. Finally, the classifier built with $\mathbf{B}^{24}_{\mathrm{nn}}$ appears less sensitive to thermal fluctuations than its counterpart using $\mathbf{B}^{6}_{\mathrm{cut}}$.

In order to examine the performance of our classifiers in extrapolation conditions with changing density (beyond their trained domain at a density of $\rho_{0}$), we generate 100 synthetic samples for each of the four crystal structures at densities $\rho\in[0.5\rho_{0},1.5\rho_{0}]$. These samples correspond to simulation cells at different lattice parameters containing 864 atoms, each atom being shifted from its crystalline position by Gaussian noise with $\sigma=0.1$ Å. This setup allows building artificial configurations far from the domain of validity of the interatomic potential, spanning a wide range of material densities.

Figure 5 compares the results between a-CNA, IDS, PTM, and the two classifiers based on $\mathbf{B}^{6}_{\mathrm{cut}}$ and $\mathbf{B}^{24}_{\mathrm{nn}}$ descriptors. Here, $\rho_{0}$ is taken as the reference density for all measurements and classifiers since we aim at comparing the results w.r.t. the samples present in the database. For the three metallic systems, a-CNA and PTM algorithms provide very good results and correctly predict 100 % of the crystal structures in the entire range of density spanned. The CNA-based classifiers are not adapted for the diamond structure, and we use the IDS classifier instead. Both IDS and PTM predict 100% of c-DIA structure for any density, proving that they are not sensitive to volumetric strain. Thus, the three tested standard classifiers can perform well in the presence of large volumetric deformation, at least in the presence of moderate thermal noise. The $\mathbf{B}^{6}_{\mathrm{cut}}$ based classifier performs well for the densities near $\rho_{0}$. However, due to the fixed cutoff, it largely fails at predicting correct crystal structures when density changes and may even predict a wrong crystal structure, e.g. BCC instead of FCC/HCP. Using the classifier trained on the $\mathbf{B}^{24}_{\mathrm{nn}}$ descriptor allows for the full restitution of the correct crystal structures in BCC, FCC, HCP, and c-IDA systems, over the whole range of densities. Together with the low sensitivity to thermal noise, these results prove its applicability to various thermodynamic conditions, i.e. when large hydrostatic pressure and high temperature are involved.

In order to explore the sensitivity of the classifiers to large deformations, we design a database with NPT trajectories for the four tested structures. For each material, the trajectories are performed at ambient pressure and at $\frac{2}{3}T_{m}$ for 200 ps. Along each trajectory, 200 snapshots are extracted for subsequent application of large deformations. All configurations are different from those of the learning database, and simulations were carried out with different seeds for temperature initialization.

In order to explore independently diagonal and deviatoric deformations, we consider two subsets of structures, where ($F_{11}$, $F_{33}$) and ($F_{12}$, $F_{13}$) are applied. Each component of the deformation tensor $F_{ij}$ is drawn from a uniform distribution in the intervals [0.7, 1.3] and [-0.3, 0.3] for longitudinal and deviatoric strains components respectively.

The 200 couples of diagonal deformation tensor components ($F_{11}$, $F_{33}$) are assigned to the 200 snapshots of each crystal structure, by applying the corresponding deformation gradient tensor ${\mathbf{F}}$ to the simulation cell while rescaling atomic positions. Finally, the different classifiers a-CNA, PTM, $\mathbf{B}_{\mathrm{cut}}^{6}$ and $\mathbf{B}_{\mathrm{nn}}^{24}$ are used to analyse the deformed samples. The same process using deviatoric deformation tensor components ($F_{12}$, $F_{13}$) has been employed. Results for each crystal structure are displayed in Figure 6.

Similar trends are observed for BCC, FCC, HCP, and c-DIA, where the performance of the classifiers is roughly independent of the crystalline structure. For diagonal deformations associated with compression and tension of the simulation cell, the CNA and PTM analysis are limited to small or very small deformation only (lower than 5%), and their accuracy is strongly reduced beyond this point. This effect is likely attributed to the combined effects of deformation and temperature. On the other hand, the performances of the $\mathbf{B}_{\mathrm{cut}}^{6}$ and $\mathbf{B}_{\mathrm{nn}}^{24}$ classifiers remain robust up to deformation as large as 20% (30% in the best cases). We note that $\mathbf{B}_{\mathrm{nn}}^{24}$ even shows an extended range of accuracy compared to $\mathbf{B}_{\mathrm{cut}}^{6}$, and both classifiers exhibit some anisotropic response to diagonal deformation. Concerning deviatoric deformations the performances of our classifiers are even better than CNA and PTM, with correct structural assignment up to 30% deformation. Once again we attribute the high accuracy of our classifiers to their capabilities to handle temperature and deformation effect conjointly. This is highly valuable in the context of in situ classification of large-scale simulations of materials at extreme conditions. In the following, we demonstrate the applications of our most robust classifier constructed with the $\mathbf{B}^{24}_{\mathrm{nn}}$ descriptor for the analysis of structures challenging to perform with traditional methods.