Assembly of particle strings via isotropic potentials

Owing to the chemical and structural versatility of their building blocks, colloidal materials can be designed to assemble into a variety of microstructures. One motif that may enable functionality is chains or strings of colloids (so-called “colloidomers”¹). When percolated, these structures can serve as a template for conductivity, facilitating photonic and electronic transfer along the connected interparticle network.^2,3 Mechanically stable particle networks with a high surface-to-volume ratio are also interesting morphologies for nanoporous catalysts,^4,5 and open network structures may be advantageous for applications that require materials that can dynamically reconfigure in response to a stimulus.

Self-assembling colloidal strings have typically been designed by choosing building blocks with anisotropic interactions commensurate with the targeted morphology. For instance, thin chains of particles have been assembled in silico from hard spheres, each decorated with two collinear attractive patches.^6–9 Such directionally specific interactions have been engineered in practice by grafting appropriate functional groups (e.g., DNA) to the surface of colloids.^1,10,11 Analogous physics can be achieved via short-range anisotropic dipolar interactions, which have been shown to promote linear chain growth between charged gold nanoparticles.¹² Strong dipole interactions have been argued to be the driving force behind the self-organization of other nanoparticles into one-dimensional chains^13–15 and three-dimensional percolated fractal chain networks.¹⁶ 

Anisotropic colloid shape can similarly be tailored to obtain flexible colloidal chains or stringy structured fluids.^17,18 For instance, Sacanna and co-workers employed Fischer’s lock and key recognition mechanism between a homogeneous sphere and a sphere with a cavity to assemble compact clusters as well as more complex and flexible colloidal polymers. Particle and interaction anisotropy can also be induced by the assembly process itself. Under certain experimental conditions, spherical and uniformly grafted nanoparticles in a homopolymer matrix self-assembled into linear chains.¹⁹ The short-ranged depletion attraction in these systems is argued to be counterbalanced by the entropy of distortion which arises when the grafted brushes on two nanoparticles compress due to steric constraints upon approach. This can lead to an anisotropic distribution of the local graft density,^19,20 imparting an amphiphilic character to the nanoparticles. Such anisotropic assembly of uniformly grafted nanoparticles has also been predicted via theory as well as simulations where the ligands are modeled explicitly.^21–27 

A few attempts have been made to design an isotropic pair potential that causes a single-component fluid of particles to self-assemble into stringlike structures. To this end, Rechtsman et al.²⁸ proposed a complex “five-finger potential” containing five repeating attractive wells at intervals set by the particle diameter that are separated by repulsive barriers that inhibit formation of compact objects. In two dimensions, simpler potentials, though still possessing competing attractions and repulsions, have been shown to generate stringy structures. For example, coarse-graining multicomponent simulations of grafted nanoparticles revealed several single-component isotropic pairwise potentials that promote self-assembly of distinct morphologies: dispersed particles, long strings, and a percolated network.²⁶ Moreover, a class of potentials characterized by a single attractive well followed by a repulsive barrier furnished by a piecewise function of linear components has displayed different microstructures, ranging from monomers to aggregates to short strings to a labyrinthine chain network, as a function of area fraction and range of the repulsion.²⁹ 

Common to the above studies is the presence of competing interactions, such as a short-range attractive and long-range repulsive (abbreviated SALR) potential, also known to promote the self-assembly of more compact particle clusters. The repulsive interactions in such potentials naturally limit the aggregation that would otherwise be promoted by the attractions. Self-assembly of clusters and strings requires growth that is self-limited—for clusters with respect to their overall size and for strings with respect to the dimension normal to growth. A few studies further reinforce potential connections between compact and elongated cluster morphologies. For instance, it was shown in both experiment³⁰ and simulation^31,32 that both clusters and thick ramified structures are possible when systems possess competing SALR interactions. Similarly, in simulation, an SALR potential was shown to produce percolating states with a mixture of transient filamentous and spherical aggregates when the packing fraction exceeded 0.148.³³ 

Motivated by the investigations described above, the aim of this paper is twofold. The first is to use methods of inverse design^34,35 [specifically, a recent strategy³⁶ based on relative entropy (RE) coarse-graining^37,38] to discover isotropic potentials that promote self-assembly of one-monomer wide chains or compact clusters, in a one-component system of spherical particles. The second aim is to identify, on the basis of the designed interactions, a simpler model pair potential that favors assembly of these and related structures as a function of the length scales of the competing attractive and repulsive interactions.

The balance of this paper is structured as follows. Section II outlines the relative entropy based method we adopt for inverse design and presents details of the molecular simulations. The pair potentials and structures resulting from the inverse design for both compact clusters and strings are described in Sec. III. Motivated by the qualitative forms of the designed interactions, Sec. IV introduces a simpler related pair potential and explores the various morphologies that it favors as a function of its parameters. Conclusions and possible directions for future research are presented in Sec. V.

Relative entropy (RE) coarse-graining,^37,38 also known as likelihood maximization in probability and statistics, is used in this work to obtain isotropic potentials capable of self-assembling particles into different target structures. Commonly applied to obtain a reduced dimensionality description of complex molecules for simulation, RE course-graining has more recently been used to design isotropic pair interactions that lead to self-organization of a rich variety of equilibrium structures including fluidic clusters,^36,39 porous mesophases,^36,40 and crystalline lattices.^36,41,42

In brief, the RE course-graining protocol considers a target ensemble of particle configurations that collectively exhibits a desired structural motif (e.g., strings or compact clusters), discussed below. The optimized interactions are those that maximize the overlap of the probability distribution for configurations at equilibrium with that of the target ensemble. Here, we consider an isotropic pair potential, $Ur|θ$⁠, with m tunable parameters, θ = [θ₁, θ₂, …, θ_m]. θ is the corresponding amplitude associated with each knot in an Akima spline.⁴¹ According to RE coarse-graining, the parameters are updated in an iterative manner via

where $β=(kBT)−1$⁠, k_B is the Boltzmann constant, T is the temperature, α is the learning rate, g(r|θ^k) is the radial distribution function of the system in the kth iterative step of the optimization, and g_tgt(r) is the radial distribution function of the target ensemble. In practice, g(r|θ^k) is obtained from the equilibrium particle configurations of a molecular simulation using U(r|θ^k).³⁶ A rigorous mathematical derivation of the above update scheme is reviewed in Refs. 36, 41, and 43. The outcome of a successful optimization is a thermally nondimensionalized interaction βU_opt(r) that results in an equilibrium structure that closely mimics that of the target ensemble.

The first step in the inverse design protocol described above is the construction of an ensemble of target configurations from which g_tgt(r) can be computed. Target distributions can be specified in any way that yields an ensemble of desired configurations with convergent statistics. Here, we employ equilibrium statistical mechanics via canonical-ensemble molecular dynamics simulations with N spherical particles of diameter σ in a periodically replicated cubic cell of side length L [i.e., packing fraction η = Nπσ³/(6L³)] at temperature T. Dimensionless (generally many-body) interparticle potentials βV, given below, are selected for each target ensemble to yield configurations characteristic of the desired morphology.

Target ensembles of monodisperse, compact clusters of size N_tgt = 2 (dimer), 4 (tetramer), and 8 (octamer) are generated via molecular dynamics simulations in the canonical ensemble at packing fraction η = 0.025. The following interactions are chosen to mimic the desired target morphology. Particles of diameter σ interact with hard-sphere-like repulsions represented via a Weeks-Chandler-Andersen (WCA) potential V_WCA(r),

where βε_w = 5. Particles are assigned to a particular cluster, the compactness of which is enforced by applying an additional finitely extensible nonlinear elastic (FENE) spring potential V_FENE between each pair of particles in the cluster

with k = 30k_BT/σ² and r₀ = 1.5σ (N_tgt = 2) or 5.0σ (N_tgt = 4, 8). Additionally, a minimum distance of separation between the clusters is ensured by introducing an isotropic Yukawa repulsion between particles in different clusters

where βε_y = 30, 10, and 0 for N_tgt = 2, 4, and 8, respectively, κ = 0.5σ⁻¹, and r_cut = L/2.5. After ensuring equilibration, the target radial distribution function g_tgt(r) for N_tgt = 2, 4, 8 sized cluster fluids is computed.

For strings, we created the target ensemble from simulations of linear particle chains of size N_tgt = 10. Monomers interact via the repulsive WCA potential of Eq. (2) with βε_w = 1, and adjacent beads interact via the FENE spring potential V_FENE of Eq. (3) with k = 30k_BT/σ² and r₀ = 1.5σ. Nonbonded monomers also interact with the Yukawa potential V_Yukawa of Eq. (4) with βε_y = 30, κ = 0.5σ⁻¹, and r_cut = L/3.0. The system of linear chains is allowed to evolve via molecular dynamics and, upon equilibration, g_tgt(r) is calculated at different packing fractions η = 0.05, 0.1, and 0.15.

Molecular dynamics simulations for the target structures and the optimization are performed in the canonical ensemble with a periodically replicated cubic simulation cell.

The software package HOOMD-blue 2.3.4^44,45 is used to generate the target configurations of particles of mass m. A time step of $dt=0.005σ2m/kBT$ is adopted, and the Nosè-Hoover thermostat with a time constant of τ = 50dt is employed. For target compact clusters of N_tgt = 2, 4, and 8 at a packing fraction of η = 0.025, N = 384 particles are used in a periodic cell of size L = 20σ. Target structures for strings are generated using N = 320, 420, and 630 particles in cubic cells of dimension L = 15σ, 13σ, and 13σ for different packing fractions of η = 0.05, 0.10, and 0.15.

Simulations within the iterative RE optimization as well as for the forward runs with the optimized potentials were performed in GROMACS 5.1.2 with a time step of $dt=0.001σ2m/kBT$⁠. Constant temperature was maintained using the velocity-rescale thermostat with a time constant τ = 50dt. The phase diagram in Sec. IV was generated using N = 968 particles in a box size L = 15σ (i.e., η = 0.15), collecting configurations for 3 × 10⁷ time steps. In order to extract improved cluster statistics as well as to check for finite-size effects, simulations with the larger box size of L = 30σ with 7744 particles were also performed for select state points. For these same state points, simulations were also conducted with different initial configurations (randomly placed particles, compact aggregates, and linear rod-like clusters) to test the robustness and sensitivity of the final structures as well as to ensure equilibrium is attained. The simulation snapshots were rendered using the Visual Molecular Dynamics⁴⁶ software.

The structures obtained from the simulations described above are characterized by the cluster size distribution (CSD), the distribution of the number of nearest neighbors, the effective dimension of the resulting aggregates d_f, and a percolation analysis. The CSD quantifies the fraction of clusters that contain n particles, where a particle is considered to belong to a cluster if its center is within a prescribed cut-off distance r_cut from the center of at least one other particle in the same cluster. The smallest range of attraction for the optimized potentials in this work is ∼1.1σ, which makes it a natural choice for r_cut. For both cases of spherical clusters and strings, various cutoffs [1.05σ − 1.25σ] yield nonperceptible changes in the CSD and other structural properties.

To distinguish between thick ramified structures and “thin” chains of one-monomer diameter width, the number of nearest neighbors of each particle is evaluated. Single strands of strings have predominantly two bonds per particle. Additionally, to characterize the anisotropy of the resulting objects, the effective dimension d_f is determined via $Rg∼n1/df$⁠, where n is the number of particles in a cluster and R_g is the radius of gyration. The latter for a cluster of size n is defined as

where r_i and R_CM are the coordinates of the ith particle and the center of mass of the cluster of n particles, respectively. The average is performed over all clusters of size n. The effective dimension is expected to be bounded by two limits: ∼3 for compact, homogeneously spherical aggregates to ∼1 for linear objects. For stringlike objects, we sometimes find that a single value for d_f does not satisfactorily fit the data for every aggregate size n, in which case we segregate the data that visually appear to have different slopes in the R_g vs log(n) plot and compute distinct values for d_f for the different regions of n.

Some of the structures discussed below are found to be percolating. Percolation is defined if a cluster spans the length of the box in at least one direction such that, under periodic boundary conditions, the cluster wraps around the box and connects to itself. Accounting for periodicity, a percolating cluster is thus infinitely long and hence the corresponding d_f or R_g is not computed for percolated structures. If at least 50% of the configurations conform to the said definition, the resulting morphology is deemed to be percolating.⁸ 

Given the connection between clusters and stringlike structures described in the Introduction, it is instructive to use inverse design to discover pair potentials that favor these morphologies. As described in Sec. II, we first use RE optimization here to determine isotropic pair potentials that lead to self-assembly of compact N_tgt-mer clusters, namely, dimers, tetramers, and octamers (target sizes of N_tgt = 2, 4, 8) at a packing fraction of η = 0.025, where the choice of η is motivated by prior work.³⁹ In each case, the RE optimization discovers an interaction βU_opt(r) capable of successfully assembling the target morphology, despite not quantitatively reproducing g_tgt(r). A detailed comparison is provided in the supplementary material (Fig. S1). The region of depletion in g(r) due to the repulsion between separate clusters is less pronounced in the equilibrium assembled structures relative to the target ensemble for N_tgt = 2 and 4. This is a consequence of the presence of particle exchange in the pairwise self-assembling simulations but not the target simulations. The monodispersity of the target clusters is ensured by enforcing real physical bonds between the constituent particles in the clusters and strings. This mechanism disallows particle exchange between the clusters. The optimized ensembles are a consequence of the equilibrium self-assembly of particles with no such strict enforcement, and thus, the exchange of particles between the optimized clusters yields a weakened depletion zone. Compared to the target ensemble, the optimized structures are also polydisperse in nature as is quantified below and reflected in the associated radial distribution function. For instance, for the 4-mer clusters, the optimized g(r) has a notable peak around r/σ = 1.73 that is absent in the target. This peak is due to 5-mers that have locally crystallized into a trigonal bipyramid; see Fig. S1(b). For the 8-mers, the smearing out of the intracluster peaks in the optimized g(r) can also, similarly, be attributed to polydispersity. Thus, it is a combination of equilibrium constituent particle exchange and cluster size polydispersity that consequently results in the mismatch between the radial distribution functions of the structures self-assembled from the optimized pair interactions and those of the target structures.

Figure 1(a) shows the optimized pair potentials βU_opt(r) that form such compact clusters. The two main features of the optimized potentials are the attractive well beginning at r = σ followed by a repulsive barrier. In accord with prior work,³⁹ the magnitude of the attraction increases and the peak of the repulsive barrier shifts to larger separations with increasing aggregate size. Using the designed pairwise interaction, strong clustering emerges. The corresponding CSDs in Fig. 1(b) demonstrate that the target cluster sizes are reproduced for all three cases of N_tgt = 2, 4, and 8. There is mild polydispersity with respect to the aggregation number, though most clusters are within one particle of the targeted size. The self-assembled aggregates are well-separated and behave as an equilibrium cluster fluid. See a snapshot for N_tgt = 8 in Fig. 1(c); corresponding movies (M1.mp4 and M2.mp4) are available in the supplementary material. The compactness of the largest clusters (N_tgt = 8) is quantified by the corresponding effective dimension d_f, which is given by the inverse of the slope of Fig. 1(d) and is approximately 2.9.

The N_tgt = 2 case of self-assembling dimers shares some complexities with the problem of string formation in that growth must be limited to a single direction. A pair of particles must associate attractively at r = σ, but formation of triangles, where an incoming particle bonds between the dimer pair at the point of closest approach to both centers (i.e., r = σ), must be suppressed. Indeed, a dimer can be considered as both the smallest cluster and the smallest string. Looking forward to self-assembly of strings, we might anticipate that the characteristics of the potential optimized for forming dimers (a net-repulsive potential with a relatively narrow attractive well) are favorable for string formation more generally.

As described in Sec. II, the target ensemble of configurations for particle strings considered here comprises chains of particles of size N_tgt = 10. In general, the isotropic pair potentials resulting from the corresponding RE optimizations βU_opt(r) successfully self-assemble fluids of stringy particles. The radial distribution functions g(r) of particles interacting via the optimized potentials capture the salient features of the target structure g_tgt(r), though—as shown in the supplementary material (Fig. S2)—the depleted region present between r = σ and r = 2σ is less prominent in the former compared to the latter. This is consistent with that of the clusters where the depletion was muted due to the equilibrium exchange of constituent particles. Additionally, the strings formed by the optimized pair potential also branch out occasionally leading to further weakening of the depletion region.

Figure 2 shows βU_opt(r) at three different packing fractions η = 0.05, 0.1, and 0.15. Analogous to the potential optimized for forming dimers in Sec. III A, the attractive well is relatively narrow for all η. The repulsive barrier, on the other hand, is more sensitive to the packing fraction, increasing in range and magnitude as η is reduced. For η = 0.05, the pair potential is more complex due to the emergence of secondary features on the scale of a monomer diameter σ. The repulsive hump for η = 0.05 is followed by three secondary attractive minima at r = 2σ, 3σ, and 4σ, respectively. This potential with alternating attractions and repulsions is qualitatively reminiscent of the “five-finger potential” proposed in Ref. 28 to form colloidal strings. These features are progressively muted as the packing fraction of the target structure and the optimization is increased. For η = 0.1, there is a slight hint of a dip at r = 2σ which is further reduced for η = 0.15 where the repulsion terminates at r = 2.25σ.

Using the optimized pair potentials, linear stringlike structures are observed to form at all three packing fractions. Figures 2(b) and 2(c) show representative simulation snapshots in equilibrium for η = 0.05 and 0.15, where stringlike objects are visually apparent (M3.mp4 and M4.mp4 in the supplementary material). More quantitatively, R_g as a function of string size at three different packing fractions are reported in Figs. 3(a) and 3(b). The corresponding effective dimension for the chainlike structures at η = 0.05 and 0.10 is approximately 1.1. The fractal analysis at η = 0.15 shows two power laws, d_f ∼ 1.20 for n ≤ 8 and 1.80 for n > 8, implying that the shorter chains are more linear, while the longer strings are more curvilinear. Thus, the resulting optimized aggregates span from rod-like to chain-like, with the linearity of the chains increasing as η decreases. Figures 2(b) and 2(c) highlight two illustrative examples where it is apparent that the selected chain at η = 0.05 is more linear than that at η = 0.15.

The chains, at all three volume fractions, predominantly have two nearest neighbours P(N_nn) ≈ 2; see Fig. 3(c). It may be somewhat surprising that there is such a large percentage of particles (≈20%) with three neighbors, particularly for η = 0.05, where d_f indicates that the objects are nearly linear. By considering the size n and N_nn for every aggregate, we discern that a fraction of aggregates are actually compact clusters. For example, compact tetramers are characterized by N_nn = 3. Indeed, we find that when the compact clusters are removed from the calculation of P(N_nn), the value of P(N_nn = 3) is significantly reduced; see Fig. S3 in the supplementary material. The percentage of aggregates that are compact clusters at η = 0.05, 0.10, and 0.15 are 28.1%, 20.1%, and 6.7%, respectively.

Finally, unlike the previous case of compact clusters where distinct peaks are formed at the desired target cluster size from the optimized interactions, here we note no such size-specific assembly for the strings. The CSDs for the optimized potentials are shown in Fig. 3(d), where polydispersity is obviously high and increases with packing fraction. Thus, while our optimization procedure illustrates that isotropic pair potentials can readily assemble monomer-wide stringlike particle structures, they are limited in their ability to control their length.

Self-assembly is the result of an interplay of energetic and entropic contributions that determine which types of structures minimize the free energy of a given system. Despite this inherent complexity, we can gain insight into the propensity for a given potential to form either stringlike or compact objects by considering the energy for a test particle to approach a dimer, where “end-attachment” to the dimer gives rise to a short string and “middle-attachment” yields a compact triangle. A slice of the potential energy landscape seen by the test particle relative to an ideal dimer is shown in Fig. 4 as a heat map for two illustrative cases: (a) the potential optimized for compact tetramers at η = 0.025 and (b) the potential optimized for strings at η = 0.15. The corresponding heat maps for the remaining cluster and string cases are demonstrated in Figs. S4 and S5 in the supplementary material. For the compact cluster-forming potential, the energy is lowest when the test particle bonds to both particles of the dimer simultaneously due to the deep attractive well at r = σ, promoting middle-attachment. Because the potential has been optimized to form compact tetramers, the repulsive corona surrounding the pair of particles in Fig. 4(a) penalizes the middle- or end-attachment of additional particles necessary to form larger compact clusters or strings.

By contrast, when the potential optimized for string formation is used, end-attachment to the dimer is more energetically favorable than middle-attachment, fostering chain growth. End-attachment is favored for this particular interaction because the energy of the attractive well at r = σ is greater than at r = 2σ, in part due to the relatively narrow repulsive barrier. Unlike for clusters, the chain-forming potential is net repulsive so that the lowest energy position for the test particle is to not attach to the dimer at all. However, in a bulk system, sufficiently high pressures induce particle association. Note that the preceding preference for end-attachment over middle-attachment will continue to be present as particles are added to the chain; that is, there is no energetic mechanism for controlling the length of the chain inherent to the potential itself.

While the above analysis is limited in that potentially important effects (e.g., the impact of surrounding aggregates on the energetics and the role of entropy) are omitted, this simplified model lends insights into how the length scales of the attractive well and the repulsive barrier might influence self-assembly. In particular, the position of the repulsive barrier controls the size of the compact clusters, in keeping with prior work.³⁹ Furthermore, the relatively narrow repulsive barrier in the string-forming potential promotes end attachment.

Comparing the pair potentials optimized for compact clusters to those optimized for chains, the common features are an attractive well at short separations and an outer repulsive barrier with a shorter range than that characteristic of the Yukawa potential routinely used in SALR potentials. As noted in Sec. III C, the location of the repulsive barrier (which is controlled by the widths of the attractive well and the repulsive barrier) is one of the key parameters in promoting either end or middle attachment.

Motivated by these basic features shown in Figs. 1 and 2, we propose a simple and tunable pair potential form that favors various equilibrium structures, from disordered monomeric fluid to thin particle strings to compact particle clusters. By “simple,” we imply a pair potential that (1) has a minimal number of features and parameters that are essential to capture different morphologies computationally and (2) can be expressed in a closed-form. Specifically, we consider a sum of a steep WCA repulsion at a short distance followed by two half harmonic potentials that mimic the attractive well and the repulsive barrier,

The three sub-potentials are defined as

where βε_wca = 1.5, α = 12, w₁ = σ + δ_A, and w₂ = σ + δ_A + δ_R. The remaining adjustable parameters (βε_R, δ_A, δ_R) are determined by fitting the above form to the optimized potential that promotes self-assembly of strings at η = 0.05. This particular potential is used as the reference because strings are generally more challenging to self-assemble than compact clusters from an isotropic potential, and η = 0.05 is the packing fraction at which the designed potential resulted in strings with the smallest effective dimension. To simplify the reference, the optimized potential is truncated beyond r = 2σ and vertically shifted so that U(r = 2σ) = 0. As shown in Fig. 5, the resulting fit approximates the short-ranged (r ≤ 2σ) reference potential well with βε_R = 6.76, δ_A = 0.5, and δ_R = 0.5. To avoid discontinuities in the force profile, the above potential is weakly smoothed using a successive two-point averaging scheme where, beyond r = σ, U(r_i) is twice replaced by an average of U(r_i−1) and U(r_i+1).

In this section, we explore the effects of tuning the ranges of attraction (δ_A) and repulsion (δ_R) in the model pair potential introduced in Sec. IV, while holding βε_R = 6.76 constant at a packing fraction of η = 0.15 (i.e., the lowest value of η for which we observed the percolated string morphologies described below). Based on the discussion in Sec. III C, we expect that modifying δ_A and δ_R will bias the potential toward favoring assembly of either strings or compact clusters. Figure 6(a) shows four possible potentials where δ_A = 0.2 and δ_R = 0.5, 0.7, 1.0, 2.0. For the family of potentials where only δ_R is varied while δ_A = 0.2 is constant, we observe four broad classes of structures in molecular dynamics simulations, shown in Figs. 6(b)–6(e), and with corresponding CSDs and nearest neighbour distributions, shown in Figs. 7(a) and 7(b), respectively.

1. Monomers (M): A fluid of particles that remain in a dispersed state results when δ_A = 0.2 and δ_R = 0.3. Quantitatively, a state is defined to be monomeric if the fraction of monomers (cluster size n = 1) in the CSD exceeds 50%. The CSD for this state point [Fig. 7(a)] shows that 22% of the aggregates are dimers, while 65% are monomers. Correspondingly, the nearest neighbour histogram [Fig. 7(b)] shows that a significant majority of the particles have zero or one nearest neighbor.

2. Strings (S and SP): Single-stranded stringy particle assemblies are obtained at values for {δ_A, δ_R} of {0.2, 0.5} and {0.2, 0.7}. We distinguish between shorter stringlike objects (S) and percolated networks of strings (SP) on the basis of the percolation analysis described in Sec. II D. Snapshots of short strings (δ_R = 0.5) and percolated chains (δ_R = 0.7) are depicted in Figs. 6(b) and 6(c), respectively; see the supplementary material for a movie of the simulation (M5.mp4). Stringlike structures, regardless of whether percolated, have predominantly two bonds per particle. The primary difference upon transitioning from shorter strings to a percolated porous network is the growth in the number of junctions, imparting greater connectivity to the branches. Thus, the corresponding nearest neighbour distribution for the percolated network shows a higher value for N_nn = 3 as compared to that of strings; see Fig. 7(b). For short strings, we confirm that the aggregates are elongated on the basis of visual inspection and the d_f. As observed for the pair potential optimized for strings at η = 0.15 with {δ_A, δ_R} values of {0.2, 0.5}, we observed two distinct regimes of d_f as a function of aggregate size n. For n ≤ 11, d_f = 1.35, and for larger aggregates, d_f = 1.9. More generally, we identify structures as stringlike on the basis of both a peak at P(N_nn) = 2 and a d_f in the range of 1–1.5 for small (n ≲ 10) objects and a d_f around 1.9–2.0 for the longer chains.

Compared to the results presented in Sec. III B, there are markedly fewer compact aggregates in coexistence with the strings: 1% of aggregates for the short strings and 0.5% of aggregates for the percolated strings are compact. However, in keeping with the inversely designed potentials, there is no apparent size-specificity for the strings. The CSDs in Fig. 7(a) show a monotonic trend of decreasing probability of larger aggregates for the unpercolated structures. To the best of our knowledge, this is the first demonstration of the self-assembly and stabilization of spherical particles into a three dimensional open, porous network of single-stranded chains via an isotropic potential at nonzero temperature.

1. Clusters (C): Well-defined compact clusters are observed [Fig. 6(d)] using the proposed potential at {δ_A, δ_R} = {0.2, 1.0}; see the supplementary material for a movie of the simulation (M6.mp4). Clusters are identified by a sharp maximum in the CSD at a cluster size n > 1 and by an effective dimension of ≲ 3 owing to their compact and isotropic shape. At this state point, the clusters are tetramers; see the prominent peaks at P(n) = 4 and P(N_nn) = 3 in Figs. 7(a) and 7(b). The corresponding d_f is 2.5. Instead of a fluid of clusters as studied in Sec. III A, the self-assembled clusters crystallize under these conditions onto a lattice.

2. Cylindrical Spirals (CS): Elongated, cylindrical spirals of colloids are observed at {δ_A, δ_R} = {0.2, 2.0}. These are multistranded, percolated networks of anisotropic aggregates. A special case of these structures with three helical chains is commonly referred to as Bernal spirals^30–32,47 which have six nearest neighbours, and a representative snapshot is shown in Fig. 6(e). A movie of the simulation is available in the supplementary material (M7.mp4). Accordingly, Bernal spirals are identified as percolated structures with a peak in P(N_nn) = 6.

In Fig. 8, we further explore how the morphologies identified above emerge in this model as a function of δ_R and δ_A. For the conditions studied, note that the quantity δ_A + δ_R appears to be the primary determinant of which self-assembled structures are observed. When δ_A + δ_R ≲ 0.5, the particles form a fluid of well-dispersed monomers (M). For progressively larger δ_A + δ_R, short single-stranded strings (S) of particles form followed by interconnected percolating networks of strings (SP). For sufficiently low attractive ranges (δ_A ≤ 0.2), the physical bonds are labile and percolated strings are fluidic, continually breaking and reforming flexible uniaxial structures. However, for δ_A > δ_R, the increased attraction impedes dissociation of the particles especially at the junctions. As a result, though the chains still fluctuate in local order, larger length-scale motions are suppressed.

When δ_R + δ_A ≳ 1, the space-spanning stringy particle network morphs into a crystalline arrangement of compact particle clusters (C). Consistent with Fig. 1(a), the clusters grow with increasing δ_A. For δ_A ≲ 0.6, the clusters are tetramers, but for larger values of δ_A, aggregation numbers of 5–7 are observed. The corresponding d_f values range from 2.5 to 2.6. Increasing δ_R + δ_A further compels the spherical clusters to coalesce, eventually resulting in kinetically arrested percolated networks of cylindrical structures (CS), of which the Bernal spiral is a special case indicated in Fig. 8. We note that, at these particular state points, the total energy per particle displays a very slow drift with time on a logarithmic scale which is a typical signature of nonequilibrium structures that are dynamically arrested.

As detailed above, the interplay between the length scales of the attraction and repulsion in this model pair potential results in a rich variety of self-assembled structures. The analysis in Sec. III C suggested that comparing the energy for a test particle to bond to either end of a dimer vs the mid-point is a helpful, though simplistic, predictor of whether the clusters formed by a given potential will be stringy or compact, respectively. Here, we also find that this analysis helps understand why δ_R + δ_A is an important parameter in determining the observed morphologies. Specifically, we compare the energetics of end-attachment (U_end ≡ U(r = σ) + U(r = 2σ)) and middle-attachment (U_mid ≡ U(r = σ) + U(r = σ)) of a test particle to an isolated dimer. Regions of phase space where end-attachment is preferred (U_end ≤ U_mid) are shaded in Fig. 8. Note that when δ_A + δ_R ≲ 1.0, end-attachment to a dimer is preferred; when δ_A + δ_R is larger, middle-attachment is favored. Interestingly, δ_A + δ_R ≈ 1.0 also approximately corresponds to the crossover between percolated string networks and compact clusters observed in the simulations. The corresponding heat maps for four morphologies are individually demonstrated in Fig. S6 in the supplementary material. Thus, though based on a simplistic energetic analysis, we can gain insights into how the length scales of the short-range attractions and longer-range repulsions in a pair potential can favor the formation of stringlike vs compact self-assembled structures.

We used an inverse design strategy based on RE optimization to determine and study isotropic potentials capable of driving one-component systems of particles to self-assemble into compact versus linear stringlike clusters. Simulations using particles interacting via the optimized potentials demonstrated spontaneous formation of the targeted morphologies, though successful design of specific aggregation numbers was only achievable for compact clusters, a limitation that might be expected for isotropic potentials. The simplicity of the optimized potentials for these structures is remarkable given that prior computational efforts to arrive at monomer-wide stringlike clusters in three dimensions have employed directional bonding or anisotropy of the colloidal building blocks to control assembly.

Motivated by the RE optimized potentials, a universal potential with a simple functional form is proposed that is capable of assembling a rich variety of complex architectures: monomeric fluid, fluid of short chain-like structures, percolated networks of strings, crystalline assemblies of compact clusters, and percolated thick cylindrical structures including Bernal spirals. The proposed model potential is a combination of short-range attraction at contact, which can be realized by polymer-mediated depletion in chemistry-matched systems (for smaller values of δ_A) and a medium-ranged repulsive barrier which approximately mimics that of suspensions of noncharged brush-grafted nanoparticles.⁴⁸ Polymer depletants that are responsive to external stimuli (e.g., pH,⁴⁹ temperature,⁵⁰ light,⁵¹ and other fields⁵²) represent another interesting avenue to tune δ_A to switch between different morphologies. We recognize that it might be challenging to design colloidal systems with interactions similar to the isotropic potential presented here to realize self-assembled colloidal strings in the laboratory. An alternative strategy would be to extend the recently introduced multicomponent RE optimization strategy⁵³ to design multisite (e.g., patchy) particles that can assemble into such morphologies, an approach that is complementary to that recently introduced to design patchy particles for assembly into a variety of two-dimensional crystal structures.⁵⁴ More generally, our results provide qualitative insights into the rich morphological phase diagrams that can potentially be realized in colloidal systems with (approximately) isotropic interactions with competitive repulsive and attractive components.

See supplementary material for additional details regarding the simulation movies. The supplementary material also contains the radial distribution functions for the clusters and strings obtained via direct design as well as optimization and the corresponding heat maps.

The authors thank Sanket Kadulkar and Michael P. Howard for valuable discussions and feedback. This research was primarily supported by the National Science Foundation through the Center for Dynamics and Control of Materials: an NSF MRSEC under Cooperative Agreement No. DMR-1720595 as well as the Welch Foundation (Grant No. F-1696). We acknowledge the Texas Advanced Computing Center (TACC) at The University of Texas at Austin for providing HPC resources.