2.1. Process

For each .cif file provided as an input by the user, Polynator traverses the process outlined in Fig. 3. We will proceed by briefly explaining each step.

Figure 3 Simplified flow chart of the working principle. Multiple crystal structures and multiple atom arrangements per crystal structure can be processed in one run, not reflected here for simplicity.

Step 1. The input .cif file is read and the relevant atom vectors are generated according to symmetry instructions.

Step 2. Atom arrangements to be fitted can be classified into coordination environments and molecules. A coordination environment may be generated around a selected central atom (or dummy atom coordinates) following one of several methods. The most basic method is to include all atoms whose distance from the central atom is within a user-specified range. Alternatively, the user may assign atomic radii to each atom. It is furthermore possible to specify a maximal coordination number; excess ligands with the highest distance from the central atom are excluded first. Coordination environments can also be constructed by a simple algorithm, finding gaps in consecutive center–ligand distances. Last but not least, an algorithm based on the Voronoi polyhedron around the central atom is available to the user. It considers an atom a ligand if the solid angle subtended by its associated Voronoi face from the perspective of the central atom is greater than a given threshold (20° by default). This last algorithm is loosely based on the definition by O'Keeffe (1979) for the coordination number and very similar to the method used by ChemEnv (Waroquiers et al., 2020). Since there is no single unbiased method to determine which atoms should be counted as ligands, Polynator leaves the choice to the user.

Molecules are constructed on the basis of their connectivity according to user-specified or default atomic radii. Subsections (e.g. ethane units, icosahedra, porphyrin rings, fullerenes) can automatically be extracted from molecules or infinite frameworks by a number of predefined or user-specified topological strategies. This does not require any pre-defined bonds, just a regular .cif file and the name of the desired topological shape.

Step 3. A model prototype with the appropriate number of vertices for a given atom arrangement is constructed. Model vertices are stored not as Cartesian coordinates but as a compressed set of parameters directing their construction and manipulation. These parameters do not affect single vertices, but entire symmetry-aligned groups of them, referred to as belts. These belts are orthogonal to and are stacked along a model axis. For example, as illustrated in Fig. 4, an anti­cuboctahedron (or rigid triangular orthobicupola) is split into three belts. The information associated with each belt is listed in Table 2, namely the number of vertices it contains, the basic cylindrical coordinates for the entire belt (width, height, torsion angle) and any free parameters (in this case only the same scaling parameter for each belt). The individual vertices are constructed by rotating them around the model axis by multiples of 2π/n, where n is the number of vertices in that belt. Fig. 5 illustrates the deformations allowed by the free parameters of a dynamic twisted trigonal prism.

Table 2 Widths, heights, torsion angles and free parameters for the construction of an anticuboctahedron Belt n vertices Width coefficient Height coefficient Torsion angle Free parameters 1 3 0 sc1 2 6 1 0 π/6 sc1 3 3 0 sc1

Figure 4 Anticuboctahedron (also known as a trigonal orthobicupola), split into belts along its 6 axis, which is used as the model axis. The individual model vertices within a belt are labeled a–f.

Figure 5 Twisted trigonal prism with three free parameters, namely h1 for its height, w1 for its width and φ1 for the torsion angle between the bases. In general, a model may have any number of parameters.

Step 4. Identifying the optimal pairing scheme of model vertices with their real counterparts found in the crystal structure is sometimes easy to do for humans but deceptively difficult for computer programs. The brute force approach of checking every possible permutation of model–real vertex pairs is the only way to guarantee the optimal result. However, the number of required evaluations grows factorially with the number of vertices, quickly rendering this approach infeasible for all but the lowest coordination numbers. As explained by Waroquiers et al. (2020), the number of evaluations can be reduced slightly by taking into account the symmetry of the model. Unfortunately, this does not ultimately solve the problem of unmanageable numbers of vertex permutations. The CoSyM program takes advantage of the topology of the vertex figure (Alon & Tuvi-Arad, 2018) to achieve a significant reduction in the number of permutations. Although this approach is very well suited for molecular structures due to their rather predictable bond lengths and comparatively low connectivities, it can run into problems when dealing with more irregular coordination polyhedra. Instead, depending on the model, Polynator chooses one out of several strategies to find the best fit without checking many combinations. When fitting a coordination polyhedron or cage, it is assumed that the real ligands are spread out more or less equally on a spherical shell around their centroid. The atom arrangement is sequentially scanned along the directions of 92 vectors pointing from the centroid to equally spread out locations on a unit sphere surface. It is then estimated how well the scanned atom arrangement agrees with the model for each of these directions. To do this, for each scanning step, the ligand vectors are first filled into the model belts according to their dot product with the scanning direction vector. The approximate model axis is then constructed as the eigenvector with the largest eigenvalue of a covariance matrix formed from the centroids of all belts. The fit quality is estimated according to the distance of each belt centroid from the model axis and the difference in dot product with the model axis between each vertex and the respective belt centroid. A small number of assignment schemes with the best estimated agreement qualify for full evaluation.

When handling strongly prolate models, such as the bi­phenyl skeleton, a very simple strategy of finding the belt assignment scheme is employed instead. The model axis, or the `longest' axis of the atom arrangement, is constructed as the eigenvector with the largest eigenvalue of a covariance matrix formed from all atom vectors. The vectors are ranked according to their dot product with that axis and filled into the belts accordingly. For multi-belt oblate models, such as the porphyrin skeleton, the atom vectors are simply ranked by their distance from the centroid. The performance and accuracy of the assignments is further improved by the possibility for models to inherit their assignment from a less general version of themselves, provided that version is a good fit. For example, a trigonal antiprism may effectively copy the assignment of a previously fitted regular hexagon or octahedron. However, when copying from an octahedron, four different orientations must still be fitted. Though none of these strategies are technically guaranteed to find the best assignment scheme, they appear to be very reliable as long as the atom arrangement resembles the model to some degree.

Step 5. After the atom vectors have been assigned to belts, they have to be paired with the model vertices within their belt. This is done by rotating the model so that one of its vertices has the same torsion angle around the model axis as one of the atom vectors in the same belt. This torsion angle is set to 0° and the clockwise angles of the other vectors are measured with respect to it. The vectors within each belt are then sorted according to their angle and paired up according to this order. The agreement between model and real atoms is estimated according to the individual angular deviations between paired vectors, weighted by the distance from the model axis. The final pairing scheme is determined by repeating this process for each atom vector and selecting for the best agreement.

Step 6. Once the vectors are paired up, the model orientation can be adjusted to optimally fit the real coordination environment. This task is carried out by a Kabsch algorithm (Kabsch, 1976).

Step 7. Subsequently, the values of free shape parameters defined individually for each model are determined via least-squares fits. In the majority of cases, this can be done analytically by simply measuring averages of vector components (height, width, magnitude). Otherwise, the optimal value must be determined iteratively. Since the optimal solutions for shape and orientation are generally interdependent, the sequence of steps 6 and 7 has to be executed several times until a convergence criterion is reached (the δ value ceases to improve by more than 10⁻⁵ between cycles).

3.1. Shape recognition and classification

Since Polynator is capable of automatically recognizing polyhedra and quantifying their similarity to a wide range of models, a very basic use case is the identification of `exotic' polyhedra, such as the elongated gyrobifastigium. This is found as a coordination polyhedron of the Th site in the α-ThSi₂ structure type, but has not been identified as such in the publications we are aware of. This elongated gyrobifastigium is distinct from the regular gyrobifastigium, which can be found in the same structure type (Alvarez, 2017), as depicted in Fig. 6. The high-pressure CaFeTi₂O₆ phase, as described by Leinenweber & Parise (1995), contains both the regular gyrobifastigium and another somewhat uncommon Johnson polyhedron, the sphenocorona. Polynator also natively includes over 40 model polyhedra with more than 12 vertices, many of which are represented in intermetallic phases (Zintl & Haucke, 1938, 1940; Kasper et al., 1965; Ban & Sikirica, 1965), zeolites (Broussard & Shoemaker, 1960) or boron clusters (Popov et al., 2015).

Figure 6 Both the regular gyrobifastigium and the elongated gyrobifastigium are found in the crystal structure of the ThSi2 type.

Another simple application is finding the best description of a given `irregular' coordination polyhedron. For example, in the scheelite-type Bi[VO₄] compound (Mariathasan et al., 1986), bis­muth is coordinated by eight oxygen atoms in an arrangement that could be understood as a distorted version of either a triangular dodecahedron, a bicapped trigonal prism or a tetragonal antiprism. As illustrated in Fig. 7, when comparing the fit for rigid versions of all three models, it becomes clear that the triangular dodecahedron is the best description, although not by a very large margin. These differences become more apparent when comparing dynamic versions of these models, as the distortion of the optimal triangular dodecahedron amounts to less than a quarter of the distortion of the optimal versions of the other candidates. This comparison between dynamic models is also arguably more in line with the intuitive understanding of a given polyhedron. For example, most coordination polyhedra commonly identified as triangular dodecahedra are not as prolate as the equilateral version, but more or less spherical.

Figure 7 Coordination environment of bis­muth in Bi[VO4], viewed from three different perspectives, corresponding to the optimal orientation of an equilateral bicapped trigonal prism (left), an equilateral tetragonal antiprism (center) and an equilateral triangular dodecahedron (right). Distortion values quantify the fit of the aforementioned rigid equilateral models, as well as dynamic models which share the basic structure of the rigid models, but are otherwise only constrained by the respective point group.

3.2. Measuring characteristic parameters

Polynator can be used as a way to obtain exact values of parameters which can otherwise be somewhat challenging to measure correctly. A good example is the dihedral angle between the phenyl rings in a bi­phenyl unit, as pictured in Fig. 8. Though it is easy to get an approximate value by calculating the dihedral angle between an atom in the first ring and an atom in the second ring, using the bridging atoms as the hinge axis, this approach falls short if more precision is desired. Polynator provides a clean value for this angle as a by-product when finding the closest-matching ideal molecule. The model for a bi­phenyl skeleton has three free parameters: one for the radius of the phenyl rings, one for the length of the bridging bond and one for their dihedral angle. These values are accessible in the output files.

Figure 8 Overlay of real molecule (black) and fitted model (blue) for bi­phenyl skeletons cut from two different molecules. The C24H30N4O4 molecule (top) (Holý et al., 2001) has somewhat bulky substituents on all ortho positions, hence the rather large dihedral angle 2φ. The bi­phenyl skeleton (bottom) is part of a fluorene molecule (Burns & Iball, 1955) with point group m. Though the latter molecule is not perfectly flat, its 2φ value is still exactly 0°. The considerable distortion is due to the skewed phenyl rings.

3.3. Quantification of distortions

Polynator is suited just as well for more intricate investigations, such as uncovering trends and patterns in the coordination behavior of a given structure type. Since any number of input files, atom arrangements and models can be processed at once and the process is quite fast (typically 10–20 models per second), large numbers of similar structures can be compared almost effortlessly. This is also facilitated by .csv output files listing important results in table form and the topological package allowing for the automatic extraction of specific structural units from large molecules or frameworks.

Because of its rich structural diversity, the perovskite family lends itself well to the analysis of its coordination environments. A good example is provided by the eponymous Ca(TiO₃) compound which, on heating to 1600 K, undergoes two phase transitions on its way from the Pbnm room-temperature modification via I4/mcm to the Pm3m aristotype structure. Fig. 9 shows that the distortion of the cuboctahedral coordination polyhedron of the Ca site decreases with increasing temperature. This is true with respect to the scaleable but otherwise rigid Archimedean model, but also with respect to dynamic models which represent lower-symmetry versions of the cuboctahedron (see Fig. 10). These are automatically deformed within the constraints of their point group in order to achieve the best fitting version. In the example at hand, the 4/mmm model starts out with a slightly better fit than the Archimedean model, but this difference all but disappears with the first-order phase transition to the I4/mcm structure. This implies that the distortion of the cuboctahedron is at first in part due to an elongation or compression along one of its fourfold axes, which then vanishes following the phase transition. The 42m model fits best from the start and seems to improve roughly at the same rate as the higher-symmetry models. However, after the phase transition, its deviation from the real coordination environment falls to zero, indicating that it matches the point group of the crystallographic site. It follows that a model having all the degrees of freedom permitted by the 1 point group would not be a better fit for the cuboctahedra in the I4/mcm structure than the 4/mmm model, as the remaining distortion of the latter is only due to the constraint of inversion symmetry imposed on it. The ratio of distortions relative to these point groups observed in Ca(TiO₃) turns out to be overall typical behavior for perovskites. To illustrate this, distorted cuboctahedron models with the point groups 43m, m3, 4/mmm, 3m, 422, 4mm, 42m, 4/m, mmm (two models), 32, 3m and 3 were fitted to the cuboctahedra found in a large selection of structures from four major tilted perovskite structure types and their elpasolite counterparts. The results are displayed in Fig. 11 and Table 3.

Figure 9 Temperature-dependent distortions of CaO12 cuboctahedra in Ca(TiO3) with respect to three model polyhedra which conform to different point groups. Experimental data according to a single sample powder X-ray diffraction measurement series (Yashima & Ali, 2009).

Figure 10 Tree diagram of distorted cuboctahedra. Higher positions in the diagram indicate higher symmetry. Note that there are two differently distorted cuboctahedra with the point group 42m. Model 42m (1) is a combination of a bidisphenoid and a square, whereas 42m (2) can be described as a combination of an anticuboid and a disphenoid. Each is viewed along its 4 axis. Similarly, model mmm (1) is a combination of a cuboid and a rhombus, whereas mmm (2) is a combination of a diagonal orthobicupola and a rectangle.

Figure 11 Various models fitted to the cuboctahedra found in four perovskite structure types and their respective elpasolite counterparts. Each structure type is represented by an array of crystal structures belonging to it. Each individual distortion δ was first divided by the distortion of the Archimedean cuboctahedron model fitted to the same coordination environment. Thus the arithmetic mean of all values generated was then taken for each combination of a model and a structure type in order to obtain a data point shown in the diagram. See the supporting information for a discussion of the methodology applied in this plot.

Several observations can be made. Firstly, while the 43m cuboctahedron model fits the cuboctahedra in the Fm3m elpasolite with untilted octahedra exactly, it fits only slightly better than the Archimedean m3m cuboctahedron model for each of the tilted elpasolite structure types investigated. Therefore, in most cases, the distortion caused by the elpasolite-type superstructure is much smaller than the distortion caused by the octahedral tilting. The difference is even more pronounced with lower space-group symmetries. This result is to be expected, as the distortion by tilting increases with decreasing space-group symmetry, while the difference in octahedron sizes should stay roughly the same on average. The maximal hettohedra with point groups 4/mmm and 3m do not fit noticeably better than the Archimedean cuboctahedron, except in the case of the Pbnm and P2₁/n structure types, where the 4/mmm model shows some relative improvement. This suggests that there is generally no significant elongation or compression along the three- and fourfold (pseudo-)axes for the higher-symmetry structure types, as both would improve these fits. However, the Pbnm and P2₁/n structure types, which represent the a⁻a⁻c⁺ tilting scheme (Glazer, 1972), adopt an elongation along one of the fourfold pseudo-axes, more precisely one of the axes along which the tilting has alternating rotational directions. Note that the elongation does not stem from an elongation of the octahedra, as compounds with octahedral distortions δ > 3.2 were filtered out and the majority of octahedra in the admitted structures were better described as trigonal antiprisms than as tetragonal bipyramids.

The overall best fits are, in all cases, achieved by the 32 and 42m (1) models. These perfectly match the site symmetry for the I4/mcm (a⁰a⁰c⁻) and R3c (a⁻a⁻a⁻) structure types, respectively. With the elpasolite counterparts I4/m and R3, the site symmetry is reduced to 4 and 3, respectively. This goes along with a significant deterioration of the fit for these models, which far outweighs the concomitant relative improvement of the 43m model fit.

Somewhat interestingly, the 32 model fits better than the 42m model for the Ibmm and related I2/m structure types. This is despite the fact that the cuboctahedral site in the Ibmm structure type has mm2 symmetry, which is a maximal subgroup of 42m, but not closely related to 32.