Crystal structure map for materials classification and modeling

ABSTRACT

For classifying and modeling properties of crystalline materials in terms of structure, a three-step workflow with (1) generation of structure feature vectors, (2) evaluation of distances among the feature vectors as a measure of similarity in structure, and (3) mapping of each structure in a low-dimensional space with principal components using dimension reduction is proposed. The obtained distance and resulting principal components are useful for classifying similar crystal structures for a given set of materials systems and for constructing descriptors in machine-learning analyses of properties. The eigenvector of the principal components indicates which part of the original structure feature vector is contained as important information. Examples are demonstrated for classification and property modeling of Al${}_2$O${}_3$ polymorphs including amorphous structures and of the alloy configurations of Si-doped LaFe${}_{13}$ compounds.

GRAPHICAL ABSTRACT

KEYWORDS:

• Structure classification
• Multidimensional scaling
• DFT calculations

Disclaimer

As a service to authors and researchers we are providing this version of an accepted manuscript (AM). Copyediting, typesetting, and review of the resulting proofs will be undertaken on this manuscript before final publication of the Version of Record (VoR). During production and pre-press, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal relate to these versions also.

1. Introduction

Today, a wealth of materials information is being generated, collected, and accumulated, experimentally, theoretically, and computationally. One of the current hot issues in materials research is materials discovery assisted by data-science approach with the information [1]. For materials exploration and design, it is important to know in advance where each existing materials system is positioned in a map and more basically how the coordinate system of the map is defined for a given set of materials group. This raises a novel concept “materials map” [2]. Since the materials space is super high dimensional, the key to the materials map for a concrete target materials system is to extract appropriate features which define the low-dimensional axes of the map like the latitude and longitude of the ordinary geological maps.

Atomic configuration in space simply called structure is the most fundamental attribute of molecular and condensed matter systems, especially crystalline materials, and determines the Hamiltonian for electronic states together with their composition under the Born-Oppenheimer approximation [3]. Our understanding of materials property is often based on its correlation with the structure. Usually, the huge sum of structure types exist in nature and can be classified by using symmetry in point group or space group (SG), that is ingeniously connected to the electronic states and related electronic properties.

For example, the electronic state of a condensed matter is often described by analyzing the symmetry of local structure around the target atom in terms of orbital splitting and hybridization. The electronic band structure of a crystalline system obeys its translational symmetry with the definition of the Brillouin zone on the basis of the Bloch theorem. In general, however, the structure dependence of the obtained electronic state and resulting properties is not so easily understood because of high-dimensionality in the structure space. With the desired structure map, the low-dimensional coordinates should classify huge related types of structure and the corresponding axes should provide structure descriptors for materials modeling following.

A common language is needed to start the analysis of the huge number of structures for a given general materials group. It should contain all the necessary but unambiguous information of the structure we investigate. In crystalline materials, the components of structural information are uniquely described with the crystallographic information framework (CIF) defined by the dictionary CoreCIF [4] established by International Union of Crystallography (IUCr) [5]. Most of the crystallographic databases such as ICSD [6], COD [7], and AtomWork [8] adopt the CIF format requisite for structure data standardization. Some crystallographic visualization tools also accept the CIF format data [9, 10].

Recently, accelerated by automated high-throughput computations with first-principles density-functional-theory (DFT) electronic structure calculation methods, computational materials databases and archives containing optimized, predicted, and designed crystal structure information compatible with CIF have been constructed like Materials Project [11], NOMAD [12], OQMD [13], AFLOW [14] and etc., forming big data on crystal structures continuously growing with the experimental structure databases. To understand materials properties associated with structure for further development and deployment, it may be indispensable and worthwhile to extract important features in structure for a given target group of materials. Generally, crystal structure consists of high-dimensional information and a dimensional reduction technique is required for the purpose. In this study, a workflow to draw such important features by using a dimension-reduction technique is proposed and demonstrated for some examples.

Obtained low-dimentional features of crystal structure can be used for classifying the structures (clustering as a data-science term) among the related materials systems and also for machine-learning modeling of structure-associated electronic properties of materials. In addition, a simple technique to create possible configurational structures by substituting atomic sites of a materials system with different elements in all combinational degrees of freedom is proposed for applications to alloy and impurity problems with supercell models.

2. Methods

A three-step procedure starting from a set of the CIF files of the target crystal systems for mapping the mutual position of each crystal structure in a low-dimensional space is described in the following sub-sections. Moreover, a generation method of structures with possible required configurations for a given set of compositions is explained in connection with the three-step workflow. The whole structure of the workflow is shown in Fig. 1.

Figure 1. Three-step workflow described in the present study. Capital italic name expresses each detailed procedure, CIF2ESC: conversion from a CIF file to electronic-structure-calculation data with output of fingerprint as structure feature for post structure analyses; DISTANCE: evaluation of distances between structure feature vectors; DISTANCE_DISTRIBUTION: calculation of distance distribution; DIAGNOSIS: extraction of inequivalent structures; CMDS: classical multidimensional scaling; PCEV: calculation of principal component eigenvectors; PROJECTION: projection of an arbitrary structure on the dimension-reduced map; ATLS: generation of atom-type list; FINDSYM${}^{*}$: space-group finder ${}^{*}$: Ref. [35, 36]. Open and solid (and dotted) arrows denote data transfers between the steps and inside each step, respectively. Dotted arrow is applicable if necessary.

2.1. STEP-1: Feature vector of crystal structure (abbreviated as CIF2ESC)

Crystal structure is the starting point of the present analysis and we assume that the detailed information for a set of materials systems is given with the CIF format. To consider a feature for a general structure, one should define a vector that is simple, tractable, translationally and rotationally invariant to treat, and often possibly highly dimensional to distinguish. Coordination number is the simplest information of local structure around a particular atomic site in a crystal. A cutoff distance to determine the nearest neighbor (NN) is rather subjective and not well defined for general structures. Radial distribution function (RDF) contains complete pair-wise structure information within a given radius and is directly measurable in extended x-ray absorption fine structure (EXAFS) experiments. An extension of RDF to more than pair correlations may be possible like cluster expansion [15].

Oganov and Valle [16] have proposed $\mathit{F}$-fingerprint as a structure feature defined for an $\mathit{AB}$ element pair of a system as

(1) $F_{\mathit{AB}}(\mathit{R})=\frac{V}{4\pi R^2{N}_A{N}_B}{\displaystyle \sum _{A_i\in \text{cell}}}{\displaystyle \sum _{B_j}}\mathit{\delta }(\mathit{R}-R_{\mathit{ij}})-\mathrm{1,}$(1)

where $A_i$ is the $\mathit{i}$-th site of element $\mathit{A}$, $R_{\mathit{ij}}$ is the distance between $A_i$ and $B_j$, $N_A$ is the number of sites of $\mathit{A}$ in cell, and $\mathit{V}$ is the volume of unit cell. In Eq (1(1) $F_{\mathit{AB}}(\mathit{R})=\frac{V}{4\pi R^2{N}_A{N}_B}{\displaystyle \sum _{A_i\in \text{cell}}}{\displaystyle \sum _{B_j}}\mathit{\delta }(\mathit{R}-R_{\mathit{ij}})-\mathrm{1,}$(1) ), the first term is nothing but RDF and the second term $(-1)$ is introduced so that the $\mathit{F}$-fingerprint becomes short-ranged with an asymptotic behavior of $\mathit{F}(\mathit{\infty })=0$. The $\mathit{\delta }$ function is smoothened with the Gaussian function of an appropriate broadening width $\mathit{\sigma }$ in a practical numerical process as

(2) $\mathit{\delta }(\mathit{R})\approx \frac{1}{\sqrt{2\mathit{\pi }}\sigma }\mathit{exp}\left[-\frac{R^2}{2{\sigma }^2}\right]\mathrm{.}$(2)

The $\mathit{F}$-fingerprint in Eq (1(1) $F_{\mathit{AB}}(\mathit{R})=\frac{V}{4\pi R^2{N}_A{N}_B}{\displaystyle \sum _{A_i\in \text{cell}}}{\displaystyle \sum _{B_j}}\mathit{\delta }(\mathit{R}-R_{\mathit{ij}})-\mathrm{1,}$(1) ) is discretized as $F_{\mathit{AB}}(R_l)\mathit{}(\mathit{l}=\mathrm{1,}\cdots \mathrm{,}\mathit{m})$ with an interval $\mathit{\Delta }$ up to an assumed maximum radius $R_m\equiv R_{\text{max}}$ and by connecting all independent RDFs for $\mathit{p}$ element pairs in series, a $\mathit{pm}(\equiv \mathit{n})$-dimensional ($\mathit{n}$-D) vector is finally obtained in a vector form. The parameters included in the feature vector, $R_{\text{max}}$, $\mathit{\sigma }$, and $\mathit{\Delta }$, should be determined by judging the evaluation of structure distance according to user’s target system and purpose. Generally speaking, $\mathit{\sigma }$ corresponds to the resolution of RDF for distinguishing given structures. Thus, typical values of $\mathit{\sigma }$ may be about one tenths of Å taken from the NN and next nearest neighbor (NNN) distance in the compounds. Accordingly, the discretizing interval $\mathit{\Delta }$ should be smaller to some extent than $\mathit{\sigma }$ to make the resolution meaningful. The cutoff radius $R_{\text{max}}$ should be selected about the cell size of the system. The RDF information beyond it might be redundant because of the periodicity in a model system.

The $\mathit{F}$-fingerprint has been successfully used in our previous works [17, 18, 19] for accelerating crystal structure prediction with the Bayesian optimization technique. Especially, some evaluations was made for investigating the efficiency of the crystal structure search with different choices of the parameters included in the $\mathit{F}$-fingerprint [18]. Extensions of the $\mathit{F}$-fingerprint representing higher-order crystal structure information might be also applicable in the present scheme by including angular degrees of freedom as local many-site correlations [20, 21].

2.2. STEP-2: Structure distance as a measure of similarity (DISTANCE and DIAGNOSIS)

Although similarity of crystal structure is not uniquely defined, let us consider feature-vector distance as a quantity to represent the similarity. The degree of the similarity can be considered in such a semi-quantitative way that two structures are equivalent if a structure distance $\mathit{D}$ satisfies $\mathit{D}<D_{\text{th}}$ with a given threshold $D_{\text{th}}$, and a pair of two structures 1 and 2 with distance $\mathit{D}(\mathrm{1,2})$ is closer in structural similarity than another pair of 1 and 3 with distance $\mathit{D}(\mathrm{1,3})$ if the distances satisfy $\mathit{D}(\mathrm{1,2})<\mathit{D}(\mathrm{1,3})$. In any case, $D_{\text{t}\mathit{h}}$ can be determined by inspecting the distance distribution and the error bars in the structure data (typically significant figures of lattice constants and atomic coordinates).

Given two $\mathit{n}$-D vectors $\mathit{X}$ and $\mathit{Y}$ representing the feature of crystal structures, the simplest distance between them is the Euclidean distance defined as

(3) $d_{\text{E}}(\mathit{X}\mathrm{,}\mathit{Y})={\left[{\displaystyle \sum _{\mathit{l}=1}^n}{(X_l-Y_l)}^2\right]}^{1/2}\mathrm{.}$(3)

For the present purpose, the Pearson correlation coefficient (PCC) can be another choice as a measure of structure similarity as

(4) $\mathit{r}(\mathit{X}\mathrm{,}\mathit{Y})=\frac{{\displaystyle {\sum }_{l=1}^n(X_l-\overline{X})(Y_l-\overline{Y})}}{{\left({\displaystyle {\sum }_{l=1}^n{(X_l-\overline{X})}^2{\displaystyle {\sum }_{l=1}^n{(Y_l-\overline{Y})}^2}}\right)}^{1/2}}\mathrm{,}$(4)

where $\overline{X}$ denotes the average of the vector components. Since PCC is distributed within the range of $(-1\le \mathit{r}\le +1)$, it can be transformed into a kind of distance as

(5) $d_{\text{C}}(\mathit{X}\mathrm{,}\mathit{Y})=\frac{1}{2}(1-\mathit{r})\mathrm{,}$(5)

known as the Cosine distance. In the present case with the vectors given in Eq (1(1) $F_{\mathit{AB}}(\mathit{R})=\frac{V}{4\pi R^2{N}_A{N}_B}{\displaystyle \sum _{A_i\in \text{cell}}}{\displaystyle \sum _{B_j}}\mathit{\delta }(\mathit{R}-R_{\mathit{ij}})-\mathrm{1,}$(1) ), the Cosine distance between structures 1 and 2 is explicitly computed as

(6) $d_{\text{C}}(\mathrm{1,2})=\frac{1}{2}\left(1-\frac{{\mathit{F}}_1\cdot {\mathit{F}}_2}{\left|{\mathit{F}}_1\right|\left|{\mathit{F}}_2\right|}\right)\mathrm{,}$(6)

where the inner product is defined as

(7) ${\mathit{F}}_1\cdot {\mathit{F}}_2={\displaystyle \sum _{\mathit{AB}}}{\displaystyle \sum _l}{F}_{\mathrm{1,}\mathit{AB}}(\mathit{l})F_{\mathrm{2,}\mathit{AB}}(\mathit{l})w_{\mathit{AB}}$(7)

with a weight proportional to the pair probability

(8) $w_{\mathit{AB}}=\frac{N_A{N}_B}{{\displaystyle {\sum }_{\text{cell}}{N}_A{N}_B}}\mathrm{.}$(8)

This means that the present scheme is limited to equal composition systems when comparing more than one distance. In the same way as Eq (6(6) $d_{\text{C}}(\mathrm{1,2})=\frac{1}{2}\left(1-\frac{{\mathit{F}}_1\cdot {\mathit{F}}_2}{\left|{\mathit{F}}_1\right|\left|{\mathit{F}}_2\right|}\right)\mathrm{,}$(6) ), the Euclidean distance between structures 1 and 2 is obtained as

(9) $d_{\text{E}}(\mathrm{1,2})={\left[{\displaystyle \sum _{\mathit{AB}}}{\displaystyle \sum _l}{\left[{\mathit{F}}_{\mathrm{1,}\mathit{AB}}(\mathit{l})-F_{\mathrm{2,}\mathit{AB}}(\mathit{l})\right]}^2{w}_{\mathit{AB}}\right]}^{1/2}\mathrm{.}$(9)

In high-throughput DFT calculations, a huge set of similar crystal structures are generated for a given composition. Some of them are essentially equivalent but are often judged to be different due to numerical errors or insufficient structural optimizations. The distance obtained above can be used for listing up crystallographically inequivalent structures among the structure set in a reasonable manner with a given threshold distance. This procedure may reduce the total computing time of DFT calculations and can be combined with a successive process for primitive-cell and symmetry determination. Without this, it is tempting to think of a large number of possible atomic configurations in an impurity or alloying problem that may contain many equivalent structures as we shall see in Sect. 3.2.

2.3. STEP-3(1): Structure mapping (CMDS)

For a given Euclidean or non-Euclidean distance matrix $(\mathit{N}\times \mathit{N})\mathit{}\mathit{D}$ constructed from $\mathit{N}$ structure sets of the $\mathit{n}$-D feature vectors $(\mathit{N}\times \mathit{n})\mathit{}\mathit{X}$, dimension-reduced coordinates of the structures $(\mathit{N}\times \mathit{k})\mathit{}{X}_{+}$ approximated in a low $\mathit{k}(\ll \mathit{n})$-D space can be obtained by using the multidimensional scaling (MDS) [22,23,24]. (For clarity, the size of matrix is indicated like $(\mathit{N}\times \mathit{k})$ when it first appears.) The MDS scheme is often called “dimension reduction” in data mapping and is equivalent to principal component analysis (PCA) when Euclidean distance is applied. In MDS, the pair-wise Euclidean distance matrix $(\mathit{N}\times \mathit{N})\mathit{}{D}^X$ calculated using $X_{+}$ is the closest approximation to D. The algorithm of MDS is quite simple as follows, starting with squared distance matrix $(\mathit{N}\times \mathit{N})\mathit{}{D}^{(2)}$ of which the components are squared values of the original distance matrix components. First, the squared distance matrix is transformed by the so-called double centering $(\mathit{N}\times \mathit{N})\mathit{}\mathit{J}$ into a matrix $(\mathit{N}\times \mathit{N})\mathit{}\mathit{B}$ as

(10) $\mathit{B}=-\frac{1}{2}\mathit{J}{D}^{(2)}\mathit{J}\mathit{}\text{with}\mathit{}\mathit{J}=\mathit{I}-\frac{1}{N}11^T\mathrm{,}$(10)

where $\mathit{I}$ and $1$ are identity matrix and identity column vector, respectively, and the superscript $\mathit{T}$ stands for matrix transpose. Next, the matrix $\mathit{B}$ is eigenvalue decomposed as

(11) $\mathit{B}=\mathit{Q\Lambda }{Q}^T\mathrm{.}$(11)

Here, the matrix $(\mathit{N}\times \mathit{N})\mathit{}\mathit{\Lambda }$ contains the eigenvalues in the diagonal elements with zeros in the off-diagonal ones and $(\mathit{N}\times \mathit{N})\mathit{}\mathit{Q}$ is the corresponding eigenvector matrix. Finally, by taking the first $\mathit{k}$ positive-eigenvalue part $(\mathit{k}\times \mathit{k})\mathit{}{\Lambda }_{+}$ and $(\mathit{N}\times \mathit{k})\mathit{}{Q}_{+}$ in the order from the positive-largest one, the principal components (PCs) $(\mathit{N}\times \mathit{k})\mathit{}{X}_{+}$ approximately representing the structure coordinates can be obtained as

(12) $X_{+}=Q_{+}{\Lambda }_{+}^{1/2}\mathrm{,}$(12)

where $(\mathit{k}\times \mathit{k})\mathit{}{\Lambda }_{+}^{1/2}$ is the matrix of which diagonal elements are the square root of the eigenvalues as ${\Lambda }_{+}^{1/2}{\Lambda }_{+}^{1/2}={\Lambda }_{+}$. Then, the matrix B is approximated in the form of scalar product as

(13) $\mathit{B}\approx B_{+}=X_{+}{X}_{+}^T\mathrm{.}$(13)

As indicated in Eq (12(12) $X_{+}=Q_{+}{\Lambda }_{+}^{1/2}\mathrm{,}$(12) ), the deviation of the coordinates of each structure is proportional to the square-root of the eigenvalue, since the eigenvector $\mathit{Q}$ in Eq (11(11) $\mathit{B}=\mathit{Q\Lambda }{Q}^T\mathrm{.}$(11) ) is usually taken in a unitary form. Thus, the size of the eigenvalue determines the importance of each PC obtained by the dimension reduction.

As Sibson has pointed out [25], the $\mathit{k}$ value that we may choose can be determined by judging the magnitude of positive eigenvalues since the sum of the eigenvalues in ${\Lambda }_{+}$ approximate the sum of all eigenvalues in $\mathit{\Lambda }$. This point is practically quite important and crucial when determining the number of PCs and evaluating the performance of feature vectors and distances we adopted [26]. To see how much original information remains after the dimension reduction, it is convenient to check the so-called proportion defined by

(14) $P_k\equiv \frac{{\displaystyle {\sum }_{i=1}^k{\Lambda }_{\mathit{ii}}}}{{\displaystyle {\sum }_{i=1}^N{\Lambda }_{\mathit{ii}}}}\mathrm{.}$(14)

Note that the summation in the numerator of Eq (14(14) $P_k\equiv \frac{{\displaystyle {\sum }_{i=1}^k{\Lambda }_{\mathit{ii}}}}{{\displaystyle {\sum }_{i=1}^N{\Lambda }_{\mathit{ii}}}}\mathrm{.}$(14) ) is taken in the descending order starting from the positive largest eigenvalue.

2.4. STEP-3(2): Projection on map and its inverse operation (PCEV and PROJECTION)

Let us describe a way to get the eigenvector of PCs in the original $\mathit{n}$-D space and to make the projection of a given new structure on the obtained low-dimensional map above, by considering the relation to PCA [27]. In PCA, the data covariance matrix $(\mathit{n}\times \mathit{n})\mathit{}\mathit{S}=X^T\mathit{X}$ is the target quantity and is maximized with a unitary matrix $(\mathit{n}\times \mathit{n})\mathit{}\mathit{u}$ and the Lagrange multiplier $\mathit{\lambda }$ as

(15) $\frac{\delta }{\delta u^T}\left[u^T\mathit{Su}+\text{λ}(1-u^T\mathit{u})\right]=\mathrm{0,}$(15)

leading to an eigenvalue problem of

(16) $\mathit{Su}=\mathit{u\Lambda }\mathrm{.}$(16)

Note that the coordinate matrix $(\mathit{N}\times \mathit{n})\mathit{}\mathit{X}$ defined in this paper follows Borg-Groenen’s prescription [23] in a transpose form of Bishop’s definition [27]. By combining Eqs (12(12) $X_{+}=Q_{+}{\Lambda }_{+}^{1/2}\mathrm{,}$(12) ), (13(13) $\mathit{B}\approx B_{+}=X_{+}{X}_{+}^T\mathrm{.}$(13) ), and (16(16) $\mathit{Su}=\mathit{u\Lambda }\mathrm{.}$(16) ), the PC eigenvector $(\mathit{n}\times \mathit{k})\mathit{}{u}_{+}$ can be computed in a normalized form as

(17) $u_{+}=\left(X^T-1{\overline{X}}^{\mathit{T}}\right)Q_{+}{\Lambda }_{+}^{-\frac{1}{2}}=\left(X^T-1{\overline{X}}^{\mathit{T}}\right)X_{+}{\Lambda }_{+}^{-1}$(17)

with the center coordinate of the original feature vectors as

(18) ${\left\{\overline{X}\right\}}_{1\mathit{l}}=\frac{1}{N}{\displaystyle \sum _{\mathit{i}=1}^N}{\left\{\mathit{X}\right\}}_{\mathit{il}}\mathit{}(\mathit{l}=\mathrm{1,}\cdots \mathrm{,}\mathit{n})\mathrm{.}$(18)

The original coordinate matrix $\mathit{X}$ which is used to calculate the distance is not necessarily centered and is enforced to be so in Eq (17(17) $u_{+}=\left(X^T-1{\overline{X}}^{\mathit{T}}\right)Q_{+}{\Lambda }_{+}^{-\frac{1}{2}}=\left(X^T-1{\overline{X}}^{\mathit{T}}\right)X_{+}{\Lambda }_{+}^{-1}$(17) ). By using the obtained PC eigenvector $u_{+}$, the projection of any arbitrary structure feature vector $(1\times \mathit{n})\mathit{Y}$ in the original dimension space on the reduced $\mathit{k}$-D map can be made as

(19) $Y_{+}=\left(\mathit{Y}-\overline{X}\right)u_{+}\mathrm{.}$(19)

With help of the minimum-error formulation of PCA [27], an inverse projection, that is the transformation of an arbitrary reduced coordinate on the map $(1\times \mathit{k})\mathit{}{Z}_{+}$ into a feature vector $(1\times \mathit{n})\mathit{}\mathit{Z}$ in the original structure space, can be made as

(20) $\mathit{Z}=Z_{+}{u}_{+}^T+\overline{X}\mathrm{.}$(20)

Since a certain amount of information is lost due to the dimension reduction, the inverse projection can only be accomplished approximately. The projection and inverse projection procedures described above are reasonably made if the distance matrix is calculated as the Euclidean distance since MDS and PCA are rigorously equivalent. Even in non-Euclidean distance case, the procedures may possibly be carried out.

2.5. Generation of structures with all possible configurations (ATLS)

In addition to stoichiometric ordered materials, alloyed and doped impurity systems are often seen as targets in materials research. In such systems, counting all the possible configurations is crucially indispensable for optimizing materials properties. In electronic structure calculations for alloy systems, coherent-potential approximation (CPA) [28, 29] is usually adopted by taking single-site average of scattering path operators within the Korringa-Kohn-Rostoker (KKR) Green’s function formalism [31, 32]. This approximation works quite well for randomly distributed alloy systems. Similar idea with supercell models has been proposed with special quasi-random structures (SQS’s) method, where the site configuration is optimized by making multiple site-correlation functions best close to those of random structure [33, 34].

However, partially ordered or site-preference features are often seen and essential to govern the materials properties in real materials. Listing up all the possible configurational combinations is quite useful to investigate such local-site-order associated properties for a given composition with specific elements included. Such a tool is developed as ATLS. For a example, let us consider a $(2\times 2\times 2)$ body-centered supercell model containing eight atoms (four Fe and four Co) with BCC-base lattice structure. All the possible combinations are ${}_8{\text{C}}_4=$70 ways. With given lattice constants, lattice centering types, atomic composition, and corresponding atomic coordinates, SG and related structural properties can be obtained by using FINDSYM software [35, 36]. By using the workflow described above, only four structures are found to be crystallographically inequivalent (32 $\mathit{R}\overline{3}\mathit{m}$ + 12 $\mathit{Cmmm}$ + 24 $\mathit{I}{4}_1\mathit{amd}$ + 2 $\mathit{Pm}\overline{3}\mathit{m}$) out of the 70 configurations. In the present work, crystallographically equivalent structures mean that their $\mathit{F}$-fingerprints are equivalent within the reasonable threshold and have the same SG. Interestingly, the $\mathit{Cmmm}$ and $\mathit{I}{4}_1\mathit{amd}$ structures extracted above have the exactly same $\mathit{F}$-fingerprint within minor numerical error bars and cannot be distinguished only by the pair-wise information. This is a drawback of the present choice of the fingerprint and structure features beyond the pair-wise correlation should be included if necessary. Actually, it is not the case shown in the example of La(Fe,Si)${}_{13}$ (see Sect. 3.2). In the case of a $(2\times 2\times 2)$ primitive supercell model containing 16 atoms (eight Fe and eight Co) with BCC-base lattice structure, there are ${}_{16}{\text{C}}_8=12870$ combinations and 51 crystallographically inequivalent configurations. Similar to the case above, some configurations have the exactly same fingerprint, resulting in 32 independent pair-wise fingerprints. Some interesting situations may happen if the mother system contains more than one crystallographically inequivalent atomic sites, leading to a site-preference problem. Such a case will be given below in Section 3.2.

2.6. DFT calculations

All the DFT electronic structure calculations in the present study are performed with the all-electron full-potential linearlized augmented-plane-wave (FLAPW) method [37, 38, 39] implimented in the HiLAPW program [40]. The exchange and correlation are incorporated within the Perdew-Burke-Ernzerhof form [41] of the generalized gradient approximation to DFT. The energy cutoffs for wavefunction and potential expansions are 20 Ry and 160 Ry, respectively. A set of $\mathit{\Gamma }$-centered $\mathit{k}$-point mesh with the interval of approximately 0.1 Bohr${}^{-1}$ is used to sample the Brillouin zone with the improved tetrahedron integration scheme [42] in all the self-consistent-field DFT calculations.

3. Results and discussion

3.1. Classification of Al${}_2$O${}_3$ polymorphs

There often exist several stable and meta-stable types of crystal structure with the same composition, called polymorph, especially in oxide systems. Silica SiO${}_2$ and alumina Al${}_2$O${}_3$ are such most typical cases seen in materials science and engineering. For example, there are 15 different crystal structures of Al${}_2$O${}_3$ found in the Materials Project site [11]. One of them, mp-1245063 seems to be an amorphous model structure with 100 atoms in a unit cell (SG: $\mathit{P}1$). In addition, another amorphous model structure including 120 atoms in a unit cell was generated by a first-principles molecular-dynamics simulation with a quenching scheme from a liquid state [43]. Structural analyses on Al${}_2$O${}_3$ for several different phases have been carried out in relation to their properties [44,45,46].

A map for the Al${}_2$O${}_3$ structures is created by using the three-step procedure proposed above, and especially it’s quite interesting to look at the positions of two amorphous model structures on the map of the polymorph structures as shown below. (One structure, mp-985587 is a model for an artificial 2D film, which will be excluded in the following analysis.)

Calculated DFT total energies and energy gaps are listed in Table 1 together with the corresponding values given in the Materials Project [11]. Here, the same lattice constants and atomic coordinates as the previous works [11, 43] were assumed in the present study. The present DFT total energies show almost equivalent values to the Materials Project data, while the energy gaps are slightly underestimated as a general tendency. Here, the highest energy value of the valence bands and the lowest energy value of the conduction bands at the calculated $\mathit{k}$-mesh points are used to evaluate the energy gaps.

Table 1. List of Al${}_2$O${}_3$ structures in the order of calculated total energy $\mathit{E}$ (meV/atom) relative to the most stable $\mathit{R}\overline{3}\mathit{c}$ structure, and the energy gap $E_G$ (eV). Figures in parentheses are values given in the Materials Project site [11]. Color scheme is used for distinguishing different structures in Figs 2 and 4.

ref.	color	SG	$\mathit{E}$	$E_G$	remarks
mp-1143	black	$\mathit{R}\overline{3}\mathit{c}$	0 (0)	5.596 (6.044)	most stable
mp-7048	maroon	$\mathit{C}2/\mathit{m}$	+17 (+10)	4.248 (4.605)	most stable phase of Ga${}_2$O${}_3$
mp-2254	blue	$\mathit{Pna}21$	+18 (+17)	4.592 (4.830)
mp-776475	pink	$\mathit{Ia}\overline{3}$	+43 (+30)	4.982 (5.216)	most stable phase of In${}_2$O${}_3$
mp-32591	green	$\mathit{Cm}$	+63 (+57)	2.854 (3.075)
mp-754624	teal	$\mathit{R}\overline{3}$	+71 (+66)	5.383 (5.797)
mp-754401	magenta	$\mathit{Cmc}21$	+74 (+46)	4.022 (4.235)
mp-755175	purple	$\mathit{P}\overline{3}1\mathit{c}$	+79 (+72)	5.227 (5.494)
mp-1938	lime	$\mathit{Pbcn}$	+92 (+94)	5.047 (5.313)
mp-755483	olive	$\mathit{R}3\mathit{m}$	+93 (+91)	3.863 (4.098)
mp-754531	yellow	$\mathit{P}21/\mathit{c}$	+98 (+77)	4.141 (4.374)
Momida${}^{\#}$	red (thick)	$\mathit{P}1$	+179 (—)	3.725 (—)	amorphous model
mp-985587	—	$\mathit{P}321$	+210 (+210)	4.686 (4.807)	artificial bilayer film
mp-1245063	red (broken)	$\mathit{P}1$	+217 (+204)	3.137 (3.326)	amorphous model
mp-642363	cyan	$\mathit{Cmcm}$	+266 (+281)	4.009 (4.249)	elastically unstable
mp-638765	navy	$\mathit{P}\overline{1}$	+682 (+809)	1.102 (1.136)

^#: Ref. [43]

To calculate the $\mathit{F}$-fingerprint for all the systems except for the 2D film model mp-985587, parameters $\mathit{\sigma }=0.2$ Å and $\mathit{\Delta }=0.1$ Å were set according to the arguments given in Sect. 2.1. Actually, clear dips are commonly seen just beyond the first NN O-Al pair shell in the $\mathit{F}$-fingerprint of all the systems as shown in Fig. 2. Two amorphous structures are quite similar each other in the $\mathit{F}$-fingerprint while some recognizable differences can be seen in the O-Al and Al-Al pair regions beyond the NN distance. They are also similar to the crystalline polymorphs within the first NN shell of the O-O, O-Al, and Al-Al pairs and turn to be remarkably broadened beyond. Such features were found and discussed in detail by Momida et al. [43]. Distances are calculated using Euclidean and Cosine distance schemes up to the radius cutoff of $R_{\text{max}}=10$ Å and $R_{\text{max}}=5$ Å and the eigenvalues and eigenvectors of the double-centered squared distance matrix are computed to get the structure coordinates as described in Sects. 2.2 and 2.3.

Figure 2. $\mathit{F}$-fingerprint of all the Al${}_2$O${}_3$ systems listed in Table 1 except for the 2D film model. Radius cutoff $R_{\mathit{max}}=10$ Å is assumed and the fingerprint of dimensions between 0 and 100, 101 and 200, and 201 and 300 corresponds to O-O, O-Al, and Al-Al pairs, respectively. Line colors follow color scheme in Table 1 for distinguishing different structures.

In Fig. 3, calculated eigenvalues are drawn for two radial cutoffs $R_{\text{max}}=10$ Å and $R_{\text{max}}=5$ Å with Euclidean and Cosine distances (Eqs (9(9) $d_{\text{E}}(\mathrm{1,2})={\left[{\displaystyle \sum _{\mathit{AB}}}{\displaystyle \sum _l}{\left[{\mathit{F}}_{\mathrm{1,}\mathit{AB}}(\mathit{l})-F_{\mathrm{2,}\mathit{AB}}(\mathit{l})\right]}^2{w}_{\mathit{AB}}\right]}^{1/2}\mathrm{.}$(9) ) and (6(6) $d_{\text{C}}(\mathrm{1,2})=\frac{1}{2}\left(1-\frac{{\mathit{F}}_1\cdot {\mathit{F}}_2}{\left|{\mathit{F}}_1\right|\left|{\mathit{F}}_2\right|}\right)\mathrm{,}$(6) ), respectively). Typical lattice constants of the systems we consider are up to about 10 Å and the cutoff $R_{\text{max}}=10$ Å should naívely be appropriate as the first choice. To see effects of the assumed cutoff on structure analyses, a study with $R_{\text{max}}=5$ Å is also tested. Looking at the Euclidean-distance eigenvalues in Fig. 3, $R_{\text{max}}=10$ Å gives relatively larger eigenvalues than $R_{\text{max}}=5$ Å. However, the proportion with $R_{\text{max}}=10$ Å ($P_5=0.91$ and $P_8=0.98$) is not so different from that with $R_{\text{max}}=5$ Å ($P_5=0.84$ and $P_8=0.95$). In the case of Cosine distance, the behavior of eigenvalues closely resembles the Euclidean one but information estimation only by the proportion is a bit ambiguous due to the existence of negative eigenvalues. It should be done by testing the precision in models constructed with the obtained PCs as descriptors, as we shall see below.

Figure 3. Eigenvalues calculated from the fingerprint of Al${}_2$O${}_3$ polymorphs using (a) $R_{\text{max}}=10$ Å and (b) $R_{\text{max}}=5$ Å with Euclidean (blue) and Cosine (red) distances.

Figure 4 shows coordinate map of all the Al${}_2$O${}_3$ structures on five PCs ($C_1$ through $C_5$ as the first five components of $X_{+}$ given in Eq (12(12) $X_{+}=Q_{+}{\Lambda }_{+}^{1/2}\mathrm{,}$(12) )) with different radial cutoff and type of distance. The distribution of the structures in the maps with different cutoff and type of distance is generally quite similar. Most typically, the three lowest-energy structures (SG: $\mathit{R}\overline{3}\mathit{c}$ (black), $\mathit{C}2/\mathit{m}$ (maroon), and $\mathit{Pna}21$ (blue)) form line-shaped clustering on the $C_1$-$C_2$ dimension. Relatively low-energy structures up to about 100 meV/atom surround the lowest-energy cluster forming a basin except for $\mathit{Ia}\overline{3}$ structure (pink) that is the most stable for In${}_2$O${}_3$, far outside along the PC $C_1$. Two amorphous structures are located near the edge of the basin close to the highest-energy structure ($\mathit{P}\overline{1}$). The $\mathit{Pbcn}$ structure is situated close to the lowest-energy one ($\mathit{R}\overline{3}\mathit{c}$) along $C_1$ through $C_3$ axes and found to be distinguishable along $C_4$ or $C_5$. It is also interesting that the stable structures of Ga${}_2$O${}_3$ ($\mathit{C}2/\mathit{m}$: maroon) and In${}_2$O${}_3$ ($\mathit{Ia}\overline{3}$: pink) reveal relatively low energy while their positions in the map are far especially along $C_1$. Since the PC $C_1$ dominantly contains information of the first NN O-Al pairs as shown below, the coordination number and the length variation in the NN O-Al bonds by local symmetry seem to affect the distance (dissimilarity) in the structure map.

Figure 4. Map of Al${}_2$O${}_3$ structures with the five principal components (PCs) with (a) Euclidean-distance $R_{\text{max}}=10$ Å, (b) Cosine-distance $R_{\text{max}}=10$ Å, (c) Euclidean-distance $R_{\text{max}}=5$ Å, and (d) Cosine-distance $R_{\text{max}}=5$ Å. Dot colors follow color scheme listed in Table 1. Square dots in each figure denote the three lowest energy structures (SG: $\mathit{R}\overline{3}\mathit{c}$ (black), $\mathit{C}2/\mathit{m}$ (maroon), and $\mathit{Pna}21$ (blue)), and small and large red dots express amorphous model structures (SG: $\mathit{P}1$) by the Materials Project [11] and Momida [43], respectively.

To investigate the information and its meaning immanent in each PC and the corresponding coordinates, the eigenvector of the obtained PCs is calculated using Eq (17(17) $u_{+}=\left(X^T-1{\overline{X}}^{\mathit{T}}\right)Q_{+}{\Lambda }_{+}^{-\frac{1}{2}}=\left(X^T-1{\overline{X}}^{\mathit{T}}\right)X_{+}{\Lambda }_{+}^{-1}$(17) ) and plotted in Fig. 5. By searching the largest amplitude of each PC eigenvector in Fig. 5, one can find which radial part of the $\mathit{F}$-fingerprint, namely RDF gives the most important contribution to the PC. The fist PC $C_1$ has the largest amplitude at the first NN position of the O-Al pair, irrespective to the choice of the radial cutoff and type of distance. The second PC $C_2$ shows the largest contribution from the second and third NNs of O-Al pair and almost equally the first NN of O-O pair. From these findings, it seems that the characteristics in the energy landscape seen in Fig. 4 may be quantitatively explained by the O-Al and O-O neighboring environmental features. However, the PCs beyond $C_1$ and $C_2$ are also significantly important to reproduce the energy more precisely, as shown below. In this context, the largest amplitudes for $C_3$, $C_4$, and $C_5$ can be seen for the Al-Al RDF.

Figure 5. Eigenvector of the five principal components (PCs) with (a) Euclidean-distance $R_{\text{max}}=10$ Å, (b) Cosine-distance $R_{\text{max}}=10$ Å, (c) Euclidean-distance $R_{\text{max}}=5$ Å, and (d) Cosine-distance $R_{\text{max}}=5$ Å, calculated by using Eq (17(17) $u_{+}=\left(X^T-1{\overline{X}}^{\mathit{T}}\right)Q_{+}{\Lambda }_{+}^{-\frac{1}{2}}=\left(X^T-1{\overline{X}}^{\mathit{T}}\right)X_{+}{\Lambda }_{+}^{-1}$(17) ) in Al${}_2$O${}_3$ polymorphs. Red, lime, blue, maroon, and green lines denote $C_1$, $C_2$, $C_3$, $C_4$, and $C_5$ PC, respectively. Horizontal axis denotes the dimension of the original fingerprint space (see Fig. 2). In (a) and (b), the dimensions between 1 and 100, 101 and 200, and 201 and 300 express the O-O, O-Al, and Al-Al pair radius, respectively, and in (c) and (d), the dimensions between 1 and 50, 51 and 100, and 1001 and 150 do the O-O, O-Al, and Al-Al pair radius, respectively.

A sparse modeling technique is useful to search the important descriptors among a given set of ones. In the present study, we prepared 21 descriptors of the zero-th (constant term) through second orders of the obtained five PCs ($C_1$—$C_5$), and detected and removed their collinearity and multi-collinearity (linear dependence), and then used exhaustive search or genetic algorithm in linear regression analysis. The best set of the descriptors for a fixed number of descriptors is selected by the leave-one-out (LOO) scheme, a kind of cross validation techniques [47]. Results of the regression analysis for the total energy using the PCs as the descriptors are shown in Fig. 6. $R^2$ is the coefficient of determination widely used for expressing the fitting performance [48], $Q^2$ is the coefficient of determination for the prediction residual sum of squares (often called PRESS) showing a measure of generalization capability [49] and RMSE denotes the root mean square error in the regression. The explicit equations of these quantities are given and discussed in our previous study [47]. At glance, good modeling is achieved as the number of descriptors is increased in all the cases and any over-fitting behaviors those are often observed in $Q^2$ are not seen, possibly because of the use of orthogonal PCs as the descriptors. Although the products of the PCs used may not necessarily be orthogonal to the others, any collinearity or multi-collinearity tendency has not been detected up to the second order of the PCs. As a more rigorous test for investigating any existing over-fitting, a double cross-validation technique is adopted. Here, one sample data are completely hidden during the modeling and then the best model trained by the LOO method with the remaining data is tested for checking the generalization capability using the hold-out data. Double cross-validation results for the total-energy model constructed from the Euclidean distance with $R_{\text{m}\mathit{ax}}=5$ Å are shown in Fig. 7. As the number of descriptors is increased, the standard deviation for the prediction error of the hold-out data shows first a decline, then minimum ($\sim 0.1$ eV) at the number of descriptors of eight and an ascent beyond, indicating the over-fitting behavior with large numbers of descriptors. The explicit model form with the number of descriptors of eight is

Figure 6. Regression analysis of the total energy with descriptors up to the second order of the five PCs in Al${}_2$O${}_3$ polymorph structures. (a) Euclidean-distance $R_{\text{max}}=10$ Å, (b) Cosine-distance $R_{\text{max}}=10$ Å, (c) Euclidean-distance $R_{\text{max}}=5$ Å, and (d) Cosine-distance $R_{\text{max}}=5$ Å. Blue, red, and green dots and lines denote the coefficient of determination ($R^2$), its leave-one-out coefficient ($Q^2$), and root mean square error (RMSE), respectively.

Figure 7. Double cross-validation test of the total-energy regression model constructed with the descriptors from the Euclidean distance and $R_{\text{max}}=5$ Å in Al${}_2$O${}_3$ polymorph structures. Dot colors follow color scheme listed in Table 1. Red open circle denotes one of the amorphous structure (mp-1245063). The model of the number of descriptors of eight shows the smallest standard deviation ($\sim$0.1 eV).

$E^{(8)}=0.049C_1+0.122C_3+0.025C_1{C}_2$

$+0.123C_1{C}_5+0.110C_3^2+0.638C_3{C}_5$

(21) $+0.228C_4^2+0.358C_5^2\mathrm{.}$(21)

All of the PCs, $C_1$ through $C_5$ are selected in the model. Detailed energy landscape is actually very complicated than just looking at the map (Fig. 4). The double cross-validation test tells that the line-shaped dip seen in the $C_1$-$C_2$ map is rather shallow in the basin landscape and not well precisely reflected in the regression model, though higher-energy structures than 0.1 eV including the amorphous models are well distinguished in the energy model.

It is found, in the same way as the total-energy regression, that the energy-gap data given in Table 1 can be also modeled with the same structural descriptors, namely $C_1$ through $C_5$ and their products up to the second order plus the constant term, at the similar level of precision.

3.2. Site preference in La(Fe,Si)${}_{13}$

As a second application of the proposed workflow, we consider the classification of alloy-type configurations for analyzing the site preference of doped Si atoms in LaFe${}_{13-\mathit{x}}$Si${}_x$. This compound is well known as a magnetocaloric materials in applications for magnetic refrigeration [50]. The crystal structure of LaFe${}_{13-\mathit{x}}$Si${}_x$ is of cubic NaZn${}_{13}$ type (SG: $\mathit{Fm}\overline{3}\mathit{c}$) for $1\le \mathit{x}\le 2.8$ [51, 52, 53] and drawn in Fig. 8 with VESTA [9]. There are two kinds of crystallographic sites for Fe: Fe1 ($8\mathit{b}$ Wyckoff position) and Fe2 ($96\mathit{i}$), and doped Si atoms share the $96\mathit{i}$ site with Fe2. Configurations of Fe2 and Si are discussed as “coloring problem” [53]. Let us consider below a primitive cell La${}_2$Fe${}_{26-\mathit{y}}$Si${}_y$ ($\mathit{y}=2\mathit{x}$) of the system with fixed experimental lattice constants and atomic positions taken from the database COD-1528546 [7].

Figure 8. Crystal structure of La${}_2$Fe${}_{26-\mathit{y}}$Si${}_y$. Green, orange, and brown balls denote La, Fe1 and Fe2 atomic sites, respectively. Partial occupation of Si at the Fe2 sites is depicted with blue. The drawing is made with VESTA [9].

Looking at the local coordinates around Fe, the Fe1 ($8\mathit{b}$) site is coordinated by 12 Fe2 at the distance of 2.46 Å, which is the similar NN distance to that of pure BCC Fe. This means that the Fe1 site is closely packed with a high coordination number. Around the Fe2 ($96\mathit{i}$) site, on the other hand, there are two Fe2 at 2.45 Å, one Fe1 at 2.46 Å, and seven Fe2 in a range of 2.50–2.69 Å, showing less densely packed than Fe1. Although the atomic radius of Si is 1.17 Å slightly smaller than that of Fe of 1.24 Å [54], the site preference of Si at the Fe sites is not so easily understandable at first glance.

In the case of La${}_2$Fe${}_{25}$Si${}_1$, there are only two possible inequivalent configurations with one Si atom either at the $8\mathit{b}$ or $96\mathit{i}$ site (Case 1). It turns to be a bit complicated situation in La${}_2$Fe${}_{24}$Si${}_2$. There is one configuration with two Si atoms fully occupying at the $8\mathit{b}$ site (Case 2-1). If two Si atoms are placed one by one at the $8\mathit{b}$ and $96\mathit{i}$ sites (Case 2-2), total number of possible configurations is 48 (=${}_2$C${}_1$$\times$${}_{24}$C${}_1$). By analyzing the structural similarity as described above, it is reduced to two crystallographically inequivalent structures. If two Si atoms are accommodated at the $96\mathit{i}$ sites (Case 2-3), there are 276 (=${}_{24}$C${}_2$) possible configurations. And then, one can easily obtain 12 crystallographically inequivalent structures by the distance diagnosis. Similarly, in La${}_2$Fe${}_{23}$Si${}_3$ there are one crystallographically inequivalent structure out of 24 possible configurations for two Si at $8\mathit{b}$ and one Si at $96\mathit{i}$ (Case 3-1), 17 crystallographically inequivalent structures out of 552 (=${}_2$C${}_1$$\times$${}_{24}$C${}_2$) configurations for one Si at $8\mathit{b}$ and two Si at $96\mathit{i}$ (Case 3-2), and 50 inequivalent structures out of 2024 (=${}_{24}$C${}_3$) configurations for three Si at $96\mathit{i}$ (Case 3-3). In all the cases of this sub-section, the crystallographically inequivalent structures correspond to the distinct pair-wise fingerprints.

Heats of formation per atom are evaluated from total energies by first-principles DFT calculations with the FLAPW method for the configuration ($\mathit{i}$) in La${}_2$Fe${}_{26-\mathit{y}}$Si${}_y$ ($\mathit{y}=0-3$) and most-stable elemental crystal systems, La ($\mathit{Fm}\overline{3}\mathit{m}$), ferromagnetic Fe ($\mathit{Im}\overline{3}\mathit{m}$), and Si ($\mathit{Fd}\overline{3}\mathit{m}$) as

(22) $\begin{array}{ll}\mathit{H}(\mathit{y};\mathit{i})\hfill & =\frac{1}{28}\{\mathit{E}[\text{L}{a}_2\text{F}{e}_{26-\mathit{y}}\text{S}{i}_y;\mathit{i}]-2\mathit{E}[\text{L}\mathit{a}]\hfill \\ \hfill & -(26-\mathit{y})\mathit{E}[\text{F}\mathit{e}]-\mathit{yE}[\text{S}\mathit{i}]\}\hfill \end{array}$(22)

and plotted in Fig. 9 together with spin magnetic moments per Fe atom. To see the configurational energy for a given $\mathit{y}$, the heats of formation relative to the lowest energy value at $\mathit{y}$ that correspond to the configurational energy defined as

Figure 9. Heats of formation $\mathit{H}$ in eV/atom (Eq (22(22) $\begin{array}{ll}\mathit{H}(\mathit{y};\mathit{i})\hfill & =\frac{1}{28}\{\mathit{E}[\text{L}{a}_2\text{F}{e}_{26-\mathit{y}}\text{S}{i}_y;\mathit{i}]-2\mathit{E}[\text{L}\mathit{a}]\hfill \\ \hfill & -(26-\mathit{y})\mathit{E}[\text{F}\mathit{e}]-\mathit{yE}[\text{S}\mathit{i}]\}\hfill \end{array}$(22) )) and $\mathit{\Delta H}$ in meV/atom (Eq (23(23) $\mathit{\Delta H}(\mathit{y};\mathit{i})=\mathit{H}(\mathit{y};\mathit{i})-{\text{min}}_j\mathit{H}(\mathit{y};\mathit{j})$(23) )), and spin magnetic moments in ${\mu }_{\text{B}}$ per Fe atom calculated for La${}_2$Fe${}_{26-\mathit{y}}$Si${}_y$. Red dots, green dots, and blue circles denote values for the configurations with Si occupying at only $8\mathit{b}$, at $8\mathit{b}$ and $96\mathit{i}$, and at only $96\mathit{i}$, respectively. Black dots indicate the corresponding values for La${}_2$Fe${}_{26}(\mathit{y}=0)$.

(23) $\mathit{\Delta H}(\mathit{y};\mathit{i})=\mathit{H}(\mathit{y};\mathit{i})-{\text{min}}_j\mathit{H}(\mathit{y};\mathit{j})$(23)

are also drawn. The calculated heat of formation for non-doped La${}_2$Fe${}_{26}$ is $\mathit{H}(\mathit{y}=0)=+60$ meV/atom, indicating unstable to the constituent elemental systems consistent with no intermetallic-compound phase in the experimental phase diagram [55]. As the composition of Si $\mathit{y}$ is increased in La${}_2$Fe${}_{26-\mathit{y}}$Si${}_y$, the systems show a nearly linear stability to $\mathit{y}$ with negligibly small configurational variation at the energy scale of eV/atom within the same $\mathit{y}$. Minor in magnitude but interesting behaviors in the configurational energy $\mathit{\Delta H}$ at the 10meV/atom energy scale are discussed below. Calculated spin magnetic moment per Fe atom is monotonically decreased as $\mathit{y}$ is increased, leading to reduction in total magnetization per cell more than that simply by the decrease in the number of Fe atoms. Interestingly, the spin magnetic moment shows certain dependence on the configuration within the same Si concentration possibly due to differences in the local hybridization between Fe and Si atoms.

The configurational energy at $\mathit{y}=1$ in Fig. 9 tells that Si prefers to be accommodated at the Fe2 ($96\mathit{i}$) site by 11 meV/atom compared with the Fe1 ($8\mathit{b}$) site, being consistent with the experimental data (COD-1528546) [7]. Since the multiplicity of the Fe2 site (24) is significantly larger than that of the Fe1 site (2), the entropic effect must additionally contribute to the site-preference of Si at the Fe2 site. A rough estimation of the difference in the configurational entropy gives $\mathit{\Delta S}=k_{\text{B}}(\mathit{ln}24-\mathit{ln}2)=0.214$ meV/K-cell, implying the same order contribution as the energy difference to the stability at temperatures of several hundreds K. By inspecting the configurational energy $\mathit{\Delta H}$ more carefully in Fig. 9, clustering features can be seen for Si occupying at the $96\mathit{i}$ sites (most remarkably in blue circles of $\mathit{y}=2$ and $\mathit{y}=3$). Since the clustering features are possibly related to the origin of the configurational stability, structure map for all the inequivalent configurations may be useful for understanding it with the benefit of the prescriptions explained in Sec. 2.

Figure 10 shows the $\mathit{F}$-fingerprints and eigenvector of three PCs on the Si-Si and Fe-Si dimensions of La${}_2$Fe${}_{24}$Si${}_2$ (Case 2-3). The $\mathit{F}$-fingerprints on the Fe-Si dimension show very minor differences compared with those on the Si-Si one and the other $\mathit{F}$-fingerprints are much weakly distinguishable (not included in the figure). Among the 12 inequivalent configurations, four structures have two Si atoms located at the NN $96\mathit{i}$ sites and are in relatively high energy (+7.9 — +20.7 meV/atom) above the lowest energy in the configurational energies $\mathit{\Delta H}$ at $\mathit{y}=2$ in Fig. 9. The rest of the configurations have no NN Si pair and their energies are below +2.6 meV/atom. The PC $C_1$ of the structure map of the Euclidean distance in Fig. 11 clearly distinguish the four high-energy configurations and the rest of eight ones. The proportion of $\mathit{k}=3$ is $P_3=0.77$. The eigenvector of $C_1$ shown in Fig. 10(b) actually has the largest amplitude around the Si-Si NN distance. The PCs $C_2$ and $C_3$ have the largest amplitude of the eigenvectors around approximately 4 Å and 6 Å, respectively, and are related to the different stability within the $C_1$-clustering structures. The three structures (blue, cyan, and teal) with the NN Si-Si pair are quite similar in the map up to the PC $C_3$ and the PCs beyond must determine the minor energy difference. It is concluded that the interaction between the Si atoms are repulsive in La${}_2$Fe${}_{24}$Si${}_2$ and they may tend to be distributed apart.

Figure 10. (a) $\mathit{F}$-fingerprints and (b) Euclidean eigenvectors of three principal components of Si-Si (0-80) and Fe-Si (80-160) dimensions of La${}_2$Fe${}_{24}$Si${}_2$ (Case 2-3). Line colors used for the $\mathit{F}$-fingerprints in (a) correspond to those of dots in Fig 11. Red, lime, and blue lines in (b) denote the eigenvectors of $C_1$, $C_2$, and $C_3$, respectively.

Figure 11. Structure map of La${}_2$Fe${}_{24}$Si${}_2$ (Case 2-3) obtained by dimension reduction of the Euclidean distance. Dot colors denote the different configurations corresponding to the $\mathit{F}$-fingerprint shown in Fig 10(a).

Some selected $\mathit{F}$-fingerprints showing distinctive behavior and eigenvectors of three PCs on the Si-Si and Fe-Si dimensions of La${}_2$Fe${}_{23}$Si${}_3$ (Case 3-3) are drawn in Fig. 12. The corresponding structure map of the Euclidean distance with $\mathit{k}=3$ is shown in Fig. 13, resulting in the proportion $P_3=0.80$. The eigenvectors of the three PCs contain the structure information around the essentially same regions of the feature vector dimension as the previous case of La${}_2$Fe${}_{24}$Si${}_2$ (Case 2-3). In this case, there are four kinds of the $\mathit{F}$-fingerprint peak height at the NN Si-Si distance. The structures (magenta and lime) possessing the highest RDF at the NN distance has Si triangular trimers, forming a clustering with two more structures (small black) in the structure map. The structures (cyan, yellow, and green) with the second highest RDF have V-shape Si trimers, those (maroon, teal, red, blue, and black) with the third highest RDF do Si dimers, and those (purple and navy) have no NN Si pair. These locally distinct Si structures are clearly classified in the structure map along the PC $C_1$, also being in accordance with the general repulsive behavior between the Si atoms in the configurational stability in La(Fe,Si)${}_{13}$.

Figure 12. (a) Selected $\mathit{F}$-fingerprints and (b) eigenvectors of three principal components of Si-Si (0-80) and Fe-Si (80-160) dimensions of La${}_2$Fe${}_{23}$Si${}_3$ (Case 3-3). Red, lime, and blue lines in (b) denote the eigenvectors of $C_1$, $C_2$, and $C_3$, respectively.

Figure 13. Structure map of La${}_2$Fe${}_{23}$Si${}_3$ (Case 3-3) obtained by dimension reduction of the Euclidean distance. Dot colors denote the different configurations corresponding to the $\mathit{F}$-fingerprint shown in Fig. 12(a), where the $\mathit{F}$-fingerprints of the small black dots here are not included.

4. Summary

The three-step procedure to draw crystal structure map for a given set of materials systems with the same composition is proposed and demonstrated to classify the structure in a low-dimensional space and to model the materials properties. Combined with the scheme for counting all the possible combinations for alloy systems within supercell models, configurational energy preference can be discussed with the multiplicity of each configuration as the entropic effect. Further applications of the present methods with first-principles DFT calculations and machine-learning analysis are highly expected. In addition, regarding future prospects the present procedure may be applicable to any vector-type information such as spectral and image data for a given set of related materials systems. Map obtained in a low-dimensional space classifies the existing materials systems in terms of the reduced features and may possibly assist to explore desired candidate systems by using the projection and inverse-projection techniques as shown in Eqs (19(19) $Y_{+}=\left(\mathit{Y}-\overline{X}\right)u_{+}\mathrm{.}$(19) ) and (20(20) $\mathit{Z}=Z_{+}{u}_{+}^T+\overline{X}\mathrm{.}$(20) ).

Disclosure staement

The author declares no known competing financial interests.

Statement of Novelty

Crystal structure map constructed from the feature of structures by a dimension-reduction technique for classification and modeling for a group of related materials such as polymorphs, alloys, and doped systems. (30 words)

Acknowledgements

The author would like to thank Hiroyoshi Momida for providing the structure data of the amorphous Al${}_2$O${}_3$ model, Yukihiro Makino, Tetsuya Fukushima, and Kazunori Sato for invaluable discussion in the study of La(Fe,Si)${}_{13}$, and Hitoshi Fujii, Tomoki Yamashita, Hiori Kino, and Takashi Miyake for suggestive comments on data-science techniques.

Data availability

The raw data required to reproduce these findings are available by making an e-mail request to the corresponding author.

Additional information

Funding

This work is supported partly by JST-CREST (Grant No. JPMJCR22O2).

Funding

This work was supported by the Japan Science and Technology Corporation [JPMJCR22O2].