A theory of entropic bonding.

Entropy alone can self-assemble hard nanoparticles into colloidal crystals of remarkable complexity whose structures are the same as atomic and molecular crystals, but with larger lattice spacings. Molecular simulation is a powerful tool used extensively to study the self-assembly of ordered phases from disordered fluid phases of atoms, molecules, or nanoparticles. However, it is not yet possible to predict colloidal crystal structures a priori from particle shape as we can for atomic crystals from electronic valency. Here, we present such a first-principles theory. By calculating and minimizing excluded volume within the framework of statistical mechanics, we describe the directional entropic forces that collectively emerge between hard shapes, in familiar terms used to describe chemical bonds. We validate our theory by demonstrating that it predicts thermodynamically preferred structures for four families of hard polyhedra that match, in every instance, previous simulation results. The success of this first-principles approach to entropic colloidal crystal structure prediction furthers fundamental understanding of both entropically driven crystallization and conceptual pictures of bonding in matter.

In his 1704 treatise, Opticks, Sir Isaac Newton wrote of an “attractive” force that holds particles of matter together. By observing that matter stays together in a host of different everyday situations, he inferred the presence of what we now refer to as chemical bonds, one of the most fundamental concepts in science. Chemical bonds originate in electronic interactions. Together with thermodynamics, chemical bonds allow us to predict and understand materials structure and properties, including the way in which atoms arrange to form crystals, from the ionically bonded simple lattice of table salt, to the covalently bonded structure of diamond, to multicomponent alloys with thousands of atoms in a unit cell held together by delocalized metallic bonding.

Like atoms, nanoparticles (NPs) can self-assemble into crystal structures of extraordinary complexity and diversity (1–6). What “force” would Newton infer from observing these colloidal crystals? From his writings, he would undoubtedly infer an attraction—a bond—between particles as for atoms, but of a different origin and at larger length, time, and energy scales. DNA linkers, ligand–ligand pairs, and protein–antibody pairs are all examples of bonding elements used to link NPs together into colloidal crystals via self-assembly. Because the self-assembled structures are often isostructural to known atomic crystals, NPs are loosely considered the colloidal analogs of atoms.

But what would Newton make of colloidal crystals self-assembled, without bonding elements, from so-called “hard” particles whose only attribute, like billiard balls, is the inability to overlap? Onsager (7), in 1949, was the first to recognize the importance of shape for equilibrium phases of hard colloidal particles (8, 9). Frenkel (10, 11) further elucidated the connection between shape and entropy. Today, we know that many reported hard-particle colloidal crystals are, remarkably, isostructural to known atomic crystals with covalent, ionic, and metallic bonding (12–19). For example, all 14 Bravais lattice types (12), Frank–Kasper (12) and clathrate structures (16), open host–guest crystals (19), quasicrystals (12, 13), and a Bergman-like crystal with a 432-particle unit cell (16) self-assembled from disordered fluid phases in computer simulations of hard, otherwise noninteracting, polyhedra. Experiments likewise observed colloidal crystallization in systems of hard particles, although the diversity of structures reported, to date, is limited by challenges in making shapes devoid of any energetic interactions (17, 20, 21).

Because hard-particle crystals have no potential energy, they are stabilized solely by entropy maximization, which, counterintuitively, is achieved when the ordered crystal has more available microstates (possible particle configurations) than the disordered fluid (7, 10, 11, 22–24). In the absence of explicit bonding elements, entropy maximization leads to the emergence of effective, attractive entropic forces with directionality that arises from particle shape and crowding. This directionality creates local valence and, ultimately, long-range order (25, 26). Such directional entropic forces (DEFs) are related to the way in which the system distributes increasingly limited local free volume upon crowding. Free volume refers to the unoccupied space into which a particle may move within a fixed arrangement of particles. A large pocket of free volume surrounding a pair of particles produces weak to no DEFs between the pair, whereas small regions of free volume produce strong DEFs. In this regard, DEFs are reminiscent of chemical bonding, where low and high local electron density correspond to weak and strong bonding, respectively. The extent to which the same complexity and diversity of structure is possible from chemical bonding and from entropic “bonding” is unexpected and demonstrates that emergent phenomena such as crystallization are agnostic to the origin of these bonds, provided the system is subject to the laws of thermodynamics. We assert that, if entropic bonds are to be considered a useful conceptual picture like chemical bonds, then entropic colloidal crystal structure prediction using ab initio methods should be possible, in the same way that atomic crystal structures are routinely predicted by (numerically) solving Schrödinger’s equation (27, 28).

Here we report a theory of entropic bonding that predicts thermodynamically preferred colloidal crystal structures directly from particle shape and a mean-field description of DEFs. Using statistical mechanics, we develop a framework for calculating entropic bonds between particle shapes that does not involve computer simulations. We test our theory on several challenging examples to show that it correctly selects the thermodynamically preferred crystal structures reported in simulations. For brevity, we limit our discussion to single-component crystals composed of hard polyhedral particles, although generalization to multicomponent crystals of arbitrarily shaped hard or patchy particles (29, 30) is foreseeable. We briefly discuss how our theory can be extended to include depletion interactions. Lastly, by comparing mathematical similarities between our statistical mechanical theory and quantum theory, we posit that theoretical ideas from quantum mechanics might be useful, if only by analogy, for understanding colloidal systems and vice versa. The success of our theory not only provides a way to predict colloidal self-assembly but also presents a potentially unifying framework that bridges our understanding of crystallization processes across scales.

Results and Discussion

Entropic crystallization builds upon two fundamental tenets of statistical mechanics (10, 11). First, thermodynamic systems evolve to the most probable state. Second, the probability of any microstate depends on the energy of that microstate. Because hard particles do not interact other than through excluded volume, there is no potential energy, and thus all microstates are equally probable. The most probable thermodynamic state is then the one with the largest number (Ω) of microstates, or, equivalently, the maximum entropy S via Boltzmann’s expression .

The starting point of our microscopic theory of entropic bonding is that hard particles optimize the sharing of free volume with neighbors to maximize the entropy of the system. Depending on a given particle configuration, the local free volume can vary considerably in both size and shape. Optimization—for example, via Monte Carlo or hard-particle molecular dynamics simulation (31)—selects configurations where the size and shape distribution of local free volume maximizes the entropy of the system. Instead of free volume, we can consider the excluded volume around a particle. Excluded volume refers to space that is inaccessible to a particle, due to the presence of another particle. As such, it is dependent on the relative orientation of the particles. While irrelevant for spheres, relative orientation is important in systems of hard shapes and, importantly, can be larger or smaller depending on how two shapes align (10, 11). A decrease in local excluded volume (e.g., of two polyhedra orienting to align their facets) creates free volume nearby. Minimizing a system’s excluded volume and maximizing its free volume are equivalent, and both maximize the entropy (10, 24). A first-principles approach to developing a microscopic theory of entropic bonding would naturally start with one or the other of these volumetric quantities. We choose to consider local excluded volume. In what follows, we present the key steps in the development of entropic bond theory. First, we introduce the ansatz of pseudoparticles (pPs) as a fictitious, mathematical construct that will enable us to quantify local excluded volume. Using this ansatz, we next derive an effective interaction between particles and pPs using mean-field theory that will ultimately represent the effective entropic attraction between NPs arising from crowding. We then derive an eigenvalue equation whose solution, when combined with the mean-field result, gives the free energy of a hard-particle colloidal crystal structure (as the energy of the optimized pP density) without need for computer simulation or thermodynamic integration (32, 33). Mathematical details may be found in .

Mean-Field Derivation of Effective Interaction

We consider a system of N hard anisotropic NPs—or, simply, particles—occupying a volume V (Fig. 1). To quantify the local excluded volume for a configuration of particles, we start by uniformly filling all empty space with pPs (Fig. 1). By definition, the fictitious pPs are significantly smaller than the NPs, interact with each other with ideal gas interactions, cannot overlap NPs, and interact with NPs via an interaction U, which we now derive. We quantify pP density using statistical mechanics in the ensemble, where the general partition function for the proxy system of N NPs and pPs is

Fig. 1.

Schematic of pP ansatz. (A) A cubic crystal of hard cubes. (B) Uniform filling of all unoccupied space with fictitious pPs (gray). (C) Energy minimization of the proxy system of pPs and cubes in B creates regions of high and low pP density that stabilize NP positions and orientations. Inset shows top-down view.

Here , k is the Boltzmann constant, and U is the total interaction energy of the system; represents a configurational integral over all possible positions and orientations of the ith particle, and represents the configurational integral over all pPs. For ease of notation, we define and , where n refers to particles and q refers to pPs. We can express the probability of a single reference particle having orientation α at position r by integrating over all pPs and over particles in Eq. . This integration giveswhere represents the configurational integrals over all pPs and particles, and we define . The interaction for the reference particle is separated from the interactions for all remaining particles and pPs, . can be decomposed into

Up to this point, the derivation is general insofar as two-body interactions are concerned. and represent the summation over all pairwise repulsive and attractive interactions, respectively, between the reference particle and other particles. For hard particles, , and is when particles overlap, and zero otherwise.

The quantity is the mean-field interaction potential of the reference particle with the equilibrium distribution of all pPs. By definition, is a function only of the reference particle’s position and orientation, and can be moved outside the integrals in Eq. . Plugging Eq. into Eq. giveswhere the primes in Eq. indicate that the integrals are performed over only configurations without overlapping particles. We next recognize that the configurational integrals in Eq. give the partition function of the system less one particle, , and thus . To further simplify , we imagine that the reference particle is removed, and the system is compressed to remove the void left by the missing particle. This assumption allows us to approximate as (34). Here, is the volume contributed by the wth differential volume element of the removed particle, given configuration γ. Each wth element represents a radial layer moving away from the center of a reference particle. C serves as a normalization constant and will be absorbed into the normalization factor for . Plugging into Eq. , it is straightforward () to derive for the reference particle (dropping the implicit γ dependency for notation brevity),

where μ is the chemical potential of a particle, and is the pressure exerted on the wth volume element (). We give an alternative derivation of Eq. by explicitly maximizing the entropy using a lattice model, in .

A derivation similar to the derivation of Eq. from Eq. can be also performed for a reference pP, givingwhere μ is the chemical potential of a pP, and other terms have the same meaning as in Eq. . To develop a functional form for the pP–NP interaction potential U, we take Eq. , rearrange for π, plug π into Eq. , and solve for U (). For a given configuration γ, setting w to be the distance r between a pP and an NP (in units of NP in-sphere diameter σ) yields, finally, Eq. for the pP–NP interaction U,

where U represents a hard-core repulsion preventing pPs from overlapping NPs, which arises in the pP equation analog of Eq. . The mean-field interaction between NPs and pPs describes, in effect, the cohesive energy holding an NP in place within a given configuration of particles (Fig. 2).

Fig. 2.

Functional form of pP–NP potential. The pP–NP interaction becomes stronger for increasing μ, indicating high pP localization.

From Mean Field to Eigenvalue Equation

Eq. can be inverted to give , the spatially dependent, mean-field density distribution of the pPs. This quantity serves as a proxy for the excluded volume. To calculate the optimal pP distribution for a given configuration of particles, we write the convection–diffusion equation for a system of N particles and pPs. This equation is also known as the Smoluchowski equation, a special case of the Fokker–Planck equation in statistical mechanics describing diffusion in the presence of a driving force. The functional form of this equation, given a set of N particle orientations and positions and the interaction potential U in Eq. , is

Solving Eq. by separation of variables transforms the operator in the Smoluchowski equation into a Hermitian, Schrödinger-like Hamiltonian operator (35, 36). Taking the equilibrium, steady-state limit immediately results in an eigenvalue equation for ρ,

The eigenvalues E obtained from solving Eq. can be interpreted as inverse equilibration times for pPs to move through the system—that is, relaxation frequencies associated with the dynamical modes of the pP density distribution. The lowest eigenvalue corresponds to the lowest-frequency mode, or longest equilibration time (i.e., pPs are more localized), corresponding to the equilibrium pP distribution.

Because the operator acting on ρ in Eq. represents the Hamiltonian for the proxy system of N particles and pPs, the eigenvalues E may also be interpreted as energies corresponding to the various dynamical modes of the pP density distribution for a given configuration of particles. These energies in the proxy system are equivalent to free energies in the original system of hard particles. Comparison of the computed eigenvalues for NPs arranged in different crystal structures should predict the thermodynamically preferred crystal as the structure with the lowest eigenvalue of all tested crystal structures. We test this assertion in the next section.

Solving the Eigenvalue Equation

Because Eq. is an eigenvalue equation, we can solve it using numerical methods found in quantum mechanical (QM) codes (e.g., ref. 37) that solve the Schrodinger equation (and its variants, such as the Kohn–Sham equation) for energies of atomic crystal structures. We first write the pP probability density as a series of “shape harmonics” in a reference frame conducive to the particle shape, similar to the spherical harmonics used in QM (). The resulting shape harmonics are functions of two angular components. We call these shape harmonics “shape orbitals” in analogy with the correspondence between spherical harmonics and atomic orbitals. The lowest shape orbitals for a cube, tetrahedron, dodecahedron, and hexagonal prism are shown in Fig. 3 , respectively. We then take a linear combination of shape orbitals [LCSO, analogous to a linear combination of atomic orbitals, LCAO, used in QM codes (28)] as a trial solution to solve for the optimal configuration of particles. The solution entails optimizing orbital overlap for a given arrangement of NPs, considering all possible relative orientations of NPs. Optimizing orbital overlap minimizes the energy of the proxy system and gives the equilibrium distribution of pP density. Because of the relationship between pPs and excluded volume, this optimization equivalently maximizes entropy in the original system of hard particles. Details of the calculation are given in Materials and Methods.

Solutions of Eq. are expected to produce pP distributions with regions of high and low density that correspond, respectively, to low and high local excluded volume around the NPs (Fig. 1). The pPs in high-density regions of the pP distribution can be considered localized, and the pP chemical potential μ provides a physical interpretation of localization. A small μ implies easy placement of new pPs, expected in the dilute NP limit where entropic forces are weak and isotropic. In this limit, the local excluded volume of an NP is maximal and equal to the NP circumsphere volume, as reflected in weak pP localization. Conversely, μ increases as NPs become crowded, and it becomes harder to accommodate another pP. As a result, the pP density distribution develops strong localization, reflected in deepening energy wells for high μ (Fig. 2).

Below, we solve the eigenvalue equation for a single pair of particles within a crowded system of N particles to define the entropic bond. Subsequently, we show how the eigenvalue equation can be applied to crystal structure prediction.

Entropic bonding

For a pair of polyhedral NPs, solving Eq. predicts the most stable configuration in a crowded system of polyhedra as the configuration favoring face–face alignment, in agreement with previous simulations (38). The predicted “bonding orbitals” between face-aligned polyhedra are visualized in Fig. 3 . Additionally, we find energy minima as a function of center-to-center separation between the two polyhedra (Fig. 3 ). These minima correspond to “bond strengths” for face-aligned tetrahedra, cubes, and hexagonal prisms ranging between 1.5 and 1.6 , with bond lengths of 1.1σ, where σ is the polyhedron in-sphere diameter. Dodecahedra exhibit a stronger bond of 2.3 due to additional orbital overlaps between small facets. The calculations employed in Fig. 3 provide a quantitative description of an entropic bond in terms of pP local density and are extensible to particles of arbitrary shape.

Fig. 3.

Shape orbitals and entropic bonds. Shape and bonding orbitals for (A) cube, (B) tetrahedron, (C) dodecahedron, and (D) hexagonal prism. (1) Visualization of the shape orbitals defined by angular expansions of the local pP density distribution. (2) The theoretically predicted lowest-energy pairwise configurations (both position and orientation). (3) The bonding energy as a function of center-to-center separation distance at the predicted most stable orientation for a pair of polyhedra. All polyhedra prefer face–face alignment with bond lengths of 1.1σ, where σ is the in-sphere diameter of each respective shape. Cubes, tetrahedra, and hexagonal prisms exhibit bond energies of 1.5 kT to 1.6 kT. For the pair of dodecahedra, secondary stabilization (from local facets near the main face–face contact) lowers the bond energy to 2.3 kT. The LCSO approximation employs the lowest shape orbital for each shape.

Crystal structure prediction with entropic bond theory

To illustrate the application of entropic bond theory to colloidal crystal prediction, we begin with a simple example and compute the optimal arrangement of hard cubes at particle volume fraction η = 0.55 in three lattices—body-centered (BCC), face-centered (FCC), and simple cubic (SC). Fig. 4 shows the computed lattice energy of formation, indicating that cubes are more stable when arranged in an SC lattice relative to BCC and FCC, as expected (14, 39, 40). Visualization of the bonding orbitals (Fig. 4 ) helps explain why. Cubes in an FCC arrangement suffer a deficit of two face–face bonding orbitals per cube that cannot be offset by gains in edge–edge bonding with neighbors. Similarly, the center cube in the BCC lattice gains 12 edge–edge bonding contacts but loses all six of its face–face bonds. Since edge–edge bonding is higher in energy than face–face bonding, both FCC and BCC lose out to SC.*

Fig. 4.

Crystal predictions for cubes. (A) Energy per cube within each respective lattice. Visualization of cubes with bonding orbital in (B) SC, (C) FCC, and (D) BCC. The LCSO approximation employs the first shape orbital for the cube. Error bar determination is described in Materials and Methods.

Entropic bond theory can predict solid–solid transitions as particle shape is modified. Consider the truncated tetrahedron shape family (43). Simulation studies reported the self-assembly of a series of crystal structures, shown in Fig. 5, as the corners of a regular tetrahedron are truncated (full truncation produces a regular octahedron) (43). Fig. 5 shows a direct comparison between the lowest-energy crystal structure predicted using our framework (colors in Fig. 5 match the corresponding unit cells in Fig. 5) with the transitions reported in Damasceno et al. (43) shown in black lines, indicating excellent agreement between theory and simulation.

Fig. 5.

Crystal predictions for shape-driven solid–solid transition of truncated tetrahedra. (A) Snapshot of thermodynamically preferred structure with relative orientations between truncated tetrahedra predicted from theory. Colors correspond to the same lattice in C. (B) Lattice energy per particle used for phase diagram. We additionally computed the SC and rotator crystal to show that other structures are not stable within this regime for the vertex-truncated tetrahedron shape family. (C) Phase diagram for vertex-truncated tetrahedron shape family. Black lines indicate transitions between respective regions found in simulations (43). With increasing truncation, the stable phase transitions from QCA to diamond to β -Sn to BCC. The LCSO approximation employs the first shape orbital for each polyhedron.

Our theory can also predict solid–solid transitions for a given shape as a function of particle volume fraction. We consider a system of regular trigonal bipyramids (TBPs) for which a thermodynamic transition from a dodecagonal quasicrystal approximant (QCA) to a triclinic dimer lattice with increasing volume fraction was reported in simulations (44). Representative snapshots of TBPs within the QCA and dimer lattices are shown in Fig. 6 , respectively. Entropic bond theory predicts the reported transition (Fig. 6, black horizontal line in vertical color bar, far right) at a particle volume fraction of 0.78 ± 0.01 as compared to the reported value of 0.79 ± 0.008 (44). Visualization of bonding orbitals for a TBP within the approximant and dimer phases indicates stark differences in pP distributions (Fig. 6). To understand which bonding orbital pattern leads to the lower energy above and below the transition, we examine the bonding energy of a reference TBP and a neighboring TBP, at varying particle volume fractions. The bonding energy curves reveal a distinct change in shape at the transition volume fraction that drives the shift from the dimer lattice to QCA (Fig. 6).

Fig. 6.

Crystal predictions for complex lattices. (A–D) Lattice prediction of solid–solid transition for TBPs. Snapshot of (A) 41-particle unit cell of QCA at η = 0.79 and (B) seven two-particle unit cells of the trigonal dimer packing at η = 0.79. (C) Calculated bonding energy of a reference TBP and a neighboring TBP a distance r away, at varying particle volume fractions η. At the previously reported (44) thermodynamic transition from the QCA to the dimer phase, undergoes a sharp transition from a single-well to a double-well function. The vertical color bar (far right) shows the lower-energy structure for TBPs determined from lattice energy calculations as a function of η (colors correspond to A and B unit cells). The reported transition from QCA to dimer phase obtained from Monte Carlo simulation is show as a black line (44). (D) Visualization of the entropic bonding orbitals around a reference TBP in QCA (Left) and dimer (Right) phases just below and above the solid–solid transition. The LCSO approximation employs the first two shape orbitals for every TBP. (E and F) Local bonding environment prediction for a complex and multilayered (oF244) lattice (16). (E) Representative snapshot of oF244 lattice. (F) Layer-specific (L1, L2, L3) energy per PBP in the oF244 lattice. (G) While each PBP has similar energy within each layer, visualization of the entropic bonding orbital at the same energy contour reveals that the local bonding environments differ significantly among L1 (Left), L2 (Center), and L3 (Right). Top and Bottom in G show views from the top and side, respectively, of the PBP. The LCSO calculation employs the first shape orbital for PBP.

Complex atomic crystals have multiple Wyckoff positions, meaning that different positions within the crystal possess different local environments. The same is true for colloidal crystals, and entropic bond theory predicts this. We consider a system of hard pentagonal bipyramids (PBPs) that, in simulation, self-assembled from a fluid to a complex, multilayered crystal with a 244-particle unit cell (oF244) (16). Fig. 6 shows a snapshot of the PBPs within the oF244 lattice, colored according to their embedding layer L1, L2, or L3 within the crystal (Fig. 6). The energy in the system is partitioned equally among the layers by adopting a different pP density distribution. Visualizations of the bonding orbitals show strong pP localization on the PBP faces in L1 (Fig. 6 , Left), while pPs localize in a ring around the center of each PBP in L2 (Fig. 6 , Center) and are much less localized in L3 (Fig. 6 , Right). These results elucidate how diversity in local bonding environments can promote crystal complexity.

Conclusion

We have presented a theory of entropic bonding that can be used to predict colloidal crystal structures of hard particles a priori. We formulated the theory using a mathematical ansatz in which fictitious pPs are used to approximate excluded volume in a system of NPs. A mean-field derivation of a pP–NP interaction was combined with a convection–diffusion equation, resulting in an eigenvalue equation whose solution for a given colloidal crystal structure gives the free energy of the structure as the energy of the pP distribution. By comparing solutions among many structures, we can select for the one with maximum entropy. We showed that our framework predicts an effective entropic interaction between hard polyhedra that we quantify in terms of regions of high pP density and describe as an entropic bond. Our theory correctly predicts the thermodynamically preferred colloidal crystal structures for four different families of hard polyhedra as a function of shape or particle volume fraction. It also shows how different crystal structures—and different local environments within a crystal structure—facilitate entropy maximization through visualization of bonding orbitals.

Entropic bond theory is different from, but complementary to, two conventional approaches for determining the thermodynamically preferred colloidal crystal structure of hard NPs—molecular simulations and Frenkel–Ladd free energy calculations (14, 32, 33, 40). Molecular simulation finds maximum entropy solutions (provided kinetics allows) by sampling millions of microstates in pursuit of the most probable state. Frenkel–Ladd calculations are used to compute free energies for given test crystals, but require simulation: Monte Carlo simulations are employed to compute an average free energy as a function of a given thermodynamic integration coupling parameter. Our theory requires no simulation. Rather, it employs an iterative, self-consistent solver akin to those used in electronic structure calculations to solve an eigenvalue equation for the entropy of a crystal structure.

Entropic bond theory shares some features with chemical bonding theory, where significant theoretical developments over many decades have enabled calculation of ground-state energies by optimizing the spatial distribution of electron density throughout an atomic crystal (37, 45–47). In entropic bond theory, calculation of entropy in a colloidal crystal is performed by optimizing the spatially varying density of a distribution of fictitious pPs that serve as a proxy for excluded volume. The pP density optimization is realized by maximizing the overlap of pP shape orbitals, which stabilizes the distance and orientation between neighboring particles. By extension, solving for the optimal spatial distribution of pPs to find the most probable crystal structure of colloidal particles is akin to solving for the optimal spatial distribution of electrons to find the ground-state crystal lattice of atoms. In other words, pPs act, mathematically speaking, as entropic analogs to electrons for entropically driven self-assembly. Whereas quantum mechanics is the language of chemical bonding, statistical mechanics is the language of entropic bonding. Table 1 provides an overview of these analogies.

Table 1.

Chemical vs. entropic bonding

	Chemical bonding	Entropic bonding
Core	Atomic nucleus	Hard shape
Mediator	Electrons	pPs
Attraction	Nuclear–electron	Shape–pP
Repulsion	Nuclear–nuclear	Shape–shape
	Electron–electron	—
Orbitals	Atomic (s, p, d, f)	Shape
Bonding	Orbital hybridization	Orbital overlap
Governing equation	Eigenvalue equation for electron density	Eigenvalue equation for pP density

The development of entropic bond theory presented here relies on an unusual ansatz where excluded volume is approximated by fictitious pPs. This mathematical device provides a quantifiable proxy for excluded volume, which facilitates the development of the theory. Explicit molecular dynamics simulations (48–50) of a system of NPs and pPs interacting via the derived effective potential (Eq. —while not necessary for crystal structure prediction—provide additional validation of the reasonableness of this ansatz ().

Our treatment of entropic forces in systems of hard particles is different from a theoretical treatment of depletion interactions, which are also entropic in origin. Although it may be tempting to interpret pPs as depletants, this would be incorrect unless one considers a purely repulsive U (negligible μ; Fig. 2). For that case only, pPs mimic “hard-sphere” depletants (24, 51), and our theory provides a unique dynamical interpretation of the classical Asakura–Oosawa depletion interaction (51). Critical to our theory of entropic bonding is the attractive nature of the NP–pP interaction U, which serves as a proxy for the effective entropic interaction between NPs arising from crowding in the absence of depletants. As considered here, pPs are not depletant-like as they are attractive to NPs, and become increasingly attractive as the particle packing fraction (and thus crowding) increases (high μ; Fig. 2), driving crystallization. An extension of entropic bonding theory to include depletion simply requires the addition of depletant particles as a second species of pP, one that interacts with NPs solely via an excluded volume interaction (). Such a treatment of depletion interactions can be considered an “implicit depletion model,” as implemented in ref. 52.

While predictive, the theory presented here suffers from the same limitations faced by electronic structure theory—namely, the need to guess lattices for comparison. For ground-state atomic crystals, where electronic structure theory is most useful, unit cells are comparatively small and simple. Only at nonzero temperatures do atomic crystal lattices typically become complex. For colloidal crystals, the same is true when comparing putative densest packings [“ground states” (53–55)] against the complex crystal structures observed away from infinite pressure. It is the latter we are most interested in, and guessing lattices may prove unwieldy.

Regardless, that there exists a first-principles theory based solely on entropy that, like electronic structure theory, is capable of predicting preferred crystal structures is, in itself, important and profound. For one, it confirms that generalization of the term “bond” to describe the DEFs (25, 56) arising in crowded systems of hard polyhedra is mathematically appropriate and physically meaningful (26). Second, it helps explain why similar crystal structures can be observed on length scales orders of magnitude apart by demonstrating that crystallization and self-assembly are agnostic to the origin of the forces between building blocks.

Materials and Methods

See for additional details.

Eigenvalue Solver

To solve Eq. , we set up a framework that determines the lowest eigenvalue for a given set of particle configurations. We first guess a trial solution for the pP density ρ via a similar LCAO approach employed in quantum mechanics. Here, since we combine shape orbitals (), we term the approach LCSO.where ρ is the guessed solution, is the sth shape orbital for the ith particle, and N is the total number of particles in the system. Eq. is then substituted for ρ in Eq. to give

Eq. can be converted into a series of N linear, homogeneous equations where each ith equation is defined by multiplying Eq. by . The resulting set of equations is solved using linear algebra. We then iteratively vary the set of C until the minimum eigenvalue is obtained for a given configuration (center-to-center distance and relative orientation). The calculation is then repeated across different center-to-center distances and relative orientations until the lowest eigenvalue overall is obtained. We employ SE propagation analysis to determine the numerical uncertainties in our implementation of an iterative eigenvalue solver for Eq. .