Comparing crystal structures with symmetry and geometry

Measuring the similarity between two arbitrary crystal structures is a common challenge in crystallography and materials science. Although there are an infinite number of ways to mathematically relate two crystal structures, only a few are physically meaningful. Here we introduce both a geometry-based and a symmetry-adapted similarity metric to compare crystal structures. Using crystal symmetry and combinatorial optimization we describe an algorithm to arrive at the structural relationship that minimizes these similarity metrics across all possible maps between any pair of crystal structures. The approach makes it possible to (i) identify pairs of crystal structures that are identical, (ii) quantitatively measure the similarity between crystal structures, and (iii) find and rank structural transformation pathways between any pair of crystal structures. We discuss the advantages of using the symmetry-adapted cost metric over the geometric cost. Finally, we show that all known structural transformation pathways between common crystal structures are recovered with the mapping algorithm. The methodology presented in this study will be of value to efforts that seek to catalogue crystal structures, identify structural transformation pathways or prune large first-principles datasets used to parameterize on-lattice Hamiltonians.


INTRODUCTION
The crystal structure of a material determines much of its functional behavior, including its electronic, ionic, and mechanical properties. The number of possible crystal structures is infinite, and a large number of crystal structures have already been characterized 1 . This number will undoubtedly continue to grow as researchers expand into the vast, high-dimensional composition spaces of multiprinciple element alloys and functional compounds.
Although great effort has been expended to cataloging 2-4 and categorizing unique crystal structures and their chemical manifestations, challenges remain to realizing a universal classification scheme. Such a scheme is becoming increasingly desirable with the growth of first-principles databases, such as the Materials Project 5 , AFLOW 6 , OQMD 7 , NOMAD 8 , and the Materials Cloud 9 . Essential to a successful classification scheme is a robust method with which to measure the structural similarity of two crystals. Unfortunately, the comparison of two crystal structures is made challenging due to their periodic nature, which allows for an infinite number of mathematical representations. Symmetry can be used to help classify crystal structures. However, while there are only 230 crystallographic space groups, two crystal structures belonging to the same space group can nevertheless be very different from the point of view of local connectivity, unit cell size and number of distinct sites within the unit cell.
Pathways that connect different crystal structures are also of great importance. Many crystalline materials undergo structural transformations as a function of temperature and pressure [10][11][12][13][14][15][16][17][18] . Only a limited number of continuous pathways connecting common crystal structures have been identified [19][20][21][22][23][24][25][26][27] . Knowledge of these pathways, however, is often invaluable to understand and manipulate phase transformations in shape-memory alloys and actuator applications. The identification of pathways requires an enumeration of possible mappings of one crystal onto another 24 and a robust definition of a distance metric between crystal structures to enable a ranking of such mappings.
Several methods and accompanying software packages have been developed to compare and classify crystal structures. These approaches tend to rely on structure standardization, space group and Wyckoff position identification and comparisons of crystal invariants [28][29][30][31] . Other methods rely on comparisons of pair distribution functions 32 . Some of the well known or recently released software packages include STRUCTURE-TIDY 28 , CRY-COM 30 , CMPZ 33 , XTAL-COMP 34 , COMPSTRU 31 , SPAP 35 , the pymatgen Structure Matcher 36 , and AFLOW-XtalFinder 37 . While they all have their merits, none determine the similarity between crystal structures using combinatorial mapping approaches with a rigorous quantification of strain and displacement costs that account for symmetry. Furthermore, there are no software tools generally available with which to systematically enumerate and rank mappings between pairs of crystal structures.
In this paper, we describe an algorithm implemented in the CASM software package 15,38,39 to enumerate mappings between pairs of crystal structures. As part of the mapping approach, we introduce a geometric and symmetry-breaking distance metric that measures the similarity between two crystal structures. The tools allow users to answer three questions. Given two input crystal structures: • Are they identical up to a symmetry-preserving scaling factor?
• If the two crystal structures are not identical, how similar are they?
• What are ranked pathways between the two crystal structures and what is the orientation relationship between these crystal structures for a particular path?
Answering these questions requires (i) a robust way to map one crystal onto another and (ii) metrics of similarity between any pair of crystal structures. In the next sections we describe the mathematical underpinning of the crystal mapping tool implemented in the CASM software package.

RESULTS
The objective of this study is to obtain a relationship between the spatial coordinates of atoms in two crystals, C 1 and C 2 . The approach relies on finding a distortion of crystal C 2 such that, after deformation, the atoms of C 2 coincide exactly with sites of C 1 . In this context, we refer to C 1 as the parent crystal and C 2 as the child crystal. There are many possible ways to deform C 2 such that it maps onto C 1 . We denote by M the set of all possible ways to relate C 2 to C 1 and refer to a particular element M m as a mapping of C 2 onto C 1 . Different mappings M m will require different degrees of deformation.
Crystal structures are represented by lattice vectors and a list of coordinates that correspond to the sites in the unit cell. We may compactly specify the lattice vectors as columns of a 3 × 3 matrix The fractional ( f ! ) and cartesian ( r ! ) coordinates of an arbitrary point in space are related to each other as r ! ¼ L f ! . A particular mapping of C 2 onto C 1 comprises two separate operations. First, in the lattice-mapping step, C 2 is rigidly rotated and homogeneously deformed until a shared lattice is obtained that can be used to define both structures simultaneously. Second, in the atomic mapping step, each atom of C 2 is assigned a companion in C 1 and displaced so that it is coincident with its companion. In the next section we describe how the strain and displacement fields that relate these crystal structures may be obtained.

Lattice mapping
We will first consider the mapping of the lattice of C 2 , denoted L 2 , onto the lattice of C 1 , denoted L 1 . Our lattice-mapping approach is similar to that of Trinkle et al. 24 . In the simplest case, the two sets of lattice vectors are related via matrix multiplication according to The deformation gradient tensor, F, which relates the two lattices, can be uniquely decomposed as the product of a reorientation matrix, Q, and the left stretch tensor, V, which is symmetric and has positive determinant. We define a strain tensor in terms of the left stretch tensor as where I is the 3 × 3 identity matrix. The strain tensor B (sometimes called the Biot strain) is also symmetric. At small deformation, the left stretch tensor is close in value to the identity matrix, and B approaches zero. We note that the application of a rotation followed by a left stretch tensor to the child lattice until it maps onto the parent lattice is equivalent to an application of a right stretch tensor followed by a rotation to the parent lattice. The matrix L 1 is one of many matrices that generates the same lattice of crystal sites in the parent crystal. Any other matrix L ðNÞ 1 ¼ L 1 N generates the same lattice as L 1 if N is an integer matrix having determinant ± 1. Such matrices are referred to as unimodular matrices. A more general equation relating L 2 to L 1 is thus Given L 1 and L 2 , there are many possible lattice mappings, corresponding to different choices of N, and each results in a distinct left stretch tensor and reorientation matrix. Thus, B (N) and Q (N) can both be considered as functions of the integer matrix N. This is illustrated in Fig. 1 in two dimensions where the parent C 1 has a rectangular unit cell and where the child crystal C 2 is simply a rotated version of C 1 . Two possible mappings are shown. In the first mapping, the short axis of C 2 is mapped onto the short axis of C 1 . In this case, the mapping involves only a rigid rotation with a zero strain cost. A unit sphere in C 1 remains a unit sphere in C 2 . In the second mapping, the long axis of C 2 is mapped onto the short axis of C 1 . This mapping requires not only a rotation, but also a deformation, as illustrated by the transformation of the circle in C 1 into an ellipsoid in C 2 . Right-multiplying L 1 by a different integer matrix, T 1 , whose determinant has a magnitude greater than one would specify a supercell of C 1 . For example, a determinant-three matrix would specify a supercell having a volume three times that of the primitive cell, L 1 , and contain three times as many basis atoms. By considering supercells of C 1 with lattice vectors S 1 = L 1 T 1 , it becomes possible to identify mappings in cases when the number of basis atoms in the primitive cell of C 2 is an integer multiple of the number of basis atoms in the primitive cell of C 1 . The integer matrix T 1 with determinant greater than one is referred to as the supercell transformation matrix. The lattice-mapping relationship of Eq. (4) can be generalized by allowing for mappings onto supercells of C 1 according to The integer volume of the supercell, defined as the ratio of the volume of the supercell S 1 to the volume of L 1 is equal to detðT 1 Þ. The problem of finding lattice deformations that link two periodic crystals can thus be reduced by Eq. (5) to counting over all possible integer and unimodular transformation matrices and recording the deformation and rotation tensors that link the two structures. This process, while computationally challenging, will provide an exhaustive list of deformations connecting the two periodic unit cells. However, as detailed in the next section, the symmetry of the underlying structures can be exploited to greatly simplify this process.
Accounting for symmetry during lattice mappings The point group symmetry operations that leave the parent or child structures unchanged result in mappings between their lattices that are equivalent. In this section we describe how 68°2 2°R eorientation Fig. 1 Multiplicity of lattice mappings. Depiction of two possible orientation relationships between C 1 and C 2 . In one case, C 2 is rotated by −68 ∘ without an accompanying strain, in the other case, C 2 is rotated by 22 ∘ and also strained. Bottom illustration depicts the definition of the strain cost both before and after removal of rigid rotation. The strain cost is minimized at a unique value when the rigid rotation is removed. J.C. Thomas et al. symmetrically equivalent maps can be avoided while enumerating all possible lattice deformations that relate two structures.
First, we only count over symmetrically distinct supercell lattices of the parent crystal at each integer supercell volume. For a particular supercell, S 1 , there are many possible representations of its lattice vectors, each of which is generated by a particular unimodular matrix N, as is the case for the primitive cell lattice vectors. To enumerate distinct supercells, we employ a widely used approach based on Hermite Normal Forms of integer matrices 40 .
Once the symmetrically unique supercells of the parent lattice have been identified, we can proceed with Eq. (5) by counting over all unimodular matrices and recording the deformation tensor that links the parent and child lattices: We can reduce the number of unimodular matrices N that need to be considered by only using those that generate symmetrically distinct transformations. To determine the symmetrically unique unimodular matrices, we use a convenient identity relationship between different representations of the same symmetry operation. Let us denote the crystal point group symmetry operations that leave the supercell and child lattices unchanged as G ð1Þ and G ð2Þ , respectively. These symmetry operations may be represented in their cartesian form, which we denote as G, or in their fractional form, which we denote as J . The cartesian and fractional representations are related according to: Cartesian symmetry operations are applied to the lattice by left multiplication while the fractional symmetry operations are multiplied to the right of lattice vectors. The cartesian symmetry representation is unitary such that its inverse is equal to its transpose (i.e., G À1 ¼ G T ). The fractional symmetry representation and its inverse are both unimodular integer matrices. Simultaneously applying the cartesian symmetry operation and the corresponding inverse fractional operation leaves the underlying lattice unchanged. For example, Invariance relationships of this type for S 1 and L 2 can be inserted into Eq. (6) to yield G ð1Þ S 1 ðJ ð1Þ Þ À1 N ¼ F ðNÞ G ð2Þ L 2 ðJ ð2Þ Þ À1 : Upon rearranging, this simplifies to where M ¼ ðJ ð1Þ Þ À1 NJ ð2Þ and F ðMÞ ¼ ðG ð1Þ Þ À1 F ðNÞ G ð2Þ are symmetrically equivalent unimodular and deformation matrices. Counting over such deformation tensors that are related under the symmetries of the parent or child structures may be avoided by only using one of the symmetrically related unimodular matrices. Thus for each distinct supercell of C 1 we only count over symmetrically unique unimodular matrices when recording the set of deformation tensors that relate C 1 and C 2 .

Atomic maps
If the primitive cell of C 1 or C 2 contains more than one basis atom, the atoms in C 2 may need to be repositioned inhomogeneously to coincide with sites of C 1 . This deformation of the atomic positions can be represented by a displacement vector for each atom in the unit cell of C 2 . We denote the site positions in the unit cell of C 1 as f r ! ð1Þ 1 ; ; r ! ð1Þ nA g and the atom positions in the unit cell of C 2 as f r ! ð2Þ 1 ; ; r ! ð2Þ nA g. All positions are in Cartesian coordinates. By convention, we measure displacements by first deforming C 2 via a particular lattice mapping, making its lattice sites coincident with lattice sites of C 1 . Then, we rigidly translate the basis of C 2 by a vector t ! to bring the basis atoms of C 2 into closer registration with the sites of C 1 . An atomic displacement is then the vector that translates a particular basis atom j of C 2 to the position of atom i of C 1 . Given a particular rigid translation t ! , the displacement vector that translates atom j of where the operation bÁc S1 selects the shortest possible displacement vector with respect to periodic boundary conditions by considering all possible lattice translations of the lattice S 1 . Equation (12) demonstrates how the position of an atom in C 2 must first be reoriented and deformed, according to the lattice mapping, before relating it to an atom in C 1 via translation. The i;j indicates that it is the displacement measured in the strain state of C 1 .
There are n A ! ways to form (i, j) pairs between sites of C 1 and C 2 . We express a particular one-to-one assignment of atoms in C 2 to sites in C 1 as the permutationp, such that ði;p½iÞ is an assignment of atomp½i in C 2 to site i in C 1 .

Geometric mapping cost
As shown in the previous sections, two arbitrary crystal structures can be related to each other geometrically in a multitude of ways. The set of strains relating the lattices of the two structures are infinite, and for each possible strain map there are several ways to displace the atoms to bring the two structures into perfect coincidence. In this section and the next section we describe how geometry and symmetry may be used to define useful "cost metrics" that can capture the complexity of every enumerated transformation pathway.
Equations (3) and (5) describe the relationship between two lattices through a symmetric left stretch tensor (V) and the Biot strain tensor (B). These tensors describe changes in both the volume and shape of crystal C 2 . Volume change alone is a poor measure of crystal structure similarity since volume depends strongly on chemistry, temperature, and other considerations. The volumetric dependence of the left stretch tensor is given simply by det V ¼ V L2 =V L1 . Consequently, we can renormalize the left stretch tensor to make it independent of volume changes, yielding the volume-normalized stretch tensor and the corresponding volume-normalized strain tensor B ¼Ṽ À I: We use the volume-normalized strain to define cost functions of lattice deformation in order to isolate shape-change effects from volume-change effects.
We now introduce a lattice-deformation cost function that measures the magnitude of lattice strain. In general, the Biot strain tensor in Eq. (14) depends on the orientation of the coordinate system used to define the lattice vectors of C 1 and C 2 . However, the lattice-deformation cost function should be invariant to changes in coordinate system. A common metric of tensor magnitude that is invariant to coordinate system is the Frobenius norm: To demonstrate that the Frobenius norm of B is a physically meaningful metric of deformation, we consider a uniform sphere of radius ρ as it is deformed by the volume-normalized stretch tensorṼ. The displacement of a point on the sphere is and the mean-squared displacement over the surface of the sphere is the integral taken over the domain of solid angle in three dimensions. Expanding the integrand term by term and integrating yields avg δ ρ The spherical radius ρ can be set to a physically meaningful parameter, such as the Wigner-Seitz radius of the child crystal, to obtain a per-atom metric of lattice deformation. However, Eq. (18) can also be normalized by ρ 2 to obtain a scale-invariant metric of deformation, which is how we shall specify the latticedeformation cost function: Measuring the degree of atomic distortion is simpler than measuring the lattice strain. We measure the quality of a particular atomic assignmentp with the mean-squared atomic displacement: As in the case of the lattice-deformation cost function, we normalize the atomic-deformation cost by the squared Wigner-Seitz radius to obtain a scale-invariant cost. The cost function Eq. (21) is a function of both the rigid translation t ! and atomic assignmentp, and we will show that it can be minimized with respect to these two parameters. It is also implicitly, via Eq. (12), a function of the integer matrices N and T 1 , which determine the lattice mapping. The lattice-deformation cost function c L M m ð Þ and the atomicdeformation cost function c A M m ð Þ measure qualitatively distinct aspects of crystal deformation. We define mappings to be "optimal" if they simultaneously minimize both the lattice and atomic-deformation costs. Specifically, we define the following scalar total cost function: which is a weighted average of the atomic and lattice-deformation cost functions specified by a lattice weight parameter w L . The presence of two cost functions (lattice and atomic costs) leads to the well-known problem in multi-objective optimization where there is usually no unambiguous 'optimal' mapping (The exception to this is when C 1 and C 2 are structurally identical, in which case some mappings have both zero lattice deformation and zero atomic deformation). By introducing the total cost function that is related to both these costs, the minima is welldefined. Because the atom-and lattice-deformation cost functions defined in Eqs. (20), (21) are both scaled to remove volumetric dependence and are normalized per atom, we anticipate that a 50/50 weighted average, corresponding to w L = 0.5 should be appropriate for many applications. Values of w L > 0.5 increase the relative penalty of lattice deformation (thus allowing atomic displacements to be relatively larger) while values of w L < 0.5 decrease the relative penalty of lattice deformation (thus biasing atomic displacements to remain small).
Symmetry-adapted mapping cost Crystal structures that differ by more than just their volumes can still be very similar. For instance, the hexagonal close-packed (hcp) structure may be deformed to have any c a ratio without qualitatively changing the crystal. Hcp structures that differ with respect to their c a ratio will have a nonzero geometric cost despite being crystallographically identical. To enable the easy comparison of such structures we define a "symmetry-adapted" mapping cost that is invariant under deformations that preserve the symmetry of the parent crystal.
The symmetry-adapted cost is obtained by removing any symmetry-preserving deformations or displacements from the mapping field. The symmetry-preserving Biot strain tensor, B sym , can be obtained by averaging the effect of each point group operation on the calculated stretch tensor obtained from Eq. (5): where G 2 is the point group of the child structure (C 2 ) and N G2 is the number of operations in the point group. Formally this corresponds to applying the Reynolds operator of the point group to the Biot strain tensor. Obtaining the average over all possible symmetrically equivalent Biot strain tensors leaves us with only that part of the deformation that is invariant under symmetry. Examples of symmetry-preserving deformations include simple hydrostatic deformations of the unit cells and changes to the c a ratio in hcp. The symmetry-adapted mapping score is computed using the symmetry-breaking strain tensor defined as The symmetry-adapted lattice cost is then calculated as: The symmetry preserving score relative to the parent may be calculated by utilizing the stretch tensor relative to the parent and the point group operations of the parent.
The symmetry-adapted atomic cost can be obtained in a manner similar to that described for the lattice-deformation cost. First, we obtain the collective atomic displacements that are invariant under the factor group of the parent structure by again applying the Reynolds operator. Subsequently, the displacements Eq. (12) that map the atoms of the child onto the atoms of the parent can be projected on to the symmetry-invariant collective displacements of the parent and removed to estimate the symmetry-breaking displacements. To perform these steps it is convenient to collect atomic displacements relative to the parent crystal atomic positions as a vector containing 3n entries These atomic shuffles can be transformed by symmetry operations in the factor group of the parent crystal according to where G 1 is now a 3n × 3n matrix representation of a symmetry operation in the factor group of the parent structure that acts on the unrolled displacement vector D ! . The 3n × 3n matrix representations can be constructed from the 3 × 3 Cartesian representation and the permutation representation of the factor group operations.
Symmetry-invariant modes can be constructed by applying the Reynolds operator to the set of basis vectors (denoted I, the 3n × 3n identity matrix) that span the full 3n × 3n space 15,41 : where 〈D〉 is the symmetry-averaged basis vectors. Each column of the 〈D〉 matrix corresponds to a collective displacement field that is invariant to the symmetry of the parent crystal. Not all columns are independent from each other and a set of spanning vectors may be identified by performing singular value decomposition (SVD) of the matrix The dimensionality (r) of the symmetry-invariant space is equal to the number of nonzero singular values (or diagonal elements) in Σ, while a set of orthogonal vectors that span this space can be derived from the r column vectors of W corresponding to the nonzero singular values. This space of symmetry-preserving distortions can be expanded into larger supercells, and their symmetry-breaking deformation fields can then be obtained as: where D ! sym is a matrix of symmetry-invariant collective atomic modes calculated from Eq. (29) and expanded so that it accounts for the total number of atoms in the supercell of C 1 that is being considered. D ! i sym refers to the ith collective displacement mode in this matrix.
The symmetry-breaking cost is then obtained as: A summary of the mapping algorithm The mapping algorithm described above is schematically illustrated in Fig. 2. The figure shows the steps involved in mapping a two dimensional honeycomb arrangement of atoms (labeled C 2 ) on to a parent square lattice (labeled C 1 ). First the algorithm counts over all symmetrically distinct supercells of the parent square lattice that contains two atoms. The first lattice map, labeled F 1 , involves a large shear and tetragonal deformation of the child lattice. The second lattice map (F 2 ) relates the child structure to an equivalent but geometrically different supercell. The lattice undergoes a tetragonal deformation and a smaller shear distortion. Finally, F 3 relates the child structure to a supercell of the square lattice that is different from those depicted in the first two deformations. Having enumerated some of the strain deformations, the algorithm then deforms the child structure and rigidly translates the atoms such that at least one of the atoms is in perfect coincidence with the lattice points of the parent structure. Figure 2 illustrates the effect of applying the homogeneous strain distortion and rigid translation to the child structure across all three lattice maps. The other atom then needs to be displaced along the highlighted vector for each of the three maps to bring it into perfect coincidence with the parent structure. Having recovered the parent crystal structure from the honeycomb child structure, the algorithm estimates the scalar costs associated with each deformation and displacement through Eqs. (20), (21), (31), (25).

Finding optimal maps
In the previous sections we outlined how a set of mappings M linking any two arbitrary crystals can be generated. First a set of strains linking the crystals is generated by counting over all symmetrically distinct lattices. For each strain map, a set of n! atomic maps can be generated by assigning and displacing every atom in the child so that it comes into perfect coincidence with its counterpart in the parent. We also defined two ways of quantifying the extent of distortion in each of these maps through cost functions. The first, which we refer to as the geometric cost calculates a scalar quantity for the strain and atomic deformations in each map M m . The second, called symmetry-adapted cost, uses the symmetry of the underlying parent crystal structure to calculate these costs over only the symmetry-breaking parts of the strain and displacement fields. In both cases we define a total cost function that combines the individual contributions from the strain and atomic distortions to calculate a single scalar mapping cost for every mapping in M. In this section we discuss optimization algorithms that can be leveraged to quickly and efficiently obtain the map with the lowest total mapping cost and an algorithm to systematically traverse the multitude of mappings in increasing order of their mapping costs The lattice-deformation cost function for a particular supercell transformation T 1 can be optimized simply by counting over unimodular matrices N in Eq. (5). Some of the simplest examples of N matrices describe operations such as a swap or permutation of lattice vectors, or the replacement of a lattice vector by its negative direction. Although N must always have a determinant of ±1, there is no limit on the magnitude of its coefficients. However, if N includes coefficients of very large magnitude, the lattice vectors S 0 1 ¼ S 1 N that define the lattice become increasingly elongated, and favorable mappings become unlikely. To minimize the range of values of N that must be considered, we first obtain the reduced cells of S 1 and L 2 using, for example, the Lenstra-Lenstra-Lovász algorithm 42 . The reduced cells are maximally compact, such that their lattice vectors have the minimal possible lengths and the angles between lattice vectors are as close to 90 ∘ as possible. The reduced cell is not unique, although conventions do exist, such as the Niggli reduced cell, that provide a unique definition, to within a rigid rotation. If both L (2) and S (1) are in reduced form, we need only consider lattice-equivalent transforms that have small entries, e.g., {−1, 0, 1}, or {−2, …, 2}. Once a range for the elements of N has been decided upon, it is straightforward to count first over distinct supercells {T 1 } and then over unimodular matrices {N} within the specified range and record the lattice-deformation cost for each (T 1 , N). We find it useful for many cases to simply retain all such mapping relationships, sorted by their weighted lattice-deformation cost or to find the best or K-best lattice-mapping relationships. The algorithm remains unchanged when using the symmetry-adapted lattice costs as defined by eq. (25).
A total mapping of C 2 onto C 1 comprises the ordered tuple M m ¼ ðT 1 ; N; t ! ;pÞ, where T 1 and N are integer matrices specifying the lattice mapping, t ! is a rigid translation of crystal C 2 relative to C 1 , andp is a permutation specifying the atomic mapping. The optimal atomic mapping for a given lattice mapping (i.e., N and T 1 ) is one that minimizes Eq. (21) with respect to both t ! andp. The rigid translation that minimizes Eq.
(21) for a given permutation operation,p, is simply the translation that yields zero average displacement of the atoms in C 2 (if nonzero, the average displacement amounts to an overall 'drift' of C 2 relative to C 1 ). The permutation operationp, in contrast, is a discrete variable-the minimization of Eq. (21) with respect top for a particular choice of t ! is a discrete combinatorial optimization problem. Fortunately, since Eq. (21) is a simple linear combination of contributions from ði;p½iÞ pairs, this optimization is classified as a linear assignment problem, for which well-known solution methods exist 43 . Because the atomic-deformation cost function in Eq. (21) depends on both the discrete assignmentp and the continuous translation vector t ! , we propose an iterative procedure to optimize the atomic-deformation cost. In the first step, we identify potential rigid translations of C 2 that are likely to be near the optimal translation. These are translations that bring an atom in C 2 into coincidence with a chemically compatible site of C 1 . It is preferable at this step to pick an atom from among the minority species. For example, if C 1 and C 2 each have two Cu atoms and four Au atoms per unit cell, we would pick as candidate translation vectors t ! 0 κ those that map one Cu atom of C 2 onto each of the two Cu atoms in C 1 . In general, we denote the number of minority species in C 2 as m s Next, we find the optimal atomic assignment for each attempted translation t ! 0 κ , κ = 1, …, n ms where n ms is the number of minority species. To do so, we construct the atomic assignment cost matrix, Cð t ! 0 κ Þ, for each rigid translation. The (i, j) element of the cost matrix is the pairwise cost contribution of assigning site j in C 2 to site i in C 1 for the rigid translation t ! 0 κ . The elements of the cost matrix are where an element of the cost matrix is infinite (or at least very large) when the species at sites i and j, labeled s i and s j , respectively, are chemically incompatible. This ensures that chemically incompatible assignments are excluded from the optimal assignment. Once the cost matrix is constructed according to eq. (32), the optimal assignmentp ? κ is the solution tô The minimization of Eq. (33) can be solved by the Hungarian algorithm 44 (used in this work) or via linear programming methods 43 .
Once the optimal assignment is found for an attempted translation, t ! 0 κ , we update the translation by subtracting the average displacement (i.e., drift) due to that assignmentp ? κ according to i;p ? Fig. 2 Schematic of the mapping algorithm. A summary of the steps involved in generating a map between a parent structure with a square lattice and a child structure with the atoms arranged in a honeycomb pattern. Three lattice and atomic maps that relate these structures are illustrated, with the final step being the evaluation of the lattice and atomic cost functions.
Next, we calculate the atomic-deformation cost, eq. (21), for the assignment/translation pair ðp ? κ ; t ! κ Þ and record the result. This must be performed for each attempted translation κ = 1, …, n ms . The optimal atomic assignmentp ? κ for the imposed lattice map is the one with the lowest atomic-deformation cost.
In certain extreme cases, it may be necessary to perform additional optimization steps. When calculating the revised atomic displacements due to the augmented translation t ! κ , it should be verified that the average total displacement (i.e., drift) is indeed equal to zero, i.e., 1 Violations of Eq. (35) may occur due to the periodic boundary considerations imposed in Eq. (12). If an outlier atom has a very large atomic displacement under translation t ! 0 κ , it may indicate that atomp½i in C 2 is nearly midway between two periodic images of site i in C 1 . In such a case, a small change in rigid translation may affect a large change in the displacement of the atom (e.g., the revised translation swaps the firstand second-nearest periodic image of atom i). If a violation of Eq. (35) is encountered, it indicates that the optimal atomic assignment may be unstable relative to small rigid translations in the vicinity of t ! 0 κ . As such, the augmented translation t ! κ may be recalculated iteratively until Eq. (35) is satisfied, with the augmented translation at each iteration added to the list of plausible translations to be considered for solution of the linear assignment problem.
Minimizing the symmetry-adapted atomic cost requires additional steps. Unlike the geometric-cost function of eq. (21), the symmetry-adapted atomic cost can only be calculated once a particular assignment of atoms has been chosen. As a result the Hungarian algorithm cannot be directly applied to obtain the minimum cost assignment. The algorithm can, however, be adapted to systematically search over suboptimal maps to arrive at the assignment of atoms that minimizes the symmetry-adapted cost. To do this, we adopt Murty's algorithm 45 to generate a queue of suboptimal maps that are arranged in increasing magnitude of their geometric cost (Eq. (21)). We then search this queue up to some maximum depth that is user specified and recalculate the symmetry-adapted cost for each of these elements. The element with the minimal symmetry-adapted cost is then recorded as the optimal map. This algorithm is guaranteed to provide the optimal map if all n! assignments are enumerated. However in practice we have found that the assignment that minimizes the symmetryadapted cost is often found by searching over the first 5-10 suboptimal assignments. A detailed algorithm outlining a methodology to systematically enumerate all possible suboptimal maps is outlined in the supplementary information.

Case studies
We illustrate the mapping algorithm by first comparing the geometric-cost functions against the symmetry-adapted cost function for two parent structure prototypes. Next we demonstrate how well-known structural transformation pathways can be identified by systematically analyzing all possible mappings and their associated cost function values.
The hexagonal close-packed (hcp) crystal structure is usually characterized by two lattice parameters that are conventionally denoted a and c. The distance between the atoms in a particular triangular layer of hcp is equal to a while the spacing between the layers is half the c lattice parameter. The underlying crystal structure remains qualitatively unchanged across all values of a and c. It is desirable to use a mapping function that is insensitive to such changes in the crystal structure. Figure 3 compares the difference between the geometric and symmetry-adapted lattice costs for a child structure with a wide-range of c a ratios that is mapped on to a parent crystal with a c a % 1:6. The geometric cost rises rapidly at values that deviate from the parents' lattice parameters. On the other hand the symmetry-adapted cost remains exactly zero across the full range of values. Figure 4a shows the crystal structure of an orthorhombic phase, denoted O 0 , that is frequently formed in titanium-rich alloys 46 . Similar to hcp, the O 0 crystal structure consists of triangular layers; however, they do not adopt the ideal stacking as in hcp as is illustrated in Fig. 4a. Furthermore, the unit cell of O 0 is slightly strained relative to that of hcp. As a result, the O 0 structure has a lower symmetry than hcp. We ignore the lattice strain for the purpose of demonstrating our algorithm. Figure 4b compares the geometric cost against the symmetry-adapted cost when mapping O 0 onto a reference O 0 structure as a function of the shuffle amplitude of alternating close-packed layers of O 0 (i.e., along the arrows in Fig. 4a). The geometric cost rapidly increases with the shuffle amplitude until it reaches a maximum before it starts decreasing. The decrease in geometric cost is due to the periodicity of the parent crystal structure. In contrast to the periodically varying geometric cost, the symmetry-adapted cost remains identically zero. The shuffle of alternating close-packed planes as shown in Fig. 4 do not result in further symmetry breaking of O 0 and thus do not incur a cost penalty. We note that as the layers are translated relative to each other in O 0 , a particular shuffle amplitude will result in the formation of perfect hcp. It is not possible to easily identify shuffles that increase the symmetry  of the underlying crystal with just the symmetry operations of the low-symmetry structure.
There is a large number of structural relationships between any pair of crystals. As an illustration, we explore mapping relationships between well-known crystal structures such as bcc, fcc, hcp, and ω. The ω phase can be described as alternating triangular and honeycomb layers, with the atoms of the triangular layers residing above and below the centers of the hexagons of the honeycomb layers. Figure 5 plots the lattice and atomic displacement costs for different mappings of fcc, hcp, and ω onto bcc. Three pathways connecting bcc with hcp, fcc, and ω are highlighted in Fig. 5b. Also included are different mappings of bcc onto bcc, which are denoted with a circle.
Several of the mappings with small mapping scores in Fig. 5 correspond to the end points of well-known pathways that convert one crystal structure into another. Among them are the Bain path that connects bcc to fcc 19,26,27 , the Burgers path that connects bcc to hcp 20,26,27 and a path that converts bcc into ω 21,22,27 . Each pathway is linked to a particular mapping relationship. The Bain path, illustrated in Fig. 6, converts bcc into fcc through a tetragonal tensile strain of the bcc unit cell vectors along one of its cubic axes. The mapping corresponding to this path only has a lattice cost, but no atomic displacement cost. The resulting fcc crystal can be converted to another bcc crystal through a Bain path along a different fcc cubic axis. This leads to the mapping relationship between two bcc crystals denoted by the circle with a lattice cost of ≈0.1 in Fig. 5.
The Burgers path, which connects bcc to hcp, connects end states with a mapping that has both a lattice cost and an atomic displacement cost. As illustrated in Fig. 6, the Burgers path couples a compressive tetragonal distortion of bcc along one of its cubic axes, followed by a shuffle of (110) planes as referenced to the original bcc lattice. While the lattice cost of the Burgers path is smaller than that of the Bain path, it does have a sizable atomic displacement cost.
A third well-known path, due to Silcock 21 and de Fontaine 22,47 , connects bcc to ω. Early descriptions of the bcc to ω path invoked a transverse shuffle of ð112Þ planes along the (111) direction 21,48 . Later, de Fontaine pointed out a more intuitive description in terms of longitudinal shuffles of (111) planes 22 . This is illustrated in Fig. 7, where every two out of three (111) planes of bcc collapse to form a honeycomb layer, thereby forming the ω crystal structure consisting of alternating honeycomb layers and triangular layers. The supercell of bcc must have an integer volume of three and has a hexagonal unit cell, with its A ! and B ! vectors spanning one family of (111) planes and its C ! axis perpendicular to the same family of (111) planes of bcc (i.e., A ! is parallel to ½101, B ! is parallel to ½011 and C ! is parallel to [111]). As is evident in Fig. 5, the mapping between bcc and ω corresponding to the Silcock-de Fontaine path has a very small lattice cost but a large atomic displacement cost due to the collapse of adjacent (111) planes. Figure 8 explores the structural relationship between hcp and common crystal structure prototypes such as bcc, fcc, ω, and a family of topologically close-packed phases. The Burgers transformation relating bcc and hcp requires a large atomic shuffle and a small lattice strain. The transformation pathway to the fcc structure with an atomic cost of ≈0.07 can be viewed as successively applying the Burgers and Bain transformation to hcp.
Structural relationships between hcp and the ω crystal structure have been extensively explored by Trinkle et al. 23,24 . Our mapping algorithm shows that two of these transformation pathways, labeled TAO-1 and Silcock, have the smallest cost. The Silcock pathway primarily involves atomic shuffles within the basal plane such that the (0001) planes of hcp are transformed into the ð1120Þ planes of ω. The TAO-1 transformation pathway, which is predicted to have the lowest barrier in pure titanium 23 , transforms the basal plane of hcp to ð0111Þ planes of ω. Both transformation pathways require significant shuffles but are accompanied by small strain distortions.  The Bain and Burgers paths. The bcc crystal can be converted to fcc through a tetragonal strain known as the Bain path 19 or it can be converted to hcp along the Burgers path, which combines a compressive tetragonal strain with an atomic shuffle 20 . Figure 8 also shows mapping costs for transformation pathways that link hcp to a family of topologically close-packed structures. These pathways, as described by Natarajan and Van der Ven 25 , involve a transformation of a subset of the hcp close-packed triangular layers into kagome nets and a dissociation of other close-packed triangular layers into three new, less dense triangular layers. The formation of a kagome net from a triangular lattice can be realized with in-plane shuffles while the dissociation of triangular lattices into three new triangular lattices occurs with out-of-plane shuffles. A large number of topologically closepacked phases that include C14, C15, and C36 Laves phases, CaCu 5 , Be 3 Nb, and Co 7 Gd 2 can be generated from hcp through a combination of in-plane and out-of-plane shuffles 25 .

DISCUSSION
We have described an approach to enumerate mappings between any pair of crystal structures. The approach first identifies symmetrically distinct lattice maps and accompanying strains that deform the lattice of the child structure onto that of the parent. For each lattice map, the basis atoms within the child structure are displaced until they coincide perfectly with the parent structure. As there are an infinite number of ways to mathematically link two crystal structures, we introduced two cost metrics that serve as a measure of the distance along a map. The geometric-cost metric measures the size of the strain and displacement fields required to bring the two structures into coincidence. The symmetry-adapted cost metric measures only that part of the distortion field that breaks the symmetry of the parent crystal structure. Finally we outlined optimization and search algorithms that can be employed to obtain the structure mapping with lowest cost metric and a family of suboptimal structural transformations. The mapping algorithm can be used to connect ordered crystals to unordered host structures or other ordered phases at the same composition. The algorithm has been implemented within the CASM software package 15,38,39 . The implementation can be used to generate maps between any pair of crystal structures and to assign a cost to each map.
The pathway generation scheme outlined here has several similarities to that introduced by Trinkle et al. 24 . For example, the enumeration over possible lattice mappings using unimodular matrices has much in common with the approach of Trinkle et al. 24 . There are, however, several crucial differences. First, the approach of this work explicitly leverages the symmetry of the underlying parent and child crystal structures to prune the set of possible maps and to discern a symmetry-adapted mapping score. As outlined in section Accounting for symmetry during lattice mappings, strain deformations that are identical under the point group operations of the two crystal structures are screened using eq. (10). Trinkle's algorithm relies on comparisons of the eigenstrains of the deformation tensors to prune elements of the map set M that represent the same deformation field but are numerically different. Second, our approach employs the Hungarian algorithm (and Murty's algorithm) to systematically enumerate all feasible assignments of atoms between the child and parent structures when generating displacement fields. In contrast, a brute force approach that only counts over "small" displacements of each atom is used in 24 . Finally, we propose a way of characterizing the complexity of the transformation pathway using a symmetry-adapted cost metric. This cost metric accounts for only symmetry-breaking deformation fields and enables the comparison of qualitatively similar crystal structures that might have very different lattice parameters and atomic coordinates. This is an especially Fig. 7 The de Fontaine path between bcc and omega. The ω crystal structure can be formed from bcc by collapsing two out of every three bcc (111) planes 22 . powerful feature that, to our knowledge, has not yet been exploited in any previously developed crystal mapping algorithms.
Our systematic enumeration of pathways between common crystal structures as summarized in Figs. 5 and 8 has revealed that all known structural transformations can be associated with pathways that have the smallest or close to the smallest mapping costs. The mapping costs enable the systematic exploration of physically meaningful structural transformation pathways. An analysis of transformation pathways relating several commonly occurring crystal structure prototypes 49 using the algorithm outlined in this study has revealed seven new structural transformation pathways. Using this algorithm to link all known crystal structure prototypes will undoubtedly generate a rich database of phase transformation pathways that will provide insights about a wide variety of structural phase transformations.
The algorithms described here and its implementation within CASM should prove invaluable in organizing high-throughput firstprinciples databases by enabling a precise cataloguing of crystal structures and the identification of relaxations to new structures due to dynamical instabilities. Similarly, the algorithms will benefit the field of alloy theory 39,50,51 , where cluster expansion Hamiltonians are trained to first-principles energies of a large number of chemical decorations over the sites of a particular parent crystal structure [52][53][54][55][56][57][58][59][60][61][62][63][64] . In many instances, a subset of orderings on a parent crystal exhibit dynamical instabilities and relax to a decoration on a different parent crystal structure 59,65 . These relaxations contaminate the training dataset of a cluster expansion and must be detected and eliminated from the training pool. A robust mapping and distance metric is necessary to establish whether a particular ordering can still be considered a decoration on the original parent crystal structure. The mapping tool can also be used to identify symmetric equivalence of deformed crystals, as for example, the crystallographic models of sheared bicrystals used to calculate generalized stacking fault energies 66 .

Mapping crystal structures
Transformation pathways connecting various crystal structures were generated with the Clusters Approach to Statistical Mechanics (CASM) code. Figures  4 and 3 were generated by employing the symmetry-adapted and geometriccost metrics to quantify the complexity of the transformation pathways. The mapping costs shown in Figs. 5 and 8 were generated by searching for symmetrically distinct transformation pathways connecting hcp and bcc to several commonly occurring crystal structures. For each set of parent and child structures, we considered all symmetrically distinct supercells of the child containing at most 12 atoms. Lattice mappings relating the child superstructure to a corresponding superstructure of the parent were generated by counting over all unimodular matrices containing entries with a maximum magnitude of 2. We then recorded the five best atomic maps for each lattice map with the optimality measure being given by the geometric cost. The symmetry-adapted cost metrics were calculated for each of these sets of lattice and atomic maps using Eqs. (24), (31). The symmetry-adapted lattice costs were calculated relative to the parent structures.

DATA AVAILABILITY
The crystal structures used in this study are provided in the Supplementary Information.

CODE AVAILABILITY
The algorithms described in this study have been implemented in the Clusters Approach to Statistical Mechanics (CASM) open source code that is available for download from https://github.com/prisms-center/CASMcode. The mapping algorithms are seamlessly integrated into the workflow of CASM and are leveraged when populating internal configurational databases from electronic structure calculations or when importing arbitrary structures through the "casm import" command.