Entropically engineered formation of fivefold and icosahedral twinned clusters of colloidal shapes

Fivefold and icosahedral symmetries induced by multiply twinned crystal structures have been studied extensively for their role in influencing the shape of synthetic nanoparticles, and solution chemistry or geometric confinement are widely considered to be essential. Here we report the purely entropy-driven formation of fivefold and icosahedral twinned clusters of particles in molecular simulation without geometric confinement or chemistry. Hard truncated tetrahedra self-assemble into cubic or hexagonal diamond colloidal crystals depending on the amount of edge and vertex truncation. By engineering particle shape to achieve a negligible entropy difference between the two diamond phases, we show that the formation of the multiply twinned clusters is easily induced. The twinned clusters are entropically stabilized within a dense fluid by a strong fluid-crystal interfacial tension arising from strong entropic bonding. Our findings provide a strategy for engineering twinning behavior in colloidal systems with and without explicit bonding elements between particles.

Sangmin Lee 1,3 & Sharon C. Glotzer 1,2 Fivefold and icosahedral symmetries induced by multiply twinned crystal structures have been studied extensively for their role in influencing the shape of synthetic nanoparticles, and solution chemistry or geometric confinement are widely considered to be essential. Here we report the purely entropydriven formation of fivefold and icosahedral twinned clusters of particles in molecular simulation without geometric confinement or chemistry. Hard truncated tetrahedra self-assemble into cubic or hexagonal diamond colloidal crystals depending on the amount of edge and vertex truncation. By engineering particle shape to achieve a negligible entropy difference between the two diamond phases, we show that the formation of the multiply twinned clusters is easily induced. The twinned clusters are entropically stabilized within a dense fluid by a strong fluid-crystal interfacial tension arising from strong entropic bonding. Our findings provide a strategy for engineering twinning behavior in colloidal systems with and without explicit bonding elements between particles.
Twinning arises from a type of grain boundary called a twin boundary, where two separate crystals sharing the same lattice plane intergrow with a certain symmetry. Fivefold and icosahedral symmetries, which are generally known to be incompatible with long-range order, can be induced in relatively large atomic crystal clusters (from several nanometers to a few microns) via multiple twinning [1][2][3][4] . For instance, materials that form face-centered cubic (fcc) 2,3,5-7 or cubic diamond 8,9 crystals can form twin boundaries toward the (111) or its equivalent lattice directions, promoting the formation of a ∼ 70 twin angle that can produce multiply twinned structures with fivefold or icosahedral symmetry. Twinning has been widely used to obtain interesting properties in synthetic nanomaterials, such as enhancing the mechanical properties of nanowires 10 and increasing oxidation resistance 11 , and to synthesize noble metal nanostructures with decahedral or icosahedral shape 3,12 , which are useful as catalysts 13 . Many studies have been conducted to control the growth mechanism 12,14,15 and stability of multiply twinned structures, e.g., tailoring solution chemistry 5,6,[16][17][18] and using multicomponent materials 6,11 . However, the multiplicity of variables in multiply twinned crystals makes it difficult to know which variables are responsible for, or have the most influence on, the formation of fivefold and icosahedral twins.
Compared to most systems, hard particle systems are regarded as relatively simple because the phase behavior is driven solely by entropy maximization 19,20 . The diversity and complexity of entropydriven phase behavior, however, are enormously rich. The selfassembly of complex crystals and quasicrystals 20,21 , as well as complex assembly pathways such as multistep crystallization, have been reported in hard particle systems 22 . The spontaneous formation of icosahedral twins has been observed in hard sphere systems both in experiments and in computer simulations, but only under spherical confinement 4,7,15,23 , which produces an artificial "surface tension" to overcome the internal strain within the cluster arising from local particle packing. That, and the negligible free energy difference between fcc and hexagonal close-packed (hcp), two competing crystal structures of hard spheres 24 , conspire to produce twins. However, the formation and stabilization of fivefold and icosahedral twinned clusters without geometrical confinement or in any other hard particle systems or have yet to be reported. We hypothesized that artificial surface tension provided by spherical confinement could be achieved naturally, without any confinement, through the judicious choice of particle shape.
Here, we show the purely entropy-driven, confinement-free assembly of fivefold and icosahedral twinned clusters from multiply twinned diamond crystals in equilibrium with a dense fluid phase. Monte Carlo (MC) simulations show that hard truncated tetrahedrons (TTs) self-assemble into either cubic and hexagonal diamond crystals depending on the amount of edge and vertex truncation (Fig. 1). We tuned the TT shape to have a negligible free energy difference between the two diamond crystal phases, and found the formation of twin boundaries is easily induced in a fluid phase (Fig. 2). This finding demonstrates that the stability of twin boundaries in a fluid can be controlled by particle shape design, even entropically. Using this strategy, we induce the formation of fivefold and icosahedral twins through seed-assisted growth of a single-crystalline seed of cubic diamond. Since the formation of twin boundaries is easily induced toward multiple directions during crystal growth, multiply twinned structures are easily formed. We show that, through an error-and-repair process, a multiply twinned structure with defects transforms into a fivefold twinned cluster in a dense fluid phase (Fig. 3). We also show the formation of an icosahedral twin in a dense fluid upon further growth of the fivefold twin. The icosahedral twin of hard TTs is entropically stabilized within the fluid without spherical confinement (Fig. 3), unlike the icosahedral cluster of hard spheres that rapidly destabilizes and falls apart (Fig. 4). We show that icosahedral clusters of hard TTs have twice the fluid-solid interfacial free-energy (or entropy) compared to icosahedral clusters of hard spheres as a natural consequence of stronger entropic bonding 25 in the former system. Our study isolates the essential role of surface tension in stabilizing the inherently strained icosahedral twinned cluster in a dense fluid, and suggests approaches for engineering the surface tension.

Results
Cubic diamond is different from hexagonal diamond (or lonsdaleite) in the conformation of local tetrahedral bonds 26  diamond has three staggered bonds and one eclipsed bond. For a hard particle with a tetrahedron shape, the two types of chemical bond conformations can be mapped to the staggered and eclipsed entropic bond conformations associated with face-to-face contacts. (Fig. 1a). An appropriate degree of tip truncation of hard tetrahedra allows them to exclusively form staggered contacts, resulting in the self-assembly of the cubic diamond phase 27,28 (Fig. 1b). However, the self-assembly of the hexagonal diamond phase, where the two bond conformation types are mixed, has not been reported in hard particle systems. We hypothesized that, through judicious truncation of a regular tetrahedron, we can find an appropriate shape that introduces both staggered and eclipsed contacts and favors hexagonal diamond (Fig. 1c), and, intermediate between the cubic diamond-forming and hexagonal diamond forming shapes, one that produces both diamond structures with negligible difference in entropy. Controlling the relative stability of cubic and hexagonal diamond will give us a route to control  control of diamond crystals via particle shape design. a Perparticle Helmholtz free-energy (F=Nk B T) plot, in units of k B T, of cubic (red) and hexagonal (blue) diamond crystals of hard TTs as a function of vertex and edge truncation parameters at constant particle volume fraction (ϕ = 0:62). b Phase diagram of hard TTs in the shape space determined by the free-energy calculation. c Free-energy difference between cubic and hexagonal diamond (4F = F H À F C ) in shape space. d, When the TT is designed ða = 1:20, c = 2:16Þ to have negligible 4F ∼ + 0:007, the initial cubic diamond (left) single-crystalline cluster forms multiple stacking faults (right) in fluid, suggesting that the free energy loss from the twin boundaries is small. Each TT is represented by a tiny sphere at the TT center of mass. Red and blue spheres represent cubic and hexagonal diamonds, respectively. e When the TT is designed (a = 1:22, c = 2:16) to have 4F ∼ À 0:028, the initial cubic diamond cluster (left) completely transforms into single-crystalline hexagonal diamond (right) cluster in a relatively short simulation time. f The change of the number ratio of particles in cubic diamond over time for three different 4F = + 0:138, + 0:007 and À0:028 (black X markers in c). Orange and cyan circles indicate each snapshot in d and e, respectively. g, Local volume fraction distribution plots calculated from the last snapshot from d (upper) and e (lower). Two peaks in the distribution show that the system is in coexistence between fluid and solid. twinning and stacking behaviors of diamond crystals ( Fig. 1b-d). Cubic and hexagonal diamond share an equivalent plane along the ð111Þ and ð0001Þ directions (Fig. 1b, c); fcc and hcp share the same equivalent plane. If the entropy difference between the two crystals can be engineered to be sufficiently small, the crystal will be prone to form twin boundaries. The twinning of cubic diamond occurs towards ð111Þ and its equivalent planes such as ð 111Þ, ð1 11Þ and ð11 1Þ; thus when twinning occurs in multiple directions, two twin planes meet with a specific angle, cos À1 ð1=3Þ = 70:5 , which is close to 72 , the same angle that promotes the formation of multiply twinned structures with fivefold symmetry 2 (Fig. 1e). The structure of the multiply twinned crystals made by TTs is easily analyzed by classifying particles based on the number of staggered and eclipsed contacts (Supplementary Figs. 1, 2), following the notation S n E m , where n + m = 4. For instance, the cubic diamond crystal consists of S 4 E 0 particles (Fig. 1a, b), and the hexagonal diamond crystal consists of S 3 E 1 particles (Fig. 1a, c). In a fivefold twin, five cubic diamond clusters with a tetrahedron shape made by S 4 E 0 particles together form a pentagonal bipyramid super-cluster with five twin planes of S 3 E 1 particles between the faces of the tetrahedron clusters (Fig. 1e). At the point where the five twin planes meet, S 2 E 2 particles form a column where small pentagonal bipyramids (N = 5) are linearly stacked. An icosahedral twin (Fig. 1f) comprises twenty cubic diamond clusters (S 4 E 0 particles), each with a tetrahedron shape (Fig. 1g). Thirty twin planes (S 3 E 1 ) exist between the faces of the tetrahedron clusters (Fig. 1h), and there twelve fivefold columns made by S 2 E 2 particles where five tetrahedron clusters meet ( Fig. 1i). At the center, there is a dodecahedron super-cluster ðN = 100Þ comprising three different shells with icosahedral symmetry: 60 S 2 E 2 particles comprising the outer shell form a rhombicosidodecahedron, 20 S 4 E 0 particles comprising the intermediate shell form a dodecahedron and 20 S 1 E 3 particles comprising the smallest shell form a dodecahedron. This hierarchical structure is equivalent to the supercluster of water ðH 2 OÞ 100 when connecting oxygen atoms 29 .
A competition between staggered and eclipsed contracts arises from many-body interactions responsible for the entropic forces producing the effective attraction between neighboring particles, and preferences to form a certain type of face-to-face TT contact-which are determined by the strength of this attraction-can be controlled by varying the amount of edge and vertex truncation 28 . We know that regular tetrahedra prefer to form a dodecagonal quasicrystal 21 in which every particle has eclipsed alignments (S 0 E 4 ). Truncation introduces additional facets, weakening the strength of the entropic bond between two primary eclipsed faces, and thereby weakening the preference for eclipsed alignment. The entropic bonding between primary eclipsed faces is weakened with increasing truncation because less and less free volume is gained by maintaining this alignment (though this alignment still maximizes free volume and thus system entropy). Beyond a certain amount of truncation, however, the system has more entropy (more free volume) if the neighbors align in a staggered way, and thus the staggered configuration becomes preferred over the eclipsed configuration (S 4 E 0 ) 27 . Motivated by this, we searched the shape space 28 between S 0 E 4 and S 4 E 0 to find a shape intermediate between the two that favors the hexagonal diamond FCC. e-g The change of interfacial profiles of the (111) direction of (e) cubic diamond and (g) FCC. f-h Fluid-solid interfacial free energy γ of (f) cubic diamond and (h) FCC. For all three lattice directions, the hard TT system has more than twice the value of γ of the hard sphere system, indicating its stronger fluid-solid surface tension.
phase ðS 3 E 1 Þ. To quantify the relative thermodynamic stability of cubic and hexagonal diamond as a function of TT shape, we constructed free energy surfaces in TT shape space varying the edge (1:14 ≤ a ≤ 1:30) and vertex (2:10 ≤ c ≤ 2:30) truncation at constant volume fraction (Fig. 2a), using the Frenkel-Ladd free energy calculation method 24,30 ("Methods"). We found a region where hexagonal diamond is thermodynamically more stable than cubic diamond, as shown in the phase diagram generated from the bottom view of the free energy landscape (Fig. 2b). We also confirmed those shapes self-assemble hexagonal diamond (Supplementary Fig. 3).
To assess the stability of the twin boundaries depending on the TT shape, we chose three systems with different relative stabilities between the two diamond phases, ΔF = F H À F C = 0:138,0:007, and À 0:028 (Fig. 2c). For each system, we prepared a cubic diamond cluster fully surrounded by a fluid phase at coexistence and equilibrated the system ("Methods"). We first distinguished particles in fluid and crystal based on the local volume fraction ("Methods"), which shows two peaks in distribution plots that indicate the two coexisting phases (Fig. 2g). Then, we colored particles in the crystal phase based on their face-to-face contact type ðS n E m Þ and analyzed the crystal structures. When the free energy difference is large enough ΔF ∼ + 0:138 ð Þ , the cubic diamond cluster is stable without any notable changes in the crystal structure (Fig. 2f). However, when ΔF is negligible (ΔF ∼ 0:007), the initial cubic diamond cluster repeatedly developed stacking faults (Fig. 2d) during the simulation, rather than stabilizing either diamond phase. This shows that the formation of twin boundaries between the two diamond phases is easily induced when the thermodynamic stability of the two crystals is comparable. On the other hand, when we further increased the stability of hexagonal diamond (ΔF = À0:028), a complete phase transformation into a single-crystalline hexagonal diamond was observed (Fig. 2e). The phase transition behavior of the three systems was quantitatively compared by the change of the number ratio of S 4 E 0 particles in a solid phase (N cubic =N solid ) (Fig. 2f). We confirmed a more gradual decrease of the N cubic =N solid in the ΔF = 0:007 system compared to that of the ΔF = À0:028 system, due to frequent formation of stacking faults.
Based on these findings, we hypothesized that the formation of a multiply twinned structure might be easily induced during the growth of the cubic diamond when ΔF is negligible because of the easy formation of twin boundaries in multiple directions, which kinetically promotes the formation of~72 twin angles. To confirm the hypothesis, we conducted seed-assisted growth simulations ("Methods") with a system ða = 1:20, c = 2:16Þ where ΔF is negligible ðΔF = 0:007Þ. We placed a small spherical seed of cubic diamond crystal (N = 500) in a fluid and allowed the seed to grow at the coexistence volume fraction ϕ = 0:58 where the resulting crystal cluster is fully surrounded by a fluid phase in equilibrium. In the early stage of growth, we observed the formation of twin boundaries in multiple directions (Fig. 3a, b). Through subsequent multiple twinning events, a fivefold center was formed (Fig. 3c). Although the fivefold centers form a new local structure, S 2 E 2 ( Supplementary Fig. 2)-a pentagonal column -it has a similar local structure to that of hexagonal diamond, resulting in minimal entropy loss. Note that an error or defect can easily occur during the growth process because the formation of twin boundaries independently in multiple directions can cause a mismatch in the stacking sequence between two adjacent directions (Fig. 3d). However, interestingly, we observed a self-repair process during growth, where a growing direction with stacking errors re-melts until the boundary with its adjacent plane matches (Fig. 3e, f). Once the boundary matches, the fivefold twin crystal fully grows and is stabilized in fluid with a truncated pentagonal dipyramid shape (Fig. 3g, h). This error-andrepair mechanism of fivefold twinned crystals occurs via particle migration on the crystal surface, differently from the error-and-repair mechanism recently reported, in which particle rearrangement occurs throughout the fivefold twinned crystal during growth 14 , and differently from that reported recently for a dodecagonal quasicrystal of hard tetrahedra 31 . The observation of this new mechanism demonstrates that there exist multiple types of error-and-repair mechanisms for the formation of fivefold twinned clusters as well as for entropically stabilized colloidal crystals.
Next, we investigated the formation of an icosahedral twinned cluster within a fluid phase. Experiments have shown that an icosahedral twinned cluster can be formed from additional multiple twinning from a decahedral cluster 12,32 . Simulations have shown it can be formed in a purely entropy-driven system provided there is artificial confinement of the cluster. Our simulation results show that this is possible also in a purely entropy-driven system, but without confinement. We conducted a growth simulation of a seed with a fivefold twinned structure ðN = 1020Þ in a system where ΔF is negligible (ΔF ∼ 0:007 at a = 1:20, c = 2:16) ("Methods"), and where the solid cluster is fully surrounded by fluid throughout the simulation. We observed the additional twinning from the surface of the seed, which eventually grew into an icosahedral twin (Fig. 3i-l), a process reminiscent of a simultaneous successive twinning of a palladium icosahedral cluster 32 . The fluid-solid coexistence of our ico-twin was confirmed by the local volume fraction distribution (Fig. 3n) showing two peaks. We classified particles into three types: fluid, interface and crystal. After removing particles at the interface, we confirmed that the icosahedral cluster exposes twenty (111) planes of the cubic diamond crystal toward the fluid (Fig. 3m). This implies that the interfacial entropy difference between fluid and the (111) plane plays an important role in the stability of the icosahedral cluster 17 .
The stability of the icosahedral cluster in the fluid can be checked by tracking the change in reduced system pressure (P * = Pv 0 =k B T) during a simulation run (Fig. 3o). In the early stage ( ≤ 50 × 10 5 MC steps), we observed a rapid decrease in pressure due to the growth of the icosahedral cluster. After the growth (>50 × 10 5 MC steps), the pressure stabilized at a constant value (P * ∼ 12:6) with small fluctuations, indicating that the ico-twin structure had stabilized in the fluid. From the simulation snapshots and diffraction patterns at 100 × 10 5 and 175 × 10 5 MC steps, we confirmed that the icosahedral symmetry of the solid cluster was maintained (Fig. 3o, insets).
Structurally, an icosahedral twin has internal and surface strain 33 , thus an internal or external effect to compensate for the strain energy is required to stabilize the cluster. For instance, we know that hard sphere systems require spherical confinement to entropically stabilize an icosahedral cluster 4,7,23 , otherwise, the icosahedral cluster rapidly destabilizes (Fig. 4a). This indicates that the fluid-solid interfacial tension of hard spheres is not strong enough to overcome the entropy loss from the strain within the icosahedral cluster. On the other hand, the icosahedral cluster of hard TTs is stable in a fluid without geometric confinement, suggesting that the fluid-solid interfacial tension is enhanced by the shape of the particle. Figure 4b shows the crystal structure of the (111) plane in the clusters of hard TTs and hard spheres, which is the largest surface of the icosahedral cluster exposed to the fluid. The (111) surface of the hard TT crystal is flat because the faces of the TTs are exposed toward the plane. On the other hand, the (111) surface of the hard sphere crystal is rough due to the spherical shape of the particles, implying that the roughness of the surface can affect the surface tension of the crystal.
To quantify the difference in the surface effect, we calculated the fluid-solid interfacial free energy of hard TTs and hard spheres using the capillary fluctuation method 34 ("Methods"). We exposed a specific lattice plane of a crystal toward a fluid in a thin slab (z-axis in Fig. 4c, d) and monitored the fluctuation of the interface in equilibrium. The interfacial profile normalized by the circumscribed radius of the particle, ðh À h 0 Þ=r circ , shows that the fluctuations of the (111) plane of the hard sphere system are larger than that of the hard truncated tetrahedra system (Fig. 4e and Supplementary Movies 1, 2). From 200 samples of the interfacial profiles obtained in equilibrium, we calculated the interfacial stiffness (e γ) of the crystal planes (Supplementary Fig. 6 and "Methods") and measured the fluid-solid interfacial free energy ðγÞ of several major crystal orientations (Supplementary Table 2). The interfacial stiffness of the (111) plane of the fcc crystal of hard spheres, e γ HS ð111Þ = 0:684k B Tσ À2 , is 56% lower than that of the cubic diamond crystal of hard truncated tetrahedra, e γ TT ð111Þ = 1:547k B Tσ À2 ("Methods"). In addition, we confirmed that the fluid-solid interfacial free energy difference (surface tension) of the (111) plane of the hard TT crystal, γ TT ð111Þ = 1:347k B Tσ À2 , is ∼ 2:4 times higher than that of the hard sphere crystal, γ HS ð111Þ = 0:560k B Tσ À2 (Supplementary Table 2). This indicates that a strong fluid-solid interfacial tension of the hard TT system stabilizes the icosahedral twinned cluster in fluid without geometric confinement.
The dependence of the fluid-solid interfacial tension on particle shape is indirectly indicated by the change in local volume fraction distributions arising from TT shape change ( Supplementary Fig. 5). The volume fraction of fluid and crystal at fluid-solid coexistence depends on the fluid-solid surface tension; thus, the local volume fraction change due to TT truncation suggests a change in fluid-solid interfacial energy from the TT shape change. To confirm this hypothesis, we compared three different hard particle systems with very different local volume fraction distributions at fluid-solid coexistence ( Supplementary Fig. 7): (1) Fig. 7). Next, we conducted capillary fluctuation simulations for the three systems to calculate the interfacial stiffness of each system. We observed that the TBPs with the largest local volume fraction difference (4ϕ = 0:116) show the smallest fluctuations of the interface, resulting in the strongest interfacial stiffness (e γ ∼ 7:21k B Tσ À2 ) ( Supplementary Fig. 7b, c). On the other hand, the HSs with the smallest local volume fraction difference (4ϕ = 0:05) show the largest fluctuations of the interface, resulting in the weakest interfacial stiffness (e γ ∼ 0:74k B Tσ À2 ) (Supplementary Fig. 7h, i). This comparison provides clear evidence that particle shape determines the volume fraction of fluid and solid at coexistence, which in turn determines the interfacial energy.
In summary, we studied the entropy-driven assembly and stabilization of fivefold and icosahedral twinned clusters in a one-component fluid of hard truncated tetrahedra (TT). We showed that the thermodynamic stability of the cubic and hexagonal diamond phases can be entropically controlled by designing the TT shape. If the particle shape is designed to have a negligible free energy difference between the two diamond crystals, the formation of twin boundaries is easily induced. This strategy can be used to induce the formation of fivefold and icosahedral twinned clusters in fluid via seed-assisted growth. We found that the formation of a fivefold cluster follows an error-andrepair mechanism that removes mismatches of the stacking sequence of adjacent twin planes through particle rearrangement at the cluster surface. We showed the formation of an icosahedral cluster by additional growth of twinned structures from a fivefold twinned seed. We showed that the icosahedral cluster of TTs can be entropically stabilized in a fluid, which is not possible for hard spheres without confinement. The capillary fluctuation method showed that a hard TT system has a much higher fluid-solid interfacial free energy (here, entropy) difference than that of a hard sphere system, directly implicating fluid-solid interfacial tension in stabilizing the TT ico-twin.
Our findings provide a quantitative understanding of the formation and stabilization of fivefold and icosahedral twinned clusters in hard particle systems. Importantly, we showed that the twinning behavior and the interfacial free energy difference may be entropically engineered by particle shape design. Experimentally, colloidal tetrahedral particles are synthesizable 26,35-38 with tunable shape 37 . Non-complementary DNA could be used to make non-interacting particles that might closely approximate those studied here. Fivefold and icosahedral twinned clusters should be also attainable with attractive TT particle shapes, provided the explicit attraction favors face-to-face alignment. Such interparticle attraction could be realized through organic ligands 36,38 via van der Waals forces or self-complementary DNA 26 , both of which are commonly used to drive nanoparticle assembly.

Particle geometry
The truncated tetrahedron (TT) shape is a member of the spherical triangle invariant 323 family, with truncation parameters a,b,c ð Þ according to previous convention 39 . In this study, b = 1:0 for all systems, and a and b vary in a range: 1:14 ≤ a ≤ 1:30 and 2:10 ≤ c ≤ 2:30.

Identification of staggered and eclipsed pairs
Each TT in a crystal phase has four nearest neighbors that form face-toface contacts. We identified the type of each face-to-face contact (staggered or eclipsed), following a similar identification protocol described in ref. 28. Briefly, for each face of the particle, we assigned three vectors that are defined from the face center to each tip of the face (Supplementary Fig. 1a). Then, for each pair contact (i,j), we can define three vectors for i particle and three other vectors for j particle ( Supplementary Fig. 1b, c). We calculated every pair combination between the i and j vectors, which is nine, and found the minimum angle (θ). All the θ from every pair contact accumulated from the entire system give an angle distribution. For instance, Supplementary Fig. 1d was calculated from hexagonal diamond crystal. The angle distribution shows a clear bimodal distribution, and the height of the peak at θ 1 and θ 2 are around a 1:3 ratio, indicating that the angles around θ 1 and θ 2 come from eclipsed pairs and staggered pairs, respectively. This allows us to select the range of the pair contact angle to identify the eclipsed pair (θ 1 ± 10 ) and the staggered pair (θ 2 ± 10 ).

Monte Carlo simulation
Simulations were performed by the hard particle Monte Carlo (HPMC) implemented in the HOOMD-blue simulation package 40,41 , which is available at https://github.com/glotzerlab/hoomd-blue. The system size for the self-assembly simulations (Fig. 2) is N = 2000, and all simulations were performed under periodic box condition. The translational and rotational movement of particles is decided to have around a 20% acceptance ratio for every simulation. The self-assembly simulations were initialized from a dilute isotropic fluid ϕ = Nv 0 =V <0:01 and compressed until the desired thermodynamic condition (ϕ or reduced pressure P * = Pv 0 =k B T) was reached. Here, v 0 and V are the volume of a particle and the volume of the simulation box, respectively. The unit length of simulation is defined by σ, and the particle volume of a TT is set to v 0 = 1:0σ 3 regardless of its truncation amount. After initialization, each run was continued in the isochoric ensemble (NVT) at constant particle volume fraction ϕ until equilibration was reached. For the most systems, crystallization occurs within 1:0 × 10 8 Monte Carlo steps.

Stability of diamond crystals in a dense fluid
We checked the stability of a cubic diamond crystal in a dense fluid for three different systems: (1) a = 1:16, c = 2:16 at ϕ = 0:57, (2) a = 1:20, c = 2:16 at ϕ = 0:575 (Fig. 2d), and (3) a = 1:24, c = 2:16 at ϕ = 0:58 (Fig. 2e). All three systems are in a fluid-solid coexistence state, where a single-crystalline cluster is fully surrounded by a dense fluid phase. System size is N = 20,000 for all three systems. For the MC simulation, we first placed a spherical cubic diamond single-crystalline cluster (N ∼ 8000) at the center of a large simulation box ðϕ = 0:01Þ and put other particles around the solid cluster without overlapping each other. Then, holding the particles of the crystal phase (no movements), we rapidly compressed the system to the target particle volume fraction. Then, we released the particles of the crystal and equilibrated the whole system.

Free-energy calculation
The reduced Helmholtz free energy per particle F=Nk B T of the cubic and hexagonal diamond phases were calculated using the Frenkel-Ladd method 24,30 at a constant volume. The details of the calculation method can be found in ref. 30. Briefly, we first constructed an Einstein lattice (i.e., reference state) of each phase, for which the free energy is analytically solvable. Then, we applied strong translational and rotational harmonic potentials between the hard TTs and the Einstein lattice to tether each particle to its reference lattice site, which makes the hard TT system a reference state. Then we gradually released the harmonic potential until the strength of the potential became nearly zero, yielding the real system without correlation to the reference state. This process allows us to set up a continuous and reversible path between the reference state and a real state. Thus, by integrating the potential energy along the path, we can calculate the free energy (or equivalently, entropy) difference between the reference state and the real system. Because that the free energy of the reference crystal is analytically solvable, we can calculate the absolute free energy of the real crystal.
Using this method, we calculated the free energy of each phase with different shape parameters ða,cÞ to construct a free energy surface (Fig. 2a) in the shape space: 1:14 ≤ a ≤ 1:30 with 0.02 interval and 2:10 ≤ c ≤ 2:30 with 0.02 interval, for a total of 99 systems with different shape parameters ða,cÞ. The system size is N = 512 for the cubic diamond and N = 992 for the hexagonal diamond. For each phase, we calculated the free energy of the 99 systems at ϕ = 0:64, 0:66, 0:68, 0:70 and 0:72, and for each ϕ, the free energy surface was constructed by interpolating the 99 data points (Supplementary Fig. 4). The free energy of the crystal at ϕ = 0:62 (Fig. 2a) was estimated by extrapolating the free energy surface plots of the five particle volume fractions studied.

Seed-assisted growth MC simulation
We conducted seed-assisted growth simulation at a fluid-solid coexistence state (ϕ = 0:58 for a = 1:20, c = 2:16 system). For the growth of the fivefold twinned cluster, we used a spherical cluster of a cubic diamond crystal (N = 500) as a seed, and the total system size including the seed is N = 20,000. At a very dilute condition (ϕ = 0:01), we placed the seed at the center of simulation box and placed other particles around the seed without overlap. Fixing the seed particles in place, we compressed the simulation box to the target particle volume fraction. Once it reached the target particle volume fraction, we released the seed particles and equilibrated the system. For the growth of the icosahedral twinned cluster, we used a fivefold twinned cluster (N ∼ 1000) as a seed, and the total system size including the seed is N = 20,000. The compression and equilibration protocol is the same in both cases.

Local volume fraction calculation
Seed-assisted growth simulations for fivefold and icosahedral twinned clusters were performed in at coexistence between a fluid phase and a crystal phase, where the solid clusters are fully surrounded by fluid. To distinguish the two phases, we calculated the local particle volume fraction ϕ loc , which is defined as the particle volume fraction around a particle within a certain radius, r cut . The value of r cut in this study was taken to be 3 times the distance from a particle center to its nearest neighbor particle center, in order to properly average the local environment. TTs are identified as belonging to the fluid region if ϕ loc ≤ 0:60 and to the crystal region if ϕ loc >0:60. We utilized the freud python library for this calculation 42 .

Stability check of an icosahedral twinned cluster of hard spheres
We performed MC simulations of a hard sphere system to check the stability of an icosahedral twinned cluster of hard spheres at fluid-solid coexistence (ϕ = 0:515). We first constructed an icosahedral cluster (N = 5971) with an ideal structure, following the method described in ref. 4. At a very dilute condition (ϕ = 0:01), we placed the ico-twin cluster at the center of simulation box and placed other hard spheres around the cluster without overlap. Temporarily "freezing" the icotwin particles in place, we compressed the simulation box to the target particle volume fraction. Once it reached the target particle volume fraction, we released the ico-twin cluster particles and equilibrated the system. To check if the ico-twin cluster is stable in the coexistence phase, we monitored the bond-order diagram and the morphology of the ico-twin cluster during the simulation. We observed the icosahedral structure rapidly destabilize within 5 × 10 5 MC steps (Fig. 4a).

Capillary fluctuation method
The fluid-solid interfacial free energy of hard TT (a = 1:20, c = 2:16) and hard sphere systems was calculated by the capillary fluctuation method 34 . We prepared systems in a fluid-solid coexistence state (ϕ = 0:582 for hard TTs and ϕ = 0:515 for hard spheres) containing a crystal, where the crystal of each system exposes different orientations toward a fluid phase: ð111Þ½ 110, ð100Þ½001, ð100Þ½ 110 and ð110Þ½001. Here, the crystal orientation is defined as ðijkÞ½lmn, where ðijkÞ is a crystal plane toward the fluid phase (z-direction in Fig. 4c, d) and ½lmn is a crystal plane perpendicular to the fluid phase along the short direction (y-direction in Fig. 4c, d). The system size is N ∼ 40,000 for every system, and the shortest direction of the simulation box contains 2-3 unit cells (y-direction). At coexistence, we distinguished fluid and solid phases using the local volume fraction ϕ loc for the hard TT system (ϕ loc >0:6 for crystal phase) and the average bond-orientational order parameter 42,43 q 6 for the hard sphere system ( q 6 >0:3 for crystal phase).
For the calculation of the fluid-solid interfacial free energy, we followed the calculation process described in ref. 34. Briefly, we obtained the interfacial profile hðxÞ as a function of MC steps ( Fig. 4e-g), which is then normalized by the circumscribed sphere radius of a particle. The values of hðxÞ were measured at a discrete set of points x n = nΔ, where n = 1, . . . ,N and Δ = L=N, where L is the length of crystal along x-direction. The Fourier modes h q are defined as h q = 1 N P N n = 1 hðx n Þe iqx n , where wave number q = 2πk L with k = 1, . . . ,N. From the equipartition theorem, the size of the capillary fluctuation modes is related to the interfacial stiffness e γ as follows 44 : where A is the area of the fluid-crystal interface. The e γ for different crystal orientations were calculated and averaged for various q and Δx ( Supplementary Fig. 6). The interfacial stiffness is related to the interfacial free energy γ by the formula: where θ is the angle between the average orientation of the interface and the local normal of the interface. Based on this relationship, the interfacial free energy can be obtained by comparing the following two equations: γ n ð Þ=γ 0 = 1 + A 1 ϵ 1 + A 2 ϵ 2 ð3Þ e γ n ð Þ=γ 0 = 1 where n is the unit vector normal to the interfacial plane. The values of A 1 , A 2 , B 1 and B 2 for each interface direction is listed in Supplementary  Table 1. Then, the interfacial free energy γ for each crystal plane was Article https://doi.org/10.1038/s41467-022-34891-5 obtained using the coefficients (Supplementary Table 1) and the interfacial stiffness ( Supplementary Fig. 6). The calculation results are listed in Supplementary Table 2.

Data availability
All data needed to evaluate the conclusions in this manuscript are present in the main text or the supplementary materials. Additional data are available upon request to the corresponding author.