Predicting kinetics of polymorphic transformations from structure mapping and coordination analysis
To extend rational materials design and discovery into the space of metastable polymorphs, rapid and reliable assessment of their lifetimes is essential. Motivated by the early work of Buerger (1951), here we investigate the routes to predict kinetics of polymorphic transformations using solely crystallographic arguments. As part of this investigation we developed a general algorithm to map crystal structures onto each other and construct transformation pathways between them. The developed algorithm reproduces reliably known transformation pathways and reveals the critical role that the dissociation of chemical bonds along the pathway plays in choosing the best (low-energy barrier) mapping. By utilizing our structure-mapping algorithm we were able to quantitatively demonstrate the intuitive expectation that the minimal dissociation of chemical bonds along the pathway, or in ionic systems, the condition of coordination of atoms along the pathway not decreasing below the coordination in the end compounds, represents the requirement for fast transformation kinetics. We also demonstrate the effectiveness of the structure-mapping algorithm in combination with the coordination analysis in studying transformations between polymorphs for which a detailed atomic-level picture is presently elusive.
Strong dependence of physical properties on crystal structure offers compelling opportunities for discovering new functional materials among metastable polymorphs. The search, however, critically requires: (i) accurate and efficient structure and property prediction algorithms to screen the potential energy surface (PES) for polymorphs with attractive properties, and (ii) rapid and reliable assessment of the kinetics of transformations to lower energy states to help identify long-lived polymorphs among the multitude of possible low-energy structures. The former problem has attracted the attention of the scientific community due to what was initially thought to be insurmountable challenges associated with predicting the structure of periodic systems. As a result, a range of structure prediction methods have been developed Woodley and Catlow (2008); Oganov (2010); Evrenk (2014) and successfully combined with property predictions Hautier et al. (2011); Lyakhov and Oganov (2011); Botti et al. (2012); Li et al. (2014); Peng et al. (2015); Jain et al. (2016); Needs and Pickard (2016); Kvashnin et al. (2017). In the context of metastability and polymorphism, identification of bonafide metastable structures among tens if not hundreds of low-energy structures that typically result from structure predictionsZwijnenburg et al. (2010); Zwijnenburg and Bromley (2011) still remains a challengeSun et al. (2016). The recent discovery of a correlation between volumes of configuration space occupied by different PES local minima (their basins of attraction) and experimentally realized metastable structures De et al. (2014); Stevanović (2016); Jones and Stevanović (2017) offers some help in narrowing down the list of candidate metastable structures. However, when targeting metastable forms of matter, knowledge of the kinetics of transformations to lower energy states remains invaluable.
Efforts in developing approaches to predict the minimal energy transformation pathways and associated structures and energies of transitions states, to assess the kinetics of polymorphic transformations have also been (and still are) under development. The most common is probably the generalized solid-state nudged elastic band (ssNEB) formalism Caspersen and Carter (2005); Sheppard et al. (2012); Qian et al. (2013), a technique for locating saddle-points on the potential energy surface to determine the lowest-energy transformation pathways (trajectories of atoms) between different crystal structures, and to compute the corresponding energy barriers. The ssNEB is a generalization of the original NEB method Henkelman and Jónsson (2000); Henkelman et al. (2000) developed for non-periodic systems. The main challenge in applying the ssNEB method is choosing the best initial pathway that will then be modified by the algorithm to the closest minimum energy path.
Despite the impact that these developments have had on determining mechanisms and predicting kinetics, they are not suitable for the screening of large number of structures for metastable polymorphs. A method that could examine large sets of polymorphs, even if offering only a qualitative assessment of the kinetics of polymorphic transformations, would greatly benefit metastable polymorph discovery. One such qualitative assessment can be found in the work of Buerger Buerger (1951), who classified polymorphic transformations as slow or rapid using mainly crystallographic arguments. In short, the existence of a diffusionless (displacive or dilatational) transformation between two structures is, according to Buerger, a signature of fast kinetics. One transformation that disobeys this classification is diamond to graphite for which a relatively simple, diffusionless pathway can be constructed Khaliullin et al. (2011); Xiao and Henkelman (2012), but diamond is a known long-lived allotrope of carbon at ambient conditions. To account for this, Buerger devised an additional class of transformations he called bond type, covalent to metallic in this case, which he argued should be slow irrespective of geometry.
Motivated by the Buerger’s work, we investigate in this paper the routes to generalize his classification and convert it from a collection of heuristics into a computational approach with the goal of accelerating the assessment of the kinetics of polymorphic transformations. The key component of our work is an algorithm to map crystal structures onto each other and construct diffusionless transformations between them (see Fig. 1). The mapping algorithm relies on two basic principles: (a) minimization of the total Euclidian distance atoms need to travel between the end structures, and (b) minimization of the change in coordination of atoms in the first coordination shell (number of atoms in the first shell) along the map. In other words, the goal is to find an optimal mapping that is diffusionless in nature and which minimizes the dissociation of chemical bonds.
We will show that the developed algorithm successfully reproduces known transformation pathways in well studied systems including to , diamond to graphite, CsCl-type to rocksalt, etc. It reveals the critical role of minimizing the change in coordination in searching for the best mapping. Beyond simple and well studied systems, we applied our mapping algorithm to polymorphic transformations between various SnO high-pressure phases. Our results show that if only a qualitative assessment of the kinetics of polymorphic transformations is needed, the condition of minimal dissociation of chemical bonds along the pathway, can be used as a signature of fast kinetics. Note that, in partially ionic systems, the condition of minimal dissociation, can be formulated as the coordination of atoms in the first shell not decreasing below the coordination in the end structures. Lastly, if a more accurate estimate of the transition pathway and the associated energy profile is necessary, our mapping algorithm is demonstrated to provide a good starting point for subsequent solid-state nudged elastic band (ssNEB) calculations.
Similarly to other methods Capillas et al. (2007); Caspersen and Carter (2005); Sheppard et al. (2012); Qian et al. (2013), in our work we also draw conclusions about the kinetics of polymorphic transformations from the assumed collective (concerted) motion of atoms. While phase transitions in solids occur mainly via nucleation and growth, we adopt the view that the features of collective transformations can provide useful insights into the overall kinetics and can serve as the starting point for more realistic modeling of nucleation processes.
Ii Structure mapping algorithm
The developed structure mapping algorithm consists of two steps as illustrated in Fig. 1. In the first step the algorithm searches for the most compatible representations of the two structures, that is the optimal mapping of their unit cells. This is done in the following way. The two primitive unit cells are brought to the same – least common multiple – number of atoms by constructing all symmetry inequivalent supercells. The enumeration of the symmetry inequivalent supercells given the number of atoms is done using Hart-Forcade theory Hart and Forcade (2008). The most compatible representations are then defined by the pair of supercells that minimizes the strain between them, or, in other words, the two supercells with the largest spatial (volumetric) overlap. This condition can be formulated as searching for the pair of unit cells that minimizes the metric defined as the weighted sum of the absolute differences in cell parameters and the total surface areas (S) of the two cells:
where represent positive weights of quantities that are introduced solely to make the numerical values the same order of magnitude. The search is accomplished by transforming the cells to the corresponding (unique) reduced cell according to the formulation of Niggli-Santora-Gruber Niggli (1928); Santoro and Mighell (1970); Krivy and Gruber (1976), which allows implicitly exploring all permutations of the unit vectors and all isometric transformations of the two cells (rigid rotations and reflections) to bring them to the positions of largest spatial overlap, as shown in Fig. 1.
In the second step, the atoms are placed back inside the two cells. The optimal atom-to-atom mapping is then searched for by performing the following operations on the two sets of atomic positions: (i) all symmetry operations of the parent Bravais lattices, (ii) translations of the origin of coordinate frame to each atom site, and (iii) permutations of indices of chemically identical atoms. Every choice of symmetry operation, position of the origin, and the permutation of atom indices defines a single mapping between the two end structures and a pathway in configuration space that connects the atoms with same indices and continuously deforms the unit cell.
Out of many possible atom-to-atom mappings our algorithm selects as the optimal mapping the solution that yields minimal dissociation of chemical bonds, i.e., the coordination of atoms along the pathway does not decrease below the coordinations in the end structures. If such a mapping cannot be found, then the one that minimizes the sum of Euclidian distances between the corresponding atoms in the two structures is chosen. Finally, if more than one solution is found with the appropriate change in coordination then the sum of distances between the atoms is used to rank them and narrow down the choice. To achieve this, the described operations are performed using Hart-Forcade theory Hart and Forcade (2008) that allows enumeration of the symmetry inequivalent atom sites. The Hungarian Algorithm of Kuhn and Munkres Kuhn (1955) is used to find the optimal permutation of atom indices (job scheduling problem) in polynomial time.
Ideally, the globally best atom-to-atom mapping would always be found in the unit cell pairing with highest overlap, but this is not the case in general. Therefore, in the first step, we form a list of the top unit cell pairing (e.g. =10) and perform the second step on all of them. In practice this has been found to yield the desired globally optimal atom-to-atom mapping.
A similar approach was developed previously by Sadeghi and Goedecker Sadeghi et al. (2013) for the purpose of measuring configuration space distances between non-periodic systems. Also, Lonie and Zurek Lonie and Zurek (2012) developed a search algorithm designed to identify identical (duplicate) periodic structures which involves mapping of the unit cells. Another class of recently developed approaches utilize the descriptor/feature based fingerprinting to quantify similarity of different periodic structures by comparing selected set of features (not atom-by-atom). These include the work of Yang et al. Yang et al. (2014) and Zhu et al. Zhu et al. (2016). Our mapping of the unit cells can be viewed as a generalization of the ideas of Lonie and Zurek Lonie and Zurek (2012) to the case where input structures are presumed to be different and where the goal is to discover the optimal alignment of the two structures. Concerning the atom-to-atom mapping, we extended the algorithm of Sadeghi and Goedecker Sadeghi et al. (2013) to periodic systems. Another important distinguishing feature of our approach is the new objective function, which includes the physical principle of the minimal dissociation of chemical bonds.
We tested our structure mapping algorithm on a number of well-studied polymorphic transformations in elemental and simple binary systems, some of which are listed in Table 1. The space groups of the lowest symmetry structures that occur along the pathway are also shown. The transformation pathways from Table 1 reproduce well the maximal symmetry transition paths for the reconstructive phase transitions derived by Capillas et al. Capillas et al. (2007). Further, our algorithm often finds the pathway that preserves the largest common subgroup of the two structures. This is a consequence of the supercell sizes chosen to accommodate the least common multiple number of atoms. If this condition is relaxed and larger supercells with compatible numbers of atoms are considered, our algorithm would also find lower symmetry pathways, some of which have been discussed previously for the transformations from Table 1. The details of the transformations from Table 1 are provided in the supplementary materials.
|FCC (Fmm)||I4/mmm||BCC (Imm)|
|BCC (Imm)||Cmcm||HCP (P6/mmc)|
|diamond (Fdm)||C2/m||graphite (P6/mmc)|
|CsCl-type (Pmm)||Rm||rocksalt (Fmm)|
|rocksalt (Fmm)||Cmc2||wurtzite (P6mc)|
|zincblende (F3m)||Imm2||rocksalt (Fmm)|
Iii Change in the coordination of atoms and its relevance to the magnitude of kinetic barriers
Another important finding that emerged during the development of our algorithm is that the two conditions, minimal distance between the corresponding atoms in the two structures and minimal change in coordination, do not always coincide. In this section we analyze relatively simple transformations for which the coordination of atoms, i.e, the number of atoms in the first coordination shell, is relatively straightforward to define. A more through discussion on the definition of the first shell coordination is presented in Section V.
The test cases from Table 1 revealed the critical significance of the change in coordination. For example, minimizing the distance alone does not necessarily result in the best pathway in the case of to transformation as there are multiple possible mappings with nearly degenerate distance values, but different coordination of atoms along the pathway. The transformations between the ground state and high temperature tin-sulfide polymorphs, Pnma and Cmcm, also exhibit similar features (see supplementary materials). The condition itself has been referred to previously in case-by-case studies (e.g. Refs. Catti (2001); Miao and Lambrecht (2003)) and has a relatively simple physical interpretation. It reflects the requirement of minimizing the number of broken bonds in going from one end structure to the other, which is expected to be correlated with the magnitude of the kinetic barriers. In ionic and covalent systems, because of the directionality of chemical bonds, this condition can be implemented just by requiring that the first-shell coordination of atoms along the pathway does not decrease below the coordination in the end structures.
We illustrate the correlation between the change in the first shell coordination and the potential energy profile by analyzing the rocksalt (Fmm) wurtzite (P6mc) transformation in binary ionic systems GaN, ZnO, MgO, and CdO. In the upper part of Fig. 2 five snapshots along the pathway are shown viewed from the wurtzite -axis. In the middle part, potential energy profiles along the pathway are shown. They are calculated using density functional theory (DFT). The energy axis is relative to each ground state structure, wurtzite for GaN and ZnO and rocksalt for MgO and CdO. For every chemical composition the end structures are fully relaxed (volume, cell shape and atomic positions), while the total energies of the snapshots along the pathway are computed without any relaxations. For each system the pathways are discretized into 100 equally spaced successive snapshots shown in Fig. 2. We employ a relatively standard DFT-GGA numerical setup described in Section VII.
As evident from Fig. 2 for all four systems calculated energy profiles exhibit relatively small barriers for transformation from the higher energy into the corresponding ground state state structure, which corresponds well to the known, fast kinetics of this transformation Decremps et al. (2002); Limpijumnong and Lambrecht (2001). The highest energy points along the pathways are all below 100 meV/atom, which is usually considered low activation energy Trinkle et al. (2003). Because the calculated energy profiles can be thought of as providing the upper bounds for the true activation energies, this result shows that irrespective of the chemistry, this particular transformation is associated with relatively low energy barriers. In addition, the change in coordination of the atoms from the 6-fold coordinated rocksalt phase to 4-fold coordinated is evidently monotonic, i.e. only the absolutely necessary number of bonds dissociate along the pathway, which is consistent with the above discussion. Similar results are obtained for other transformations from Table 1 including: BCC HCP in metallic titanium for which we find a nearly barrierless energy profile in line with previous findings Trinkle et al. (2003); CsCl-type rocksalt in CsCl, with the monotonic change in coordinations and the upper bound for the activation energy of 13 meV/atom. Moreover, for the zincblende rocksalt transformation in SiC our algorithm confirms previously discussed energetic preference of the orthorhombic (Imm2) pathway relative to the R3m mechanism because of the monotonic change in the coordination of atoms from four to sixCatti (2001); Miao and Lambrecht (2003). The calculated energy profiles peak at 112 and 250 meV for the pathways passing through the Imm2 and R3m, respectively, in agreement with previous studies.
Iv Diamond graphite transformation in elemental carbon
Our results also confirm the exception to the rule, the cubic diamond hexagonal graphite transformation for which our algorithm finds the same transformation pathway as the one discussed in Refs. Khaliullin et al. (2011); Xiao and Henkelman (2012) that is accompanied by a monotonic change in coordination of atoms (see Fig. 3) and yet, the maximal energy point along the pathway is 470 meV/atom above the diamond structure. Similar to other systems studied herein, the calculations were performed by discretizing the pathway constructed by our algorithm (interpolation). Fig. 3 also illustrates the energy profile in grey, which resulted from the ssNEB calculation that used our mapping algorithm pathway as the initial pathway. The energy profiles are calculated using the optB86 exchange-correlation functional that includes contributions coming from Van der Waalls interactions KlimeÅ¡ et al. (2010) as implemented in the VASP software package Kresse and Furthmüller (1996). As expected, the calculated ssNEB energy barrier is lower, but the qualitative picture does not change. The transformation is slow despite the monotonic change in coordination reproduced by both sets of calculations.
So, why is elemental carbon different? As correctly observed by Buerger, the diamond to graphite transformation requires a change in chemical bonding. In the diamond structure the four carbon hybrid atomic orbitals form four strong bonds per carbon atom. These are replaced in the graphite structure by three hybrids, which form three bonds per carbon, and a orbital, which forms one weaker, resonating bond per carbon. Relatively weak Van der Waals interactions between the layers further stabilize the graphite structure. If one counts the number of chemical bonds per atom rather than the geometric coordination, then both diamond and graphite have four bonds per carbon atom and the transformation from diamond to graphite would imply dissociation and formation of one bond per C. Consistent with the previous discussion, the intermediate decrease in the number of chemical bonds to three from four in the end structures would lead to a high barrier.
In ionic systems on the other hand, significant contribution to the energy differences between different atomic configurations comes form purely electrostatic interactions Stevanović et al. (2010, 2011). Hence, the increase in energy along the pathway is influenced by the changes in the charge distribution and to a lesser extent is due to vanishing overlaps of atomic overlaps. The evidence of the remaining charge transfer along the pathways are the non-vanishing band-gaps for all ionic systems studied here. Therefore, the argument here is that the geometric coordination of atoms along the pathway is more appropriate when trying to understand the kinetics of polymorphic transformations in ionic systems, while in covalent systems it needs to be replaced by a chemical bonding analysis. In both cases however, the condition of minimal dissociation of chemical bonds serves as a signature of rapid transformations.
V SnO polymorphs
To further validate the previous discussion we extend our study to SnO, a partially ionic system for which a number of polymorphs have been realized under pressure. With increasing pressure the structures appear in the following sequence: P4/mnm Pnnm Pbcn Pa3 Pbca Fmm Shieh et al. (2006); Das and Jayaraman (2014). Upon releasing the pressure, however, all phases either relax back to the ground state rutile (P4/mnm) structure following the same sequence or to a phase mixture between rutile and the Pbcn structure (-PbO structure type) Das and Jayaraman (2014); Liu (1978). So, the only phase that survives at ambient conditions is P4/mnm, occasionally in combination with small amounts of Pbcn. A previous study Stevanović (2016) has shown that these two SnO structures, P4/mnm and Pbcn, have the “largest” local minima, that is, they occupy larger regions of configuration space than any other. Here, we extend this result by investigating transition pathways between different SnO polymorphs.
In Fig. 4, a chart illustrating the crystal structures of all six SnO polymorphs is shown with thick arrows connecting structures for which
our structure mapping and coordination analysis suggest fast polymorphic transformations. The arrows point in the direction of lowering
total energy. Interestingly, the highest pressure Fmm phase is connected to all other phases by a fast polymorphic transformation. The DFT
calculated energy profiles all exhibit energy barriers lower than 30 meV/atom.
We also find another classe of transformations between the SnO polymorphs, those marked in Fig. 4 by dashed arrows. For these transformations our classification based on the coordination analysis is very sensitive to how the coordination of atoms is computed. Namely, along the pathways the symmetry of the end structures is usually broken, and the distances between the atoms and their first neighbors are not all the same. As a result, the calculated number of neighbors will in some cases depend on the cutoff radii and/or tolerances on these distances. For binary ionic systems we define the first coordination shell as consisting of the atoms of the type other than the central atom that are all separated from the central atom within a certain tolerance factor. In other words, the tolerance factor represents the thickness of a spherical shell around the central atom. The first shell is then defined as the number of atoms of the other type that are found inside this shell. The counting starts from the closest atom and stops if either the tolerance or the atom of the same type as the central one is reached, whichever comes first. We check the stability of all our results by varying the tolerance factor from 0.1 to 0.5 Å.
The dashed arrows in Fig. 4 represent the transformation pathways that would either be classified as slow for the smallest tolerance factor (0.1 Å) and would change to fast upon increasing the tolerance, or for which the results would be ambiguous, i.e., the trends in the coordination of different types of atoms would be different (the opposite). The former physically means that along the pathway some chemical bonds are strained more than others, but the atoms remain relatively close to each other. The latter means that the atoms of the same type as the central atom approached closer to the central atom and intermixed with its first coordination shell and in that way decreased the evaluated coordination number. Based on the subsequent ssNEB calculations we argue that it is more appropriate to classify these types of transformations as likely rapid for the following reason. Namely, our procedure is based solely on geometry and does not allow atoms to relax to more stable configurations as, for example, the ssNEB method would. Therefore, because of the energy minimization with respect to atomic positions within the ssNEB the atoms will have the opportunity to re-bond during the atomic relaxations, which would then lower the barrier and imply fast transformation based on the coordination analysis. Based on this discussion, all dashed-line transformations from Fig. 4 should be classified as rapid.
Finally, the coordination analysis along the Pbca P4/mnm and Pbca Pnnm pathways clearly shows the dissociation of chemical bonds and consequently, we do not connect these structures with arrows. Unlike those marked by the dashed arrows, the change in coordination along Pbca P4/mnm and Pbca Pnnm does not depend on the details of how is the first shell coordination of atoms evaluated.
To better illustrate previous discussion we show in Fig. 5 calculated energy profiles and the average coordination of atoms along the pathways together with the ssNEB results for a selected set of transformations from Fig. 4. We consider Fmm P4/mnm (thick arrow), Pbcn P4/mnm (dashed arrow), Pa3 P4/mnm (dashed arrow), and Pbca P4/mnm (no arrow). The first transformation, which would be classified as fast according to our coordination analysis, clearly indicates the existence of a low energy barrier. The highest energy point along the pathway is only 22 meV/atom above the high-pressure Fmm phase. The transformation is nearly barrierless in the ssNEB showing a very similar change in the coordination of atoms. The second and third transformations that are denoted by dashed arrows in Fig. 4 both fall in the category of undetermined based on the initial coordination analysis. Namely, the coordination numbers for Sn and O follow opposite trends. While the coordination of Sn grows along the pathway from six to seven the average number of first shell Sn atoms surrounding oxygen drops from 3 to 2.5. The reason for this is already mentioned intermixing of O atoms within the first shell of other O atoms, which decreases their first shell coordination number (as defined here) and increases Sn coordination numbers. Minimization of energy within the ssNEB formalism would allow for some rearrangements of atoms and lowering of the energy barriers along these pathways. Indeed, the energy barriers calculated from the direct interpolation along our pathways are about 216 and 124 meV/atom for Pbcn P4/mnm and Pa3 P4/mnm, respectively. The subsequent ssNEB calculations lower the energy barriers in both cases, but while for the Pbcn P4/mnm the symmetry of the initial pathway remains the same with the final barrier of 115 meV/atom, the ssNEB result for the Pa3 P4/mnm relaxes from the initial pathway to the one that has the Pbcn structure as the intermediate. The ssNEB energy barrier for the Pa3 Pbcn transition is 32 meV/atom and is in qualitative agreement with the direct interpolation one of 97 meV/atom. Given that the Pa3 structure is the next in the sequence of increasing pressure after the Pbcn, the ssNEB result from Fig. 5 explains the appearance of the Pbcn phase mixed with the ground-state rutile phase (P4/mnm) upon releasing pressure. Namely, the barrier of 32 meV/atom implies Pa3 will transform to Pbcn fairly quickly, while the Pbcn phase will transform to P4/mnm at a slower rate. Interestingly, just from the coordination analysis and the calculated energy profiles one could conceive the Pa3 Pbcn P4/mnm route instead of the direct Pa3 P4/mnm.
Similarly, the ssNEB result for the Pbca P4/mnm transformation (see Fig. 5(d)) is clearly consistent with the conclusions that could be drawn form the coordination analysis. Namely, from the coordination of atoms this transformation would undoubtedly be considered slow, and the lowering of energy from Pbca would likely proceed along the pathways marked by the arrows in Fig. 4, that is, either through the Pa3 and then Pbcn, or directly to the Pbcn structure. As previously discussed both of these transformations can be considered rapid. This is exactly what the ssNEB predicts would happen. The Pbca would, according to ssNEB, transform rapidly into the Pbcn (barrier of 44 meV/atom) and would then continue along the Pbcn P4/mnm path.
As the SnO results demonstrate, structure mapping and careful coordination analysis can offer qualitative guidance and accelarate the classification of the polymorphic transformation into rapid and slow. While replacing the ssNEB results with just the energy profiles calculated from the interpolation of our pathways might be tempting, one needs to remember that the barriers calculated in this way only represent the upper bounds for the true activation energy. The problem occurs if the upper bound for the activation energy is relatively large. We argue here that under these circumstances careful coordination analysis can still be useful in providing qualitative assessment as discussed for the transitions denoted by dashed lines in Fig. 4. Of course, the ssNEB in these cases would provide the ultimate answer.
We further tested our mapping algorithm and the qualitative conclusions based on coordination analysis using SiO as another case study (see Supplementary Materials for details). The resulting classification based on the coordination of atoms shows fast transformations between - and -quartz, - and -cristobalite, and - and -tridymite. All other transformations between them and also including high-pressure moganite and stishovite would be classified as slow (coordination does not depend on the tolerance factor). This entirely geometric- and coordination of atoms- based classification reproduces qualitatively well the available knowledge of the SiO polymorphs and the kinetics of polymorphic transformations between them.
In conclusion, with the main goal of accelerating the assessment of the kinetics of polymorphic transformations we developed a general algorithm to map crystal structures onto each other and construct transformation pathways between them. The algorithm itself is based on the physical principles of minimizing the distance atoms need travel between the two structures (diffusionless transformations) and minimizing the change in the first shell coordination of atoms along the pathway (minimal dissociations of chemical bonds). The developed algorithm reliably reproduces known transformation pathways and reveals the critical role the change in coordination of atoms plays in identifying the correct mapping. Moreover, we find that the change in coordination along the pathway can also be used as a proxy for the kinetics of polymorphic transformations in partially ionic systems. Namely, the condition of the first-shell coordination of atoms not decreasing below the end structures, i.e., minimal dissociation of chemical bonds along the pathway, can be used as a signature of rapid transformations. In covalent systems the geometric coordination of atoms needs to be replaced by the proper bonding analysis. This allows qualitative, quick and reliable classification of polymorphic transformations into rapid and slow just from the coordination analysis. For more quantitative assessments integration of our algorithm with ssNEB calculations is shown to provide a robust computational procedure for predicting kinetics of polymorphic transformations. Ultimately, the methods presented here, in combination with structure and property predictions, offer a route to identifying novel, realizable and long-lived functional materials among metastable polymorphs.
All calculations were performed by employing relatively standard DFT computations explained in details elsewhere Stevanović et al. (2012). In short, the PBE form of the exchange-correlation functional Perdew et al. (1996) was used within the projector augmented wave (PAW) method Blöchl (1994) as implemented in the VASP code Kresse and Furthmüller (1996). In case of elemental carbon, the VASP code implementation of the optB86 exchange-correlation functional KlimeÅ¡ et al. (2010), which includes contributions arising from van der Waalls interactions, was employed. To generate energy profiles along the mapping (pathway) all initial and final structures were fully relaxed (volume, cell shape and atomic coordinates), whereas only the static calculations were performed for the intermediate structures along the pathway.
The solid state nudged elastic band (ssNEB) calculations were performed using the implementation in the Transition State Tools for VASP (VTST) code developed by the Henkelman group at UT Austin.Sheppard et al. (2012) Each ssNEB calculation was initialized with an approximately 20 image band determined by taking an interpolation of the entire desired polymorphic transformation generated by the mapping algorithm discussed in this paper. The initial pathway is relaxed until the forces on all images were less than 0.01 eV/Å. If any intermediates appeared along the ssNEB, these were relaxed using a standard geometry relaxation, and the ssNEB path was split, taking the intermediate as a new end point for one ssNEB, and a starting point for another. This was done so that each ssNEB calculation would only have a single local maximum. Images were either added or removed from the chain so that each intermediate ssNEB had 10-12 images. Once the NEB images were relaxed, the climbing image NEB method was used to better identify the true saddle point.Henkelman et al. (2000)
Acknowledgements.This work was supported as part of 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. The research was performed using computational resources sponsored by the Department of Energy’s Office of Energy Efficiency and Renewable Energy and located at the National Renewable Energy Laboratory.
- Complete information about all SnO polymorphs, their relative energies as well as the energy profiles for the transformations between them are provided in Supplementary Information.
- S. M. Woodley and R. Catlow, Nat Mater 7, 937 (2008).
- A. Oganov, Modern methods of crystal structure prediction (Wiley-VCH John Wiley distributor, Weinheim Chichester, 2010).
- S. Evrenk, Prediction and calculation of crystal structures: methods and applications (Springer, Cham Switzerland, 2014).
- G. Hautier, A. Jain, S. P. Ong, B. Kang, C. Moore, R. Doe, and G. Ceder, Chemistry of Materials 23, 3495 (2011).
- A. O. Lyakhov and A. R. Oganov, Phys. Rev. B 84, 092103 (2011).
- S. Botti, J. A. Flores-Livas, M. Amsler, S. Goedecker, and M. A. L. Marques, Phys. Rev. B 86, 121204 (2012).
- Y. Li, J. Hao, H. Liu, Y. Li, and Y. Ma, J. of Chem. Phys. 140, 174712 (2014).
- H. Peng, P. F. Ndione, D. S. Ginley, A. Zakutayev, and S. Lany, Phys. Rev. X 5, 021016 (2015).
- A. Jain, Y. Shin, and K. A. Persson, Nature Reviews Materials 1, 15004 EP (2016).
- R. J. Needs and C. J. Pickard, APL Materials 4, 053210 (2016).
- A. G. Kvashnin, A. R. Oganov, A. I. Samtsevich, and Z. Allahyari, J. Phys. Chem. Lett. 8, 755 (2017).
- M. A. Zwijnenburg, F. Illas, and S. T. Bromley, Phys. Rev. Lett. 104, 175503 (2010).
- M. A. Zwijnenburg and S. T. Bromley, Phys. Rev. B 83, 024104 (2011).
- W. Sun, S. T. Dacek, S. P. Ong, G. Hautier, A. Jain, W. D. Richards, A. C. Gamst, K. A. Persson, and G. Ceder, Science Advances 2 (2016), 10.1126/sciadv.1600225.
- S. De, B. Schaefer, A. Sadeghi, M. Sicher, D. G. Kanhere, and S. Goedecker, Phys. Rev. Lett. 112, 083401 (2014).
- V. Stevanović, Phys. Rev. Lett. 116, 075503 (2016).
- E. Jones and V. Stevanović, arXiv preprint , arXiv:1708.09026 (2017).
- K. J. Caspersen and E. A. Carter, Proc. Natl. Acad. Sci. U.S.A. 102, 6738 (2005).
- D. Sheppard, P. Xiao, W. Chemelewski, D. D. Johnson, and G. Henkelman, The Journal of Chemical Physics 136, 074103 (2012).
- G.-R. Qian, X. Dong, X.-F. Zhou, Y. Tian, A. R. Oganov, and H.-T. Wang, Comput. Phys. Commun. 184, 2111 (2013).
- G. Henkelman and H. Jónsson, The Journal of Chemical Physics 113, 9978 (2000).
- G. Henkelman, B. P. Uberuaga, and H. Jónsson, The Journal of Chemical Physics 113, 9901 (2000).
- M. J. Buerger, in Phase transformations in solids, edited by R. Smoluchowski, J. E. Mayer, and W. A. Weyl (John Wiley and Sons, New York, 1951) pp. 183–211.
- R. Z. Khaliullin, H. Eshet, T. D. Kühne, J. Behler, and M. Parrinello, Nat. Mater. 10, 693 (2011).
- P. Xiao and G. Henkelman, The Journal of Chemical Physics 137, 101101 (2012).
- C. Capillas, J. M. Perez-Mato, and M. I. Aroyo, J. Phys.: Condens. Matter 19, 275203 (2007).
- G. L. W. Hart and R. W. Forcade, Phys. Rev. B 77, 224115 (2008).
- P. Niggli, Handbuch der Experimentalphysik 7, part 1 (1928).
- A. Santoro and A. D. Mighell, Acta Cryst. A26, 124 (1970).
- I. Krivy and B. Gruber, Acta Cryst. A32, 297 (1976).
- H. W. Kuhn, Naval Research Logistics Quarterly 2, 83 (1955).
- A. Sadeghi, S. A. Ghasemi, B. Schaefer, S. Mohr, M. A. Lill, and S. Goedecker, J. Chem. Phys. 139, 184118 (2013).
- D. C. Lonie and E. Zurek, Comput. Phys. Comm. 183, 690 (2012).
- L. Yang, S. Dacek, and G. Ceder, Phys. Rev. B 90, 054102 (2014).
- L. Zhu, M. Amsler, T. Fuhrer, B. Schaefer, S. Faraji, S. Rostami, S. A. Ghasemi, A. Sadeghi, M. Grauzinyte, C. Wolverton, and S. Goedecker, J. Chem. Phys. 144, 034203 (2016).
- M. Catti, Phys. Rev. Lett. 87, 035504 (2001).
- M. S. Miao and W. R. L. Lambrecht, Phys. Rev. B 68, 092103 (2003).
- F. Decremps, J. Pellicer-Porres, F. Datchi, J. P. Itié, A. Polian, F. Baudelet, and J. Z. Jiang, Applied Physics Letters 81, 4820 (2002).
- S. Limpijumnong and W. R. L. Lambrecht, Phys. Rev. Lett. 86, 91 (2001).
- D. R. Trinkle, R. G. Hennig, S. G. Srinivasan, D. M. Hatch, M. D. Jones, H. T. Stokes, R. C. Albers, and J. W. Wilkins, Phys. Rev. Lett. 91, 025701 (2003).
- J. KlimeÅ¡, D. R. Bowler, and A. Michaelides, Journal of Physics: Condensed Matter 22, 022201 (2010).
- G. Kresse and J. Furthmüller, Comput. Mater. Sci. 6, 15 (1996).
- V. Stevanović, M. d’Avezac, and A. Zunger, Phys. Rev. Lett. 105, 075501 (2010).
- V. Stevanović, M. d’Avezac, and A. Zunger, Journal of the American Chemical Society 133, 11649 (2011).
- S. R. Shieh, A. Kubo, T. S. Duffy, V. B. Prakapenka, and G. Shen, Phys. Rev. B 73, 014105 (2006).
- S. Das and V. Jayaraman, Prog. Mater Sci. 66, 112 (2014).
- L.-G. Liu, Science 199, 422 (1978).
- Complete information about all SnO polymorphs, their relative energies as well as the energy profiles for the transformations between them are provided in Supplementary Information.
- V. Stevanović, S. Lany, X. Zhang, and A. Zunger, Phys. Rev. B 85, 115104 (2012).
- J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996).
- P. E. Blöchl, Phys. Rev. B 50, 17953 (1994).