Identity crisis in alchemical space drives the entropic colloidal glass transition

A universally accepted explanation for why liquids sometimes vitrify rather than crystallize remains hotly pursued, despite the ubiquity of glass in our everyday lives, the utilization of the glass transition in innumerable modern technologies, and nearly a century of theoretical and experimental investigation. Among the most compelling hypothesized mechanisms underlying glass formation is the development in the fluid phase of local structures that somehow prevent crystallization. Here, we explore that mechanism in the case of hard particle glasses by examining the glass transition in an extended alchemical (here, shape) space; that is, a space where particle shape is treated as a thermodynamic variable. We investigate simple systems of hard polyhedra, with no interactions aside from volume exclusion, and show via Monte Carlo simulation that glass formation in these systems arises from a multiplicity of competing local motifs, each of which is prevalent in—and predictable from—nearby ordered structures in alchemical space.


A. Simulation protocols
To verify the self-assembly behavior reported by Klotsa et al. 1 and to generate our own data, we sampled one-component systems composed of particles in the spheric triangle invariant 323 family, as detailed in the main text. Simulations of 4,096 particles were run for about 100 million MC sweeps or until self-assembly was observed. At all state points corresponding to α a < 0.3 and α c > 0.4 for which self-assembly was not observed, we also simulated smaller systems of 2,624 particles for about 70-100 million MC sweeps. At one state point, (α a , α c , φ) = (0, 0.6, 0.6), self-assembly was observed around 120 million MC sweeps. To reach each density, we initialized our system in a sparse cubic array inside a cubic box, randomized the system by running isochoric Monte Carlo sampling for 10,000 MC sweeps, and progressively rescaled box vectors by a scale factor of 0.9995 until the target density was reached. After every rescaling step, isochoric Monte Carlo sampling proceeded until all particle-particle overlaps were eliminated. During the compression process, the translational trial move size was identical to its value during equilibrium sampling, but the rotational trial move size was larger in order to facilitate the removal of overlaps during fast compressions to high densities. During equilibrium sampling, rotational and translational trial move sizes were constant across all systems studied, and were chosen to most efficiently structurally relax a typical system.
To perform the "doping simulations" detailed in the main text, we employed a similar compression and equilibration strategy, with a few differences to accommodate the larger size and aspect ratio of the dimer dopants: we thermalized the system prior to compression for 1 million MC sweeps rather than 10,000 MC sweeps, and switched rotation move size during compression to its smaller equilibration value if compression was proceeding slowly.

B. Environment matching for crystallinity determination
We determined fractions of particles in a particular crystalline environment in our simulations via an "environment matching" scheme, which we have made available as part of freud 2 , our group's open-source analysis package. We define particle i's environment as the set of vectors { r im }, where r im points from the center of particle i to the center of particle m, and m is an index over i's M nearest neighbors. We then determine whether i's environment "matches" the environment of any of its K nearest neighbors: for particle k, where k is a member of the set of the K nearest neighbors of i, we determine k's environment { r km }, where m loops over k's M nearest neighbors. We then attempt to find a one-to-one mapping between the vectors of i's environment and the vectors of k's environment such that | r im − r km | < t for every mapping pair (m, m ) for some threshold t. If that mapping is found, i's environment matches k's environment. We cluster particles according to common environment, searching over the K nearest neighbors of every particle, and label a particle as "crystalline" if it belongs to a cluster of size s > s * , where s * is some cut-off. In almost all cases, we chose t = 0.3r * , where r * is a length scale we chose to be 1.7, a typical nearest-neighbor distance in our systems. We occasionally set t = 0.2r * for finer-grained analysis of systems that were suspected to assemble into the high-pressure Lithium structure. For tracking crystallization during alchemical Monte Carlo runs, we chose the size cut-off s * = 5, while for determining crystallization time, we chose the size cut-off s * = 1. Environment parameters K and M were tailored according to the crystal structure whose growth or presence we wanted to detect: for fcc, we set {K, M } = {12, 12}, for bcc, we set {K, M } = {14, 14}, for diamond, we set {K, M } = {16, 4}, for high-pressure Lithium, we set {K, M } = {80, 11}, and for γ-brass, we set K = 200, and found M for each particle as the number of neighbors that fell within a cut-off distance of r = 1.7. Then, only particles with identical values of M were tested for similarity. This was necessary due to the spread in nearest neighbor number among particles in the more complicated γ-brass structure.

C. Dynamical characterization
To establish the glass-forming ability of disordered systems in this shape space, we chose particle shapes that remain disordered on the long time scales of our simulations at the densities investigated to establish phase behavior, and re-sampled one-component systems of those shapes in a broader density range. We used Monte Carlo methods identical to those described in detail above. We equilibrated our systems for approximately 50 million MC sweeps, and then collected dynamical data for 100 million MC sweeps.
For each system, we calculated the following order parameters at logarithmic timescales: the mean-squared displacement ∆r 2 (t) of all particles in the system, the self-part of the intermediate scattering function F s ( k, t), computed for the k-value associated with the first peak of the static structure factor, the non-Gaussian parameter α(t) 3 , and the self-part of the four-point susceptibility χ SS 4 (t) 4,5 . These quantities are defined as follows: In all definitions above, ∆ r j (t) ≡ r j (t) − r j (0). In the final expression, H is the Heaviside step function: H(x) = 1 for x > 0 and 0 otherwise. a is a length-scale associated with the self-overlap of any particle in the system; in this paper we took a to be the inscribing sphere radius of the particle shape for each system. In all cases, angle brackets indicate ensemble averages. We determined that relaxation in most systems is complete by about 10 million MC sweeps; thus, we broke each 100 million MC sweep trajectory into 10 windows and took appropriate ensemble averages over these windows. Error bars were determined through either error propagation or jackknife resampling; see Supplementary Method I D for more detail.
We computed F s ( k, t) at every pertinent lagtime t by averaging over computed values of F s ( k, t) for 10 randomly generated vectors with magnitude k in a similar manner to that described elsewhere 6 . We did this to speed up our calculations, as F s ( k, t) = F s (k, t) in an isotropic medium. We computed F s ( k, t) for the k value corresponding approximately to the location of the first peak of the static structure factor of each system. Supplementary Fig. 1 shows the static structure factors S(k) for the disordered systems (α a , α c ) = (0, 0.5) and (α a , α c ) = (0.2, 0.5). We found S( k) ≡ 1 N N j,l=1 e i k·( r l − rj ) via the squared FFT of the number density ρ( r) of the system, Gaussian-blurred for smoothness, since S( k) = 1 N d rρ( r)e i k· r 2 . We then found S(k) by assuming that the system is isotropic, and spherically averaging S( k) using a channel-sharing method 7 . The static structure factors we show here were calculated for the first frame of each trajectory only. They are given as functions of kσ, where σ is a length scale that characterizes the particle size of each system: where v p is the particle volume (1 in all cases).
In glass-forming systems generally, ∆r 2 (t) and F s ( k, t) increasingly display three regimes as density increases or temperature decreases: a regime at short timescales in which particles move without colliding with any others, a caging regime at intermediate timescales in which particles are caged by their neighbors and relaxation slows, and a regime at long timescales in which particles escape the confines of their cages and eventually diffuse through the system. α(t) gives a measure of the degree to which the distribution of particle displacements in the system deviates from a Gaussian distribution. It typically has a peak in glassforming systems at times of large dynamical heterogeneity 8,9 , when some particle motions are cooperative and therefore a subset of particle displacements is higher than that given by the expected Gaussian distribution. χ SS 4 (t) gives a direct measure of the dynamical heterogeneity of the system, as it is the scaled variance of the 2-point self-correlation function Q S (t): χ 4 grows from zero as heterogeneity in the dynamics of the system increases over a time window t, and decreases back to zero at long times in the dense fluid.

D. Error determination of dynamical order parameters
We determined error bars for each dynamical order parameter with a variety of techniques. We detail them here.
In practice we calculated ∆r 2 (t) in the following manner: where r jm (t) is the position of particle j at time t in frame m, and there are M total frames making up the ensemble average. Then we computed the variance of this average as where σ 2 m (t) is the variance of the mean squared displacement taken across all particles in frame m. This is error propagation via a Taylor series expansion of the total error, where we assume that measurements taken in each frame m are independent of each other. This is a reasonable assumption, since we took care to ensemble average over frames separated by times greater than the autocorrelation time.
We determined error associated with the self-part of the intermediate scattering function in a similar manner: where k l is a point randomly chosen on the surface of a sphere of radius k via the Marsaglia 10 method during iteration l, and there are L total iterations. Then the variance, again by Taylor series error propagation, is: where σ 2 m (t) is the variance of the mean self-intermediate scattering function taken across L values of k l . Error associated with χ SS 4 (t) could not be dealt with in an identical manner, since it is a function of the clearly correlated variables Q S (t) and Q 2 S (t) : where w jm (t) = 1 if | r jm (t) − r jm (0)| < a and 0 otherwise, for some length a. χ SS 4 (t) is itself a (scaled) sample variance of Q S (t) across M samples. The sample variance distribution is well-characterized (see, eg., http://mathworld.wolfram.com/SampleVarianceDistribution.html) and we can use its established (scaled) variance as the variance of χ SS 4 (t): where µ Q 4 (t) is the fourth moment of the distribution of Q S (t) across M frames, and µ Q 2 (t) is its second moment.
To find the error associated with α(t), a function of the higher-order sample kurtosis of average particle displacements, we resorted to using a jackknife resampling technique 11,12 . We computed α(t) as: We then used the jackknife variance estimate: where α (m) (t) is α(t) computed with the m-th simulation frame removed.

E. Crystal stability tests
To test the stability of various crystal structures at any location in shape space, we replaced the particles of a well-equilibrated simulation snapshot of the relevant crystal at φ = 0.62 with particles of the desired shape, while leaving particle positions and orientations unchanged. We then sampled in the isochoric ensemble to eliminate particle overlaps, isotropically enlarging the simulation box by a small scale factor every 10,000 MC sweeps if overlaps still existed. We subsequently compressed the system to some initial density between φ = 0.62 and φ = 0.66 in the manner detailed in Supplementary Method I A, and allowed it to slowly melt. During the melting process, system densities were decreased in increments of ∆φ = 0.01 every 10 million MC sweeps.
Step size tuning was performed during the melting process to maintain sampling efficiency.
Pressure was calculated using a volume perturbation technique 13-15 that extrapolates pressure in hard particle isochoric simulations through evaluations of the volume scaling needed to cause particle-pair overlaps throughout the system. Its implementation in HOOMD-blue is detailed elsewhere 16,17 . Melting of the crystal structures was determined by eye and corroborated by melting equations of state; for each melting event, pressure exhibited behavior characteristic of a phase transition, and often showed non-monotonic behavior in the form of Mayer-Wood loops 18,19 . Three melting replicates were run for each state point.
Equations of state for each disordered, super-compressed fluid, shown in the upper left panel of Supplementary Figs. 6 and 7, were determined from system snapshots used to calculate dynamical order parameters. At each density, pressure was calculated for 100 snapshots, taken every 1 million MC sweeps. Pressures for each snapshot were then ensemble averaged, and the equation of state shows these averages and associated standard deviations of the mean at each packing fraction. Pressure is reported in reduced units: p * ≡ βpσ 3 , where β ≡ 1/k B T , p is pressure, and σ = 1 is the characteristic length of our systems. Supplementary Fig. 2 is rather straight-forward, and shows that at densities relevant to crystallization, the disordered fluid at (α a , α c ) = (0, 0.5) contains fewer motifs associated with the diamond structure than nearby fluids that assemble into diamond, and fewer motifs associated with the dodecagonal quasicrystal than the nearby fluid that is capable of assembling into the dodecagonal quasicrystal. Thus, the disordered fluid is structurally different from nearby fluids that assemble into crystal structures, and reflects a higher competition between the face-to-face twisted (pink) motif associated with the diamond structure and the face-to-face aligned (red) motif associated with the dodecagonal quasicrystal. Note that in panel (B), at φ = 0.54, the system at (α a , α c ) = (0, 0.4) forms a coexistence between the diamond structure and the fluid, and the crystal (dot-dash) motif fractions at this density reflect that. In both panels (B) and (C), at high enough densities, the fraction of twisted face-to-face pairwise motifs shown in dashed pink triangles is high enough to promote self-assembly into the diamond structure, at which point the aligned face-toface motif shown in dot-dash red triangles is strongly suppressed in favor of the twisted face-to-face motif shown in dot-dash pink triangles. We observe some portion of the particles in the motif indicated by green triangles, associated with the fcc structure, in all systems at all densities. However, plots of P obs (f, θ) for these systems show that this is essentially just an artifact of imposing cut-offs on the misorientation angle to define our motifs-these systems do not have any special spike in probability near 58 • , the misorientation angle associated with the face-edge connected fcc motif. Supplementary Fig. 3 is more complicated, due to the presence of more types of pairwise motifs in systems in this region of shape space, but nevertheless portrays a similar story to Supplementary Fig. 2. At densities relevant to crystallization, the disordered system at (α a , α c ) = (0.2, 0.5) contains fewer motifs associated with the diamond structure than the nearby fluid that assembles into diamond, fewer motifs (almost none in this case) associated with the bcc structure than the nearby fluid that assembles into bcc, and fewer motifs associated with the fcc structure than the nearby fluids that assemble into γ-brass and fcc. (Whether the system assembles into γ-brass or fcc seems to depend on the precise cocktail of motifs in the pre-cursor fluid at appropriate densities.) We also point out that the disordered system contains about half as many motifs associated with the dodecagonal quasicrystal than the fluid at (α a , α c ) = (0, 0.6), which is capable of assembling the quasicrystal. Thus, the disordered fluid is structurally different from nearby fluids that assemble into crystal structures, and reflects a higher competition between the face-to-face twisted (pink) motif associated with the diamond structure, the edge-to-face (green) motifs associated with the fcc structure, and the face-to-face aligned (red) motif associated with the nearby quasicrystal. In panel (C), the system at (α a , α c ) = (0.25, 0.5) does assemble into the γ-brass structure at φ = 0.6, but the fraction of motifs in the pre-cursor fluid is quite similar to the fraction of motifs in the assembled structure, so it is difficult to delineate the crystalline motif fractions in the plot. In panel (D), we also note the non-negligible presence in the fcc-forming fluid of the edge-connected motifs shown as light and dark green squares, also associated with the fcc structure. These motifs are more suppressed in the disordered fluid at all densities.  in disordered or pre-cursor fluids in regions of shape space in which fluids tend to self-assemble into the crystals associated with those motifs. Thus, (i) the motifs associated with bcc are strongly suppressed in all fluids except near the α a = α c line, where vertex truncation is highest, (ii) motifs associated with fcc are more abundant in regions with higher α a , or edge truncation, (iii) the quasicrystal motif is more abundant near the (α a , α c ) = (0, 1) corner of shape space corresponding to a non-truncated tetrahedron, and (iv) the diamond motif is more abundant near the (α a , α c ) = (0, 0.35) location in shape space that corresponds to a vertex-truncated tetrahedron. As density increases, the regions in which the bcc and fcc motifs are abundant become smaller and more concentrated near the α a = α c line and at higher α a values respectively. Conversely, the regions of quasicrystal and diamond structure motif abundance grow as density increases, and they grow in directions in which the corresponding crystals still self-assemble at higher density. The behavior of all regions as density increases makes sense in the context of locally dense packing arguments in hard particle fluids 26 : vertex/edge connected motifs only appear in systems in which particle vertex/edge truncation is significant, and these connection types are suppressed as density increases because face connection is enhanced as density increases. Face connection provides higher locally dense packing of particle pairs.
Note that particles on the diagonal of the shape landscape, where α a = α c , have octahedral symmetry rather than tetrahedral symmetry. Thus, their set of possible misorientation angles is different: in particular, it is not possible for particles with this symmetry to have misorientation angles between ∼ 65 • and 90 • , so systems of these particles somewhat artificially display zero motifs associated with the dodecagonal quasicrystal or diamond structure.

C. Non-Gaussianity of MC sampling at short times
Here we provide some justification for the measured increase of the non-Gaussian parameter α(t) in our systems at short times. Monte Carlo methods can simulate a diffusive process if only local moves are made. In this case, the simulation is effectively a random walk, which becomes a Gaussian distribution as the  number of steps taken on the walk goes to infinity. As the number of steps taken goes to zero, however, the system becomes decidedly non-Gaussian, and α > 0.
To show this, we idealize the MC process, and compute α exactly in a toy model. First let's start in one dimension, and give our MC sampling method three options: particles can either remain in place with probability s, move to +L with probability m/2, or move to −L with probability m/2, with m+ s = 1. This is technically a trinomial distribution. We will follow the treatment of a random walk in Nelson's Biological Physics 20 : Let the displacement of step j be k j L, where k j = 0 with probability s, and k j = ±1 with probability m/2. Then consider all possible trajectories of N steps: To get from line 2 to line 3 above, note that k 2 N = 0s + m/2 + m/2 = m, and that x N −1 k N evaluates to zero. This is because x N −1 and k N are uncorrelated, so we can split that average of a product into a product of averages, and k N = 0. As another explanation, note that for every value of x N −1 , there are 3 contributions to the ensemble average: k N = 0 with probability s, k N = 1 with probability m/2, and k N = −1 with probability m/2. These average to zero.
We can compute that x 2 1 = x 2 0 + L 2 m = L 2 m, and then use the recursion relation to determine that: Supplementary Figure  We can also consider x 4 N : To get from line 2 to line 3, first note that k 4 N = 0s + m/2 + m/2 = m. The terms that contain odd powers of k N go to zero as per our previous argument involving x N −1 k N . Then, to evaluate x 2 N −1 k 2 N , we can again note that x N −1 and k N are uncorrelated, and k 2 N = m, so x 2 N −1 k 2 N = m x 2 N −1 . Line 4 proceeds from line 3 by the previous result: x 2 N −1 = (N − 1)mL 2 . We can compute that x 4 1 = x 4 0 + L 4 m = L 4 m, and use the recursion relation to determine that: To get from line 1 to line 2, note that the first sum can be broken into pairs, (0 + N ) + (1 + N − 1) + . . . . If N is odd, there are exactly N +1 2 of these pairs, for a total sum of N N +1 2 . If N is even, there are exactly N 2 of these pairs, and there is a left-over term N 2 that also contributes to the sum, for a total sum of N N 2 + N 2 = N 2 (N + 1). Thus, As N → ∞, the first term dominates the expression, and can be proven to be true for the Gaussian distribution in 1D. Let's set N = 1 (to mimic the short time limit in our MC simulation): As m → 0, or the probability of moving becomes increasingly unlikely, this expression diverges, and the distribution consequently gets increasingly "tail-heavy." We now move to 3D by employing a few tricks: and Thus, The non-Gaussian parameter α is correspondingly: We can note several things about the above quantity. If m = 1, and a move is always made, r 4 1 = r 2 1 2 and α is negative. The distribution resembles that of a Bernoulli distributed coin toss, which is known to have a negative excess kurtosis. If m = 1/3, α = 0. For m < 1/3, however, α > 0, and the system becomes increasingly tail-heavy. At the reasonable translation acceptance ratio m = 0.1, for example, α = 1.4. Indeed, this number is comparable to values of α at small times shown in Fig. 2 in the main text.

D. Crystal stability tests
Supplementary Figs. 6 and 7 show results of stability testing of various crystal structures as functions of particle shape, at or near glass-forming locations in shape space. Near the glass-forming location at (α a , α c ) = (0, 0.5), shown in Supplementary Fig. 6, we studied the stability of (i) the diamond structure self-assembled at (α a , α c ) = (0, 0.4) as a function of increasing α c , (ii) the dodecagonal quasicrystal selfassembled at (α a , α c ) = (0, 0.6) as a function of decreasing α c , and (iii) the diamond structure self-assembled at (α a , α c ) = (0.1, 0.5) as a function of decreasing α a . Near the glass-forming location at (α a , α c ) = (0.2, 0.5), shown in Supplementary Fig. 7, we studied the stability of (i) the diamond structure self-assembled at (α a , α c ) = (0.1, 0.5) as a function of increasing α a , (ii) the bcc structure self-assembled at (α a , α c ) = (0.2, 0.4) as a function of increasing α c , and (iii) the fcc structure self-assembled at (α a , α c ) = (0.3, 0.5) as a function of decreasing α a .
Melting plots show the solidus line, or lowest density at which systems remain fully crystallized, and the liquidus line, or lowest density at which crystals coexist with the fluid. In all cases, the densities shown are the highest found across all replicates, since our method establishes a lower bound for melting density. Highest liquidus or solidus densities across replicates thus represent the most restrictive lower bound. Each solidus line is labeled by the crystal structure that is stable above the line. In some cases, systems undergo phase transitions during the melting process to other solids. Whether the system passes through a fluid phase during that process or undergoes a solid-solid phase transition is not shown here, because our method is not rigorous enough to determine the nature of these transitions. Instead, we simply show stability lines for all observed structures.
We now describe the phases displayed in the melting plots that are not identical to those self-assembled from the fluid. The "hR6-SbSn/hR6/oC4" melting line shows the stability of three interrelated phases, whose delineation was often not clear during the melting process itself. "hR6-SbSn" is a lower (hexagonal) symmetry diamond derivative, featuring a coordination number of 4 for all particles and squashed tetrahedral local environments. Its Pearson symbol is hR6, and its space group is R3m, number 166. We use the shorthand hR6-SbSn because compounds of Sb and Sn have been found to crystallize into this structure 21 . "hR6" is a slightly distorted version of hR6-SbSn, with the same space group. Next nearest neighbor distances are shifted closer to each particle, such that there is no clear peak in the radial distribution function corresponding to four nearest neighbors, and instead the structure could be described with a larger coordination number of 14. "oC4" is a related structure of lower symmetry with coordination number 12. Its Pearson symbol is oC4, and its space group is Cmcm, number 63.
The "bcc (OO)" melting line corresponds to orientationally-ordered bcc, featuring all particles oriented in the same direction. The "bcc (OD)" melting line corresponds to orientationally-disordered bcc, featuring plastic-like particle orientations. Equations of state indicate a phase transition between these two phases, so we mark them distinctly. For more information regarding orientational phase transitions in hard particle systems, see Karas et al. 22 The "distorted/tetragonal diamond" melting line describes systems that shear and distort, occasionally managing to form a lower symmetry tetragonal diamond derivative, when melted from cubic diamond at higher α c . This signature appears in the pressure and indicates a phase transition between this strained or lower symmetry phase and (cubic) diamond at lower densities. Because we keep the box cubic, the phase transition is not clean, and strain occurs; however, "floppy box" simulations in these α c regimes, in which we allow the box aspect ratio and box shear to change randomly and independently (while keeping box volume fixed), do show clear phase transitions between a lower symmetry diamond derivative phase at high packing fraction and cubic diamond at lower packing fraction. This phase transition in these systems was first observed by Cersonsky et al. 23 The lower-density (cubic) diamond structure melts at approximately the same density in these floppy box simulations as it does in the simulations in which we keep the box fixed and cubic.
We also used floppy box Monte Carlo to investigate the melting of the other crystals shown in these figures  at identical locations in shape space; we did not observe any qualitatively different behavior in the melting lines or the stable structures exhibited by each system with decreasing density.

E. The identity crisis in the 423 family
Here we describe an identity crisis in a disordered system in the shape space composed of polyhedral particles in the spheric triangle invariant 423 family 24 . This set of convex polyhedra is formed by truncating the vertices and edges of an octahedron by sets of planes at varying radial distances from the polyhedron center. As we did for the 323 family described in the main text, we define truncation parameters α a and α c such that (α a , α c ) = (0, 0), denoting a cuboctahedron, (α a , α c ) = (0, 1), denoting a cube, (α a , α c ) = (1, 1), denoting a rhombic dodecahedron, and (α a , α c ) = (1, 0), denoting an octahedron, form the corners of this shape space. Supplementary Fig. 8A shows the full shape space, with representative polyhedra superimposed above their corresponding positions in the space. It was previously found that at certain locations in this shape space near the octahedron corner (α a , α c ) = (1, 0), systems failed to assemble into any ordered structure at densities ranging from φ = 0.50 to 0.65, while at nearby locations, bcc or occasionally highpressure Lithium formed only at high packing fractions at or above φ = 0.6 1 .
We investigated the region near the octahedral corner of this shape space using Monte Carlo methods We tested for the stability of the high-pressure Lithium structure at every investigated state point by running melting simulations in which we initialized our systems in the Lithium structure at packing fractions ranging from φ = 0.56 to 0.64, and subsequently sampled in the isochoric ensemble for 39 to 61 million MC sweeps or until melting into bcc occurred. We initialized in the Lithium structure via a compression scheme similar to that described earlier, with a few changes. We set up our particles in a high-pressure Lithium structure at very low packing fraction, then compressed quickly to the desired density with a very small translational trial move size and a much larger rotational trial move size, to ensure that the particles remained on-lattice during the compression. We ran three replicate simulations at each state point and each density. We additionally ran one corresponding simulation at each state point and density in which we tested for the stability of bcc in an identical manner. We found that, on these time scales, bcc was stable at all densities for all state points, and Lithium was stable at a subset of (higher) densities for all state points. We found that Lithium was stable at φ ≥ 0.58 for (α a , α c ) = (0.85, 0.05), (0.9, 0.05), (0.95, 0.05), and (0.95, 0.1); at φ ≥ 0.6 for (α a , α c ) = (0.85, 0.1), (0.9, 0.1), (0.9, 0.15), and (0.95, 0.15); and at φ ≥ 0.62 for (α a , α c ) = (0.85, 0.15). At other densities, all replicates of Lithium melted into bcc. The trend is clear: as edge and vertex truncation generally increase, and thus α a decreases and/or α c increases, Lithium is only stable at higher and higher densities.
These assembly and melting simulations demonstrate that high-pressure Lithium can form and remain stable in this region of shape space, especially at small edge and vertex truncations, but Lithium appears generally metastable to bcc. Supplementary Fig. 8B shows the results of pairwise motif identification for the disordered system at (α a , α c ) = (0.9, 0.1) and φ = 0.62, as well as surrounding crystal structures. Motif fractions were found by averaging across all snapshots (separated by 1 million MC sweeps) for which systems were fully crystallized, or all snapshots in the case of the disordered system. The disordered system is dominated by two motifs, colored light and dark purple, that are each prevalent in nearby bcc (at φ = 0.62) and Lithium (at φ = 0.61) structures respectively. In direct parallel to the motifs that dominate the diamond structure and the dodecagonal quasicrystal within the 323 shape family studied elsewhere in this paper, we call these the "twisted" and "aligned" motifs, respectively. The "aligned" motif, detailed in row 5 in Supplementary Fig. 8B, dominates in a nearby Lithium structure and consists of two truncated octahedra face-to-face and rotated such that their truncated vertices are aligned. The "twisted" motif, detailed in row 4 in Supplementary Fig. 8B, dominates in a nearby bcc structure and consists of two truncated octahedra face-to-face and rotated such that they have the same orientation, and thus their bonded faces are twisted 60 • with respect to each other. Supplementary Fig. 8C shows a subset of results for alchemical Monte Carlo (Alch-MC) sampling at the disordered state point (α a , α c ) = (0.9, 0.1). At densities ranging from φ = 0.52 to 0.64, we equilibrated our systems for 10 million MC sweeps, then allowed only the vertex truncation parameter α a , only the edge truncation parameter α c , or both α a and α c to fluctuate as thermodynamic variables for 29.7 -40 million MC sweeps. Due to the increased computational time required for particle-particle overlap checks after every alchemical shape move attempt in these systems, we attempted shape moves with a 25% probability after every 10 MC sweeps for all vertex and vertex+edge truncation sampling simulations, rather than following the protocol described in the main text in which we sampled shape moves with a 25% probability after every MC sweep. During the edge truncation sampling simulations, we attempted shape moves with a 25% probability after every 10 MC sweeps for 26.2 -28.5 million MC sweeps, and attempted shape moves with a 25% probability after every MC sweep for the remainder of the time. We constrained each system to only explore the area inside a square of side length ∆α = 0.1 centered at its initial position in shape space. Selfassembly on our simulation time scales occurred at φ = 0.6 when we allowed only α a to fluctuate, at φ = 0.6 and φ = 0.62 when we allowed only α c to fluctuate, and at φ = 0.6 when we allowed both α a and α c to fluctuate. In these simulations, disordered systems escaped into regions of larger truncation and assembled into bcc, with a larger fraction of their particles accordingly adopting the twisted motif. Supplementary  Figs. 22  family, with the portion of the shape space explored in this study outlined. That portion is shown in greater detail, with sample particle shapes overlaid above corresponding regions of shape space, and regions colored according to their assembled structure.
(B) Pairwise motifs in the disordered system at (αa, αc) = (0.9, 0.1) compete, and are found to dominate in nearby ordered structures in shape space. Relevant locations in shape space are outlined in black in panel A. Motifs are listed in tabular form and color-coded according to their connection type ct and θ range. (C) The disordered system escapes its region of identity crisis and crystallizes when allowed to explore its surrounding shape space via Alch-MC. The empty square in the shape space diagram to the upper right indicates system position in shape space at the start of Alch-MC sampling, and letters indicate system position after ∼ 30-40 million MC sweeps of vertex truncation (v) sampling at φ = 0.6, edge truncation (e) sampling at φ = 0.62, or both vertex and edge truncation (ve) sampling at φ = 0.6. System snapshots, particle shapes, pie charts of pairwise motif fractions, and bond-order diagrams are shown for initial and final frames of each Alch-MC simulation.

III. SUPPLEMENTARY FIGURES
B A Supplementary Figure. 9: (A) Lowest packing fraction for which crystallization was observed and (B) crystallization time at that packing fraction across the shape landscape. If multiple replicates were observed to crystallize at the minimum packing fraction, the minimum crystallization time is shown. Color bars below each panel show relevant limits and scales. Squares are left uncolored if crystallization was not observed at any density at the corresponding location in shape space, or if the location in shape space is outside the bounds studied in this paper. We define crystallization time, or so-called nucleation incubation time 25 , as the first frame after which approximately all crystalline particle fractions measured over the trajectory are greater than 0.1. For nucleation of the dodecagonal quasicrystal, we estimated the crystallization time by eye, corroborating our observations by calculating pressure over the trajectory (using volume perturbation methods discussed in the Supplementary Methods) when that data was available, and confirming that pressure begins to drop to its crystal value around the estimated crystallization time. For a more detailed map of the lowest packing fraction for self-assembly across this shape landscape, see Klotsa      The dopant dimer in this case is the face-to-face aligned motif, associated with the dodecagonal quasicrystal. At high enough dopant dimer fraction, crystallization into diamond is suppressed, and the system is dominated by the face-to-face aligned motif, shown in red in the motif fraction plots. Vertical lines through plots of crystalline particle fraction indicate nucleation incubation time, defined in the main text.
Symbols to the right of the first panel indicate the connection type of each motif shown.
Supplementary Figure. 16: Alchemical Monte Carlo simulations for (αa, αc) = (0, 0.5) at various packing fractions φ, where vertex truncation is sampled as a thermodynamic variable. Plots show the fraction of particles in each motif, the truncation parameter, and the fraction of particles with a diamond crystalline environment (whose detection is detailed in the Supplementary Methods). Motif fractions are separated by black lines and colored according to the motifs identified in Fig. 3 of the main text. Regions colored gray represent (connection type, θ) regimes that were not identified with any crystal structure. Regions colored identically represent motifs characteristic of the same crystal structure that differ only by connection type. In those cases, the motif with face connection is always shown above the motif with vertex or edge connection. In most cases, systems escape the region of disorder and crystallize into the diamond structure, dominated by the pink twisted motif.
Supplementary Figure. 17: Alchemical Monte Carlo simulations for (αa, αc) = (0, 0.5) at various packing fractions φ, where edge truncation is sampled as a thermodynamic variable. Plots show the fraction of particles in each motif, the truncation parameter, and the fraction of particles with a diamond crystalline environment (whose detection is detailed in the Supplementary Methods). Motif fractions are separated by black lines and colored according to the motifs identified in Fig. 3 of the main text. Regions colored gray represent (connection type, θ) regimes that were not identified with any crystal structure. Regions colored identically represent motifs characteristic of the same crystal structure that differ only by connection type. In those cases, the motif with face connection is always shown above the motif with vertex or edge connection. At high enough packing fractions, systems escape the region of disorder and crystallize into a diamond-like structure dominated by the pink twisted motif. This diamond-like structure contains wurtzite-like structural motifs, as described in the main text.
Supplementary Figure. 19: Alchemical Monte Carlo simulations for (αa, αc) = (0.2, 0.5) at various packing fractions φ, where vertex truncation is sampled as a thermodynamic variable. Plots show the fraction of particles in each motif, the truncation parameter, and the fraction of particles with a bcc crystalline environment (whose detection is detailed in the Supplementary Methods). Motif fractions are separated by black lines and colored according to the motifs identified in Fig. 3 of the main text. Regions colored gray represent (connection type, θ) regimes that were not identified with any crystal structure. Regions colored identically represent motifs characteristic of the same crystal structure that differ only by connection type. In those cases, the motif with face connection is always shown above the motif with vertex or edge connection. At high enough packing fractions, systems escape the region of disorder and crystallize into bcc. As packing fraction increases, this structure is increasingly dominated by the blue face-vertex connected motif shown in Fig. 3 of the main text.
Supplementary Figure. Fig. 3 of the main text. Regions colored gray represent (connection type, θ) regimes that were not identified with any crystal structure. Regions colored identically represent motifs characteristic of the same crystal structure that differ only by connection type. In those cases, the motif with face connection is always shown above the motif with vertex or edge connection. At φ = 0.58 and φ = 0.6, systems escape the region of disorder and crystallize into bcc and fcc respectively. Interestingly, both structures in this case are dominated by the green face-edge connected motif shown in Fig. 3 of the main text. This motif is typically associated with fcc only.
Supplementary Figure. Fig. 3 of the main text. Regions colored gray represent (connection type, θ) regimes that were not identified with any crystal structure. Regions colored identically represent motifs characteristic of the same crystal structure that differ only by connection type. In those cases, the motif with face connection is always shown above the motif with vertex or edge connection. In all cases, systems escape the region of disorder and crystallize into fcc (or bcc for φ = 0.64). fcc is dominated by the green faceedge connected motif shown in Fig. 3 of the main text. The bcc structure assembled at φ = 0.64 displays an increased fraction of particles in the blue face-vertex connected motif shown in Fig. 3    where vertex truncation is sampled as a thermodynamic variable. Plots show the fraction of particles in each motif, the truncation parameter, and the fraction of particles with a bcc crystalline environment (whose detection is detailed in the Supplementary Methods). At φ = 0.6, the system increases its vertex truncation and crystallizes into bcc, dominated by the light purple twisted motif.  Figure. 23: Alchemical Monte Carlo simulations for (αa, αc) = (0.9, 0.1) at various packing fractions φ, where edge truncation is sampled as a thermodynamic variable. Plots show the fraction of particles in each motif, the truncation parameter, and the fraction of particles with a bcc crystalline environment (whose detection is detailed in the Supplementary Methods). At φ = 0.6 and φ = 0.62, the system increases its edge truncation and crystallizes into bcc, dominated by the light purple twisted motif.