Order parameters (OPs), q, are mathematical constructs which aim to provide a numerical measure of the immediate local environment around an atom. The simplest OP, q_CN, is a coordination number (CN (Sprik, 1998)), (2)qCN=∑j≠iS∥pj−pi∥, obtained from counting neighbors within a cutoff radius r_cut,i around a central atom i. S is a weight that is 1 if an atom j is within a distance d_i,j < r_cut,i of atom i and S = 0 otherwise. The next level of sophistication is a distance-weighted approach such as using the Fermi function (Sprik, 1998): (3)qCN,Fermi=1exp[κ(di,j−rcut,i)]+1, where κ⁻¹ defines the transition width in which the contribution of an atom to the OP changes fastest to go from around 1 toward 0 as d_i,j increases (Sprik, 1998). Note that neither of these approaches provides information about how closely an environment resembles a given structural motif.

Bond-orientational order parameters, introduced by Steinhardt et al. (1983), can, to some extent, be used to discern different structural motifs (Mickel et al., 2013). Thermal motion and other small distortions (e.g., caused by relaxation to the ground state of a compound from an ideal initial prototype structure) can, however, yield overlap between the OP distributions so that identification of the motif type becomes difficult. Hence, there is a need to more reliably identify structural motifs with order parameters.

Peters (2009) and Zimmermann et al. (2015) have introduced order parameters based on pattern-matching ideas put forward by Shetty et al. (2002). The OPs specifically recognize body-centered cubic-like (Peters, 2009) as well as tetrahedral and octahedral environments (Zimmermann et al., 2015). The pattern-matching ansatz places a reward on local environments that are similar to the target structure, thus, resulting in a value, q_i, of (close to) 1 for perfect resemblance. When the surrounding atoms are not in a configuration resembling the perfect prototype motif, penalties force the order parameter to attain values around zero. The pattern matching is achieved by setting up a spherical coordinate system around a central atom i with a subset of neighboring atoms j and k (Figure 3). This allows the determination of the remaining neighbors’ (m, …) polar angles, θ, and azimuth angles, φ. If a neighbor is not located at angles that are commensurate with the expected positions in the underlying structure motif, the decline in reward for being away from the perfect position follows a Gaussian function. Conversely, a full reward is given if any expected remaining position is exactly assumed. The procedure provides rotationally invariant order parameters because it is applied to each neighbor j being used as the North pole position and any remaining neighbor k for defining the prime meridian (cf., Figure 3); all possible combinations are then averaged. Below, we provide the definitions of the OPs that we use in this work.

The tetrahedral order parameter, q_tet, is given by (Zimmermann et al., 2015): (4)qtet=1Nngh(Nngh−1)(Nngh−2)∑j≠kNnghexp−(θk−109.47∘)22Δθ2∑m≠j,kNnghcos2(1.5φ)exp−(θm−109.47∘)22Δθ2, where N_ngh denotes the number of neighbors bonded to the central atom i in a motif, θ_k is the polar angle formed between the bonds of neighboring atoms j and k with their mutually bonded atom i, φ is the azimuth angle between bond i − m with the plane spanned by i, j, and k, and Δθ = 12° a parameter controlling the reward loss for increasingly non-ideal positions, which was optimized to distinguish tetrahedral and octahedral environments in NaCl wurtzite and conventional rocksalt (Zimmermann et al., 2015). Note that we use a slightly different variant of q_tet as in the original formulation [cos²(1.5 φ) instead of cos(3 φ)] to avoid negative values.

The octahedral order parameter, q_oct, is given by (Zimmermann et al., 2015): (5)qoct=1Nngh[3+(Nngh−2)(Nngh−3)]∑j≠kNngh3H(θk−θthr)exp−(θk−180∘)22Δθ12+∑m≠j,kNnghH(θthr−θk)H(θthr−θm)cos2(2φ)exp−(θm−90∘)22Δθ22, where H(x) denotes the Heaviside function which is 1 if the argument x > 0 and 0 otherwise, Δθ₁ = 12°, Δθ₂ = 10°, and θ_thr is a threshold angle to distinguish second neighbors that are considered to be either in a “South pole” configuration or in a “prime meridian” position; we set this threshold to 160°. Note that we change the definition of q_oct in a similar manner as we have done for q_tet to avoid negative values.

The body-centered cubic order parameter, q_bcc, is given by (Peters, 2009): (6)qbcc=1Nngh[6+(Nngh−2)(Nngh−3)]∑j≠kNngh{6H(θk−θthr)exp−(θk−180∘)22Δθ12+H(θthr−θk)∑m≠j,kNnghcos(3φ)1.6θm−90∘Δθ2exp−(θm−90∘)22Δθ12sgn(θk−90∘)}, where Δθ₁ = 12°, Δθ₂ = 19.47°, and sgn(θ) is the signum function which is −1 for θ < 0, 0 if θ = 0, and 1 if θ > 0.

The mathematical definitions of the motif-specific order parameters q_tet, q_oct, and q_bcc in equations (4)–(6) follow a mutual recipe:

The innermost sum gives the contribution of how closely neighbor m is located at its expected position with respect to polar angle, θ_m, and azimuth angle, φ. For the azimuth angle, squared cosine functions are preferably used to pinpoint locations around a circle to ensure that the OP strictly gives positive values.

For q_oct and q_bcc, there is furthermore the need to distinguish whether or not neighbor k is in (approximate) South pole position via the Heaviside function H(θ_k − θ_thr). And, q_bcc also requires incorporating the alternating pattern of subsequent neighbors m being above and below the equator in bcc via the term 1.6 × sgn(θ_k − 90°) × (θ_m − 90°)/(Δθ₂).

Finally, we use the bond-orientational order parameters (Steinhardt et al., 1983) q₄ and q₆ for identifying close-packed motifs—fcc and hcp—(Ackland and Jones, 2006): (7)qi=1Nngh∑j=1NnghYimθ(pj),φ(pj), where Y_im are the spherical harmonics of degree i. Note that (i) the angles θ and φ are here with respect to a fixed frame of reference (Jungblut et al., 2013) and that (ii), while between 0 and 1, the values of the bond-orientational order parameters are typically not close to one for any motif, which is different to the behavior of q_tet, q_oct, and q_bcc. Despite the fact that the bond-orientational order parameters q₄ and q₆ are frequently used in nucleation studies involving bcc, fcc, and hcp environments (ten Wolde et al., 1996; Peters, 2009; Jungblut et al., 2013; Limmer and Chandler, 2013), they are not highly reliable indicators to distinguish between all of these motifs when distortions (thermal noise) are introduced (Gasser et al., 2003). For this reason, we extensively use here those order parameters that were specifically designed for a given motif type (e.g., q_tet for tetrahedra).

Motif recognition is typically achieved on the basis of a threshold approach (Peters, 2009; Zimmermann et al., 2015): if a given order parameter, q_i, is larger than an appropriate threshold, q_i,thr, the coordination pattern is confirmed. Because of the design of these OPs, a threshold of 0.5 is often a reasonable a priori choice. The motif-specific order parameters q_tet, q_oct, and q_bcc should then be ideally usable as stand-alone identifiers for tetrahedral, octahedral, and bcc-like coordination environments, respectively. However, Table 1 indicates that such an approach to defining criteria for motif recognition is not effective. For example, a site in a diamond structure would be identified as having both tetrahedral and bcc-like coordination; this is not surprising because the bcc motif can be viewed as two point-symmetric tetrahedra. Therefore, slightly more complex criteria must be developed to accurately distinguish structure motifs on the basis of order parameter values. We start with following set of criteria that allows us to identify all motifs separately and unambiguously for perfect prototype structures: (8)qtet>0.5tetrahedral, (9)qoct>0.5octahedral, (10)qbcc>0.5andqtet<0.5bcc, (11)q6>0.5andqtet,qoct,qbcc<0.5fcc, (12)q6<0.5andqtet,qoct,qbcc<0.5hcp.

The fourth-order bond-orientational order parameter q₄ is hence not necessary to identify all motifs. However, as we will explain later, small modifications to these criteria are needed to more accurately distinguish non-ideal motifs. In particular, we will merge the fcc and hcp criteria into a single close-packed rule (q₆ > 0.4 and q_tet, q_oct, q_bcc < 0.4), and we will decrease the threshold of the tetrahedral order parameter to 0.3.

Validation is essential for reliable quantitative structure–property relationships (QSPRs) (Tropsha et al., 2003). We follow a three-step hierarchical approach for our structure-motif assessment on the basis of order parameters using perfect prototype structures. First, the responses of the tetrahedral and octahedral order parameters are measured when a single neighbor in the corresponding prototype structure is subjected to defined perturbations. We focus here on q_tet and q_oct because those motifs are particularly important throughout materials science and in corresponding design rules. Second, we randomly, but systematically, perturb the locations of sites in all prototype structures to mimic the effect of small distortions in equilibrium structures on all order parameters. This provides a more detailed insight into the sensitivity of order parameters to motif distortion, and it provides the necessary data for the next validation level. Third, we calculate motif-recognition likelihoods based on the histograms of the second validation level, which provides a first reasonable value of the tolerance, δ, that is later used for neighbor finding in materials from a database.

In Figure 4, we display the response of the tetrahedral (Figure 4A) and octahedral (Figure 4B) order parameter to perturbations of a single tagged neighbor in the respective perfect structure motif along the azimuth (orange) and polar angle (blue) as well as along an intermediate path (purple). The order parameters decline from 1 almost linearly to approximately 0.37 and 0.8 for polar and azimuth angular perturbations up to 27° and 18° for q_tet and q_oct, respectively. The response along the intermediate perturbation path is similar. This indicates that the OPs can be used as (linear) measures of tetrahedrality and octahedrality for non-negligible perturbations of single neighbors, making the tetrahedral OP, for example, also attractive for structure analysis of liquid and solid water phases (Tao et al., 2017). Figure 4 indicates that configurations with two atoms nearly overlapping score relatively high, which, however, does not pose a problem because those arrangements are very unlikely (i.e., two atoms will not be present at such close distances).

Next, we apply Gaussian-distributed random perturbations to all sites in the prototype structures via the polar form of Box–Muller transforms (Box and Muller, 1958) as implemented in numpy (van der Walt et al., 2011). The perturbations resemble the spatial distribution behavior of atoms in Einstein crystals (Einstein, 1906; Frenkel and Ladd, 1984; Aragones et al., 2012), the procedure of which we therefore call Einstein crystal test rig. Each lattice atom oscillates independently around its equilibrium position with a distribution width of σ_ES. We refer the reader to the Supplementary Material (Section 1.1 in Supplementary Material) for more details on the Einstein crystal calculations. Figure 5 indicates that the OPs respond to Einsteinian perturbations of magnitude σ_ES/d_id = 0.02 and a relative cutoff radius of r_cut/d_id = 1.1 with relative frequency distributions of finite width. The distributions confirm the validity of both: (i) our proposed structure motif-recognition criteria and (ii) that q_tet and q_oct are in particular well-behaved measures for the degree of a given motif. The data, however, also suggest that the ability to distinguish between fcc and hcp on the basis of q₆ will be limited. Therefore, we redefine those structure motif criteria (equations (11) and (12)) to provide a single criterion for close-packed motifs (cp = fcc + hcp): (13)q6>0.4andqtet,qoct,qbcc<0.4fcc+hcp.

In this context, the approach by Honeycutt and Andersen (1987) is worthwhile noting, which compresses information about local environments—specifically, the number of shared near neighbors of a pair of atoms and the connectivity among the shared neighbors—into a four-integer index. This index can, however, not easily be used as an automatic motif recognition tool for hcp–fcc distinction because of two reasons. First, it requires visual inspection of the index-underlying graphs because the fourth integer is an arbitrary enumeration (cf., indices 1,421 and 1,422 occurring in fcc and hcp). Second, the index is computed for pairs of atoms so that the index itself cannot be used to directly characterize the entire coordination environment of a single atom. We are currently working toward solving the hcp–fcc distinction problem via definition of additional order parameters, resulting in an order parameter feature vector that might enable distinction between hcp and fcc.

To quantitatively assess our structure-motif recognition capabilities, we systematically expand the Einstein crystal sensitivity test approach to various distortion degrees, σ_ES, and relative cutoff radii, r_cut/d_id. For this purpose, we use following basic likelihood function (Sivia, 2012), L: (14)L=∏i=1Nstr∏j=1Nmotδi,jKrNi,j+(1−δi,jKr)(1−Ni,j)Nsamp, where N_str and N_mot are the number of tested prototype structures and motifs, respectively (both are 4 here: tet, oct, bcc, and cp), N_i,j is the number of times that structure i is recognized as consisting of motif j of all N_samp samples for i and j (here: N_samp = 1,000), and δi,jKr is the Kronecker delta function, which is 1 if i = j and 0 otherwise. Since we have merged fcc and hcp to one cp motif, we average data from fcc and hcp evaluations so that each of the four regrouped motifs (tet, oct, bcc, and cp) has an equal weight on L. The likelihood function thus represents the joint probability to do both: to correctly identify true motifs and to correctly reject false motifs. Furthermore, the particular form of L as products of the separate motif and structure likelihoods more stringently requires that all motifs in all structures are correctly identified and rejected, respectively, to an acceptable degree. A mere arithmetical mean would favor compensation effects in which high recognition scores achieved with one motif and/or structure would balance entirely unsuccessful recognitions with low scores.

There exists a very localized optimal region of parameters where L is maximal (close to 1), which is found in the vicinity of small perturbations and small relative cutoff radii (Figure 6A). We expected a more broad distribution at vanishing degrees of distortion, which is, in fact, observable once the results from the bcc prototype structure are removed (Figure 6B). This is a reflection of the well-known issue that, in the bcc structure, second nearest neighbors are close to nearest neighbors (cf., Mickel et al. (2013) and Figure 5). As a result, the bcc issue sets certain limits to using a global tolerance for successful neighbor finding.

Before turning to the results from the entire database, we determine the most suitable neighbor-finding method and the optimal tolerance parameter, δ, given a reasonable test (sub)set from the MP database. Our initial test set, for which the number of different motif sites, N_i, is well known, consists of materials that are members of the structure groups listed in Table 2. Materials belonging to a given structure group were found by scanning the MP database for structures that are similar to a reference structure. Similarity between structures is hereby determined with a structure matcher algorithm implemented in pymatgen (Ong et al., 2013) using default parameters, except for turning off species matching (i.e., we match the frameworks of the structures). The references were ideal prototype structures in the case of unary materials (diamond-like, simple cubic, bcc, fcc, and hcp) and the canonical structure in the MP database for all other materials (zinc blende: mp-10695; rocksalt: mp-22862; CsCl-like: mp-22851; MgAl₂O₄-spinel: mp-3536), respectively. For the resulting 1,025 test structures, we calculate the order parameter values q_CN, q_tet, q_oct, q_bcc, and q₆ of all N_sites sites in each structure. Then, we determine the numbers of different structure motifs (N_tet, N_oct, N_bcc, and N_cp) in each material on the basis of our recognition criteria (equations (8)–(10) and (13)) as well as the number of times that a site was assigned to more than a single structure motif (N_multi).

The different neighbor-finding settings are compared by averaging the fraction, p_rec,i, of correctly recognized number of expected motifs j, in each structure i, N_j,i, (15)prec,i=Nj,i/Nsites,i, overall 1,025 test structures. Note that we put a very stringent criterion on multiple assignments: if a structure contains a single site that is assigned to different motifs, the structure has no positive contribution on the recognition fraction (i.e., p_rec,i = 0). Furthermore, we address the problem that, in spinels, the target number of tetrahedral and octahedral sites, respectively, can be overpredicted by using following functional form for prec,ispinel: (16)prec,ispinel=1−13−Ntet,i/3Nsites,i/7+23−Noct,i/3Nsites,i/7.

The results in Figure 7 indicate that all three neighbor-finding approaches yield high recognition rates (>85%) for the chosen test set if the respective tolerance parameter is small enough (≪0.1). Furthermore, the prediction quality decreases precipitously when the neighbor-finding tolerance, δ, reaches a value of around 0.1 for any of the three methods. The best performing method p¯rec=100% is the minimum distance-based approach with 0.03 ≤ δ ≤ 0.08. For all following analyses, we use this method where we, however, employ a slightly larger tolerance (δ = 0.1) because this tolerance yields a very low average coordination number of 3 for all sites in the entire Materials Project database.

Apart from being effective, the minimum distance-based neighbor-finding method is also computationally exceptionally efficient. Calculating all order parameters requires only 0.026 s per site on a compute node of NERSC’s² Edison cluster (time was averaged over all sites in the MP database). That is, analyzing a structure with 100 sites should take 2–3 s, assuming the entire procedure roughly scales linearly with the number of sites. Note (i) that we decrease the tetrahedral OP threshold to 0.3 because tests on Jahn–Teller active structures suggest 0.5 being too strict and (ii) that we also require the coordination number, q_CN, to match the ideal motif value. Thus, the final set of rules that we use for motif determination across the entire Materials Project database are: (17)qCN=4andqtet>0.3tetrahedral, (18)qCN=6andqoct>0.5octahedral, (19)qCN=8andqbcc>0.5andqtet<0.5bcc, (20)qCN=12andq6>0.4andqi<0.4fcc+hcp.

In Algorithm 1, we also provide a pseudocode implementation of our motif recognition method.

The MP database screening results are presented in Figure 8A as relative motif occurrences for each element separately, and they agree with several previous observations (Cotton and Wilkinson, 1980; Brown, 1988). For example, Li occurs similarly often as a tetrahedral site as it occurs as an octahedral site, but bcc-like Li motifs are rare and close-packed ones are absent. These motif frequencies follow trends of Li for 4-, 6-, 8-, and 12-fold coordination that were already established by Brown (1988). Furthermore, we see comparably good agreement between motif occurrence and Brown’s coordination statistics for Na, Ca, Ba, Mn, Fe, Co, B, As, and La. Examples with some differences but still overall favorable agreement include K, Sr, Pb, Ni, Cd, Sb, and Bi, whereas we obtain very different dominant motif vs coordination number distributions for Tl, Hg, Al, and In. In light of these differences, we note that our material set is far larger (and more diverse) than the set that Brown used; there is a 100-fold difference between the number of sites that we have analyzed (>15,000,000) in comparison to the previous study by Brown (1988) (14,000).

To continue the discussion of how our structure motif data relate to results from literature and common wisdom we note that (i) S and Se occur in tetrahedral coordination, (ii) S, Se, and Te can be seen in octahedra, but (iii) 8-fold coordination only occurs for Te and (iv) close-packed motifs of S, Se, and Te are not known (Cotton and Wilkinson, 1980). Silicon is known as a strong tetrahedron-former (e.g., in zeolites (Baerlocher et al., 2007)). This is confirmed when comparing both panels in Figure 8; Figure 8B gives the frequency at which a motif (here: tet) occurs relative to all motifs—recognizable and unrecognizable—with the same coordination number (here: 4). Figure 8B indicates that 86% of all 4-fold coordinated Si motifs are tetrahedra. Surprisingly, the ratio at which lithium occurs in 4-fold coordination as a tetrahedron is low (47%). By contrast, the low tetrahedral coordination fraction of Ag and Cu motifs (both 41%) was expected because these elements are frequently found in square-planar environments.

An intriguing observation for us is that the known decline in tetrahedral site preference (Navrotsky and Kleppa, 1967; Burdett et al., 1982; Rong et al., 2015) as we go from Li over Mg to Ca in spinel is reflected in the MP database as a whole. The relative occurrence of Li as a tetrahedron among all 4 recognizable motifs is is 42%, Mg 20%, and Ca 13%. Zn typically takes a place between Li and Mg in spinels (Rong et al., 2015), which is, however, different in the MP database (46%).

The unexpected motif prevalences that are observed for some elements may hint at the fact that our neighbor finding method could benefit from further study and improvement. However, the overall approach represents a fast, fully automatic method to determine coordination environment motifs over large crystal databases that is relatively trustworthy and intuitive.