Main

Clathrates are found in gas hydrate systems, where natural gas guest molecules are trapped inside ice cages1,2,3,4, or in intermetallics5,6, where metal ion guests are encapsulated by cages of group 14 elements such as Si, Ge and Sn. The length scale of the cages in these traditional clathrates is typically 1 or 2 nm, but the recent discovery of clathrate structures self-assembled from colloidal particles7,8 creates the possibility of larger cages and concomitant mesoscale phenomena relevant to biological or photonic applications. One of the most interesting properties of the molecular clathrate is the host–guest relationship, which is often used for trapping natural gas2 or storing hydrogen9, and the guest dynamics has been extensively studied3,4,9,10,11,12,13,14,15 to understand its importance in relation to the thermophysical properties and stability of the clathrate. Extending the host–guest chemistry of molecular clathrates to colloidal ones might lead to practical applications such as nanoparticle storage and the scaffolding of nanoparticle guests into a desired lattice. So far, however, none of the reported colloidal clathrates have been able to accommodate guests because of their closed cages, raising the question of whether it is possible to self-assemble colloidal clathrates with open cages of sufficient size to accommodate guest particles or macromolecules.

In this Article we report the self-assembly of seven host–guest colloidal clathrates with unit cells ranging from 84 to 364 particles in Monte Carlo (MC) simulations of hard, otherwise non-interacting truncated triangular bipyramid (TBP) shapes, in which entropy maximization alone drives assembly of the ordered crystal from the disordered fluid phase16,17. By design, truncating the TBP tips creates a cavity inside the cages that is large enough to accommodate a guest (Fig. 1), and the guest/cavity size ratio as determined by the amount of TBP truncation dictates the resulting clathrate crystal (Fig. 2), similar to that reported in gas hydrate systems2,18. It is noteworthy that, in every case, hosts and guests form two sublattices with distinctly different particle dynamics and with an unusual compartmentalization of the entropy19 such that the sublattice of guest particles contains a substantially higher entropy per particle (Fig. 3a–c). As we will show, this entropic compartmentalization thermodynamically stabilizes the equilibrium assemblies (Figs. 3d,e and 4). Five of the seven structures reported are one-component systems, in which the host and guest particles are the same species. For one of these five systems, we also introduce two different guests, with different particle shapes, to obtain traditional (binary) host–guest structures in addition to the one-component system (Fig. 2f,g). We compare these entropy-driven colloidal clathrate crystals with molecular clathrates (Table 1), in which guest molecule rotation is well known. We conclude by using entropic bonding theory (EBT) to predict the effective attraction between truncated TBPs needed to self-assemble a clathrate structure and show via simulation how to achieve such a system using DNA-functionalized nanoparticles.

Table 1 Hard colloidal TBP clathrates versus clathrate hydrates
Full size table

Results

Particle design for host–guest hard truncated TBP clathrates

In clathrate hydrates2, the size ratio of guest to cage determines the clathrate structure. Inspired by this knowledge, we engineered the shape of truncated TBP particles to tune guest/cage size ratios. Previous work8 has shown that, for regular hard TBPs, the cage interiors are almost fully filled with the body of the host particles such that there is no room for a guest particle (Fig. 1b). However, a cavity can be made by truncating the equatorial tips of the TBP (Fig. 1a,b and Extended Data Fig. 1). The truncation parameter S (0 ≤ S ≤ 1) is defined by the distance between the tip of the regular TBP (S = 0) and the centre of the truncated face (Extended Data Fig. 1). We measured the average radius of the cavity (rcavity) in the four types of cage possible for a clathrate crystal (A, 512; B, 51262; C, 51263; D, 51264, following the usual clathrate nomenclature8) from ideally constructed clathrates thermalized in hard-particle Monte Carlo (HPMC20) simulations at equilibrium (Extended Data Fig. 2 and Monte Carlo simulation section in the Methods). HPMC attempts small local particle moves (both translational and rotational) while avoiding particle overlap to mimic the Brownian motion expected for colloidal particles in water or other solvent (see the Monte Carlo simulation section in the Methods and Supplementary Table 1 for simulation details). The range of S over which a cage has a cavity large enough to contain a guest can be estimated from the size ratio α between guest and cavity (α ≡ rguest/rcavity). The size ratio α decreases monotonically with increasing S for all four cages (Fig. 1c) over the entire S range, with relative cavity sizes D > C > B > A for any given value of S. This observation indicates that the type of cage that can house a guest particle varies with S. Because the stability of any clathrate structure requires a stable combination of cage types, we hypothesized that the clathrate type obtained in self-assembly simulations would vary with S. For example, if S is large enough for the D-cage—but not other cages—to accommodate a guest, we might expect the self-assembly of a clathrate type II structure (Clath II–A0D1) with guest-free A-cages and single-guest D-cages (Fig. 1d). The A-cages are empty because the size of the empty pore is smaller than the host particle. We show in a later section that these empty cages can be filled with particles of arbitrary shape during self-assembly, provided that the size is small enough to fit inside the pore.

Fig. 1: Particle design for host–guest colloidal clathrates.
figure 1

a, TBP with truncation parameter, from left to right, of S = 0, 0.3, 0.5, 0.7 and 1.0. b, Thirty truncated TBPs form a cage-type clathrate cluster. The truncation creates a cavity (red dotted circle) at the centre of the cluster. The size of the cavity increases as S increases. c, The change of size ratio α = rguest/rcavity as a function of S for the four clathrate cages measured at a constant volume fraction, ϕ = 0.65 (Extended Data Fig. 2 and Guest-to-cavity size ratio section in the Methods). Data points are obtained from HPMC simulations; solid lines are guides to the eye. The location where the red dotted line meets each curve indicates the minimum S for the cage to have a single truncated TBP guest. d, A model structure of Clath II–A0D1 (S = 0.42), composed of guest-free A-cages and single-guest D-cages.

Source data

Full size image

Self-assembly simulations

To test our hypothesis, we conducted self-assembly HPMC simulations of hard truncated TBPs with varying truncation S. For 0 ≤ S < 0.40, a guest-free Clath I–A0B0 formed via self-assembly from a disordered fluid phase; this is the same type of clathrate found in the hard regular TBP system8, but now with an empty cavity inside the cages that increases in size with increasing S while remaining smaller than the host particle (Fig. 1c and Supplementary Fig. 1a). Above S = 0.40, we observed only host–guest clathrates. For 0.40 ≤ S ≤ 0.71, five different types of host–guest colloidal clathrate crystal self-assembled (Fig. 2a–e). The type of clathrate was confirmed by the cage identification method and the arrangement of the cages (Clathrate identification section in the Methods). For 0.40 ≤ S < 0.46, we obtained Clath II–A0D1 with empty (that is, guest-free) A-cages and single-guest D-cages (Fig. 2a). In this structure, the A-cages arrange into rhombs when viewed along the two-fold axis and, in between, the D-cages are tetrahedrally arranged. For 0.46 ≤ S < 0.49, we obtained Clath III–A0B1C1 with empty, guest-free A-cages, single-guest B-cages and single-guest C-cages (Fig. 2b). In this structure, the B-cages arrange in square and triangular tiles, which can be mapped to the Frank–Kasper σ phase21. For 0.49 ≤ S < 0.55, we observed the assembly of Clath I–A0B1 with empty, guest-free A-cages and single-guest B-cages (Fig. 2c, Supplementary Fig. 1b and Supplementary Video 1). The B-cages are arranged in square tiles, which is a unique feature of the Clath I structure. For 0.55 ≤ S < 0.57, Clath I–A1B1 was obtained, with both A-cages and B-cages occupied by single guests. Finally, for 0.57 ≤ S < 0.71, we obtained a clathrate with single-guest A-cages and tetramer-guest D-cages (Fig. 2d). Each tetramer is composed of four truncated TBPs that form a cluster with tetrahedral symmetry, reminiscent of the tetrahedral cluster of hydrogen comprising the guest in the D-cage of a clathrate hydrate system9. Here, the A-cages arrange into a rhombus, corresponding to Clath II–A1D4. For S ≥ 0.72, crystallization was not observed during the simulation time. Our findings are summarized in the phase diagrams in Fig. 2e and Supplementary Fig. 2.

Fig. 2: Self-assembly simulation results.
figure 2

a–e, MC simulation results for Clath II–A0D1 at S = 0.42 (a), Clath III–A0B1C1 at S = 0.48 (b), Clath I–A0B1 at S = 0.52 (c) and Clath II–A1D4 at S = 0.64 (d). Red and green truncated TBPs are guest particles, and blue truncated TBPs are host particles. Each clathrate can be described along an axis by a unique tiling of polygons (orange lines and green/yellow dashed circles) and a bond orientational order diagram (BOOD; insets). The BOODs are a histogram of the directions of bonds connecting nearest neighbours. Bottom left: network representation of each clathrate (Network representation of clathrate section in the Methods) for clear visualization of host particle cages. Bottom right: number ratio of cages identified in each clathrate (Clathrate identification section in the Methods). e, Phase diagram as a function of S obtained from self-assembly simulation results. Coloured dots are data points, and the four systems shown in a–d are indicated by red arrows. f,g, Self-assembly simulation results of truncated TBP clathrate with other polyhedron guests. f, A binary mixture of truncated TBPs (S = 0.52) and tetrahedra (αA ≈ 0.89) forms Clath 1–A1B1. The A-cage has a tetrahedron guest and the B-cage has a truncated TBP guest. g, A binary mixture of truncated TBPs (S = 0.52) and dodecahedra (αA = 0.97) forms Clath I–A1B1. The A-cage has a dodecahedron guest and the D-cage has a tetramer truncated TBP guest.

Source data

Full size image

From the self-assembly results, we find a critical value of α below which a truncated TBP guest enters the cavities (αc ≈ 1.2, the red dotted line in Fig. 1c), corresponding to Sc = 0.40, 0.46, 0.48 and 0.58 for cage types A, B, C and D, respectively. Thus, as hypothesized, specifying the value of S determines the type of cage that can house a guest particle, which in turn determines the type of clathrate structure that self-assembles. Because there is no potential energy in any of these systems, the clathrate structures that minimize free energy maximize entropy. Although the diversity and complexity of entropic colloidal crystals is well established16,17, the way in which these clathrates maximize entropy, which we discuss next, is not.

Entropy compartmentalization

In an isotropic fluid phase of a single-component system of hard particles, the free volume is distributed uniformly (Fig. 3a). In the clathrate structures, the free volume is not uniformly distributed among hosts and guests. For example, in the Clath I–A0B1 crystal (S = 0.56), both host and guest particles in equilibrium exhibit similar translational displacement (quantified by the mean-squared displacement \({{{{\bf{r}}}} - {{{\bf{r}}}}_0^2}\), where r is the centre-of-mass position) with comparable fluctuations about the mean, over 1 × 107 MC sweeps (Fig. 3b, upper panel), indicating a similar free volume in the two sublattices in which translational motion can occur. In contrast, guest particles exhibit much larger rotational displacements (\({{{{\bf{q}}}} - {{{\bf{q}}}}_0^2}\), where q is a unit quaternion describing orientation) than host particles (Fig. 3b, lower panel and Supplementary Video 2), indicating unequal free volume partitioning. We observe the same behaviour in all seven colloidal clathrates. Thus, the many microstates available to these self-assembled crystal phases arise from the many rotational degrees of freedom available to the guest particles. In this way, we say that the entropy is compartmentalized, and can be almost entirely attributed to the sublattice(s) of guest particles. When multiple cage types contain guests, the rotational degrees of freedom of the guests in all cages contribute to the system entropy.

Fig. 3: Entropy compartmentalization and free energy calculation.
figure 3

a, A schematic of entropy compartmentalization in a single-component system. b, Translational (upper) and rotational (bottom) displacement of the host and the guest truncated TBPs of Clath I–A0B1 (S = 0.56) at P* = 13 for 1 × 107 MC sweeps. c, Left: Helmholtz free energy of the S = 0.56 system at ϕ = 0.60 in fluid and crystal (Clath I–A0B1) phases. Right: Helmholtz free energy of host and guest sublattices in the crystal. d, Gibbs free energy of clathrates at P* = 10.0 in 0.35 ≤ S ≤ 0.75 systems. e, Phase diagram of truncated TBP clathrates (P* versus S). The melting pressure is obtained from where the free energy of the clathrate crystal crosses the free energy of the fluid as the pressure increases (Extended Data Fig. 5). The colour key is the same as in d.

Source data

Full size image

Entropy compartmentalization is confirmed in these clathrate structures by free energy calculations. Based on the Frenkel–Ladd method22,23, we developed a way to calculate the free energy of the two subsystems of host and guest particles independently by controlling each subsystem with different paths (Extended Data Figs. 3 and 4 and Free-energy calculation for the host and guest sublattices section in the Methods). For example, in the Clath I–A0B1 system (ϕ = 0.60 and S = 0.56), the per-particle Helmholtz free energy F of the guest subsystem (Fguest/Nguest ≈ 0.72kBT) has only 1/8 the per-particle free energy of the host subsystem (Fhost/Nhost ≈ 5.65kBT). Because there are far fewer guest particles than host particles and no potential energy, the guest subsystem has substantially lower free energy and thus much higher entropy than the host subsystem. The Clath I–A0B1 system has 92 host particles and 6 guest particles per unit cell (Supplementary Table 1), and thus Fguest ≈ 4.32kBT and Fhost ≈ 519.8kBT, a difference of over two orders of magnitude.

We compared the bulk free energies of the different clathrate phases to confirm the relative thermodynamic stability among the various clathrate types at a given S. We calculated the Gibbs free energy G of seven types of clathrate (six of which self-assembled and one, Clath IV–A0B1C1, did not and is metastable in all S ranges) for 0.35 ≤ S ≤ 0.71 at constant pressure (P* = 11.0) (Fig. 3d and Free-energy calculation of bulk phases section in the Methods). The calculation accurately predicts the range of S over which each of the six self-assembled clathrate crystals is stable (Fig. 3d). Melting pressure (\({P_{\rm{m}}^ \ast}\)) is another indicator to estimate the stability of a crystal in hard particle systems. \({P_{\rm{m}}^ \ast}\) can be determined from where the free energy of the fluid and crystal cross (Extended Data Fig. 5a), and the free energy of the fluid can be calculated using thermodynamic integration (Free-energy calculation of bulk phases section in the Methods). We confirmed that the values of \({P_{\rm{m}}^ \ast}\) obtained from the free-energy calculations are consistent with those we find in the simulations (Extended Data Fig. 5b). Based on the \({P_{\rm{m}}^ \ast}\) of each clathrate phase, we constructed a phase diagram (P* versus S) for 0.35 ≤ S ≤ 0.71 (Fig. 3e). We also tested the importance of the guest for the stability of the host clathrate framework by removing the guest particles from a previously self-assembled Clath I–A0B1 (S = 0.52). Upon further equilibration, the host network reorganized to sacrifice the required number of host particles to fill the empty cages, forming a new, complete host–guest structure (Extended Data Fig. 6). This indicates that the presence of the guest particle is essential to stabilize the host framework.

Three types of rotational dynamics of guest particles

The rotational dynamics of the guest particles varies depending on the relationship between the shape and symmetry of the guest particle and cavity. To investigate the rotational dynamics of the guests in each of the colloidal clathrates, we tracked the motion of the guests in equilibrium. We assigned a direction vector of unit length to the guest particle (Fig. 4a,f), tracked its position for 3 × 107 MC sweeps, and mapped the accumulated positions of the direction vector onto a unit sphere. We identified three distinct types of rotational motion of the guest particles: (1) free rotation (Fig. 4b), (2) rotation around an axis (Fig. 4d) and (3) discrete rotation (Fig. 4g). For the single truncated TBP guest of the D-cage in the Clath II–A0D1 phase (S = 0.42), we assigned a direction vector from the centre of the particle to a tip (Fig. 4a) and observed that the accumulated positions of the direction vector show a random distribution (Fig. 4b), indicating the free rotation of the guest (Supplementary Video 3). For the single truncated TBP guest of the B-cage in the Clath I–A0B1 phase (S = 0.52), the accumulated positions of the direction vector distribute along a continuous band (Fig. 4d), indicating that the guest is rotating around an axis (Supplementary Video 2). For the tetramer truncated TBP guests of the D-cage in the Clath II–A1D4 phase (S = 0.64), we assigned direction vectors from the tetramer guest centre to the centre of each truncated TBP (Fig. 4f) and observed that the accumulated positions form several clusters (Fig. 4g), indicating that the four truncated TBPs comprising the tetramer move as one, and with preferred orientations. Each cluster has a mixture of the four direction vectors (Fig. 4g, right), indicating that the tetramer guest rotates through all four discrete directions (Supplementary Video 4). We refer to this type of rotational motion as discrete rotation.

Fig. 4: Rotational motion of guests in clathrate cages.
figure 4

a, For a single truncated TBP guest, a direction vector v with unit length σ from the particle centre to the centre of a truncated face quantifies the rotational motion of the guest. b,c, Accumulated points of v of a guest particle in a D-cage of a Clath II–A0D1 crystal phase (S = 0.42) at constant pressure (P* = 10.0) from MC simulation (b) and EBT calculation (c), represented on a unit sphere (left) and spherical coordinates (r = 1.0σ) (right). d,e, Accumulated points of v of a guest particle in a B-cage of a Clath I–A0B1 crystal phase (S = 0.52) at P* = 12.5 from MC simulation (d) and the lower energy distribution as computed from EBT (e), in each case represented on a unit sphere (left) and spherical coordinates (r = 1.0σ) (right). f, For a truncated TBP tetramer guest, four direction vectors (v1, v2, v3 and v4) with unit length σ were assigned from the tetramer centre to the centre of each truncated TBP. g, Accumulated points of the direction vectors of a tetramer guest in a D-cage of a Clath II–A1D4 crystal phase (S = 0.64) at constant pressure P* = 15.0 on a unit sphere (left) and spherical coordinates (r = 1.0σ) (right). h, Distribution of prolateness Pr of the cavity of the B-cage of Clath I–A0B1 (the same system as in d) and the D-cage of Clath II–A0D1 (the same system as in b). The negative Pr of the B-cage indicates an oblate shape of the cavity, and the Pr ≈ 0 of the D-cage indicates a spherical shape of the cavity. i, Both the D-cage and the tetramer guest in Clath II–A1D4 (the same system as in g) have tetrahedral symmetry, inducing the discrete rotation of the guest. j–l, Schematic representation (left) and rotational motion (right) of the tetrahedron in cage A (free rotation) (j), the hexagonal prism in cage B (rotation around an axis) (k) and the icosahedron in cage A (discrete rotation) (l). The three cages are composed of S = 0.52 truncated TBPs.

Source data

Full size image

The shape of a guest relative to that of the cavity that contains it correlates with the three types of rotational motion. We characterized cavity shape using the prolateness Pr (Dynamics of guest particles section in Methods) and confirmed that the cavity of the D-cage is almost spherical (Pr ≈ 0), while the cavity of the B-cage is oblate (Pr ≈ −0.02) (Fig. 4h). Combined with the above findings, the spherical cavity of the D-cage does not hinder the single truncated TBP guest (~0.4 aspect ratio), resulting in free rotation (Fig. 4h, right inset). In contrast, continuous rotation around an axis and discrete rotation occur when the guest shares symmetry elements with the cavity. For example, the oblate cavity of the B-cage (symmetry group D6d) allows the oblate guest (symmetry group D3h) to align along, and rotate around, the same symmetry axis (Fig. 4h, left inset). In the case of discrete rotation, both the truncated TBP tetramer guest and the D-cage cavity have tetrahedral (T) symmetry (Fig. 4i), producing four preferred, discrete orientations within the cavity. Similar types of rotational behaviour were reported for two-dimensional (2D) colloidal crystals of polygons when the particle symmetry matches the symmetry of its local environment24.

To elucidate the correlation between guest and cage shape/symmetry, we used EBT25 to calculate the ground-state free energy for each individual orientation of a guest particle (Entropic bond calculations section in the Methods). We first validated the applicability of EBT to the host–guest clathrate systems by constructing a phase diagram for the simulated structures, finding good agreement with both simulations and Frenkel–Ladd free energy calculations (Extended Data Fig. 7). We then used EBT to compute the free energy of a guest-occupied D-cage in Clath II–A0D1 and B-cage in Clath I–A0B1 at S = 0.42 and 0.52, respectively, for all possible guest orientations. Histograms of the computed free energies from EBT show two distinct distributions for the B-cage (Extended Data Fig. 8b), but only one continuous distribution for the D-cage (Extended Data Fig. 8a). Employing the same definition of directional vector in Fig. 4a, mapping the lower free energy distribution for the B-cage to a unit sphere predicts the same rotational pattern about an axis observed in the simulations (Fig. 4e). EBT predicts free rotation for the case of the D-cage (Fig. 4c), in agreement with simulation. The excellent agreement between EBT predictions, the observed rotational dynamics and the corresponding clathrate structures reveals that thermodynamics not only drives the formation of the clathrate structures but also controls the allowed/disallowed guest dynamics for each respective clathrate cage.

The above results suggest that one can a priori select the shape of the guest particles relative to the cage shape to target specific rotational patterns, which is potentially useful for the design of stimuli-responsive sensors or signal-amplifying/propagating mesostructures. To test this idea, we replaced the truncated TBP guests with tetrahedral, hexagonal prism and icosahedral guests in S = 0.52 clathrate cages (Fig. 4j–l) and observed their rotational motion using HPMC simulation. As hypothesized, the tetrahedron guest (αA ≈ 0.89) freely rotates in the A-cage due to the symmetry mismatch between the guest (tetrahedral symmetry) and the cage (icosahedral symmetry); the hexagonal prism (αB ≈ 0.95) rotates continuously around an axis in the B-cage because both guest and cage have an oblate shape; and the icosahedron guest (αA ≈ 1.01) rotates among several discrete orientations in the A-cage due to the symmetry match between the guest (icosahedral symmetry) and the cage (icosahedral symmetry). When the guest particles are sufficiently small relative to the cage size, the shape and symmetry of cage versus particle is irrelevant and all guests rotate freely. This shape control over the rotational dynamics of the guest particles could be especially powerful in the self-assembly of host–guest clathrates from binary mixtures of shapes, as discussed next.

Self-assembly of binary host–guest colloidal clathrates

The cavities of the clathrate cages can be filled by particles different from that of the host during self-assembly, as confirmed by HPMC simulations of two binary hard-particle systems: (1) truncated TBPs with tetrahedra and (2) truncated TBPs with dodecahedra (Fig. 2f,g, Supplementary Fig. 3 and Supplementary Table 2). We chose the size of the second particle type to fit into the cavity of the A-cage in Clath I in both cases, but selected different symmetries to induce different rotational dynamics. For the first system, a fluid mixture of truncated TBPs (S = 0.52) and tetrahedra (αA ≈ 0.89) self-assembled into Clath I–A1B1 at ϕ = 0.55, with the A-cages now occupied by a single tetrahedron guest (rather than remaining empty or being filled by a host particle) and the B-cages occupied by a single truncated TBP guest (Fig. 2f). Similarly, the second system (S = 0.52) formed Clath I–A1B1 at ϕ = 0.55, with the A-cages containing a single dodecahedron guest (αA = 0.97) and the B-cages containing a single truncated TBP guest (Fig. 2g). As expected, the guest and host particles again show different rotational dynamics and entropy compartmentalization, and the second particle types with different symmetries show different rotational dynamics. Tetrahedron guests in the truncated TBP–tetrahedron system rotate freely in the A-cages (Extended Data Fig. 10a) due to the symmetry mismatch with the cage, whereas the dodecahedron guests in the truncated TBP-dodecahedron system rotate among several discrete orientations in the A-cages (Extended Data Fig. 10d) due to the symmetry match with the cage. For both systems, the single truncated TBP guest in the B-cages rotates continuously around an axis (both the truncated TBP and B-cage have an oblate shape; Extended Data Fig. 10b,d); at the same time, the host truncated TBP particles exhibit only local vibration (Extended Data Fig. 10c), as in all of the one-component self-assembled colloidal clathrates. The multiple types of rotational motion identified for the different species of guest particles in the different cages of the binary clathrates indicate that the entropy of the guest subsystem may be even further compartmentalized among different guest sublattices, by design.

Host–guest colloidal clathrate with attractive polyhedra

We anticipate that these and other open-caged colloidal clathrate crystals can be realized in the laboratory in several ways, even when purely hard particles are not available or desired. EBT calculations indicate that the key ingredient is a net attraction of at least 3 kT between the host particles to ‘freeze’ them into the cage network (Extended Data Fig. 9). In the hard-particle systems studied here, this attraction is emergent and implicit—a statistical, effective force of attraction driven by entropy maximization. However, an attractive interaction energy between particles can be explicitly encoded through the use of attractive ligands26. One way to achieve this would be to use DNA-mediated interactions to explicitly program a minimum 3-kT attraction between the facets of the host truncated TBPs. For example, mixing a batch of non-attractive, truncated TBPs (or other suitable shapes) with truncated TBPs functionalized with self-complementary DNA should produce a relatively rigid host clathrate network composed of DNA-linked particles, and a guest sublattice of unlinked particles that rotate in the cages of the network (to make them non-attractive, the guest particles could be functionalized with non-complementary DNA). Thermal annealing about the DNA melting temperature should produce net DNA base-pairing interactions between truncated TBP faces on the order of kT. Once this thermal energy scale is achieved, entropy maximization will dominate and drive the particles functionalized with non-complementary DNA to occupy the interior of the clathrate cages. A similar scheme should work for other types of patchy particle with attractive interactions between host particles and no attraction between guests or between guests and hosts.

We tested these possibilities using two different molecular dynamics (MD) simulation models (Methods): a mean-field, pairwise interaction potential using the anisotropic Lennard–Jones potential (ALJ)27 and a patchy particle model of DNA-functionalized colloidal particles7,28. Both models have successfully predicted nanoparticle assemblies7,28,29. In the ALJ model, we checked that a B-cage (S = 0.52) of attractive truncated TBPs with a purely repulsive guest TBP self-assembles at ~5-kT attraction between host particles (Supplementary Video 5). Here, the choice of 5 kT as opposed to 3 kT was chosen because EBT does not consider thermal fluctuations present in the simulations. That is, 3 kT is the minimum attractive strength in the absence of fluctuations. By selecting 5 kT we ensure that the attraction is sufficiently strong to overcome thermal fluctuations. For the patchy particle model, we constructed a Clath I–A0B1 (N = 2,315, S = 0.52) system with attractive host truncated TBPs and repulsive guest truncated TBPs. We confirmed that the clathrate is stable at \({T^ \ast /T_{\rm{m}}^ \ast \approx 0.90}\) as shown in the potential energy of the system over time (Supplementary Fig. 4). The fact that the clathrate structure can be stabilized with attractive host particles suggests that the design approach used here to select for different clathrate structures using particle shape should be useful in guiding experiments, regardless of the source of interaction (and whether explicit or implicit) between particles.

Discussion

Our study shows entropy compartmentalization in both single- and binary-component, hard-particle host–guest clathrate crystals that self-assemble from a disordered fluid phase. We observe purely entropy-driven self-assembly of five different single-component host–guest colloidal clathrate crystals by simply tuning the amount of truncation of the TBPs. We also observe purely entropy-driven self-assembly of two additional clathrate phases with some cages occupied by tetrahedra or dodecahedra, and others occupied by truncated TBPs. In all cases, the host truncated TBP particles assemble into a tightly packed, open clathrate network structure and are largely immobile, while the guest particles exhibit three different types of rotational dynamics: free rotation, continuous rotation around an axis and discrete rotation. Because rotation of the guest particles is the primary way of sampling microstates in equilibrium, we say that the open clathrate structures are stabilized by the rotational entropy of the subsystem (or subsystems) of guests. That entropy can stabilize an open colloidal crystal is not new; Mao et al. showed how an open kagome lattice of attractive triblock Janus spheres is stabilized over a close-packed structure by the additional vibrational entropy afforded to the spheres by the open structure30. In our clathrate system, the absence of explicit attractive interactions, and compartmentalization of entropy not by degrees of freedom, as in ref. 30, but by particle subsystem, makes such entropic stabilization unusual. Importantly, our description of the nature of the entropic stabilization for hard-particle colloidal clathrates is not unlike that used for molecular (for example, hydrate) clathrates3,4,12,14, where the host network of water molecules is ice-like and the guests rotate either freely or with preferred orientations (Table 1). In molecular clathrates, non-free rotation of guests is attributed to anisotropy in the electrostatic fields arising within the cages3,4,10,11,12,14,31. In colloidal clathrates, the anisotropy necessary for non-free guest motion arises simply from mismatches between the symmetries of the particle and the cage.

Unlike molecular clathrates, the pores in colloidal clathrates are not connected; thus, these colloidal crystals do not lend themselves to applications requiring the transport of guest particles within the crystal. However, they are relevant for applications involving particle storage, particle capture/release, and scaffolding (positioning) of guest particles into a desired lattice. For example, if the size of the open pores of the clathrates is designed to be tens of nanometres, the structure can be used to capture biomaterials (for example, proteins and viruses) by judiciously functionalizing the host particles. As another example, note that the open pores in Clath I–A0B1 are arranged as a body-centred cubic (bcc) lattice, so one can arrange arbitrary particles into a bcc lattice (a lattice they might never form on their own) by filling them into the open pores.

As has been demonstrated now in many articles in the literature, the ability of entropy alone to produce thermodynamically stable crystals that are isostructural to atomic and molecular crystals seems not to be limited by crystal complexity, and so perhaps we should not be surprised that purely entropic, colloidal host–guest clathrates are also possible. However, that entropy alone can produce the key ingredient of the molecular clathrate (rotating guests inside a relatively ‘frozen’ cage network) and that this rotational ability is critical to the stability of the crystal is profound. Our findings point at new solutions dictated by thermodynamics that leverage entropy maximization to drive the formation of complex structures. In the systems reported here, sublattices form in such a way that the contribution to the entropy difference from translational degrees of freedom is negligible compared to that from rotational degrees of freedom, which in turn are apportioned solely to the members of one of the sublattices. Thus, although thermodynamic entropy is a global system property and undefined on the local level, we can loosely say that the system compartmentalizes its entropy into the sublattices. Moreover, we see that multiple guest sublattices with different rotational dynamics and thus different entropies are possible through further compartmentalization dictated by the guest/cavity size ratio and symmetries of the various types of guest and cage contained within a single clathrate crystal. Once again, we see that the ‘solution’ found by thermodynamics does not depend on the origin of the forces in the system, and that, in principle, structures equivalent to those of a great many atomic and molecular crystals can be achieved with colloidal nanoparticles. Our findings also further support the long-held opinion that colloidal systems can be useful models for atomic and molecular systems by permitting the study of processes that may be hard to see on atomic and molecular scales.

Methods

Monte Carlo simulation

Simulations were performed using the HPMC20 algorithm implemented in the HOOMD-blue simulation package37,38, which is available at https://github.com/glotzerlab/hoomd-blue. The system size for all self-assembly simulations (Fig. 2) is either N = 4,000 or N = 10,000, and all simulations were performed with periodic boundary conditions. HOOMD-blue’s HPMC module logs the positions and orientations of particles as xyz coordinates and quaternions, respectively; quaternions are a simple and fast metric for post-processing of the data. The self-assembly simulations were initialized from a dilute isotropic fluid phase with volume fraction ϕ = Nv0/V < 0.1 and compressed until the desired thermodynamic target (ϕ or reduced pressure P* = Pv0/kBT) was reached. Here, v0 and V are the volume of a particle and the volume of the simulation box, respectively. The unit length is defined by σ, and the volume of a truncated TBP is set to v0 = 1.0σ3, regardless of truncation amount S. Following compression, each run was continued in the NVT ensemble at constant volume fraction ϕ or in the NPT ensemble at constant pressure P* until equilibration was reached, as determined by the appearance of a plateau in the time evolution of pressure or volume, respectively. Our experience has taught us that hard polyhedron systems typically crystallize (if they are going to crystallize) within the range of particle volume fractions 0.45 ≤ ϕ ≤ 0.65 in the NVT ensemble or within the range of reduced pressures 8 ≤ P* ≤ 20 in the NPT ensemble. For the present systems, we thus scanned those parameter ranges to find the values of ϕ and P* where colloidal clathrates formed. For most systems, crystallization occurred within 1.0 × 108 MC sweeps. As observed for the previously reported hard regular TBP colloidal clathrate system8, in most cases, the equilibrium state was that of a crystal in coexistence with a fluid phase. The simulation parameters for Fig. 2 are shown in Supplementary Table 1.

Guest-to-cavity size ratio

The guest-to-cavity size ratio α is defined as rguest/rcavity, where rguest is the radius of a circumscribed sphere of the guest particle. We occasionally use αA—the value of α for an A-cage (that is, rcavity is that of the A-cage). rcavity was obtained by the distribution of the distances from the centre of the cavity to each vertex of the cavity (Extended Data Fig. 2). To obtain the cavity size distribution of the four different clathrate cages in 0.35 ≤ S ≤ 0.71, we constructed ideal clathrate structures for Clath I–A0B0, Clath II–A0D1, Clath III–A0B1C1, Clath I–A0B1 and Clath II–A1D4, and equilibrated those structures before measuring the sizes of their cavities. After measuring the cavity size of the four clathrate cages, we plotted α versus S (Fig. 1c). Note that, if the shapes of the guest and the cavity are both aspherical, a guest can potentially fit inside a cavity even though rguest is slightly larger than rcavity because rcavity is an averaged value.

Clathrate identification

The type of clathrate can be identified by its cages and their arrangement2. We classified the cages of the truncated TBP clathrate structure into four types7 based on the number of truncated TBPs comprising the cage (A-cage, 30 truncated TBPs; B-cage, 36 truncated TBPs; C-cage, 39 truncated TBPs; D-cage, 42 truncated TBPs); these are mappable to the 512, 51262, 51263 and 51264 cages of molecular clathrates, respectively. The notation 5n6m indicates that a cage has n pentagons and m hexagons as faces. To detect the number and type of cages in the clathrate crystal, we converted the truncated TBPs into their polar tips, because the particles point their polar tips towards the cage centre. We then detected clusters of the polar tips within a cutoff and counted the number of polar tips in each cluster. The cluster sizes 30, 36, 39 and 42 then correspond to cages A, B, C and D, respectively. Each type of clathrate has its own unique arrangement of cages, which can be used to classify the type of clathrate. For clathrate type 1, the B-cages are arranged in square tiles along the four-fold symmetry axis. For clathrate type 2, the A-cages are arranged in rhombus tiles along the two-fold symmetry axis. For clathrate type 3, the B-cages are arranged in the Frank–Kasper sigma tiles along the four-fold axis.

Network representation of clathrate

To simplify visualization of the clathrates, we occasionally represent the truncated TBP host lattice by a network connecting tetramer centres (Fig. 2a–d,f,g). This representation is also shown in ref. 7.

Free-energy calculation for the host and guest sublattices

The free energies of the truncated TBP guest sublattice and the truncated TBP host sublattice were calculated along two paths (A and B) using the Frenkel–Ladd method22,23, implemented in HOOMD-blue. For each path, we selected multiple states (1, 2, 3, …) between the reference state and the equilibrium state that are connected to subsequent states by a continuous and reversible harmonic potential (Extended Data Fig. 3), which allows us to calculate the free energy difference between the reference and equilibrium states. From the reference state to the A1 state (the first state along path A), we released the guest particles and constrained the positions of the host particles as in the reference state. In this way we obtained FA1_host (=FRef_host) and FA1_guest independently (FA1_host is the free energy of the host particles at the A1 state). In the same way, we obtained FB1_host and FB1_guest (=FRef_guest) for the B1 state, by releasing only the guest particles from their positions in the reference state. For the A2 state, we partially released the host particles from A1 until the host particles in A2 showed the same translational and rotational displacement as that of the host particles in B1 (Extended Data Fig. 4), giving FA2_host ~ FB1_host. Similarly, for the B2 state, we partially released the guest particles from their positions in B1 until the guest particles in B2 showed the same translational and rotational displacements as that of A1, giving FB2_guest ~ FA1_guest. In this way, we calculated the free energy of the guest and host particle subsystems independently for each intermediate state. As we increase the number of intermediate states n, Fguest and Fhost will approach Freal_guest and Freal_host.

Free-energy calculation of bulk phases

The Gibbs free energy per particle, G/NkBT, of the bulk fluid and clathrate phases was calculated using a combination of the Frenkel–Ladd method22 and thermodynamic integration39. The approach is described in detail in ref. 23.

Dynamics of guest particles

The dynamics of the guest particles were characterized by the accumulated positions of a directional vector on the guest particle mapped onto a unit sphere (Fig. 4b–e,g, left panels). These positions were converted to spherical coordinates (r, θ, φ) and plotted as a function of θ and φ for r = 1 (Fig. 4b–e, g, right panels). For a single truncated TBP guest, the directional vector was assigned from the centre of the truncated TBP towards one of the three equatorial tips of the truncated TBP (Fig. 4a). For a truncated TBP tetramer guest, the directional vectors are a set of vectors from the centre of the tetramer towards the centre of each truncated TBP (Fig. 4f). For freely rotating guests, the distribution of accumulated positions of the directional vector is isotropic. For guests rotating around an axis, the distribution is a continuous band. For the discretely rotating guests, the distribution is localized into clusters.

The type of guest particle motion is determined, in part, by α and by the prolateness, Pr, of the cavity. The size ratio α has already been defined. Pr measures the prolateness (or the oblateness) of a cluster of points from the elements of the radius of gyration tensor40:

$${S_{mn}} = {\frac{1}{N}\mathop {\sum}\limits_{i = 1}^N {r_m^{(i)}r_n^{(i)}}}$$

where \({r_m^{(i)}}\) is the mth Cartesian coordinate of the position vector r(i) of the ith particle. Because the radius of the gyration tensor is a symmetric 3 × 3 matrix, it can be diagonalized:

$${{S} = \left[ {\begin{array}{*{20}{c}} {\lambda _x^2} & 0 & 0 \\ 0 & {\lambda _y^2} & 0 \\ 0 & 0 & {\lambda _z^2} \end{array}} \right]}$$

where the diagonal elements are chosen such that \({\lambda _x^2 \le \lambda _y^2 \le \lambda _z^2}\). These eigenvalues are used to define the Pr of the cavity as shown in ref. 40:

$${{\rm{Pr}} = \frac{{\mathop {\prod}\nolimits_{i = 1}^3 {(\lambda _i - \bar \lambda )} }}{{\bar \lambda ^3}} = 27\frac{{\det {{{\hat{\mathbf S}}}}}}{{({{{\mathrm{Tr}}}}\,{{{\mathbf{S}}}})^3}}}$$

In the limit of high prolateness, Pr = 2, and in the limit of high oblateness, Pr = −0.25.

Entropic bond calculations

Entropic bond calculations followed the same protocol as detailed in ref. 25. Briefly, we constructed each of the observed clathrate structures for truncated TBPs spanning the range of truncation values S studied in the main text. For each lattice, we solved for the energy of formation as a function of changing the truncated TBP’s degree of truncation. Construction of the phase diagram in Extended Data Fig. 7 was performed by selecting for the lowest-energy structure at each truncation value S. Energy computations for the rotational behaviours of guest particles were performed for a single cage type (either B-cage or D-cage). The guest particles were rotated so that the directional vector (defined in Fig. 4a) points at all possible angular locations on the surface of a unit sphere, each corresponding to a specific orientational configuration. The orientations of the cage particles were kept the same across all sampled guest orientations.

Molecular dynamics simulations

We used a patchy particle model7,28 that has been used previously to model DNA-functionalized nanoparticles, where a hard particle core (here, a truncated TBP) is surrounded by patches interacting with patches on other particles. The patch–patch interaction is modelled by a shifted Gaussian potential, \({V\left( r \right) = \varepsilon _{\rm{Gauss}}{{{\mathrm{exp}}}}\left[ { - 1/2\left( {r - r_0/\sigma _{\rm{Gauss}}} \right)} \right]}\), where εGauss is the depth of the potential well, σGauss is the width of the potential well, r is the centre-to-centre distance between patches, and r0 is the potential energy minimum distance. In this simulation, we used εGauss = 0.1ε, σGauss = 0.3σ and r0 = 0.45σ for the host particles, where ε and σ are the energy and length units used in MD simulations. The shifted Gaussian potential was not applied to the guest particles to make them purely repulsive. This model uses the discrete element method41 to model the hard-core TBP, and the rounding radius for the method is 0.2σ. The xyz coordinates of the patches on a TBP are listed in Supplementary Table 3.

To test the stability of the Clath I–A0B1 system using this model, we constructed an ideal structure with N = 2,315 TBPs (S = 0.52) and placed the crystal cluster at the centre of a simulation box with volume fraction VN_total/Vbox ≈ 0.01. We then stabilized the system in the NVT ensemble at \({T^ \ast /T_{\rm{m}}^ \ast \approx 0.90}\) for 4.2 × 106 MD timesteps, where T* is the dimensionless temperature (kBT/ε) and \({T_{\rm{m}}^ \ast}\) is the melting temperature of the clathrate. Each timestep (∆t) is 0.001τ, where τ is the simulation time unit. The potential energy of the system was monitored over time (Supplementary Fig. 4b) to ensure the system was thermodynamically (meta-)stable.