Inverse design of multicomponent assemblies

The synthesis and fabrication of materials with specific sub-mesoscale features or architectures is a desirable but challenging design goal. While some top-down techniques allow material patterning at these scales (e.g., lithography), they are often prohibitively slow and expensive processes for large-scale industrial applications.^1–3 Self-assembly of colloidal systems, on the other hand, presents one promising bottom-up alternative whereby particles are designed to organize into specific configurations determined by their effective interparticle interactions.^4–7 This approach, though still in its infancy, is attractive due to the possibility to systematically tune the effective interparticle interactions^8–10—including those determined by particle shape^11–14 or “patchiness”^15–20—to drive self-assembly. It is also possible to expand the available design space for assembly by employing multicomponent systems (e.g., DNA-grafted nanoparticles,^21–23 binary systems of charged colloids,^24–26 and others^27–30) to augment the diversity of possible interactions and structures.

Given the large potential parameter space available, systematic design strategies that allow one to explore how the properties of constituent particles relate to the resulting self-assembled structures are needed. To this end, one can consider forward design methods whereby the phase behavior of a system with known interactions is mapped by systematic variation of a few characteristic parameters. Alternatively, inverse design methods can be employed, where a desired particle assembly is explicitly targeted, and the parameters necessary to achieve it are found by a way of solving a constrained optimization problem.^{13,14,31–39}

A classic inverse design problem is the optimization of decision variables θ for a system of identical particles interacting via a given isotropic pairwise interaction potential, ϕ(r; θ), to stabilize a desired structure. Here, r is the distance between particle centers and θ are the parameters required for the pair potential. By minimizing either the ground state energy of the ideal configuration (relative to competing structures) or the free energy of an associated configurational ensemble at a higher temperature (relative to competing phases), researchers have successfully found isotropic interactions that stabilize a wide variety of structures and phases including, for example, two-dimensional honeycomb^40–42 and kagome^43–45 lattice assemblies as well as three-dimensional simple cubic⁴⁶ and diamond^46–48 crystals. Isotropic pair potentials that stabilize more exotic phases have also been discovered via recently introduced inverse design strategies.^49–51 An analytical reformulation⁴⁹ of the ground-state optimization problem⁴⁶ allowed for the efficient design of potentials that stabilize two-dimensional snub-square⁵² and truncated hexagonal⁵³ lattices. Furthermore, a strategy adapted from the bio-molecular coarse graining community^54–58 called relative entropy (RE) minimization^59,60—which enables “on-the-fly” optimization of the potential parameters directly within a molecular simulation—led to the discovery of isotropic interactions that promote self-assembly of a variety of two and three dimensional crystals and quasicrystals^50,51,61 as well as clustered fluids and porous mesophases.^51,62,63

While much of the earlier inverse design work focused on the wealth of structural possibilities for single-component systems, similar demonstrations in multicomponent systems are surprisingly lacking. Multicomponent systems are of particular interest given the expanded design space afforded to them through their additional self- and cross component interactions (i.e., additional degrees of freedom). Moreover, the ability to partition a desirable target structure into a larger number of components suggests it may help simplify the particle interactions required for a given structure relative to those of an equivalent single-component system. Indeed, multicomponent systems are known to stabilize entirely new and exotic structures which are inaccessible to single-component systems. For instance, it is known that moving from single-component to binary colloidal systems significantly diversifies the list of crystal structures and motifs observed in experiment.^23,64,65 Similarly, recent theoretical and computational efforts demonstrate that binary mixtures can produce a wide variety of novel structures, some of which exhibit highly specific local orderings.^66–70 

Clearly, multi-component systems hold significant promise as tailorable building blocks for novel structural architectures, but a design strategy that allows for rational design of underlying interactions remains to be demonstrated. Here, we adopt one such strategy and extend the RE optimization method for inverse design of single-component materials interacting via isotropic interparticle potentials to multicomponent systems. By independently optimizing the self- and cross interactions for species in a binary two-dimensional mixture, we demonstrate that this methodology can successfully discover pair potentials that readily self-assemble a variety of simple (intercalated square and triangular) lattices as well as more challenging target phases (stripes and highly open structures), many of which had not yet, to our knowledge, been observed to self-assemble for such systems. For select cases, we further compare the optimized binary interactions to those obtained from optimizations for the target structures assuming single-component systems. These comparisons help us to understand the trade-offs associated with using binary systems for self-assembly as well as the roles of self- versus cross interactions in stabilizing various target lattices.

The remainder of this paper is organized as follows. We elaborate on the extension of the RE optimization strategy to multicomponent systems and explain how our binary crystal target structures are chosen in Sec. II. In Sec. III, we illustrate the optimized pair potentials and demonstrate they successfully promote self-assembly of all targeted structures. We also discuss for select cases how binary interactions compare to those of analogous single-component systems. We conclude in Sec. IV by summarizing the similarities and differences between the binary and the single component interaction potentials, focusing specifically on how these observed features enable formation of the targeted structures.

The RE course graining approach is a probabilistic optimization method that addresses an inverse design problem for self-assembly by systematically tuning the interparticle interaction potential U(r|θ) in a system of particles via parameters θ ≡ {θ₁, …, θ_m} to maximize the likelihood of forming a desired structure. One advantage of the approach is that it can perform “on-the-fly” optimization of particle interactions during the course of a simulation in order to promote thermodynamic stability (and, with a judicious choice for the simulation protocol,⁵⁰ kinetic accessibility) of the assembled target phase. Although RE optimization was originally applied to single-component systems,^50,51 it can be readily extended to multicomponent materials. For a two-dimensional binary system with components A and B, we partition the energy and parameters in terms of self- and cross interactions as $u(r|θ(k,k′))$⁠, where k, k′ = A or B. Parameters $θ(k,k′)$ are then updated with a gradient ascent procedure as

where i indicates the iteration step, $g(r|θi(k,k′))$ represents the simulated radial distribution function (RDF), $gtgt(k,k′)(r)$ is the target RDF, $[∇θu(r|θ)]θ≡θi(k,k′)$ denotes the gradient of the pair potential, and $α(k,k′)$ is a learning rate chosen to ensure simulation stability and convergence. While it may be possible to choose a single value of the learning rate parameter that is effective for all interaction types, we adopt independently tuned values for the three interaction types here to attain faster convergence. We define $u(r|θ(k,k′))$ as a set of Akima splines whose knots are computed from $θ(k,k′)$ as reported in Ref. 50. Finally, the knot amplitudes are restricted to increase with decreasing values of r to ensure a monotonically decreasing (i.e., repulsive) pair potential. A derivation of the update expression presented in Eq. (1) is provided in the supplementary material.

Briefly, an RE optimization is carried out as follows. An iteration RE step starts by simulating a high temperature fluid interacting through pair-potentials $u(k,k′)(r|θ)$ and cooled slowly to a prescribed optimization temperature T*. The system is allowed to equilibrate and RDF statistics are collected. Using the differences between the calculated and target RDFs, potential parameters are updated as per Eq. (1) and this resulting potential used to simulate the high temperature fluid system where the entire process is repeated in the next iteration step. For implementation details, please see Subsection II C.

A tiling can be described by the vertices formed by the underlying polygon tiles and denoted as [n₁.n₂…], where n_i denotes the numbers of sides for the polygons that meet at each vertex.⁷¹ For instance, in the case of a square lattice, each vertex is the meeting point of 4 squares tiles; in vertex notation, this lattice is denoted by [4.4.4.4] or [4⁴] for short. Tilings consisting of k vertex types are said to be k-uniform and denoted similarly as [$n11.n21$⁠...; ...; $n1k.n2k…$]. For a crystal composed of particles at the vertices of regular polygon tilings, the number of non-equivalent origins in the crystal then corresponds to the number of k vertices necessary to create the crystal. In this work, we only consider target lattices that are 1- and 2-uniform tilings in order to guarantee a single equivalent origin for each component in the binary mixture. In addition to regular polygon tilings, we consider a tiling consisting of an octadecagonal polygon and a star polygon with 3 corners with an internal angle of 2π/9–denoted as [18².3_2π/9]. Our chosen targets are shown in Figs. 1(a)–1(f) in a formal vertex notation. Note that some realized crystals are based on a single underlying vertex (e.g., square [4⁴]) but with the components arranged so as to create new structural motifs (e.g., intercalated rows of stripes) while still conserving origin equivalency. In such a case, the actual vertex seen by the individual component may not be regular, but together with the other component still span the original vertex of the tiling.

The RE optimization algorithm described above can be interfaced with any standard simulation engine, and here we use the GROMACS 5.0.6 molecular dynamics package.^72,73 We define the component masses as m = m^(A) = m^(B) = 1 and the nearest same-component neighbor distance as σ = σ^(k,k) = 1. The energy scale is defined by $β=(kBT)−1$⁠, where k_B is the Boltzmann constant, and we take the effective temperature T* = (β)⁻¹ to be unity at the optimization temperature. Briefly, a binary system of particles of type A and B interacting via independent self- (AA and BB) and cross (AB) monotonic repulsive spline pair potentials is simulated in the canonical ensemble using a periodically replicated rectangular simulation cell with the aspect ratio chosen to accommodate the target lattice. The corresponding simulation time step is set to dt = 0.001 for all runs. The number of A and B particles, N^(A) and N^(B), respectively, is chosen to match the ideal target configuration stoichiometry N^(A) : N^(B) and kept to a combined total particle number of approximately N = N^(A) + N^(B) = 1000 for all targets. The system is initiated at a high temperature (⁠$Th*$⁠) and annealed to T* = 1 over a span of 5 × 10⁶ time steps. For triangular binary or similar compact structures, $Th*=3$⁠, while for more open or square-based structures, $Th*=1.5$ was sufficient to ensure melting. Radial distribution function $g(k,k′)(r)$ statistics are collected for an additional 10⁶ time step after cooling to T* = 1. Target RDFs are generated by tethering particles to ideal crystal positions using a quadratic restoring potential with spring constants chosen such that peaks are sharp but integrable.^50,51 Appropriate spring constants were generally in the range of 10-1000; larger values were required for the more open targets. Using the simulated structural data, the spline potential was updated per Eq. (1), and the process was iterated in this manner. In practice, $α(k,k′)=0.005–0.015$ is sufficient for all crystal target structures considered here. Convergence is typically achieved in about 80-150 iterations, and optimization is considered complete once the self-assembled crystal remains stable up to the high temperature limit $Th*$⁠. For a description of the scheme used to determine interaction cutoff radius in our systems, see the supplementary material.

Following the RE optimization protocol elaborated in Sec. II, isotropic interparticle potentials for binary mixtures that successfully self-assembled each of our nine binary crystal targets were designed. The optimized potentials and assemblies are shown for systems featuring triangular and honeycomb motifs (Fig. 2), squares and stripes (Fig. 3), and other more complex, open structures (Fig. 4). Most of the optimized potentials yield excellent particle assembly including the very open and intricate “square truncated hexagonal” (STH) or “octadecagonal star binary” lattices presented in Fig. 4.

Considering the designed interactions and the associated assembled structures shown in Figs. 2–4, two general observations can be made. First, targets featuring equivalent component sites yield comparable interactions for A and B components as should be expected based on symmetry. This can be seen, for instance, in the square binary or square stripe structures shown in Fig. 3 (top or middle, respectively), where exchanging component identities would yield near-identical configurations. Second, optimized AA, BB, and AB pair potentials become longer ranged and exhibit more features (i.e., shoulders and plateaus) as the target structures become more open and complex. Compare, for instance, potentials for square binary in Fig. 3 (top) with those of STH in Fig. 4 (middle). The former are clearly short ranged and feature a single shoulder, while STH interactions are longer ranged and exhibit multiple plateaus.

We also find evidence for two advantages of binary mixtures over single-component systems for self-assembly: (1) the optimized pair interactions of a binary system can each be simpler than an analogous single-component interaction designed to stabilize the same global (or overall) structure and (2) the expanded parameter space of a binary system can help to stabilize a richer variety of self-assembled structures than single-component systems. To illustrate the former, we consider the triangular honeycomb structure in Fig. 2 (bottom). The underlying tiling is 2-uniform as illustrated in Fig. 1(c). Therefore, the global structure of the lattice can be naturally partitioned into two simpler sub-lattices: a honeycomb lattice (red particles) and a triangular lattice with a side length of three (blue particles). By allowing two components to occupy the two distinct types of lattice sites, relatively simple interactions can be combined to favor self-assembly of a rather complex target. By contrast, when the same global structure was targeted via a single-component optimization, the best resulting interaction was not only considerably more complex, but also self-assembly of the target structure was significantly less satisfactory than for the binary mixture (see Fig. S1 of the supplementary material). As a second example, a previous study reported⁴⁵ that single-component assembly of a rectangular kagome lattice from a single-component system required an interaction possessing an abrupt attractive well within a larger repulsive profile, while the optimized binary interactions reported here for the same structure are purely repulsive and limited to a few shoulder features; see Fig. 2 (middle).

With respect to binary mixtures increasing the diversity of possible self-assembled structures relative to single-component systems, the partitioning of a desired global structure into individual sub-lattices corresponding to each species not only allows for self-assembly of intricate structures such as STH and octadecagonal star binary lattices (Fig. 4), but it also opens up the possibility of segregating components in specific local structures using an otherwise simple global lattice structure. This is seen clearly in Fig. 3 where the same underlying square lattice is partitioned to form intercalated squares (top) or stripes with either a 1:1 (middle) or a 2:1 (bottom) ratio. As expected, we find that increasing the asymmetry of the particle arrangement (as is displayed from top to bottom in Fig. 3) generally requires more complex potentials in terms of both the interaction range and the number of shoulders. A similar argument applies to the triangular binary target, where it is known that a single-component hard-core fluid favors a global triangular lattice, but achieving assemblies possessing the specific relative AA and BB ordering characteristic of the target lattice requires the binary interactions shown in Fig. 2 (top). In short, the above observations indicate that the enhanced design space of binary systems can enable the self-assembly of structures also available to single-component systems (but with simpler interactions) as well as the self-assembly of significantly more complex structures than those attainable in a single-component system.

Focusing on binary assemblies, we can also gain some insight into the mechanisms for their global stabilization by comparing the underlying component interactions to those of single-component systems optimized to stabilize similar local structures. Specifically, Eq. (1) recognizes that the binary system can comprise independent self- and cross interactions that are optimized to recreate the RDFs of the target lattice. However, the equation alone does not clarify the extent to which the binary assembly can be considered a trivial superposition of sub-structures that can be stabilized by component self-interactions (e.g., approximately equal to those expected from the analogous single-component target structure) versus a more cooperative, or coupled, assembly relying on nontrivial cross interactions. In order to help address this question, we compare the optimized self-interactions of the binary system to analogous single-component interactions that stabilize the same local structure.

We begin by considering the square binary, triangular binary, and triangular honeycomb targets shown in Fig. 3 (top) and Fig. 2 (top) and (bottom), respectively. The underlying single-component structure is that of a square lattice for square binary AA or BB components and a honeycomb lattice for the BB component of both triangular binary and triangular honeycomb structures. As such, we can plot these individual component interactions and compare them to single-component interactions known to stabilize these lattices.⁵⁰ This comparison is shown in Fig. 5 for the optimized square lattice (a) and the honeycomb lattice (b) interactions where potentials are normalized such that ϕ(r)/ϕ(1) = 1. As seen, the square binary AA(BB) interaction is largely similar, though not identical, to its single-component equivalent in range and complexity. The BB interactions for the triangular honeycomb and triangular binary [Fig. 5(b)] also approximate those of the reported honeycomb potential, though the triangular binary BB interaction deviates more strongly from the single-component results. This larger discrepancy, as we will show below, is an indication that individual component interactions need not exactly match their ideal local target lattice counterpart in order to achieve proper global assembly.

To investigate the above deviation more closely, we use one of the self-interactions (AA or BB) from the optimized binary system to carry out a single-component assembly simulation at the same temperature and box size as the binary system, where the other component has been removed from the simulation box. We then compare the resulting equilibrium assemblies to the expected perfect local lattice for that component. To this end, we choose the simpler square binary structure (Fig. 3, top) and the more elaborate STH structure (Fig. 4, middle) as two contrasting test cases. As seen in Fig. 6(a) (top right), the square lattice is the local structure for the A(B) component in the square binary target. However, as shown in the bottom right of (a), particles interacting via the AA interaction form a largely amorphous configuration. Despite this, the corresponding RDF for the extracted binary interaction (blue) shows a good match between the first two target square lattice coordination shells positions (note that, for this case, the AA pair interaction only spans the first two target coordination shells).

Similarly for the STH target structure, the A component locally forms a stretched truncated square lattice. However, particles interacting with the AA potential in a single-component simulation assemble into stripe-like phases. Nonetheless, comparing RDF peak positions shows good agreement with the target in the first three coordination shells which spans the full interaction range. Finally, looking at the STH BB component assembly in (c), we expect square clusters in a square super-lattice arrangement, but simulations yield rhomboid clusters in a triangular super-lattice orientation instead. While the peak positions in the RDFs do not overlap as closely as in Figs. 6(a) and 6(b), visual comparison of the target and the formed assembly shows that the mean inter-cluster distance is approximately the same in both cases. Together, these results demonstrate that individual component interactions cannot be expected to fully recreate the underlying local target structure in the binary target but only to ensure proper local particle positioning on average. This is why the optimized interactions need not (and generally will not) closely match corresponding single-component target interactions; the latter fail to fully stabilize the local target structures and more specific coupling between interaction types must play a key role.

To further unravel the requirements for assembly of binary structures, we performed optimizations where one of the three interaction types was fixed as a hard-core-like [Weeks-Chandler-Andersen (WCA)] potential, while the rest were optimized as usual.⁷⁴ In particular, we consider the square binary target from Fig. 3 (top) and carry out optimizations controlled for annealing schedule and iteration step so as to isolate the effects of fixing the interaction relative to the fully optimized system. Results for the square binary structure, where either AA(BB) or AB interactions are fixed, are shown in Fig. 7 (top or bottom, respectively). In the first case, fixed AA interactions resulted in stronger AB and BB interactions that helped boost global structure stability—though clearly not as efficiently as the fully optimized system (compare AA RDF, top right). This contrast is more drastic when AB (the cross coupling) is fixed, resulting in sharpened AA and BB interactions that nevertheless fail to restore original system stability (broadened peaks for all RDFs, bottom right). Fixing one of the individual interactions as hard-core-like and optimizing the others does not work at all in assembling in the more complex binary structures. For instance, when such a procedure was carried for the STH lattice, it resulted in phase separation for fixed AA or BB interaction or stripe-like configurations for fixed AB interaction (see Fig. S2 of the supplementary material). Similar results hold for the square corral structure.

Together, these results help us to highlight the individual roles of self- and cross interactions, respectively, in binary assemblies. For simple global targets like the square binary lattice, stronger AB and BB optimized interactions can compensate for a hard-sphere-like AA interaction to result in a successfully self-assembled target structure. However, for more intricate local orderings like those in the STH lattice, sharpening remaining interactions is not enough to compensate for a neutral interparticle interaction, resulting in a failed global assembly. Similarly, fixing the cross interaction places a larger load on self-self-interactions to achieve correct local positioning, but these cannot stabilize relative component ordering by definition and as such this strategy results in reduced overall global stability for simple and complex targets alike.

To summarize, the analysis of Figs. 6 and 7 strongly suggests that global binary assembly can be understood as follows: self-interactions act as a “primer” that help ensure that individual component particles assemble into the right positional order (coordination shells), while cross interactions “bind” the locally ordered particles into their correct, target positions.⁷⁵ This is why, despite the self-interactions failing to assemble the target local structure in an analogous single-component system, they nonetheless encode the relevant length scales present in the local structure. In fact, substituting the optimized single-component square lattice interaction (which necessarily contains the characteristic length scales of a square) for the self-interactions in the intercalated square binary system results in successful assembly of the target, with greatly enhanced local AA(BB) stability as seen in the RDFs in Fig. S4 of the supplementary material. The analysis also highlights the importance of the interplay among the self- and cross interactions for intricate structures, as it is through mutual coupling that the entire global structure is locked into place. This result also implies that the self-interactions that establish a necessary local structure may be much simpler than those required for an equivalent single-component system, so long as full system interactions are in place. Indeed, while the stretched truncated square local structure in the binary square corral could be achieved by the BB component interaction only spanning four coordination shells, the equivalent single-component truncated square interaction required five shells and resulted in poorer assembly (see Fig. S3 of the supplementary material).

In this work, we extended a recently introduced inverse design methodology, relative entropy optimization, to discover new isotropic interactions that favor the formation of targeted multicomponent phases. We have used this approach here to determine interactions for binary mixtures that stabilize a wide variety of two-dimensional lattice assemblies and, in doing so, have gained new insights into how adopting a multicomponent system can affect the overall prospects for self-assembly. Although the expanded parameter space for design in multicomponent systems increases the complexity of the design problem, it helps discover simpler interparticle interactions (as compared to single-component systems) to assemble a desired target phase. It also allows for designed assembly of complex phases with structures that cannot be stabilized by single-component materials.

Mechanistically, our results suggest that optimized interactions between like components in a binary mixture act as a “primer” to help ensure that such species adopt, on average, the correct positional order for the target phase. Cross interactions, in turn, act as a “binder” to further ensure that species conform to the precise local compositional order required by the target. For complex or open lattices, independent design of cross and self-interactions is required to stabilize the desired assemblies.

Physically, reported single shoulder interactions for targets such as Fig. 2 (bottom) or Fig. 3 (top) are similar to those displayed by core-corona⁷⁶ or dendritic particle systems⁷⁷ and predict structure stability at (osmotic) pressures moderately above 1 atm, assuming room temperature and particle size scales on the order of ∼10 nm. However, similar experimental analogs for the more complex interactions required for the remaining structures in this work are more challenging to envision. In future work, it will be interesting to pursue related inverse design calculations for classes of materials whose self- and cross interactions are effectively constrained in ways that can be encoded by fundamental physics or empirical mixing rules.

See supplementary material for a derivation of the RE update scheme and a description of how interaction cutoffs are determined. Additionally, we include figures for the triangular honeycomb structure assembly from a single component, self-assembly results for the square truncated hexagonal structure with fixed hard-core-like interactions, comparison of single component assembly of the B component sub-lattice in the square corral target, and an example of how using optimized single component interactions may help boost binary assembly stability.

T.M.T. acknowledges support of the Welch Foundation (No. F-1696) and the National Science Foundation (Nos. DMR-1720595 and CBET-1403768). We also acknowledge the Texas Advanced Computing Center (TACC) at the University of Texas at Austin for providing computing resources used to obtain results presented in this paper.