Crystallization of non-convex colloids: the roles of particle shape and entropy

Abstract

Particle shape and entropy play critical roles in order–disorder transitions; examples include liquid crystals and crystals formed from hard polyhedra. With few exceptions, colloidal crystals reported in the literature form from convex particles, for which the roles of particle shape and entropy are well explored. However, recent experimental work has shown that a cubic diamond lattice, elusive but long sought after for its wide and robust photonic band gap, can be assembled from non-convex particles. Here, we use simulations to explore the crystallization of colloidal diamond from non-convex tetrahedral-lobed patchy particles (TLPPs). Our results show how the entropic cost of binding, measured using umbrella sampling, can be tuned through subtle changes to the particle geometry. Geometric constraints, uniquely provided by the particle's non-convex features, select the staggered bond conformation required for cubic diamond. These constraints and the entropy associated with them can be geometrically tuned to vary the flexibility of the bonds. Depending on the particle geometry, TLPPs form liquid, diamond, or amorphous structures. Importantly, all geometrical parameters can be controlled in experiments, which should lead to larger crystals with less disorder, improving the resulting photonic band gap.

---

1 Introduction

Since the proposal over 35 years ago that periodic dielectric structures could exhibit remarkable optical properties,^1,2 including a photonic band gap,^1–3 researchers sought to self-assemble colloidal diamond.^4–6 Despite intense activity on synthesizing colloids with diamond-like bonds,^4,7,8 a solution eluded researchers until 2020, when it was demonstrated that a novel kind of colloidal particle with a non-convex shape and DNA-coated patches, could self-assemble into the desired diamond structure.⁹

Particle shape plays an important role in the self-assembly of different structures. It profoundly affects particle packing¹⁰ and the entropy of self-assembled structures.¹¹ Onsager showed that changing particle shape from spherical to rod-like can lead to nematic liquid crystalline order.¹² Later studies explored the myriad of structures, crystalline and amorphous, that self-assemble from convex hard polyhedra.^11,13 To date, most studies of crystallization have focused on particles with convex shapes.^5,11,13–16

However, the only known colloidal particles that self-assemble into a diamond crystal with a photonic band gap are non-convex.⁹ This result, unanticipated by simulations or previous experiments, suggests that using non-convex particles is an important modality to explore for the self-assembly of colloidal crystals.

While self-assembly involving non-convex particles leads to a number of interesting structures in experiments, including lock-and-key binding,^17,18 chains,^19–21 and assorted liquid crystalline phases,²² there are few examples of three-dimensional crystalline order emerging from the self-assembly of non-convex particles, and there has been only limited exploration of 3D crystallization of non-convex particles in simulations.^23–26

In this paper, we report simulations of non-convex colloidal particles that self-assemble into a cubic diamond lattice, as described by He et al.⁹ The particles consist of four overlapping tetrahedrally oriented spherical lobes with DNA-coated sticky patches nestled between the lobes (Fig. 1). The particles are designed so that DNA-coated patches on neighboring particles can bind only if their lobes interlock, as shown in Fig. 1a. The interlocking lobes enforce the staggered conformation of bonds required for a cubic diamond crystal.

Fig. 1 SEM images and 3D renders depicting the geometry of a tetrahedral lobed patchy particle (TLPP) synthesized in experiment and input into the molecular dynamics simulations, respectively. (a) Image of two bound TLPPs at the “kissing” geometry, where the patches are touching and the lobes are perfectly interlocked. (b) Image of two bound TLPPs at a larger size ratio, where the lobes are not interlocked. The TLPPs are capable of exploring a range of bond angles, off of the normal bond axis, and torsion angles, about the normal bond axis (c) and (d) images highlighting a TLPP's compression ratio d_cc/2a, the ratio of the center-to-center distance between two neighboring lobes d_cc to the diameter of the spherical lobes 2a. The compression ratio quantifies the degree to which the lobes of the TLPP overlap each other. (e) and (f) Images highlight the size ratio b/a of the TLPP, where b is the distance from the center of the TLPP to the edge of the patches. The size ratio quantifies how far the patches protrude and thus the the degree to which bound particles can rock and twist relative to each other.

The interlocking of these particles' concave features is designed to restrict the rocking and twisting of bonds between particles when their patches bind. However, perfectly interlocked particles cannot rock or twist, imposing a prohibitive entropic binding cost. Therefore, some rocking and twisting must be allowed, but how much and how can it be controlled? As we shall see, the degree of rocking and twisting can be controlled by adjusting the radial extent of the patches (concavity).

The use of patchy spheres to promote the self-assembly of cubic diamond has long been studied using simulations.^5,16,26–33 In these studies, patches on the surface of spheres are precisely arranged, shaped, and sized to enforce the valence and, if possible, the bond conformation. Simulations using patch geometries designed to enforce the staggered bond conformation propose complex features that have proved experimentally unrealizable thus far.^5,31 Alternatively, binary mixtures of patchy particles have been proposed to promote the selection of even-numbered rings, as diamond has only six-membered rings.¹⁶ Unfortunately, this does not exclusively select the staggered bond conformation required for cubic diamond. Similarly, triblock patchy particles of spheres³² and rods,³⁴ designed with patch-specific interactions or with a hierarchy of bond interaction strengths, self-assemble into mixtures of diamond-like polymorphs involving both staggered and eclipsed bond conformations.

The present study is distinguished from these previous studies in that the concave particle shape, rather than the precise size and shape of the patches, is used to impose the requisite staggered bond conformation. Moreover, the particles have been realized experimentally.⁹ In contrast to a previous study,²⁶ the attractive interaction used in our simulations is both short-ranged and narrow, which isolates and highlights the role of entropy in the self-assembly, while better modeling the experimental conditions used by He et al.⁹

In this paper, we report simulations of tetrahedral-lobed patchy particles (TLPPs) with short-range DNA-like interactions to explore the roles of particle shape and entropy. We find that the influence of these two factors is subtler and more important than one might expect. For example, the four lobes of a particle must overlap to a considerable extent for crystallization to occur. Moreover, these new simulations reveal how the entropic contribution to the binding free energy can be tuned using the radial extent of the patches, which affects the concavity of the particles. They also elucidate entropy's critical contribution to the liquid-to-crystal phase transition and determining the crystal-to-amorphous boundary, expanding on previous work exploring the role of entropy in the crystallization of convex particles.^33,35 Finally, they highlight the subtle and fundamental role of many-body interactions in stabilizing the diamond phase.

2 TLPP geometry

The tetrahedral-lobed patchy particles (TLPPs) synthesized and crystallized by He et al.⁹ comprise four spherical lobes arranged in a tetrahedron and compressed into each other, with one oil droplet at their center that is extruded and polymerized to form four patches. The resulting geometry is shown in Fig. 1. The non-convex geometry of the TLPPs is characterized by two parameters: the compression ratio and the size ratio. The compression ratio, depicted in Fig. 1c, quantifies the overlap of neighboring spherical lobes. The size ratio, depicted in Fig. 1e, quantifies the degree to which the patches extend into the concave basin formed by three lobes surrounding each patch. Together, the compression and size ratios determine the concavity of the TLPPs, whether they can bond, and whether they form the diamond structure.

A bond between two particles forms when two patches on neighboring particles touch. As shown in Fig. 1a and b, bonding can occur only if the patches extend sufficiently far out from the TLPP center, that is, only if the size ratio is greater than some minimum value. When the size ratio is at this minimum value, two TLPPs, with their lobes interlocked, as shown in Fig. 1e, have seven contacts, two contacts for each of the three interlocked lobes and one for the bound patches. We call this the “kissing” geometry, depicted in Fig. 1a. When patches extend beyond the kissing value, fewer contacts can form, which allows the two particles some freedom to rock and twist relative to each other while the patches remain bound, as in Fig. 1b.

In addition to having the proper geometry, an attractive interaction is necessary to enable binding and crystallization. In experiment, the four patches are functionalized with DNA that facilitates a short-range, attractive interaction. The strength of the DNA interaction can be tuned by changing the temperature during assembly. The four lobes of the TLPP ensure that the patches are capable of binding only when two TLPPs are rotated ∼60° with respect to each other, interlocking the lobes as illustrated in Fig. 1a and b. The lobe–lobe interlock ensures the staggered orientation of next-nearest neighbor bonds, which is essential for the selective growth of cubic diamond.⁹

The combination of directional interactions and non-convex particle geometry means that two-particle binding is dictated by a competition between the attractive DNA interaction and the entropic cost of binding. The TLPP's compression and size ratios determine how tightly two particles interlock when they bind, influencing both the entropy of binding and the disorder of the final lattice.

As noted above, the kissing geometry is the most tightly interlocked configuration possible and would lead to a “perfect” bond with no orientational disorder. However, such bonds would also completely restrict the rocking (orientational) and twisting (rotational) degrees of freedom of bound particle pairs (see Fig. 1). For particles to establish these “perfect” bonds, they must relinquish all of their entropy, a sacrifice that renders this tight fit entropically forbidden.

A similar phenomenon was noted by Sacanna et al.¹⁷ for lock and key colloids whose binding is driven by a depletion interaction between geometrically complementary features. They showed that if the lock and key are made to fit perfectly, the particles do not bind; instead, some free volume has to be left to ensure that the entropic binding cost is not too high. Similarly, we show in the following sections that retaining some orientational and rotational freedom decreases the entropic binding cost and promotes crystallization.³⁵

One consequence of loosening the lobe interlock is the introduction of orientational and rotational disorder in the bonding. In the limit opposite to the “kissing” geometry, the patches can protrude so far that they can bind and freely rotate without any lobe–lobe interlock, as illustrated in Fig. 1b. Here, there is little geometric bias to the angular conformation of the bond, so bonds have random orientations and form an amorphous structure instead of a diamond. By performing molecular dynamics simulations in HOOMD-Blue,³⁶ we vary the TLPP geometry and map out a phase diagram of the various structures that form. These include a colloidal liquid, where particles are incapable of binding, a crystalline diamond phase, where the particles bind and interlock, and an amorphous phase, where the particles bind but do not interlock or have long-range tetrahedral order. Throughout these regimes, we use umbrella sampling to numerically determine the entropic cost of binding for a pair of TLPPs. A quantitative understanding of the entropy associated with different geometries gives us the tools to manipulate the crystallization phase space by tuning the interaction strength. This command of the geometric landscape for crystal assembly should prove useful in experiments where highly tunable interactions like DNA are employed to grow larger and less disordered crystals.

3 Computational approach

A simulation TLPP is made up of four large spheres (“lobes”), which are typically about 1 μm in diameter, and four small spheres (“patches”) designed to mimic the shape of real TLPPs. The lobe centers are positioned at the corners of a tetrahedron and partially overlap to achieve the desired compression ratio. The patch centers are situated between the lobes and pushed out from the TLPP center to achieve the desired size ratio. In real TLPPs, the patches are created from a liquid droplet with a fixed volume that fills the space between the spherical lobes. Therefore, the size of the simulated patch sphere is chosen so that the interstitial patch volume remains constant. This ensures that, similar to real TLPPs, smaller size ratio patches are flatter, whereas higher size ratio patches bead up on the surface of the lobes. Each TLPP is treated as a single rigid body within HOOMD-Blue. Additionally, a ninth non-interacting sphere is placed at the center of the TLPP to serve as the center of mass.

A full simulation comprises 2400 TLPPs that start evenly distributed throughout a tall but narrow simulation box with periodic boundary conditions. The box is 4a wide, where a is the lattice constant of the perfect diamond. TLPP positions are evolved using Langevin dynamics. Additional details defining the TLPP geometry in simulation are given in ESI† Section S1. The particles are placed in a constant field applying a force of F_g = −0.025k_BT/a in the vertical (z) direction, corresponding to a nearly density-matched solution in experiments with gravitational height h_g ≈ 14 μm. At the bottom of the tube is a hard wall with seven ordered TLPPs immobilized on the surface to template and seed the diamond growth for conditions under which the crystal is thermodynamically favorable. In four parallel simulations without a template, the time to nucleate a crystal was found to vary dramatically and be much longer than in the templated case. By templating the crystal, nucleation occurs faster and more consistently, allowing us to spend more computational time observing crystal growth, necessitating only a single simulation for most TLPP geometries.

After the system has been initialized, it is left to evolve for 4 × 10⁸ time steps. By studying the diffusion of one colloidal particle and comparing it to experimental particle tracks, a simulation time step Δt was determined to be equivalent to approximately 5 μs; 1 second ≈ 2 × 10⁵ Δt. Each particle's trajectory is governed by diffusion in the presence of a gravitational field, subject to interactions with the hard wall at the bottom and particle–particle interactions described in the next section. Further simulation details are given in ESI† Section S1.

3.1 Wang–Frenkel interaction potential

Recent measurements made using total internal reflection microscopy (TIRM) by Cui et al.³⁷ have shown that the colloidal DNA interactions we want to replicate are well-described by a short-ranged Wang–Frenkel potential,³⁸ which we use for the attractive potentials in simulation. The potential has two parameters r_c and σ, which define the range of interaction, and a third parameter ε, which defines the depth of the well in k_BT:wherewhen σ < r < r_c, the potential is attractive with a minimum of depth ε, which models the range and strength of two particles' DNA brushes binding. For r < σ the potential becomes repulsive, so σ is the radius of the compressed DNA brush where it starts behaving like a hard sphere. The Wang-Frankel potential U_WF is designed so that it smoothly goes to zero as r approaches r_c. For our particles, r_c is the maximal extent of the DNA brush. For spheres of diameter 1 μm, a typical brush compression is 20 nm, with the maximal extent of the DNA interaction 10 nm beyond that. We choose ε to be 10k_BT for the DNA interaction on the patches, experimentally reached by lowering the temperature of the sample below the melting by approximately one degree. When two TLPPs come into contact, we allow the lobes to interact with a weak Wang–Frenkel potential (ε = 1k_BT) to mimic a small amount of depletion present in our experiments.

A separate regime of strong depletion and a wider range of attraction, with depletant sizes up to 40% of a TLPP lobe, was studied by Marín-Aguilar et al.²⁶ By using a shifted 96–48 Lennard-Jones potential representing a Mie interaction, they probed the phase space of crystallization when the attractive interactions have a width of 10% of their spherical particle's diameter, an order of magnitude wider than our DNA-like Wang–Frenkel potential. Here, we show that use of a narrow interaction range, like those resulting from the use of DNA in experiments, limits the crystallization phase space because the particles must overcome the steep entropic barrier associated with the concave particle shape, which would be compensated for by a longer-ranged attraction.

4 Results and discussion

4.1 Overview

In a typical simulation, we allow the particles to sediment down onto the templated substrate. Gravity creates a concentration gradient in z, biasing the crystal to nucleate in the more concentrated regions near the bottom. For favorable crystallization conditions, particle–particle interactions lead to one of three states: diamond, liquid, or amorphous. Diamond, on the bottom, and liquid, on top, phases can be seen in the simulation output of Fig. 2.

Fig. 2 (Left) Particle geometry input into the simulation. The radius of curvature of the patch is varied as a function of compression and size ratio to conserve the volume of the patch as the geometry changes. (Middle) A typical simulation output. (Right) A depiction of the orientational freedom afforded to particles that bind when the patch size ratio is larger and the lobe–lobe interlock is weaker.

When the TLPP compression and size ratios allow lobes to interlock and patches to touch, all bonds are staggered, and the cubic diamond structure forms. However, when the patch size ratio is too small, the recessed patches cannot bind, and a liquid state results. Conversely, when the patch size ratio is too high, patches bind, but lobes fail to interlock, leading to randomly oriented bonds and an amorphous structure. Between these extremes, patches bind while lobes interlock with a finite free volume, allowing particles to explore various configurations. Particles binding with finite free volume, depicted in Fig. 2, possess higher entropy compared to rigidly interlocked particles. In subsequent sections, we examine the transition from liquid to crystalline to amorphous phases as the size ratio increases, examining in some detail how these phases correlate with the entropic cost of binding.

4.2 Structure determination

We set the compression and size ratios for each simulation and then determine the final structure to be either crystal, liquid, or amorphous based on the radial distribution function g(r) and the structure factor S(k). The final structure varies as a function of the height z due to the applied gravitational field, which induces a concentration gradient. Therefore, g(r) and S(k) are calculated as a function of height using the local density, which is sufficiently constant over a dozen or more coordination shells to obtain meaningful results.

The three phases, liquid, crystalline (diamond), and amorphous, are readily distinguished by their different radial distribution functions g(r), which are shown in Fig. 3a. For comparison, we have also simulated an assembly of spherical patchy particles, compression ratio 0.0, and plot its g(r) in ESI† Section 2. Additionally, the structure factor S(k) is calculated for these structures, as shown in Fig. 3b. While S(k) and g(r) both easily identify the diamond phase, we focus on g(r) because it also distinguishes between the liquid and amorphous phases we observe in our system.

Fig. 3 Pair distribution functions g(r), structure factors S(k), and particle positions and bonds for three final structures obtained in MD simulations. (a) The pair distribution function g(r) is plotted for three different patch size ratios at compression 0.70. Size ratio 1.30 (green) has patches that stick out far enough to bind while still allowing the lobes to interlock. This geometry results in a structure with strong g(r) peaks. By comparing to the dashed lines, which are the five first nearest neighbor distances for a perfect diamond of size ratio 1.30, it is clear that this structure is a diamond. For size ratio 1.20 (red), the patches are too recessed to be able to bind; g(r) has a large liquid peak, but no longer-range correlations. Size ratio 1.40 (blue) forms bonds based on the first nearest neighbor peak but has no long-range order, forming an amorphous structure. (b) The structure factor S(k) for these structures confirms that size ratio 1.30 (green) has long-range order. Comparing the peaks to the dashed lines, which correspond to the first seven S(k) peaks for a perfect diamond, the structure can again be identified as diamond. Both size ratios 1.20 and 1.40 have no S(k) peaks, corresponding to a lack of long-range order. (c)–(e) The particle centers of mass (blue) and their bonds (gray) for the final frame of the simulations with compression 0.70 and size ratios (c) 1.20, (d) 1.30, and (e) 1.40.

The radial distribution function g(r) for diamond (Fig. 3) exhibits a large first peak at the bond length calculated for a perfect diamond crystal. Subsequent peaks correspond to the second through fifth nearest neighbors of a perfect diamond crystal. This phase is observed for intermediate values of the size ratio where the lobes loosely interlock.

For the liquid phase, the first peak in g(r) is lower, broader, and shifted to larger distances corresponding to the typical particle–particle separation in the dilute phase (Fig. 3). At larger distances, the correlation peaks diminish rapidly as particle positions become increasingly uncorrelated. This phase is observed for the smaller size ratios where no permanent bonds form.

For the amorphous phase, the first peak in g(r) is observed very near the diamond bond length, but with a diminished height compared to the crystalline phase. However, for larger distances, the correlations decay at roughly the same rate as for unbound particles. We classify these disordered structures, with short-range binding but no long-range correlations, as amorphous. We characterize and discuss the disorder of the amorphous phase in greater detail in Section 4.4.

This analysis is carried out for systems of different compression and size ratios. The results are summarized by the phase diagram shown in Fig. 4. The phase diagram shows a crystalline phase over a significant range of size ratios for compression ratios near 0.7. Increasing the compression ratio causes the range of size ratios where crystallization occurs to narrow. As noted above, unbound particles form a liquid state at the lower size ratios, while an amorphous phase is observed at the highest size ratios.

Fig. 4 Phase space of structures found in MD simulations based on TLPP geometry. At lower size ratios, the patches are too recessed, and the entropic cost of binding is too high, so the particles remain unbound in a liquid state. Once the patches protrude further, the TLPPs bind and form diamond for a large range of size ratios (at moderate compression). Once the size ratio is increased to the point where TLPPs no longer interlock, the particles bind and form an amorphous structure. The size ratio delineating the crystalline regime from the liquid and amorphous regions depends on the compression ratio. The gray line shows the ideal kissing geometry, which is the predicted absolute cutoff; to the left of this line, the diamond geometry cannot be formed.

4.3 Crystallization phase space

The phase space for crystallization as a function of size and compression ratios is shown in Fig. 4. This diagram highlights how the geometry of the concave particles influences the final state of the system. A useful reference is the kissing geometry, discussed in Section 2 and represented by the gray curve in Fig. 4.

For size ratios below the kissing geometry, the patches on the particles are too recessed to make contact and form bonds, resulting in a liquid phase. As the size ratio increases past the kissing geometry, the TLPPs remain unbound because the entropic cost of binding is too high for bonds to form (see Section 4.7 for further discussion).

At still larger size ratios, the particles can bind, forming either a crystalline or an amorphous phase.

In the crystalline phase, indicated by the blue points in Fig. 4, two particle patches can only bind if their lobes interlock. This specific arrangement ensures that the bonds adopt a staggered conformation, resulting in a cubic diamond lattice.

In the amorphous region, marked by the black points in Fig. 4 for larger size ratios, the patches extend so far that particles can bind without their lobes interlocking. Because the bond orientations are not confined to the staggered conformation, these particles do not crystallize into cubic diamond; instead, they form an amorphous structure.

When bond orientations are unrestricted, one might expect the particles to behave like convex particles with tetrahedral patches. In this scenario, tetrahedrally bound particles can adopt various bond conformations, including staggered, eclipsed, and random arrangements. Based solely on geometry, several structures are possible, including cubic diamond, hexagonal diamond, and clathrate structures. Hexagonal diamond is characterized by a 3 : 1 ratio of staggered to eclipsed conformations, while clathrate structures comprise entirely eclipsed conformations.

In previous simulations of convex patchy particles, mixtures of these structures could lead to amorphous formations.^30,39 Moreover, the structures observed could depend on kinetic factors. In some simulations, kinetic factors produced clathrate structures, even though free-energy calculations showed that diamond structures are always thermodynamically favored.³⁹

Therefore, if our system were allowed to evolve over a significantly longer timescale, it could potentially anneal into cubic diamond, hexagonal diamond, or other tetrahedrally bound phases. Even when these other phases can form through annealing, the unrestricted bond conformations will lead to competition among multiple phases within the system.^30,39

While analysis of the pair correlation function g(r) indicates that all the blue points in the phase diagram are crystalline, the crystals all differ in a few key ways. First, the degree of interlock, which is related to the entropic cost of binding, is different for different size ratios. Particles with larger size ratios have more rotational freedom in their bonds. Second, the rate of growth of the crystals differs for different geometries. The growth is relatively slow for smaller size ratios and faster for larger size ratios. Last, the entropic cost of binding varies with particle geometry, affecting all aspects of the assembly. In the following sections, we discuss and analyze all of these trends.

4.4 Torsion analysis

Understanding the degree of interlock for the various geometries requires identifying all the bonds and determining how far they are from a perfect staggered conformation. This provides insight into local disorder in the bonds of the crystalline samples and how the transition from a crystalline to an amorphous state is associated with weakening the lobe–lobe interlock.

The bonds formed during a simulation are identified by finding nearest neighbors whose patches are in contact. To quantify how tightly interlocked the bound TLPPs are, we calculate the torsion angle of each bond. The torsion angle, defined as rotation about the axis of the bond and depicted in Fig. 5, is 0° for TLPPs aligned in the eclipsed conformation and 60° for perfectly staggered TLPPs. After analyzing all of the bonds in a simulation of a given geometry, the distribution of the torsion angles is analyzed.

Fig. 5 Distributions of torsion angle between bound particles for various patch size ratios at compression 0.70 and standard deviation of torsion distributions for compression 0.70. (inset) TLPP and stick model of two particles binding with each other and three other neighbors. The bond in the middle of the two central TLPPs defines an axis about which the torsion angle θ of the bond can be measured. (a) The distribution of torsion angles for the lowest size ratio that crystallizes b/a = 1.26. The bonds that form are Gaussian distributed about 60° with a standard deviation of σ = 5°, meaning they are in the staggered conformation with relatively small angular fluctuations. (b) When the size ratio is increased to b/a = 1.30, the bonds are still peaked at 60°. However, the width of the distribution has increased to σ = 7°. The larger size ratio affords the bound particles more rotational freedom. (c) The largest size ratio b/a = 1.40 has a flatter distribution. Any conformation is allowed because the patches can bind without the interlocking lobes. The relatively small peak at 60° is due to weakly attractive depletion interactions between lobes; the staggered conformation maximizes the number of such contacts and thus lowers the free energy slightly. (d) Plotting the torsion standard deviation for compression d_cc/2a = 0.70 as a function of size ratio shows the standard deviation increases continuously with size ratio up to b/a = 1.38, where it then jumps from σ = 10° to σ = 33° at b/a = 1.39.

The distributions for three characteristic size ratios at compression ratio d_cc/2a = 0.70 are plotted in Fig. 5. Referring to the phase diagram in Fig. 4, the three torsion distributions plotted correspond to: b/a = 1.26, the lowest size ratio for which binding is observed; b/a = 1.30, a point near the middle of the crystalline regime; and b/a = 1.40, an amorphous structure assembled from the largest size ratio studied.

For all three size ratios, the torsion angle distributions peak at 60°. However, the widths of these distributions vary widely. As shown in Fig. 5a, the standard deviation of the torsion angle for the lowest crystallizing size ratio is σ = 5°. This corresponds to the least rotational freedom required to make binding entropically favorable. Comparing the two crystalline systems, the standard deviation in the torsion is increased from σ = 5° to σ = 7° as the size ratio increases from b/a = 1.26 to b/a = 1.30. The torsion standard deviation increases with the size ratio because the more extended patches can bind with less lobe interlock, leaving more free volume for the TLPPs to explore. The final amorphous size ratio shown, b/a = 1.40, differs significantly from the previous two. The distribution is flatter with a slight bias towards the staggered conformation, but all conformations from staggered to eclipsed are represented. Clearly, the TLPPs' lobes no longer restrict the bonds to form only near the staggered conformation. The patches protrude far enough at size ratio b/a = 1.40 that they can bind regardless of the orientation of the lobes.

The torsion angle distributions are analyzed for every point in the phase space and the standard deviation for all size ratios at compression d_cc/2a = 0.70, the same geometries discussed in the previous section, are plotted in Fig. 5d. The standard deviation increases continuously as the size ratio is increased from b/a = 1.26, for which σ = 5°, up to b/a = 1.38, the highest size ratio that crystallizes, for which σ = 10°, which corresponds to a looser interlock.

When the size ratio increases from b/a = 1.38 to b/a = 1.39, the standard deviation increases discontinuously from σ = 10° to σ = 33°. The discontinuous jump in the torsion angle standard deviation occurs at precisely the same point where the structure transitions from crystal to amorphous in Fig. 4, which agrees with our expectation that the amorphous structure is associated with a lack of lobe–lobe interlock and staggered conformation binding.

4.5 Kinetics of crystallization

In addition to the analysis of the final structures, it is also useful to study how each simulation reaches its final state. Does the growth of the crystals depend on the compression ratio and size ratio of the particles? We address this question by analyzing the rate of bond formation for all points in the crystalline and amorphous regions of the phase space (Fig. 4).

In Fig. 6, we plot the number of bonds per particle as a function of simulated time for the same conditions discussed in the previous section (d_cc/2a = 0.70), two that ultimately crystallize (b/a = 1.26 and b/a = 1.30) and one that ends up as an amorphous structure (b/a = 1.39). Two particles are determined to be bound when their patches are close enough together to interact via their WF attractive potentials. In all three cases, the number of bonds formed per particle increases sharply with time as particles concentrate and begin to bind. All geometries' number of bonds then plateau at approximately 10⁴ seconds. At this plateau, particles are actively binding and unbinding, but not crystallizing. After the plateau, size ratios b/a = 1.26 and b/a = 1.30 both begin to sharply rise again. Inspection of the simulation around these times shows that crystals begin to nucleate on the surface at approximately 3 × 10⁴ seconds (b/a = 1.30) and 4 × 10⁴ seconds (b/a = 1.26). This second increase in bond number is then associated with the particles crystallizing into diamond.

Fig. 6 Number of bonds per particle in simulations as a function of time. At the compression ratio of 0.70 for three size ratios, we see that the rate of bond formation increases as the size ratio increases. Size ratios 1.26 (blue) and 1.30 (green) form diamond but size ratio 1.30's more protruded patches allow it to grow faster. Size ratio 1.40 (red) quickly plateaus at a lower number of bonds, indicating an amorphous phase.

For diamond, the maximum number of bonds per particle is two since each particle has four bonds, each shared across two neighbors. Particles in the liquid phase have zero bonds, and particles that are part of crystal defects and boundaries have fewer than four bonds. Since every crystal coexists with some liquid (see Fig. 2), the number of bonds per particle is always less than the maximal value of 2.

The green and blue curves in Fig. 6 show the bond formation kinetics for two simulations that crystallize. The number of bonds per particle does not plateau by the end of the simulation, indicating that the crystals are still growing when the simulation terminates. The initial rate of bond formation, the initial slope of the curves, is greater for TLPPs with the larger size ratio b/a = 1.30 (green) compared to those with smaller size ratio b/a = 1.26 (blue). The rate of bond formation for the largest size ratio b/a = 1.40 (red), associated with an amorphous structure, is even higher. Generally, for a given compression ratio, the initial rate of bond formation is faster for larger size ratios. The increased rate of bond formation is associated with greater rotational freedom (and greater entropy), which poses fewer constraints on binding.^17,35

In contrast to the other two size ratios, the number of bonds per particle does plateau for size ratio b/a = 1.40 and it plateaus at a value of 0.6, well below two bonds per particle. The amorphous structures have far fewer bonds than the diamond lattice due to the disorder in the bond conformation. When bonds can be formed in non-staggered conformations, a particles' nearest and next nearest neighbors can impede bonding to all of its patches.

4.6 Umbrella sampling for pairs of TLPPs

While the previous simulations provide insight into the structures formed by different TLPP geometries, they do not provide quantitative information about the entropy of bound and unbound TLPPs. As noted previously, the entropic cost of binding plays an integral role in determining the equilibrium structure formed by the TLPPs. To understand the entropic cost, we use the umbrella sampling method,⁴⁰ introduced by Torrie and Valleau,⁴¹ to find the free energy of binding of two TLPPs in an MD simulation.

If an MD run samples all of phase space, the system is ergodic, and the ensemble-averaged probability distribution of center-to-center distances ξ of two interacting TLPPs equals the time-averaged probability distribution. In the case studied here, two particles diffusing and interacting typically do not sample all of phase space because of the large entropic barrier that must be overcome to reach the bound state. Umbrella sampling addresses this problem by introducing a biasing potential U_b(ξ) to ensure that the particles are able to sample the full free energy landscape despite a large entropic barrier. The biasing potential restricts a single simulation to a window of separations; multiple simulations need to be run, shifting the biasing potential to probe the full range of separations from bound to unbound states. The entire set of simulations with different biasing potentials can be analyzed to determine the free energy over a wide range of ξ using the weighted histogram analysis method (WHAM),⁴² as implemented by Grossfield.⁴³

When performing umbrella sampling simulations, it is important to choose an appropriate U_bias to ensure that the system is ergodic. For studying the free energy of two particles binding, we use a harmonic biasing potential: , where k is the spring constant, ξ is the center-to-center distance, and ξ₀ is the rest length of the spring. The two TLPPs are allowed to diffuse, rotate, and interact with each other while also being connected by a spring, as illustrated in Fig. 7a. The spring keeps the particles at separations near their rest length, giving the particles sufficient time to sample those separations. The full range of separations is sampled by performing simulations with springs of various rest lengths, spanning from separations at which the particles are bound (ξ ∼ 1 μm) to separations at which particles are unbound (ξ ∼ 2 μm). Further discussion of the spring constants used can be found in ESI† Section S3.

Fig. 7 (a) Depiction of a typical umbrella sampling simulation. Two particles have their centers of mass connected by a harmonic biasing potential (spring), while they are allowed to freely rotate, diffuse, and interact with each other. The harmonic potential biases the particles to sample separations around their rest length and overcome free energy barriers. Simulations with different rest lengths sample the range of separations from bound to unbound. (b) Free energy of two TLPPs binding for three different TLPP geometries with compression ratio d_cc/2a = 0.7 and size ratios b/a = 1.15 (blue), b/a = 1.30 (orange), and b/a = 1.40 (green). At the lowest size ratio of 1.15, the patches are geometrically incapable of binding and interacting, so there is no free energy minimum due to the attractive potential. At higher size ratios, the patches can bind, leading to a local minimum of free energy at the separation. (c) Entropic barrier as a function of size ratio for compression d_cc/2a = 0.70. (d) Free energy minimum of the potential well as a function of size ratio for compression d_cc/2a = 0.70.

The spring's rest length determines the center of the sampling region, and the spring constant determines its width. Larger spring constants are associated with sampling a smaller range of separations. When sampling particle separations close to binding, we use a spring constant of 3265k_BT/μm², which corresponds to a range of separations approximately 20 nm wide. This ensures that we effectively sample the free energy landscape where the entropic barrier and potential well are strongest, out to separations of approximately ξ = 1.37 μm. At larger separations ξ > 1.38 μm, the particles do not interact very strongly, and a weaker spring k = 326.5k_BT/μm² is used, which allows us to sample larger ranges of separations more efficiently. After performing 14 simulations, 10 with the strong spring and 4 with the weaker spring, the distributions are analyzed with WHAM, and the free energy as a function of separation between the two particles F(ξ) is determined.

4.7 Free energy of binding

The free energy of binding for TLPPs at compression d_cc/2a = 0.70 with three different size ratios are shown in Fig. 7b. Two key features characterize the free energy landscapes: (1) the entropic barrier and (2) the free energy well. The TLPPs with the lowest size ratio b/a = 1.15 have patches so recessed that they are geometrically incapable of binding. As the particles approach each other, their lobes interlock, restricting their orientational degrees of freedom and creating a rising entropic barrier. The TLPPs' entropic barrier is much higher than that found for simple spheres, whose free energy is shown in ESI† Section S4. The free energy curve has no stable position along ξ and is determined only by the entropy of the TLPPs.

Particles with size ratios b/a = 1.30 and b/a = 1.40 both have patches that protrude enough to bind. For these geometries, there is a free energy well at small separations associated with the attractive Wang–Frenkel potential. The free energy well is a local minimum and is located at a distance equal to the size ratio; size ratios of 1.30 and 1.40 have minima located at 1.3 and 1.4 μm, respectively. Outside of the well, the increase in the free energy of binding is solely due to the entropic cost associated with bringing two TLPPs together, which is shown in ESI† Section S4.

The entropic part of the free energy changes very little with the size ratio; the primary effect of changing the size ratio is moving the position of the potential minimum. Because the entropic barrier increases with decreasing particle separation, moving the potential minimum increases the entropic barrier to binding, as shown in Fig. 7c. The lower entropic barrier for higher size ratios is the reason that the rate of bond formation is higher at high size ratios, a correlation noted in the previous section. The entropic barriers for all TLPP geometries studied are shown in ESI† Section S5.

Away from the potential well, the free energy curves are quite similar. At large separations ξ > 1.8 μm, there is no potential well and no entropic cost; the particles are too far apart to interact or restrict each other's degrees of freedom.

The depth of the free energy minimum for bound particles also changes with the size ratio. Fig. 7d shows the local free energy minimum due to the attractive Wang–Frenkel potential. For size ratios b/a < 1.25, increasing the size ratio causes the depth of the potential well to decrease as expected. However, at size ratios b/a > 1.25 the depth of the well plateaus. Comparing size ratios b/a = 1.30 and b/a = 1.40 in Fig. 7b, it is evident that the attractive wells are getting wider instead of deeper as the size ratio is increased. This is due principally to the increase in bond angles that are permitted when the patches protrude farther, studied in ESI† Section S6.

4.8 Effect of concentration

As noted above, the potential well obtained from umbrella sampling is a local minimum of the free energy; the free energy is smaller at large particle separations, as seen in Fig. 7d. Thus, for two TLPPs, the 10k_BT attractive interaction used in the MD crystallization simulations is, on its own, insufficient to have dimerization be thermodynamically favored under very dilute conditions. By contrast, the crystallization experiments and simulations are conducted at high concentrations. In the crystallization simulations, the particles' chemical potentials contribute to the free energy, which, at high enough concentrations, tilts the free energy landscape to make crystalline binding energetically favorable. In our simulations and experiments, gravity induces a concentration gradient of particles which promotes the nucleation of crystals at the highest concentrations. After nucleation, the crystals grow to some terminal height which is associated with the equilibrium particle concentration, shown in ESI† Section S8.

The effect of concentration on the free energy is highly complex in this case and many-bodied, as it combines effects from multiple particles being interlocked and chains of particles interacting through specific bonds. Therefore, to illustrate how crowding alone influences the binding free energy for these interlocked particles, we perform umbrella sampling simulations on two TLPPs, interacting as described in the previous section but in a bath of other TLPPs with which they cannot bind. In this case, both particles interact with all other particles through a 1k_BT attractive Wang–Frenkel potential, both for the patches and the lobes. This weaker patch interaction facilitates the bath particles' volume exclusion but is too weak for them to bind to the two TLPPs of interest. By controlling the number of bath particles surrounding the two bound TLPPs, we can analyze the free energy of dimerization at different surrounding concentrations.

The free energy for TLPPs with compression ratio d_cc/2a = 0.70 and size ratio b/a = 1.30 is plotted at three different concentrations in Fig. 8. The highest concentration, 3.8 particles/a³, is the terminal concentration for which particles bind (see ESI† Section S8). In Fig. 8, we see that for a concentration of 3.8 particles/a³ or higher, the potential well is shifted sufficiently low to become the global minimum of free energy, meaning that the crystal is the equilibrium structure. Moreover, the entropic free energy barrier is lowered from approximately 9k_BT to 4k_BT, which improves the crystallization kinetics.

Fig. 8 Free energy of binding in a batch of TLPPs for compression d_cc/2a = 0.70 and size ratio b/a = 1.30. At the given patch–patch interaction strength of 10k_BT, the particle concentration of the liquid at the interface between the liquid and the crystal is 3.8 particle/a³. Below this concentration, the local free energy minima is positive. However, at the concentration of 3.8 particle/a³, the local free energy minima becomes the global minima. The free energy is negative, which indicates that the TLPPs prefer to bind in the thermodynamic limit.

4.9 Effect of interaction strength

Another way to lower the free energy of the bound state is to increase the strength of the attractive interaction, which pushes the minimum of the well deeper so that particles can bind at lower concentrations. This further restricts the motion of the bound particles and thus reduces the entropy, but the entropic cost in free energy is compensated by the increased binding energy. Fig. 9 shows the phase space of crystallization as a function of patch size ratio and patch–patch interaction strength for TLPPs with a compression ratio of 0.70. We see that increasing the depth of the potential from 10k_BT to 12.5k_BT causes the TLPPs to start crystallizing at a size ratio of 1.22 instead of 1.25. The extra 2.5k_BT of interaction is enough to compensate for the entropic cost of shrinking the patches, as seen by comparison to Fig. 7c where the free energy difference between size ratio 1.22 and 1.25 is approximately 1.5k_BT. In experiments, the strength of the attractive DNA-mediated interaction can be adjusted by changing the temperature, giving us fine external control over this contribution to the energy.

Fig. 9 Crystallization phase space for TLPPs of compression ratio 0.7 with varying patch–patch interaction strengths. When the patches interact more strongly, higher k_BT, the particles bind and crystallize at lower size ratios. The extra binding strength can compensate for the increased entropic cost at lower size ratios and make binding energetically favorable.

5 Conclusions and outlook

The crystallization of TLPPs into a cubic diamond lattice depends critically on their non-convex features. These non-convex features, realized here by partially overlapping tetrahedrally oriented spherical lobes, are essential for creating a physical interlock. This interlock restricts the relative orientation of the bonds between particles and ensures the staggered conformation necessary for achieving the cubic diamond structure.

The physical interlock mechanism facilitated by non-convex particles contrasts strongly with previous work that sought to realize the requisite staggered conformation using convex particles decorated with elaborate patterns of attractive patches,^5,31 which have little chance of being realized experimentally. Moreover, other schemes that seek to enforce valence but do not strictly enforce 100% staggered bond conformation inevitably result in a mixture of bond conformations and thus do not crystallize into cubic diamond.^5,16,31

It is informative to contrast how atomic and molecular systems achieve the desired staggered conformation with how colloids achieve it. In both cases, nearest neighbors must bind while the bonds of next-nearest neighbors must adopt a staggered conformation. In molecules such as ethane, the nearest-neighbor carbon atoms are bound by a covalent bond (C–C), which lowers the system's energy by allowing the electronic wave function to spread out. The staggered conformation is preferred to any other conformation because of the electrostatic repulsion between different C–H bonds. Thus, the simultaneous existence of a strong attractive bond between carbon atoms and a repulsive electrostatic interaction between C–H bonds leads to a staggered conformation. Similar considerations give rise to the staggered conformation of carbon atoms in a cubic diamond crystal. In colloids, it is difficult to simultaneously realize a strong short-range attractive interaction between colloidal particles (absent irreversible aggregation due to van der Waals attraction) and a long-range repulsion that can stabilize the staggered conformation. This is particularly true in aqueous solvents, where the Debye screening length is usually much smaller than the particle diameter. Therefore, we use DNA hybridization for short-range binding of nearest neighbors and the steric repulsion afforded by the non-convex shape of the tetrahedral lobes to maintain the staggered conformation.

In both the atomic and colloidal systems described here, the short-range interaction binding nearest neighbors is attractive—covalent bonding vs. DNA hybridization—while the next-nearest neighbor interaction that favors the staggered conformation is repulsive—electrostatic repulsion vs. steric repulsion between concave surface structures. Different systems, one quantum mechanical and the other classical, share common features but demand different solutions to achieve the same overarching goal.

A nuance of the colloidal system is the critical role played by entropy. As detailed in this article, the interlocking of lobes is associated with an entropic cost that can be prohibitive for crystallization if neighboring particles approach each other too closely. Fortunately, the entropic cost of binding can be adjusted by tuning the non-convex features of the particles, most notably by experimentally controlling how far the patches protrude (the size ratio) and also by controlling the degree of overlap of the lobes (the compression ratio).

The simulations provide an improved understanding of the role of entropy and non-convex particle geometry. They also provide a detailed guide for using non-convex features on particles to tune entropy and control orientational interactions. These advances open up new strategies for manipulating interactions between particles, enabling the self-assembly of structures never previously realized in colloidal science.

Author contributions

DJP and JPG conceived the project. JPG and MJK conducted the simulations. DJP and GMH supervised the research. All authors interpreted the results. JPG wrote the original draft, and all authors contributed to editing and finalizing the manuscript.

Data availability

HOOMD-Blue scripts are available from the following GitHub URL, https://github.com/pine-research-group/Gales-Nonconvex-Colloids-2025, and the simulation data is available from a NYU resource that is linked from that GitHub page.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

This research was supported by the Army Research Office under award numbers W911NF-17-1-0328 (JPG, MJK, and DJP) and W911NF-21-1-0011 (GMH).

Notes and references

1. E. Yablonovitch, Phys. Rev. Lett., 1987, 58, 2059–2062.
2. S. John, Phys. Rev. Lett., 1987, 58, 2486–2489.
3. K. M. Ho, C. T. Chan and C. M. Soukoulis, Phys. Rev. Lett., 1990, 65, 3152–3155.
4. V. N. Manoharan, M. T. Elsesser and D. J. Pine, Science, 2003, 301, 483–487.
5. Z. Zhang, A. S. Keys, T. Chen and S. C. Glotzer, Langmuir, 2005, 21, 11547–11551.
6. F. Smallenburg and F. Sciortino, Nat. Phys., 2013, 9, 554–558.
7. Y. Wang, Y. Wang, D. R. Breed, V. N. Manoharan, L. Feng, A. D. Hollingsworth, M. Weck and D. J. Pine, Nature, 2012, 491, 51–55.
8. Y. Wang, Y. Wang, X. Zheng, E. Ducrot, M.-G. Lee, G.-R. Yi, M. Weck and D. J. Pine, J. Am. Chem. Soc., 2015, 137, 10760–10766.
9. M. He, J. P. Gales, É. Ducrot, Z. Gong, G.-R. Yi, S. Sacanna and D. J. Pine, Nature, 2020, 585, 524–529.
10. A. Donev, I. Cisse, D. Sachs, E. A. Variano, F. H. Stillinger, R. Connelly, S. Torquato and P. M. Chaikin, Science, 2004, 303, 990–993.
11. P. F. Damasceno, M. Engel and S. C. Glotzer, Science, 2012, 337, 453–457.
12. L. Onsager, Ann. N. Y. Acad. Sci., 1949, 51, 627–659.
13. R. Ni, A. P. Gantapara, J. De Graaf, R. Van Roij and M. Dijkstra, Soft Matter, 2012, 8, 8826–8834.
14. P. N. Pusey and W. Van Megen, Nature, 1986, 320, 340–342.
15. E. Bianchi, J. Largo, P. Tartaglia, E. Zaccarelli and F. Sciortino, Phys. Rev. Lett., 2006, 97, 168301.
16. A. Neophytou, D. Chakrabarti and F. Sciortino, Proc. Natl. Acad. Sci. U. S. A., 2021, 118, e2109776118.
17. S. Sacanna, W. T. M. Irvine, P. M. Chaikin and D. J. Pine, Nature, 2010, 464, 575–578.
18. Y. Wang, Y. Wang, X. Zheng, G.-R. Yi, S. Sacanna, D. J. Pine and M. Weck, J. Am. Chem. Soc., 2014, 136, 6866–6869.
19. D. J. Ashton, R. L. Jack and N. B. Wilding, Soft Matter, 2013, 9, 9661.
20. K. V. Edmond, T. W. P. Jacobson, J. S. Oh, G.-R. Yi, A. D. Hollingsworth, S. Sacanna and D. J. Pine, Soft Matter, 2021, 17, 6176–6181.
21. J. S. Oh, S. Lee, S. C. Glotzer, G.-R. Yi and D. J. Pine, Nat. Commun., 2019, 10, 1–10.
22. C. Fernández-Rico, M. Chiappini, T. Yanagishima, H. De Sousa, D. G. A. L. Aarts, M. Dijkstra and R. P. A. Dullens, Science, 2020, 369, 950–955.
23. J. De Graaf, R. Van Roij and M. Dijkstra, Phys. Rev. Lett., 2011, 107, 155501.
24. C. Avendaño and F. A. Escobedo, Curr. Opin. Colloid Interface Sci., 2017, 30, 62–69.
25. R. L. Marson, E. G. Teich, J. Dshemuchadse, S. C. Glotzer and R. G. Larson, Soft Matter, 2019, 6288–6299.
26. S. Marn-Aguilar, F. Camerin and M. Dijkstra, J. Chem. Phys., 2022, 157, 154503.
27. F. Romano, E. Sanz and F. Sciortino, J. Chem. Phys., 2010, 132, 184501.
28. E. G. Noya, C. Vega, J. P. K. Doye and A. A. Louis, J. Chem. Phys., 2010, 132, 234511.
29. E. Bianchi, R. Blaak and C. N. Likos, Phys. Chem. Chem. Phys., 2011, 13, 6397.
30. F. Romano, E. Sanz and F. Sciortino, J. Chem. Phys., 2011, 134, 174502.
31. F. Romano and F. Sciortino, Nat. Commun., 2012, 3, 975.
32. A. B. Rao, J. Shaw, A. Neophytou, D. Morphew, F. Sciortino, R. L. Johnston and D. Chakrabarti, ACS Nano, 2020, 14, 5348–5359.
33. I. Q. Matos and F. A. Escobedo, J. Phys. Chem. B, 2023, 127, 3746–3755.
34. A. Neophytou, V. N. Manoharan and D. Chakrabarti, ACS Nano, 2021, 15, 2668–2678.
35. X. Mao, Q. Chen and S. Granick, Nat. Mater., 2013, 12, 217–222.
36. J. Glaser, T. D. Nguyen, J. A. Anderson, P. Lui, F. Spiga, J. A. Millan, D. C. Morse and S. C. Glotzer, Comput. Phys. Commun., 2015, 192, 97–107.
37. F. Cui, S. Marbach, J. A. Zheng, M. Holmes-Cerfon and D. J. Pine, Nat. Commun., 2022, 13, 2304.
38. X. Wang, S. Ramrez-Hinestrosa, J. Dobnikar and D. Frenkel, Phys. Chem. Chem. Phys., 2020, 22, 10624–10633.
39. E. G. Noya, I. Zubieta, D. J. Pine and F. Sciortino, J. Chem. Phys., 2019, 151, 094502.
40. D. Frenkel and B. Smit, Understanding molecular simulation: from algorithms to applications, Elsevier, 3rd edn, 2023.
41. G. M. Torrie and J. P. Valleau, J. Comput. Phys., 1977, 23, 187–199.
42. S. Kumar, J. M. Rosenberg, D. Bouzida, R. H. Swendsen and P. A. Kollman, J. Comput. Chem., 1992, 13, 1011–1021.
43. A. Grossfield, WHAM: The Weighted Histogram Analysis Method, version 2.0.10, https://membrane.urmc.rochester.edu/wordpress/?page_id=126.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: https://doi.org/10.1039/d5sm00158g

---