Matching Crystal Structures Atom-to-Atom
Félix Therrien Peter Graf Vladan Stevanović


The Journal of Chemical Physics








Colorado School of Mines
National Renewable Energy Laboratory

National Renewable Energy Laboratory

vstevano@mines.edu
Colorado School of Mines
National Renewable Energy Laboratory



Finding an optimal match between two different crystal structures underpins many important materials science problems, including describing solid-solid phase transitions, developing models for interface and grain boundary structures. In this work, we formulate the matching of crystals as an optimization problem where the goal is to find the alignment and the atom-to-atom map that minimize a given cost function such as the Euclidean distance between the atoms. We construct an algorithm that directly solves this problem for large finite portions of the crystals and retrieves the periodicity of the match subsequently. We demonstrate its capacity to describe transformation pathways between known polymorphs and to reproduce experimentally realized structures of semi-coherent interfaces. Additionally, from our findings we define a rigorous metric for measuring distances between crystal structures that can be used to properly quantify their geometric (Euclidean) closeness.


  Keywords
Structure matching, Polymorphism, Phase transitions, Heterojunctions



Matching Crystal Structures Atom-to-Atom
    Vladan Stevanović
    December 30, 2023
========================================






§ INTRODUCTION



Establishing an optimal match between two different crystal structures with respect to some cost function is a problem that cuts across the entire field of materials science. Perhaps the most evident example is the process of finding a suitable substrate to epitaxially grow a material <cit.>. Similarly, when studying interfaces between different phases (heterojunctions) one might be interested in the alignment and the bonding pattern between the two phases. Another important example lies in finding minimal energy pathways between different polymorphs. The initial and final structures are known, but the transformation from one to the other is not. To even begin to describe it, one needs to find the best way to map every atom of the initial structure to its counterpart in the final structure and to optimally align the structures. Once the mapping and alignment are established other methods such as the Solid State Nudge Elastic Band <cit.> can be used to determine the energetics of the transition.

In regard to interfaces, many have worked on methods to find and characterize the coincidence of lattices and orientation relationships between phases and grains. Several different approaches were developed e.g., the O-lattice theory part of the CSL/DSC Lattice Model [CSL: Coincidental Site Lattice, DSC: Displacement of one crystal Lattice with respect to the second causes a pattern Shift which is Complete] 
<cit.>, the Edge to Edge Model <cit.>, the Coincidence of Reciprocal Lattice Points (CRLP) model <cit.>, methods based on the Zur Algorithm <cit.> and the work of Jelver et al. <cit.>. While the O-lattice theory suffers from a lack of predictive capabilities, the other approaches do have the ability to predict orientation relationships, but they do not match the full structures. The Edge to Edge model only considers high density (nearly close packed) planes and directions whereas the CRLP and Zur Algorithm only match the underlying lattices  of the structures, and not the atoms inside them. Jelver et al. presented a crystal matching method that maps atoms inside a combination of the unit cells of the two structures which, as explained further in this section, has some important limitations. 


Matching is also closely related to measuring distances between crystal structures. Indeed, any definition of a distance metric requires establishing some correspondence between their atoms. For finite systems such as molecules, 
Sadeghi and Goedecker <cit.> defined a distance as the minimal l_2-norm of the vector joining the molecules in the configuration space of atoms with respect to both their relative positions (alignment) and the permutation of atomic indices [The l_2-norm in configuration space equivalent to the Forbenius norm of the position matrix]. However, the configuration space of periodic systems is, strictly speaking, ill-defined because of the infinite number of dimensions. Moreover, the permutation degeneracy in labeling atoms also poses problems. Both of these make the matching of crystal structures challenging and the definition of the distance metric between periodic structures elusive.

In structure predictions, for the purpose of identifying similar (close) structures Oganov elegantly circumvented this problem by introducing the so-called fingerprint function constructed to reflect the short-range order (coordination in various shells, etc.) and defining the distance metric between crystal structures with respect to it <cit.>. Various other fingerprint functions and representations have been proposed since <cit.>. Although they can efficiently measure the similarity between structures, those methods do not establish a one-to-one correspondence between each atom of the structures, nor do they provide the optimal alignment between them. Therefore, the question whether a Euclidean distance metric in the configuration space of atoms (l_2-norm) and its corresponding matching can be defined for periodic systems remains open.



An intuitive way to go about this problem is to take advantage of the periodicity by matching and measuring the distance between atoms inside the unit cells of the two structures. This inevitably leads to the obstacles illustrated in Fig. <ref>. Panel (a) shows two 2D crystals with different unit cells each having one atom per cell. If, for example, one wishes to minimize the distance between the atoms, a naive way to match these two structures would be to align their unit cells so that the atoms overlap. This would produce a distance of zero. Other atoms would simply be mapped based on the correspondence between the unit cells. However, this mapping would lead to the distances between the corresponding atoms diverging as one moves away from the two perfectly aligned unit cells at the center of Fig. <ref>(b). There is, however, a solution to this particular problem shown in Fig. <ref>(c) that does not suffer from this divergence and that produces equal and finite distances between all corresponding atoms. This solution yields a much shorter total distance if large portions (grains) of the two crystals are considered. While matching larger supercells and taking advantage of their Niggli-Santora-Gruber reduced cell <cit.> might seem like a solution to this particular example, choosing the size of the supercells and comparing distances between different sizes remains an issue.

One way to robustly find solutions such as the one from Fig. <ref>(c) is by disregarding the periodicity in the two structures. If it exists, the periodicity of the mapping itself can be retrieved subsequently. For this reason, as we will explain, our algorithm is constructed to use large sections of the two crystals and to minimize the total distance traveled by all the atoms. Consequently, the choice of a unit cell has no impact on the final result and the periodic unit of the transformation emerges naturally. Any other method to match two crystal structures that relies on matching some choice of the unit cell, including our previous work <cit.> as well as work of others <cit.>, will suffer from the problem illustrated in Fig. <ref>.

With the aforementioned considerations, we formulate the problem of matching crystal structures in the following way: given the positions of all the atoms in the two crystals, {a⃗_i | i=1,…,N } and {b⃗_j | i=1,…,N}, what is the best atom-to-atom mapping p_min (permutation of atomic indices) and the best alignment of the two structures (linear transformation Q_min and translation t⃗_min) that minimize a given distance (cost) function d? This is equivalent to solving the following equation when N →∞: 


    p_min, Q_min,t⃗_min = p, Q, t⃗argmin∑_i^N d(  a⃗_i ,  Qb⃗_p(i) + t⃗  ).


Formulating the problem in this way has the advantage of making the solution method easily adaptable to any given cost function, being the sum of Euclidean distances (l_2-norms) between the atoms or some other function depending on the particular problem or application. Posed as such, matching the structures is equivalent to doing a Point Set Registration (PSR)<cit.>, a well-studied process used in computer vision and pattern recognition. Our algorithm is inspired by PSR methods.

Herein, we describe our structure matching algorithm in detail and showcase its applications to phase transformations and semi-coherent interfaces. We demonstrate that it robustly reproduces known results for several well-studied polymorphic transformations. It also seamlessly reproduces and explains the experimentally observed semi-coherency and orientation relationships for the interfaces between Si and SiC, and Ni and Yttria Stabilized Zirconia (YSZ). Finally, drawing upon the results of our crystal matching algorithm we discuss and propose a proper Euclidean metric between infinitely periodic crystal structures.




§ THE ALGORITHM



We start by introducing nomenclature used in this paper. Crystal structures are represented by three objects: the unit cell matrix (C), the 3× N matrix of atomic positions (P) where N is the number of atoms in the unit cell, and the 1× N list (L) of the symbols of chemical elements occupying those positions. It is understood that the order of elements in the matrix P and the list L is the same. Finally, we combine these three into a single crystal structure object:


    A ≡{ C_A , P_A, L_A}.


It is important to note that there exist an infinite number of representations of the same structure due to the arbitrariness in the choice of the periodic unit (unit cell).






 §.§ Distance Minimization



The first step to solving equation (<ref>) in practice is to make the two structures we would like to match finite. The structures A and B are the primary input to our algorithm along with the size of the finite sections of the structures and the distance function (d) to optimize. We make A and B finite by cutting out spherical sections around the origins of the two structures ensuring that the stoichiometry is preserved in each structure and that they both have the same number of atoms. Making the sections approximately spherical is done by selecting the atoms that are the closest to a central point in a process called tiling as denoted in the flowchart of our algorithm shown in Fig. <ref>.

Next, the structures are brought to the same geometric centers and an initial random rotation Q_0 and a random translation t⃗_0 are applied to one of them (step 2, see Fig. <ref>). The structure that is rotated and translated is labeled “mapping structure” (B), to which all subsequent geometric transformations will be applied and the other one is labeled “mapped structure” (A) and it remains fixed in space. The initial translation t⃗_0 is constrained to within one unit cell of the mapped structure (A).

After this initial alignment, atoms in the two structures are mapped to each other, i.e., the permutation p_1 of atom indices in the mapping structure is chosen such that the distance function is minimized (step 3). The details of the mapping procedure are provided in section <ref>. Next, for this particular atom-to-atom map, the distance between the two structures is minimized with respect to rotations Q_1 and translations t⃗_1 using a gradient descent (step 4). It is important to note that Q does not need to be a rigid rotation. Depending on the application, it can also contain a certain amount of deformation such that 1-ε < det(Q) < 1+ε. At this point, the atom-to-atom map p_1 is not necessarily optimal, since it was established before the translation and rotation were optimized. Step 3 is therefore repeated using the new alignment Q_1 and t⃗_⃗1⃗ to obtain p_2. Then, with this new mapping, the algorithm finds Q_2 and t⃗_2 (step 4), remaps again and so on, iteratively, until the p,Q and t⃗ stop changing. This iterative procedure can be mathematically formulated as: 


    p_j  =  pargmin∑_i^N d(a⃗_i,Q_j-1b⃗_p(i) + t⃗_j-1),  
    t⃗_j, Q_j =  t⃗, Qargmin∑_i^N d(a⃗_i, Qb⃗_p_j(i) + t⃗ ),


where the index j is the iteration number of the Mapping Loop from Fig. <ref>. At the end of the Mapping Loop the algorithm has reached a local minimum. In order to find the global minimum, one needs to explore the dependence of the results on the random initialization. This we do by constructing an outer Random Minimization Loop by repeating the whole procedure (i.e., steps 2,3 and 4) a large number of times until the best local minimum stops changing.

The distance minimization part of the algorithm is a form of iterative closest point (ICP)<cit.>, a classic scheme to solve point set registration (PSR) problems that iteratively associates two sets of points and minimizes the distance between them. In our case, the data association is done by solving the assignment problem using the Khun-Munkres algorithm (see section <ref>). To the best of our knowledge, only a few ICP algorithms have used this method to associate data <cit.>. Moreover, our problem differs widely from most PSR problems for the following reasons: (1) The crystals are by definition featureless since they are infinite periodic sets of points. Macroscopic features cannot be used to reduce the number of points or to find an approximate result like it is done in many PSR algorithms including refs <cit.> (2) The two structures are three-dimensional (in the case of phase transformations), i.e., points occupy the full 3D space, they do not represent the surface of a 3D object like in most PSR problems <cit.> (3) The two structures are in many cases inherently different; they are not supposed to be similar whereas typical PSR problems aim to match different representations of a same object e.g., modeled and measured data, data from different instruments, two pieces of a broken object, etc. <cit.> (4) Every point, its mapping and the distance function to minimize are physically significant. For that reason, no point can be ignored or considered as an outlier, each atom needs to map to exactly one atom (in the case of phase transformations) and the cost function must be relevant to the problem at hand (phase transformation, interfaces, similarity measurement, etc.).
Therefore, even though the distance minimization part of our algorithm is a form of ICP, other existing PSR algorithm could not be used directly to find the optimal one-to-one mapping and alignment between crystal structures.




 §.§ Atom-to-Atom Mapping





Let us explain in more details how the permutation p of atomic indices is optimized. At a given Mapping Loop iteration, the position of each atom in both sections of the A and B structures is known. The goal is to assign each atom of the mapping structure to an atom in the mapped structure such that the sum of the distances between the pairs of corresponding atoms is minimized. As cleverly noted by Sadeghi and Goedecker <cit.>, this is exactly analogous to the assignment problem, a well-studied mathematical problem for which there exist an exact solution that can be computed in polynomial time <cit.>.

The assignment problem consists of finding an optimal way to assign agents to tasks, e.g., clients (tasks) to their taxis (agents) such that the total distance traveled by all the taxis (cost) is minimized. Once the position of each atom is known, the structure mapping problem is exactly equivalent to the assignment problem. The algorithm needs to assign each atom of the mapping structure (or task), shown in red in Fig. <ref>, to an atom in the mapped structure (or agent), shown in gray, such that the total distance (cost) is minimized. The naive route to solving this problem is to try all possible assignments of atoms, but this operation scales as N!. Instead, our algorithm uses the existing Kuhn-Munkres method <cit.> (also known as the Hungarian Algorithm) that solves the assignment problem in polynomial time. This algorithm takes as an input a cost matrix which consists of distances between each possible pair of atoms between the two structures.



Intuitively, one could think that the best way to apply the Hungarian algorithm is to map all A atoms to all B atoms. This is, in fact, problematic because there is no guarantee that the boundaries of the two spherical sections are perfectly compatible. In other words, if at the boundary of B an atom needs to be mapped in the most optimal way to an atom that is outside the boundary of A (it is not part of the finite section created at step 1), it will have to be mapped to some other atom of A regardless. This will lead to an unwanted, exaggerated, influence of the boundaries on the final result. To prevent this from happening, the mapping structure (B) is made smaller than the mapped structure (A) by making the bottom portion of the cost matrix costless (see Fig. <ref>). This means that the mapping of the outer shell of the mapping structure has no effect on the total cost and that these atoms can be considered nonexistent, they are simply placeholders. In other words, there are more agents then there are tasks; some agents will be assigned the task “do nothing.” Thereby, since there are now less atoms in the mapping structure than in the mapped structure, each atom at the boundary of the mapping structure, can find its true counterpart in the mapped structure (provided that the mapped structure is large enough).

This inevitably leads to a new problem: there is no guarantee that all the atoms close to the center of the mapped structure will actually be mapped. In other words, some atoms of A can be “skipped" by the algorithm. When studying polymorphic transformation, this can be problematic since atoms cannot disappear when going from the initial to the final structure, i.e., every atom needs to be mapped. To avoid this problem, a very high cost can be given for mapping placeholders at the outside of the mapping structure to core atoms (important atoms) inside the mapped structure (see Fig. <ref>). When an atom of the mapped structure is not mapped, it is, in fact, mapped to a placeholder in the mapping structure (it is assigned the task “do nothing"). Therefore, imposing a very high cost to mapping core atoms to placeholders will prevent those atoms from not being mapped. Or in terms of tasks and agents: important agents cannot be assigned the task “do nothing." Adding core atoms enforces the one-to-one mapping (bijection) between the core of the mapped structure and the corresponding subset of atoms in the mapping structure. The choice of core atoms (if any) and the relative size of the mapping structure compared to the mapped structures are also the parameters external to the algorithm (set by the user).

This concludes the distance minimization part of the algorithm (upper part in Fig. <ref>), which leads to the optimal alignment and the atom-to-atom mapping between two structures. The result of this stage is depicted in Fig. <ref>(a) for a simple 2D example used to illustrate various aspects of our algorithm.




 §.§ Finding Periodicity



In the previous part of the algorithm, the distance has been minimized, and the optimal mapping has been found. The resulting (p_min, Q_min, t⃗_min) is only applicable to the finite portions of the two crystals that were chosen at the Tiling step. The goal, however, is to describe the matching for the full infinite crystals, which requires finding the periodicity of the map if it exists.
 

We start from a vector field of connections, that is, the vectors that go from the mapping structure (B) to the mapped structure (A) noted  ρ⃗_i = a⃗_i - (Q_minb⃗_p_min(i) + t⃗_min). The idea is to classify equivalent connections into groups, label them, and find the unit cell of the resulting “connection crystal". To do so, the first step is to make the connections periodic. Indeed, even when the mapping p_min is periodic, the connections themselves are not necessarily periodic and, as already discussed, they can diverge in magnitude (see Fig. <ref>(a)). In the example from Fig. <ref> the volumes per atom (areas in 2D) of the structures are exactly the same which implies the existence of a solution with non-diverging connections, but, in general, if the volumes are different the divergence cannot be avoided. 

The connections can be decomposed in two components: (1) a component that accounts for the difference in volumes (stretching/compressing or strain) and (2) a non-diverging, periodic component. The magnitude of the former increases as one moves away from the center of alignment. In order to reveal the periodicity, the non-diverging component needs to be isolated from the divergent one. To do so, keeping the final atom-to-atom mapping p_min fixed, the algorithm minimizes the distance once more, but this time with respect to a linear transformation T and a translation t⃗ where det(T) is unrestricted (step 5, see Fig. <ref>).  An illustration of the resulting connection field after step 5 is presented in  Fig. <ref>(b). If the structures were infinite, minimizing the distance with respect to T would naturally eliminate the diverging component of the connections by making the volume per atom the same in both structures i.e., det(T) = det(C_A)/det(C_B). However, since, in practice, the structures are finite, the condition on det(T) is not exactly fulfilled and the connections are not yet fully periodic.     

To address this problem, the algorithm proceeds to an initial coarse classification of the connections. It is done by placing the connection vectors in different groups with respect to their norm and orientation according to a certain tolerance factor (analogous to bins when making a histogram). On Fig. <ref>, from panel (b) to panel (c), the connections are separated into two groups: blue, pointing up and orange, pointing down. The algorithm then proceeds to making the connections in each group as similar as possible  to each other (in norms and directions) by applying an additional linear transformation to the mapping structure in order to correct the finite size effects introduced at the previous step (step 5). This is done simultaneously for all classes of connections where instead of minimizing the distance function, the algorithm minimizes the class-specific standard deviation (STD) of the connections. This step is represented mathematically by the following equation:


    t⃗, T = t⃗', T'argmin∑_i^N(ρ⃗_i - ∑_j ∈Ω_iρ⃗_j/|Ω_i|)^2

where

    ρ⃗_i = a⃗_i - (T'b⃗_p_min(i) + t⃗'),


and Ω_i is the class that contains i and |Ω_i| denotes the number of elements in that class. The quantity to minimize in eq. (<ref>) is simply a standard deviation with respect to the mean of each class. The classification (step 6) and the minimization of the STD (step 7) are repeated reducing the classification tolerance iteratively until, the STD is equal to zero. This is what we call the Classification Loop on Fig. <ref>. In practice, making connections of each class exactly identical eliminates the remaining diverging component, which is confirmed by verifying that det(T) = det(C_A)/det(C_B).  After this step, the connections are perfectly periodic, and reflect periodicity in the mapping (which has remained the same) as shown in Fig. <ref>(d).

Using the classification of the connections, we can simply proceed as if we were to find the unit cell of a crystal made of connection vectors (instead of atoms). This structure can be described like any other crystal structure by D={ C_D , P_D, L_D}, but in this case L_D is a list of labels that indicates the atomic specie and the class of connection (e.g., blue or yellow on Fig. <ref>(d)). The primitive cell of that structure is the scale of the matching and also an alternative unit cell C_A' of structure A and, consequently, also determines an alternative unit cell C_B' of B.  D and T are the final results of the algorithm. They have the following properties:



    C_D    = C_A' = TC_B' 
    
        P_D   = TP_B' 
    
        P_D    = P_A' - V_D 
    
        L_D     = L_A' = L_B' ,



where A' and B' are alternate representations of A and B; in general they are not the ones that were input initially. V_D is a matrix whose columns are the connection vectors associated with each atomic position. It can easily be constructed from L_D.

A full implementation of our algorithm is available online (see section <ref>).




§ APPLICATIONS





 §.§ Solid-Solid Phase Transformations



To find transformation pathways using our algorithm we set the distance function to be the Euclidean distance between the atoms, we set ε=0 such that no amount of deformation is allowed during the distance minimization step (upper part of Fig. <ref>) and we enforce the one-to-one mapping condition. The algorithm therefore finds the transformation  for which the total distance traveled by all the atoms to go from the initial to the final structure is minimal. The “connection vectors" (arrows) represent the displacements of the atoms during the transition.  The output from the algorithm, D and T, can be used to fully describe the system at any state along the transition path. 

  

We have tested our algorithm on several well-studied transformations. Table <ref> summarizes the results. For hexagonal close-packed (HCP) to body centered cubic (BCC), Graphite to Diamond and Rocksalt to Zincblende, we find pathways that have been previously reported in literature <cit.>. The symmetries of the intermediate structures are exactly the same. For the other three transformations (face centered cubic (FCC) to BCC, Rocksalt to CsCl-type, Roscksalt to Wurtzite), we find new pathways that have not been reported yet. 

In the case of FCC to BCC and Rocksalt to CsCl-type, the newly found pathways involve a slipping process. It has the effect of reducing the total distance between the atoms. By preventing slipping from happening we find exactly the same transformation mechanisms (indicated in parentheses in Table <ref>) that have been reported before. All of the mechanisms that are discussed in literature for these transformations have been derived using periodic boundary conditions; hence, they all suffer from the problem from Fig. <ref>.  For a large number of atoms, they inevitably lead to the greater total travel distances than the same mechanisms with the added slipping process (see Fig. <ref> for FCC to BCC). As we have mentioned before, the component of the displacement of the atoms (or connections) that has the most impact on the total distance is the one associated with strain. Therefore, by minimizing distance our algorithm also minimizes strain. The slipping process appears naturally because it reduces the strain associated with the transformation. In regard to Rocksalt to Wurtzite, our new pathway does not involve a simple slipping mechanism, but a more complex process which, in turn, also leads to a shorter travel distance and consequently smaller principal strains than the path with symmetry Cmc2_1 reported in Refs. <cit.> (see Fig. <ref>). These effects are not present in the HCP to BCC, Graphite to Diamond and Rocksalt to Zincblende, for which our algorithm agrees with the mechanisms commonly discussed in the literature, because these mechanisms already minimize the strain.

These considerations show that: (1) we have reached our goal of creating an algorithm that finds the true path of minimal distance since we either find known pathways or new pathways of shorter total distance and that (2) the result from our algorithm can not only be used as a starting point for ssNEB, but it can also be interpreted directly to explain certain features of the transformation.



For example, let us analyze in more detail the FCC to BCC transformation in iron also known as the martensitic transformation. The martensitic transformation is the diffusion-less transformation of steel from the cubic face centered (FCC) austenite (γ) phase to the body-centered cubic (BCC) or body-centered tetragonal (BCT) martensite (α). For simplicity, we considered the transition of pure iron from FCC to BCC. For austenite, we used a lattice parameter of a_γ=0.3585Å and we defined the lattice parameter in martensite as a_α=√(2/3)a_γ=2.2927Å such that the closed pack directions have the same atomic density for both structures. 

Figure <ref> shows the martensitic transformation found by our algorithm. The transformation consists of a main shear of the (112)_α planes in the [111]_α direction with slip planes every six layers. Between the slip planes, the intermediate structure has the I4/mmm space group which correspond to a Bain distortion accompanied by a rotation. Transformation mechanisms that involve a rotated Bain deformation have been widely theorized <cit.>. As we mentioned, the occurrence of slip planes can be explained by the fact that they greatly reduce the strain necessary to carry out the transformation. The principal strains for our new mechanism are -5.7%, 0% and 15.5%, whereas they would be -18.4%, 15.5% and 15.5% without the occurrence of slip planes (Bain distortion). This slipping process is often used in the context of the Phenomenological Theory of Martensitic Transformation <cit.> to explain the occurrence of striations along the (112)_α; our algorithm finds it naturally by minimizing distance. 

The number of layers between the slip planes depends on the ratio between the parameters of the initial and final structure. The connection structure is composed of 6 atoms which means that the transformation occurs at a scale that corresponds to 6 primitive cells of the two end structures (they have the same number of atoms). Once again, our algorithm behaves as expected by finding a transformation that reduces the total travel distance–and thus the strain–and by being able to find transitions that occur on a larger scale. A more detailed analysis of our results for the martensitic transformation will be published elsewhere <cit.>.

In this study of the martensitic transformation, we used 1000 random initialization steps (Random Initialization Loop), a mapping structure of 180 atoms and a mapped structure of 600 atoms with the one-to-one mapping enforced. Those parameters ensure that the finite portion of the crystal is much larger than the connection cell and that the global minimum is reached. The calculation took 6 min 41 s on a 36-core Intel Xeon Gold 6154 (3.00 GHz) node and 2 h 12 min 44 s on a 4-core Intel Core i7-8550U (1.80 GHz) Laptop. Using the same parameters, the minimization for HCP to BCC, Rocksalt to CsCl-type and Rocksalt to Zincblende were done in 7 min 31 s, 2 min 47 s and 1 min 14 s respectively on the computing node. For the transition from graphite to diamond which involves a large change in specific volumes, we used a mapping structure of 160 atoms and a much larger mapped structure of 1600 atoms. This was done to ensure that there was a sufficient number of graphite layers in the mapped structure. We also increased the number of random initial steps to 3000 in order to find the global minimum. This calculation was completed in 29 h 54 min on the compute node. Similarly, for the transition from Rocksalt to Wurtzite we used a mapping structure of 600 atoms and a mapped structure of 2000 atoms with 2000 initial steps; the calculation was completed in 5 h 14 min. The limiting factor in the distance minimization is the resolution of the assignment problem, therefore the complexity of the mapping dictates the total execution time. Future development may involve the use of an implementation of the Hungarian algorithm parallelized for graphics processing units (GPUs) <cit.>.





 §.§ Semi-Coherent Interfaces



Next, we illustrate how our algorithm can also be used to find the structures of semi-coherent interfaces between different materials. In the examples that follow, we consider only the terminating planes in each structure. Therefore, in our algorithm the two structures are modeled as large disk-like 2D sections of the terminating planes (instead of spheres in 3D). For demonstration purposes, here, the plane directions and terminating layers are taken from experiment. Each connection between an atom from the mapped structure and an atom from the mapping structure represents a chemical bond. Atoms no longer have to be mapped to atoms of the same specie, they are mapped according to chemistry rules that determine which types of atoms from one structure will bond to which type of atoms from the other structure e.g., Zr atoms bond with Ni atoms. These rules need to be known in advance. Since connections now represent chemical bonds, for the distance metric in equation <ref>, we use the Lennard-Jones potential:


    d(a⃗,b⃗) = σ(r^12/||a⃗-b⃗||^12 - 2r^6/||a⃗-b⃗||^6)


where σ denotes the potential strength and r the equilibrium radius. We use this potential because it is a mathematically simple representation of the general shape of the potential between 2 atoms. In our model, atoms are bonded with at most 1 atom of the other phase. In other words, we assume that the bond with the closest neighboring atom of the other phase is the strongest and most consequential in terms of energy and alignment. Since we are only interested in the optimal alignment–we are not trying to predict the interfacial energy, the strength of the potential σ is not important and it is set to 1. Thus, the potential has only one parameter: the equilibrium radius r. It can be set based on physical or experimental arguments.
Moreover, there is no need to enforce that each and every atom of both structures form a bond. In the case of semi-coherent interfaces for example, the lattice constant of the two materials can be very different such that only a fraction of the atoms at the interfaces will form bonds. Because, in our model, an atom can form at most 1 bond (1 or 0 bond), there cannot be more bonds per unit area than there are atoms per unit area in the structure that is the least dense. Since, in general we wish to maximize the number of bonds per unit area, all the atoms of the least dense phase need to form a bond. This is done by setting the denser structure as the mapped structure and by setting the fraction of core atoms in the cost matrix to 0 such that the one-to-one mapping is not enforced. In fact, in this case, we take advantage of the fact that certain atom will naturally be “skipped" when making the mapping structure smaller than the mapped structure.
Finally, during the distance minimization step (in this case the distance is defined by equation <ref>), the structures may be slightly strained in-plane near the interface in order to maximize the bonding energy. Therefore, we usually set ε to a value between 3-8%.



We used our algorithm to find the orientation between two experimentally realized interfaces.
The first system is a solid-solid interface between Ni and yttrium-stabilized zirconia (YSZ). This interface was experimentally realized and studied by Nahor et al. <cit.>. We used the parameters from their experiment to run our simulation. For face-centered cubic Ni, we used a lattice parameter of 3.52Å and for cubic ZrO_2, we used a lattice parameter of 5.125Å (apart from its effect on the lattice parameter, the presence of yttrium was not considered in our simulation). We specified the interfacial plane (111) for both structures and used the termination (Zr for ZrO_2) specified  in Ref. <cit.>. We set the equilibrium point of the Lennard-Jones potential r to be 2.5Å, because it is close to the Ni-Ni interatomic distance of 2.49Å. The value of the equilibrium point is an estimate of the length of the bonds and therefore it determines the distance between the layers of the two phases. We find that its value does not have a strong incidence on the final result. We set ε to 0.08 such that there can be some strain in the layers close to the surface. 

We find Ni<110>(111)//ZrO_2 <110>(111) [Our algorithm does not differentiate between OR_1 and OR_2 because they are identical in-plane.] to be the optimal alignment between the two structures in accordance with the experimental observation.  Fig. <ref>(a) shows the interface viewed in the [111] projection (from above). The algorithm finds a repeating pattern of only a few unit cells in sharp contrast with the 33×33 unit cells O-lattice found using the measured orientation relationship <cit.>. This smaller cell has the advantage of directly providing the matching modes of the interfaces: 2:3 Zr-Ni which is in accordance with the “one dislocation every three Ni planes" observed by Nahor et al. This smaller cell is possible because the optimal result was found by allowing some strain in the interface layers. In fact, in the result shown in Fig. <ref>, the ratio between the area of the Zr unit cells in-plane and the Ni unit cells in-plane is increased by 6.12%. In other words, there is tension in the YSZ side and compression on the Ni side also in agreement with the observation of Nahor et al. Not only does our algorithm correctly reproduce the measured orientation relationship solely using the in-plane lattice parameters, but it also provides the matching mode and the general direction of the strain (tension or compression) in the layers near the interface. A projection of our model along the Ni[110] direction with the aforementioned in-plane strain is presented in Fig. <ref>(b) in comparison with an HRTEM micrograph of the same projection. The model is in very good agreement with the experimental result.

The second system is the solid-solid interface between Si and SiC. This interface was experimentally realized and studied by Li et al. <cit.>, once again, we used the parameters from their experiment to run our simulation. For diamond Si we used a lattice parameter of 5.43Å and for hexagonal close-packed (HCP) 6H-SiC we used a lattice parameter of 3.08Å. We set ε to 0.08. As specified in Ref. <cit.>, we matched the (110) plane of Si with the Si terminated (001) plane of SiC. We set the equilibrium point of the Lennard-Jones potential (r) to be 2.5(̊A), because it is close to the Si-Si interatomic distance of 2.35(̊A).

We find the Si<110>(110)//6H-SiC<110>(001) orientation to be the optimal alignment between the two structures in accordance with the experimental observation. Fig. <ref>(c) shows the interface viewed in the Si[110] projection (from above). Once again, our periodic pattern is in accordance with the observed 4:5 Si to SiC matching mode in the Si[111]/SiC[110] direction. In addition, the algorithm finds a 1.68% increase in the ratio between the in-plane Si(SiC) cell and the in-plane Si(Si) cell. In other words, the 6H-SiC structure is stretched and/or the Si structure is compressed at the interface to obtain that ratio. This could explain the 1.84% mismatch in the Si[001] direction and the 0.26% residual mismatch in the Si[110] direction noted by Li et al. A projection of our model along the Si[001] direction with the aforementioned in-plane strain is presented on Fig. <ref>(d) in comparison with an HRTEM micrograph of the same projection; the model is in very good agreement with the experimental result.

We obtained the result presented above using the same type of 36-core Intel Xeon computing nodes. For Ni//YSZ the minimum was obtained in 1 min 35 s using a mapping structure of 325 atoms and a mapped structure of 1300 atoms with 1000 initializations. For Si//SiC we used a mapping structure of 150 atoms and a mapped structure of 600 atoms with the same number of random initializations. The calculation was completed in 13 min 39 s. 




§ MEASURING DISTANCE BETWEEN CRYSTAL STRUCTURES



As discussed previously, defining a rigorous Euclidean distance between crystal structures or more broadly between infinitely periodic arrays of points is a challenging task. First, the infinite dimensionality of the configuration space poses problem. This difficulty can, in principle, be avoided by scaling the metric with some function of the number of atoms N. As we will show in this section, it is actually not possible as the dependence on the number of atoms involves different powers of N. Secondly, and as importantly, even a finite portion of a crystal structure is not represented by a unique point in the N-dimensional configuration space of atomic coordinates, but by several points that reflect: (a) the permutations of atomic indices that describe the same crystal structure, and (b) the variability in the choice of the N-atom section of an infinitely periodic crystal. The  definition of a distance metric between two crystal structures for any fixed N implies finding the two closest representative points of the two structures in the N-dimensional configuration space, a task that is tackled by our algorithm. In fact, once the optimal parameters (p_min, Q_min,t⃗_min) defined in equation (<ref>) are found, for a fixed N, the minimized distance may serve as a mathematical metric between periodic structures.

Let us consider the distance between two structures in the situation where the correspondence between them has already been established. Since the mapping is periodic, the two structures (A and B) in their optimal matching can be described with cells C_A and C_B which both contain m atoms and are optimally aligned, and with atomic positions {a⃗_i  |  i=1,2,…,m} and {b⃗_i  |  i=1,2,…,m} inside the cells indexed according to the optimal mapping. The shortest travel distance between the two structures with this match is: 


    d_1 = ∑_l=1^m∑_i,j,k = -n/2^n/2||a⃗_ijkl - b⃗_ijkl|| ,


where a⃗_ijkl stands for the position of an atom that belongs to the structure A. a⃗_ijkl is a periodic image of an atom with an index l located in the unit cell indexed with ijk. More precisely: 


    a⃗_ijkl = C_Ai
    j
    k + a⃗_l,


and analogously for the structure B. The total number of atoms in each structure is N=m(n+1)^3. This distance (d_1) is, by construction, the l_1,2-norm <cit.> of the 3 × N matrix formed by the connection vectors. We used the same norm when posing the matching problem for phase transitions.
So far, we used the l_1,2-norm because it represents the sum of the distances traveled by all the atoms during the transition, but one could also be interested in computing the Frobenius norm of that matrix i.e., the l_2-norm of the vector joining the two structures in configuration space. It is given by:


    d_2 = √(∑_l=1^m∑_i,j,k = -n/2^n/2||a⃗_ijkl - b⃗_ijkl||^2) .


Adding a 1/N factor inside the square root in front of the summation gives the root mean square distance (RMSD). It is important to note that, for finite N, the set of optimal parameters (p_min, Q_min,t⃗_min) is not necessarily the same for d_1 and d_2 
When they are optimized, both d_1 and d_2 fulfill the 4 requirements of a metric: 
 
 
 
  * d(A,B)≥0, 
 
  * d(A,B) = 0  A=B, 
 
  * d(A,B) = d(B,A) and 
 
  * d(A,C) ≤ d(A,B) + d(B,C). 
 
 
 The first 3 criteria follow trivially from the properties of the l_2-norm applied to the connection vectors. The fourth criterion follows from the fact that both d_1(A,C) and d_2(A,C) are defined respective to their (p_min, Q_min,t⃗_min) and represent the shortest distance between A and C.
The problem is that both d_1 and d_2 depend on N and in order to compare actual structures, we either need: (1) to derive quantities from d_1 and d_2 that are independent of N, or (2) find a way to compare distances in the limit where N→∞.
In both cases, the first step is to derive the dependence of d_1 and d_2 on N.




 §.§ Size dependence
  


The case of d_2 can be derived analytically so we will use it for demonstration purposes.
Let us first define ρ⃗_l = a⃗_l - b⃗_l and C' = (C_A - C_B) whose columns are {c⃗_ν  |  ν=1,2,3 }. From equation (<ref>), we can now write:


    d_2^2 = ∑_l=1^m∑_i,j,k = -n/2^n/2||C'i
    j
    k + ρ⃗_l ||^2 .


If the cells are exactly identical, then C'=0 and it follows that:


    d_2^2 = (n+1)^3∑_l=1^m ||ρ⃗_l||^2 = K_2N ,

where

    K_2 = 1/m∑_l=1^m ||ρ⃗_l||^2.


If the cells are not identical, C' is an invertible matrix and we can write:


    d_2^2 =   ∑_l=1^m∑_i,j,k = -n/2^n/2||C'i+ρ'_l^(1)
    j+ρ'_l^(2)
    k+ρ'_l^(3)||^2 ,


where {ρ'_l^(ν)} are the elements of C^-1ρ⃗_l. Then, using the fact that the vector norm squared is equivalent to the inner product of the vector with itself:


    d_2^2    =  ∑_l=1^m∑_i,j,k = -n/2^n/2 (Ic⃗_1 + Jc⃗_2 + Kc⃗_3)·(Ic⃗_1 + Jc⃗_2 + Kc⃗_3) 
       = T_11 + T_22 + T_33 + 2T_12 + 2T_13 + 2T_23 ,


where


    I = i+ρ'_l^(1),  J = j+ρ'_l^(2),  K = k+ρ'_l^(3) .


This gives 6 terms, that can be broken down into two cases. First:


    T_11   = ∑_l=1^m∑_i,j,k = -n/2^n/2 (i+ρ'_l^(1))^2||c⃗_1||^2
       = (n+1)^3n(n+2)m/12||c_1||^2 
       aaaaa + (n+1)^3∑_l=1^m (ρ'_l^(1))^2||c_1||^2


(the T_22 and T_33 cases are similar), and second:


    T_12 =   ∑_l=1^m∑_i,j,k = -n/2^n/2 (i+ρ'_l^(1))(j+ρ'_l^(2))c⃗_1·c⃗_2
    
    =    (n+1)^3∑_l=1^mρ'_l^(1)ρ'_l^(2)c⃗_1·c⃗_2

(the T_13 and T_23 cases are similar).
Finally, equation (<ref>) becomes:


    d_2^2 = (n+1)^3n(n+2)m/12(||c⃗_1||^2+||c⃗_2||^2+||c⃗_3||^2) 
    +  (n+1)^3∑_l=1^m∑_ν^3∑_μ^3 ρ'_l^(ν)ρ'_l^(μ)c⃗_ν·c⃗_μ .


Substituting n=(N/m)^1/3 - 1, and regrouping the constants, we get the following relation:


    d_2^2 = G_2N^5/3 + K_2N  ,


where


    G_2 =   (||c⃗_1||^2+||c⃗_2||^2+||c⃗_3||^2)/12m^2/3  ,  
    
        K_2 =   1/m∑_l=1^m∑_ν^3∑_μ^3 ρ'_l^(ν)ρ'_l^(μ)c⃗_ν·c⃗_μ
          - 1/12(||c⃗_1||^2+||c⃗_2||^2+||c⃗_3||^2)


Note that neither G_2 nor K_2 are dependant on the number of atoms N.
It is immediately evident that dividing by N to any power will not result in a size independent (scaled) metric. Therefore, the RMSD (i.e., d_2^2/N) depends on the size of the system. It cannot be applied to atoms inside a unit cell (or any finite portion of the crystal) to measure distances between periodic structures. G_2 is solely dependent on the difference between C_A and C_B; if the cells remain invariant during the transformation, G_2=0 and equation (<ref>) simplifies to equation (<ref>). We can say that G_2 is associated with the change in unit cells whereas K_2, the linear term, is associated with the displacements inside the cell.

Finding the relation for the l_1,2-norm is slightly more involved (see Appendix <ref>) but we find a similar relationship:


    d_1 = G_1N^4/3 + K_1N + 𝒪(N) ,


where G_1≥0. Once again there is no trivial way to make the distance an intensive quantity because of the presence of a non-linear term associated with the distortion in the unit cell.

Since we are only interested in d_1 and d_2 in the limit of N→∞, one might think that the leading terms in N, G_2 and G_1 respectively, could directly serve as metrics. However, G (G_1 or  G_2) does not fulfill the second criteria of a metric since there could, in principle, exist a transformation that consists of a pure reorganization of the atoms where G=0 even though the end structures are different. The only way to define a proper metric is to use all the parameters in equations (<ref>) or (<ref>) depending on the particular choice. Providing that those parameters are known, the most straightforward approach, is to compare distances in the limit. 
However, since comparing functions in the limit can be tedious and not convenient for computation, we also defined a metric function that uses both parameters (G and K), it is presented in Appendix <ref>. 




 §.§ Practical Use



We chose to study solid-solid phases transitions using d_1 because it represents the sum of the Euclidean distances travelled by all the atoms in the structure. In that context d_1 can also be seen as the true Euclidean distance between two structures and it can be used to measure distances between them. We were not able to find a general closed form for G_1, therefore it is not possible to obtain it directly from the optimal sets {C_A, P_A, L_A} and {C_B, P_B, L_B} that are found by our algorithm. However, we were able to find the general dependence of d_1 on N.  Using the optimal mapping to compute the distance d_1 at different sizes, we can show that the distance indeed grows according to equation (<ref>).



Fig. <ref> shows in blue the dependence of the total travelled distance d_1 on the number of atoms N, for 2 of the 6 transitions presented in Table <ref>. This distance is compared with the one that corresponds to the pathways previously discussed in the literature that do not involve slipping and are hence suboptimal with respect to minimizing the d_1. Panel (a) compares d_1 norms for these two pathways for the transition of ZnO from Rocksalt to Wurtzite. The red points correspond to the distance for the path with symmetry Cmc2_1 reported by Refs. <cit.>. As already noted, our new pathway (in blue) produces a shorter travel distance that grows slower with  N. 

Fitting equation (<ref>) confirms the derived dependence of d_1 on N and shows that G_1=0.05Å for the new pathway is much smaller than G_1=0.26Å for the one found in literature. The fitting curves are shown as black lines; they overlap the simulated data almost perfectly. Similarly, the red points in panel (b) correspond to the Bain deformation of Iron from FCC to BCC. Again, the new pathway has a shorter travel distance. In this case, G_1=0.07Å for the new path and 0.17Å for the Bain path. Equation (<ref>) fits both data sets very well.
Thus, in practice, one can use our algorithm to fit G_1 and K_1 and use them to compare distances in the limit or using our metric presented in Appendix <ref>. 

On the other hand, if one wishes to use d_2 to compare periodic crystal structures, our algorithm, can be used to find the optimal sets {C_A, P_A, L_A} and {C_B, P_B, L_B}  by setting the distance function as the square of the euclidean norm. Then, one can use the closed form for G_2 and K_2 provided in equation (<ref>), equation (<ref>) and equation (<ref>) to compare distances in the limit or using our metric presented in Appendix <ref>.


In sum, in this section, we showed that distances such as the RMSD cannot be used directly to compare periodic structure because they depend on the size of the system. Instead, we established the dependence on N for two metrics d_1 and d_2 and we showed how they can be used, in the limit, to quantify the similarity between structures.




§ CONCLUSIONS



In this work, we formulated the matching of two different crystal structures as an optimization problem and described our algorithmic solution to it. The methodology that we developed, inspired by the Iterative Closest Point, is constructed to work on large and finite portions of the two crystal structures rather than on some choice of a periodic unit. It consists of a sequential minimization of a given distance function with respect to the permutations of atomic indices and linear transformations (rotations and translations) of the atomic positions. The sequence is repeated iteratively until the convergence is achieved. After the optimal alignment of the structures and the optimal atom-to-atom map are found, our algorithm analyzes the result and retrieves the periodicity in the match. This last step ensures that the boundaries have no influence on the final result. 

We presented two different implementations of our algorithm tailored for their respective class of applications. First, we demonstrated our algorithm's relevance when studying phase transformations by examining six well-studied transformations. In each case, we either confirmed an existing mechanism or uncovered a new lower-strained pathway. In particular, for the martensitic transformation,  we found a new modified version of the Bain path that does not require large expansion along certain crystallographic directions. Then, we showed that, starting solely from the in-plane lattice parameters, our algorithm was capable of reproducing the features of experimental interface structures such as their orientation relationships, matching modes and strain directions for two case examples: Ni on YSZ and Si on SiC.

Finally, we analyzed and discussed a practical formulation of a rigorous distance metric between crystal structures that can be used to assess their Euclidean “closeness”.



§ DATA AVAILABILITY
 
A full implementation of our algorithm is available via Github at https://github.com/ftherrien/p2ptrans.
All the parameters necessary to reproduce the examples presented in this study are  included in this article. All other relevant data is available from the corresponding author on request.



§ ACKNOWLEDGEMENTS

This work was supported by the Center for the Next Generation of Materials by Design, an Energy Frontier Research Center funded by the U.S. Department of Energy, Office of Science, Basic Energy Sciences; and by the National Science Foundation grant No. DMR-1945010. The research was performed using computational resources sponsored by the Department of Energy's Office of Energy Efficiency and Renewable Energy, located at the National Renewable Energy Laboratory.











§ AN ADDITIONAL METRIC
 
Since using the limit might be inconvenient computationally, let's define a metric 𝒟 as such:


    𝒟(A,B) = 
     
          e-e^-K   if  G=0 
    
          e^G   otherwise


This fulfills the four criteria:

    
  * G≥0 and K≥0 if G=0, therefore 𝒟(A,B)≥0
    
  * 𝒟(A,B)=0  G=K=0  A=B
    
  * d_1,2(A,B) = d_1,2(B,A) 𝒟(A,B) = 𝒟(B,A)
    
  * 
        
  * if G(A,C),G(A,B),G(B,C) > 0:
        
        Follows from the fact that d_1 is a metric and e^G is monotone
        
  * if G(A,C),G(A,B) > 0, G(B,C) = 0 (or G(A,C),G(B,C) > 0, G(A,B) = 0):
        
        Then C_B = C_C  G(A,C)=G(A,B)=G  e^G≤ e^G + e-e^-K(B,C)   
        
  * if G(A,B),G(B,C) > 0, G(A,C) = 0:
        
        e-e^-K(A,C)≤ e^G(A,B) + e^G(B,C) since e-e^-K(A,C)≤ 1 and e^G(A,B)≥ 1
        
  * if G(A,B),G(B,C),G(A,C) = 0:
        
        Follows from the fact that d_1 is a metric and e-e^-K is monotone
        
  * if two of the three Gs are equal to 0
        
        Then C_A = C_B = C_C which is equivalent to the previous case
    




§ N-DEPENDENCE OF D1
 
Once again let's define C' = (C_A - C_B) whose columns are {c⃗_⃗ν⃗  |  ν=1,2,3} and ρ⃗_l = a⃗_l - b⃗_l, we can rewrite equation (<ref>):


    d_1 = ∑_l=1^m∑_i,j,k = -n/2^n/2||C'i
    j
    k + ρ⃗_l|| .


If the cells are exactly identical, C'=0 and it follows that:


    d_1 = (c+1)^3∑_l=1^m ||ρ⃗_l|| = G_1N ,

where

    G_1 = 1/m∑_l=1^m ||ρ⃗_l||.


Let's consider the simpler case where the vectors of C_1 and C_2 are orthogonal and where ρ⃗_l=0. One can show that aligning the 3 vectors is the optimal alignment C_1 and C_2 and therefore, the vectors of C' are also orthogonal. Equation (<ref>) becomes:     


    d_1 = ∑_i,j,k = -n/2^n/2√(i^2||c⃗_1||^2 + j^2||c⃗_2||^2 + k^2||c⃗_3||^2) .


Let's approximate the summation using the right Riemann sum:


    ∑_i=1^nf(a+i(b-a)/n)(b-a)/n
    = ∫_a^b f(x) dx + 𝒪(M(b-a)^2/n) ,


where M=max f'(x). In our case, a=-n/2, b=n/2 and:


    f(x)    = √(x^2||c⃗_1||^2 + j^2||c⃗_2||^2 + k^2||c⃗_3||^2)
    
    f'(x)    = 2 x||c⃗_1||^2/2 √(x^2||c⃗_1||^2 + j^2||c⃗_2||^2 + k^2||c⃗_3||^2)≤ ||c⃗_1|| .


Replacing in equation (<ref>), we get:


    d_1 = ∑_j,k = -n/2^n/2∫_-n/2^n/2√(x^2||c⃗_1||^2 + j^2||c⃗_2||^2 + k^2||c⃗_3||^2) dx 
    + 𝒪(||c⃗_1||n) .


Using the same argument for j and k:


    d_1 = ∭_Ω√(x^2||c⃗_1||^2 + y^2||c⃗_2||^2 + z^2||c⃗_3||^2) dx dy dz 
    + 𝒪((||c⃗_1||+||c⃗_2||+||c⃗_3||)n(n+1)^2) .


Where Ω is a cube of parameter n centered at the origin. Let's define x' = x||c⃗_1||, y' = y||c⃗_2|| and z' = z||c⃗_3||. The integral becomes:


    d_1 = 1/V∭_Ω'√(x'^2 + y'^2 + z'^2)dx' dy' dz' 
     + 𝒪(n^3) ,


where Ω' is a prism of parameter n||c⃗_1||, n||c⃗_2|| and n||c⃗_3|| centered at the origin and V = ||c⃗_1||·||c⃗_2||·||c⃗_3||. This integral can be carried out in spherical coordinates by carefully adjusting the integration limits:


    d_1 = 1/V∑_i<j<k^3∫_0^L_ϕ∫_0^L_θ∫_0^L_r  r^3sinθ dr dθ dϕ + 𝒪(n^3) ,


where:


    L_ϕ   = tan(||c⃗_j||/||c⃗_i||) ,
    
    L_θ   = arccos(1/√(1+||c⃗_i||^2/||c⃗_k||^2^2ϕ)) ,
    
    L_r    = n||c⃗_k||/2cosθ.



Integrating:


    d_1 = n^4/192 V∑_i<j<k^3 ∫_0^tan(||c⃗_j||/||c⃗_i||)||c⃗_k||^4
    ((1+||c⃗_i||^2/||c⃗_k||^2^2ϕ)^3/2-1) dϕ + 𝒪(n^3) .


Replacing n=N^1/3 - 1, and regrouping the constants:


    d_1    = G_1(N^1/3 - 1)^4 +  𝒪((N^1/3 - 1)^3) 
       = G_1(N^4/3 - 4 N + 6 N^2/3 - 4 N^1/3 + 1) + 𝒪(N) 
       = G_1 N^4/3 + 𝒪(N) ,


where 


    G_1 = 1/192 V∑_i<j<k^3 ∫_0^tan(||c⃗_j||/||c⃗_i||)||c⃗_k||^4((1+||c⃗_i||^2/||c⃗_k||^2^2ϕ)^3/2-1) dϕ


In sum, a pure reorganization of the atoms that does not change the unit cell (equation (<ref>)) makes the distance depend linearly on the size whereas a pure distortion of the cell makes the distance non-linearly dependant on the size. Considering the more general case where the cell vectors are not orthonormal would only lead to lower order terms. To highlights the two leading contributions we write:


    d_1 = G_1N^4/3 + K_1N + 𝒪(N) .










60
fxundefined [1]
 ifx#1

fnum [1]
 #1firstoftwo
 secondoftwo
 
fx [1]
 #1firstoftwo
 secondoftwo
 





noop [0]secondoftwo

ref[1]@startlink#1@href
href[1]#1@endlink
anitize@url [0]`
12`$12`&12`#12`1̂2`_12`%12
startlink[1]
endlink[0]

rl [1]href #1












@bib@innerbibempty

[Brune(2014)]brune:2014
  author author H. Brune, @noop  title Surface and Interface
  Science (publisher John Wiley & Sons, Ltd, year 2014) Chap. chapter 20, pp. pages
  421–492NoStop
[Poeppelmeier and Rondinelli(2016)]rondo_NC:2016
  author author K. R. Poeppelmeier and author J. M. Rondinelli, title title Mismatched
  lattices patched up, https://doi.org/10.1038/nchem.2477
  journal journal Nature Chemistry volume 8, pages 292 EP – (year
  2016)NoStop
[Ding et al.(2016)Ding,
  Dwaraknath, Garten, Ndione,
  Ginley, and Persson]persson_ACS:2016
  author author H. Ding, author S. S. Dwaraknath, author L. Garten,
  author P. Ndione, author D. Ginley,  and author
  K. A. Persson, title
  title Computational approach for epitaxial polymorph
  stabilization through substrate selection, booktitle booktitle ACS Applied Materials & Interfaces, 10.1021/acsami.6b01630 journal journal
  ACS Applied Materials & Interfaces volume 8,  pages 13086–13093 (year 2016)NoStop
[Henkelman, Uberuaga, and  Jónsson(2000)]henkelman_2000
  author author G. Henkelman, author B. P. Uberuaga,  and author H. Jónsson, title title A climbing
  image nudged elastic band method for finding saddle points and minimum energy
  paths, 10.1063/1.1329672 journal journal The Journal of Chemical Physics volume
  113, pages 9901–9904 (year 2000)NoStop
[Sheppard et al.(2012)Sheppard, Xiao, Chemelewski, Johnson, and Henkelman]sheppard_2012
  author author D. Sheppard, author P. Xiao,
  author W. Chemelewski, author D. D. Johnson,  and author G. Henkelman, title
  title A generalized solid-state nudged elastic band
  method, @noop  journal journal The
  Journal of chemical physics volume 136, pages 074103 (year 2012)NoStop
[Caspersen and Carter(2005)]caspersen_2005
  author author K. J. Caspersen and author E. A. Carter, title title Finding
  transition states for crystalline solid–solid phase transformations,  @noop  journal journal Proceedings of the
  National Academy of Sciences volume 102, pages 6738–6743 (year 2005)NoStop
[Qian et al.(2013)Qian,
  Dong, Zhou, Tian,
  Oganov, and Wang]qian_2013
  author author G.-R. Qian, author X. Dong, author X.-F. Zhou, author
  Y. Tian, author A. R. Oganov,  and author H.-T. Wang, title title Variable cell nudged elastic band method for studying solid–solid
  structural phase transitions, @noop  journal
  journal Computer Physics Communications volume 184, pages 2111–2118 (year
  2013)NoStop
[Note1()]Note1
  note CSL: Coincidental Site Lattice, DSC: @underline
  D@isplacement of one crystal Lattice with
  respect to the second causes a pattern @underline S@hift which is @underline C@ompleteNoStop
[Bollmann(1970)]bollmann1
  author author W. Bollmann, @noop  title Crystal Defects and
  Crystalline Interfaces, edition 1st ed. (publisher Springer-Verlag Berlin Heidelberg, year 1970)  pp. pages 83–97NoStop
[Bollmann(1974)]bollmann_1974
  author author W. Bollmann, title title O-lattice
  calculation of an f.c.c.-b.c.c. interface, 10.1002/pssa.2210210218 journal journal Physica
  Status Solidi (a) volume 21, pages
  543–550 (year 1974)NoStop
[Smith and Pond(1976)]smith_pond_1976
  author author D. A. Smith and author R. C. Pond, title title Bollmann's
  0-iattice theory; a geometrical approach to interface structure, 10.1179/imtr.1976.21.1.61 journal journal International Metals Reviews volume 21,  pages 61–74 (year 1976), http://arxiv.org/abs/https://doi.org/10.1179/imtr.1976.21.1.61
  https://doi.org/10.1179/imtr.1976.21.1.61 NoStop
[Balluffi, Brokman, and  King(1982)]balluffi_brokman_king_1982
  author author R. Balluffi, author A. Brokman,   and author A. King, title title Csl/dsc lattice model for general
  crystalcrystal boundaries and their line defects, 10.1016/0001-6160(82)90166-3 journal journal
  Acta Metallurgica volume 30, pages
  1453–1470 (year 1982)NoStop
[Zhang and Kelly(1998)]zhang_kelly_1998
  author author M.-X. Zhang and author P. Kelly,  title title Crystallography and
  morphology of widmanstätten cementite in austenite, 10.1016/s1359-6454(98)00139-6 journal journal
  Acta Materialia volume 46, pages
  4617–4628 (year 1998)NoStop
[Zhang et al.(2005)Zhang,
  Kelly, Easton, and Taylor]zhang_kelly_easton_taylor_2005
  author author M. Zhang, author P. Kelly,
  author M. Easton,  and author J. Taylor, title
  title Crystallographic study of grain refinement in
  aluminum alloys using the edge-to-edge matching model, 10.1016/j.actamat.2004.11.037 journal journal
  Acta Materialia volume 53, pages
  1427–1438 (year 2005)NoStop
[Ikuhara and Pirouz(1996)]ikuhara_pirouz_1996
  author author Y. Ikuhara and author P. Pirouz, title title Orientation
  relationship in large mismatched bicrystals and coincidence of reciprocal
  lattice points (crlp), 10.4028/www.scientific.net/msf.207-209.121 journal journal Materials Science Forum volume 207-209,  pages 121–124 (year 1996)NoStop
[Zur and McGill(1984)]zur1984
  author author A. Zur and author T. C. McGill, title title Lattice match:
  An application to heteroepitaxy, 10.1063/1.333084
  journal journal Journal of Applied Physics  volume 55, pages 378–386 (year 1984)NoStop
[Mathew et al.(2016)Mathew,
  Singh, Gabriel, Choudhary,
  Sinnott, Davydov, Tavazza,  and Hennig]mathew_2016
  author author K. Mathew, author A. K. Singh,
  author J. J. Gabriel, author K. Choudhary, author
  S. B. Sinnott, author
  A. V. Davydov, author
  F. Tavazza,  and author
  R. G. Hennig, title
  title Mpinterfaces: A materials project based python
  tool for high-throughput computational screening of interfacial systems,  @noop  journal journal Computational
  Materials Science volume 122, pages
  183–190 (year 2016)NoStop
[Jelver et al.(2017)Jelver,
  Larsen, Stradi, Stokbro,  and Jacobsen]jelver2017determination
  author author L. Jelver, author P. M. Larsen,
  author D. Stradi, author K. Stokbro,  and author K. W. Jacobsen, title title Determination of low-strain interfaces
  via geometric matching, @noop  journal journal Physical Review B volume 96, pages 085306 (year 2017)NoStop
[Sadeghi et al.(2013)Sadeghi, Ghasemi, Schaefer, Mohr, Lill, and Goedecker]sadeghi_2013
  author author A. Sadeghi, author S. A. Ghasemi, author B. Schaefer,
  author S. Mohr, author
  M. A. Lill,  and author
  S. Goedecker, title title Metrics for measuring distances in configuration
  spaces, @noop  journal journal The
  Journal of chemical physics volume 139, pages 184118 (year 2013)NoStop
[Note2()]Note2
  note The l_2-norm in configuration space equivalent to the
  Forbenius norm of the position matrixNoStop
[Oganov and Valle(2009)]oganov_JCP:2009
  author author A. R. Oganov and author M. Valle,  title title How to quantify energy
  landscapes of solids, 10.1063/1.3079326 journal journal The Journal of Chemical Physics volume 130, pages 104504 (year
  2009)NoStop
[Bartók, Kondor, and  Csányi(2013)]bartoK_2013
  author author A. P. Bartók, author R. Kondor,
   and author G. Csányi,  title title On representing chemical
  environments, @noop  journal journal
  Physical Review B volume 87, pages
  184115 (year 2013)NoStop
[Yang, Dacek, and Ceder(2014)]yang_2014
  author author L. Yang, author S. Dacek,  and  author G. Ceder, title title Proposed definition of crystal
  substructure and substructural similarity, @noop  journal journal Physical Review B volume 90, pages 054102 (year
  2014)NoStop
[De et al.(2016)De,
  Bartók, Csányi, and Ceriotti]de_2016
  author author S. De, author A. P. Bartók,
  author G. Csányi,  and  author M. Ceriotti, title title Comparing molecules and solids across
  structural and alchemical space, @noop  journal
  journal Physical Chemistry Chemical Physics volume 18, pages 13754–13769 (year
  2016)NoStop
[Zhu et al.(2016)Zhu,
  Amsler, Fuhrer, Schaefer,
  Faraji, Rostami, Ghasemi,
  Sadeghi, Grauzinyte, Wolverton et al.]zhu_2016
  author author L. Zhu, author M. Amsler,
  author T. Fuhrer, author B. Schaefer, author
  S. Faraji, author S. Rostami, author S. A.  Ghasemi, author A. Sadeghi, author M. Grauzinyte,
  author C. Wolverton,  et al., title title A fingerprint
  based metric for measuring similarities of crystalline structures,  @noop  journal journal The Journal of
  chemical physics volume 144, pages
  034203 (year 2016)NoStop
[Niggli(1928)]niggli1928handbuch
  author author P. Niggli, @noop  title Handbuch der
  Experimentalphysik, Vol. volume 7 (publisher
  Akademische Verlagsgesellschaft, year 1928) note
  part 1NoStop
[Santoro and Mighell(1970)]santoro1970determination
  author author A. t. Santoro and author A. Mighell, title title Determination
  of reduced cells, @noop  journal journal Acta Crystallographica Section A: Crystal Physics, Diffraction,
  Theoretical and General Crystallography volume 26,  pages 124–127 (year 1970)NoStop
[Krivý and Gruber(1976)]krivy1976
  author author I. Krivý and author B. Gruber, title title A unified
  algorithm for determining the reduced (niggli) cell, 10.1107/s0567739476000636 journal journal Acta
  Crystallographica Section A volume 32, pages 297–298 (year 1976)NoStop
[Stevanović et al.(2018)Stevanović, Trottier, Musgrave,
  Therrien, Holder, and Graf]stevanovic2018
  author author V. Stevanović, author R. Trottier, author C. Musgrave,
  author F. Therrien, author A. Holder,  and author
  P. Graf, title title Predicting kinetics of polymorphic transformations from
  structure mapping and coordination analysis, 10.1103/physrevmaterials.2.033802 journal journal
  Physical Review Materials volume 2 (year 2018), 10.1103/physrevmaterials.2.033802NoStop
[Lonie and Zurek(2012)]lonie_2012
  author author D. C. Lonie and author E. Zurek,  title title Identifying duplicate
  crystal structures: Xtalcomp, an open-source solution, @noop 
  journal journal Computer Physics Communications  volume 183, pages 690–697 (year 2012)NoStop
[Capillas, Perez-Mato, and  Aroyo(2007)]capillas_2007
  author author C. Capillas, author J. Perez-Mato,  and author M. Aroyo, title title Maximal symmetry
  transition paths for reconstructive phase transitions, @noop 
  journal journal Journal of Physics: Condensed
  Matter volume 19, pages 275203
  (year 2007)NoStop
[Larsen, Schiøtz, and  Schmidt(2017)]larsen2017structural
  author author P. M. Larsen, author J. Schiøtz,
   and author S. Schmidt,  @noop  title Structural analysis algorithms for
  nanomaterials (publisher Department of Physics, Technical
  University of Denmark, year 2017)NoStop
[Zhu et al.(2019)Zhu,
  Guo, Zou, Li, Yuen, Mihaylova, and Leung]zhu2019review
  author author H. Zhu, author B. Guo, author K. Zou, author
  Y. Li, author K.-V. Yuen, author L. Mihaylova,  and author H. Leung, title title
  A review of point set registration: From pairwise registration to groupwise
  registration, @noop  journal journal
  Sensors volume 19, pages 1191
  (year 2019)NoStop
[Besl and McKay(1992)]besl_1992
  author author P. J. Besl and author N. D. McKay, title title Method for
  registration of 3-D shapes, in 10.1117/12.57955 booktitle Sensor Fusion IV: Control Paradigms and Data
  Structures, Vol. volume 1611, editor edited
  by editor P. S. Schenker,  organization International Society for Optics and Photonics  (publisher SPIE, year 1992) pp. pages 586 – 606NoStop
[Pan et al.(2018)Pan,
  Yang, Liang, and Dong]pan2018iterative
  author author Y. Pan, author B. Yang, author F. Liang,  and author
  Z. Dong, title title Iterative global similarity points: A robust
  coarse-to-fine integration solution for pairwise 3d point cloud
  registration, in @noop  booktitle 2018
  International Conference on 3D Vision (3DV) (organization
  IEEE, year 2018) pp. pages
  180–189NoStop
[Bhandarkar et al.(2004)Bhandarkar, Chowdhury, Tang, Yu, and Tolllner]bhandarkar2004surface
  author author S. M. Bhandarkar, author A. S. Chowdhury, author Y. Tang,
  author J. Yu,  and author E. Tolllner, title
  title Surface matching algorithms computer aided
  reconstructive plastic surgery, in @noop  booktitle 2004 2nd IEEE International Symposium on Biomedical Imaging:
  Nano to Macro (IEEE Cat No. 04EX821) (organization IEEE,  year 2004) pp. pages 740–743NoStop
[Maiseli, Gu, and Gao(2017)]maiseli2017recent
  author author B. Maiseli, author Y. Gu,  and  author H. Gao, title title Recent developments and trends in point
  set registration methods, @noop  journal journal Journal of Visual Communication and Image Representation  volume 46, pages 95–106 (year 2017)NoStop
[Kuhn(1955)]kuhn1955
  author author H. W. Kuhn, title title The hungarian
  method for the assignment problem, 10.1002/nav.3800020109 journal journal Naval
  Research Logistics Quarterly volume 2, pages 83–97 (year 1955)NoStop
[Burgers(1934)]burgers_1934
  author author W. Burgers, title title On the process
  of transition of the cubic-body-centered modification into the
  hexagonal-close-packed modification of zirconium, @noop 
  journal journal Physica volume 1, pages 561–586 (year
  1934)NoStop
[Masuda-Jindo, Nishitani, and  Van Hung(2004)]masuda_2004
  author author K. Masuda-Jindo, author S. Nishitani,  and author V. Van Hung, title title hcp-bcc
  structural phase transformation of titanium: Analytic model calculations,  @noop  journal journal Physical Review
  B volume 70, pages 184122 (year 2004)NoStop
[Khaliullin et al.(2011)Khaliullin, Eshet, Kühne, Behler, and Parrinello]khaliullin_2011
  author author R. Z. Khaliullin, author H. Eshet,
  author T. D. Kühne, author J. Behler,  and author
  M. Parrinello, title
  title Nucleation mechanism for the direct
  graphite-to-diamond phase transition, @noop  journal
  journal Nature materials volume 10,  pages 693 (year 2011)NoStop
[Xiao and Henkelman(2012)]xiao_2012
  author author P. Xiao and author G. Henkelman, title title
  Communication: From graphite to diamond: Reaction pathways of the phase
  transition, 10.1063/1.4752249 journal
  journal The Journal of Chemical Physics volume 137, pages 101101 (year
  2012)NoStop
[Bain and Dunkirk(1924)]bain_1924
  author author E. C. Bain and author N. Dunkirk,  title title The nature of martensite,  @noop  journal journal trans. AIME  volume 70, pages 25–47 (year 1924)NoStop
[Nishiyama(2012)]nishiyama_2012
  author author Z. Nishiyama, @noop  title Martensitic
  transformation (publisher Elsevier, year
  2012)NoStop
[Watanabe, Tokonami, and  Morimoto(1977)]watanabe_1977
  author author M. Watanabe, author M. Tokonami,
   and author N. Morimoto,  title title The transition mechanism
  between the cscl-type and nacl-type structures in cscl, @noop 
  journal journal Acta Crystallographica Section A:
  Crystal Physics, Diffraction, Theoretical and General Crystallography  volume 33, pages 294–298 (year 1977)NoStop
[Blanco et al.(2000)Blanco,
  Recio, Costales, and Pandey]blanco_2000
  author author M. A. Blanco, author J. M. Recio,
  author A. Costales,  and  author R. Pandey, title title Transition path for the
  b3b1 phase transformation in
  semiconductors, 10.1103/PhysRevB.62.R10599 journal journal Phys. Rev. B volume
  62, pages R10599–R10602 (year
  2000)NoStop
[Sowa(2001)]sowa_2001
  author author H. Sowa, title title On the transition
  from the wurtzite to the nacl type, @noop  journal
  journal Acta Crystallographica Section A: Foundations of
  Crystallography volume 57, pages
  176–182 (year 2001)NoStop
[Zhang and Kelly(2005)]zhang_kelly_2005
  author author M.-X. Zhang and author P. Kelly,  title title Edge-to-edge matching model
  for predicting orientation relationships and habit planes—the
  improvements, https://doi.org/10.1016/j.scriptamat.2005.01.040 journal
  journal Scripta Materialia volume
  52, pages 963 – 968 (year 2005)NoStop
[Bowles(1951)]bowles1951crystallographic
  author author J. Bowles, title title The
  crystallographic mechanism of the martensite reaction in iron-carbon
  alloys, @noop  journal journal Acta
  Crystallographica volume 4, pages
  162–171 (year 1951)NoStop
[Bowles and Mackenzie(1954)]bowles1954crystallography
  author author J. Bowles and author J. Mackenzie, title title The
  crystallography of martensite transformations i, @noop  journal journal Acta metallurgica volume 2, pages 129–137 (year
  1954)NoStop
[Mackenzie and Bowles(1954)]mackenzie1954crystallography
  author author J. Mackenzie and author J. Bowles, title title The
  crystallography of martensite transformations ii, @noop 
  journal journal Acta Metallurgica volume 2, pages 138–147 (year
  1954)NoStop
[Wechsler, Lieberman, and  TA(1953)]Wechsler1953
  author author M. Wechsler, author D. Lieberman,
   and author R. TA, @noop
   journal journal Transactions of the American
  Institue of Mining and Metallurgical Engineers volume
  197, pages 1503 (year 1953)NoStop
[Wechsler, TA, and Lieberman(1960)]Wechsler1960
  author author M. Wechsler, author R. TA,  and  author D. Lieberman, @noop
   journal journal Transactions of the American
  Institue of Mining and Metallurgical Engineers volume
  218, pages 202 (year 1960)NoStop
[Therrien and Stevanović(2020)]therrien_2020
  author author F. Therrien and author V. Stevanović, @noop  title Unifying
  description of the martensitic phase transformation from the minimization of
  atomic displacements,  (year 2020), http://arxiv.org/abs/1912.11915 arXiv:1912.11915 [cond-mat.mtrl-sci]
  NoStop
[Date and Nagi(2016)]date2016gpu
  author author K. Date and author R. Nagi,  title title Gpu-accelerated hungarian
  algorithms for the linear assignment problem, @noop  journal journal Parallel Computing volume 57, pages 52–72 (year
  2016)NoStop
[Nahor and Kaplan(2016)]nahor_2016
  author author H. Nahor and author W. D. Kaplan, title title Structure of the
  equilibrated ni(111)-ysz(111) solid–solid interface, 10.1111/jace.14015 journal journal Journal of
  the American Ceramic Society volume 99, pages 1064–1070 (year 2016)NoStop
[Li et al.(2016)Li,
  Chen, Zang, and Feng]Li_2016
  author author L. Li, author Z. Chen, author Y. Zang,  and author
  S. Feng, title title Atomic-scale characterization of si(110)/6h-sic(0001)
  heterostructure by hrtem, https://doi.org/10.1016/j.matlet.2015.10.017 journal journal Materials Letters volume 163, pages 47 – 50 (year 2016)NoStop
[Note3()]Note3
  note Our algorithm does not differentiate between OR_1 and OR_2 because they are identical
  in-plane.Stop
[Ding et al.(2006)Ding,
  Zhou, He, and Zha]ding2006r
  author author C. Ding, author D. Zhou, author X. He,  and author
  H. Zha, title title R 1-pca: rotational invariant l 1-norm principal component
  analysis for robust subspace factorization, in @noop  booktitle Proceedings of the 23rd international conference on
  Machine learning (organization ACM, year
  2006) pp. pages 281–288NoStop
[Nie et al.(2010)Nie,
  Huang, Cai, and Ding]nie2010efficient
  author author F. Nie, author H. Huang, author X. Cai,  and author
  C. H. Ding, title title Efficient and robust feature selection via joint
  𝓁2,1-norms minimization, in http://papers.nips.cc/paper/3988-efficient-and-robust-feature-selection-via-joint-l21-norms-minimization.pdf
  booktitle Advances in Neural Information Processing
  Systems 23, editor edited by editor
  J. D. Lafferty, editor
  C. K. I. Williams, editor
  J. Shawe-Taylor, editor
  R. S. Zemel,  and editor
  A. Culotta (publisher Curran
  Associates, Inc., year 2010) pp. pages
  1813–1821NoStop



