Variational design principles for nonequilibrium colloidal assembly

The self-assembly of soft and biological matter out of equilibrium can result in novel structures and dynamical responses not constrained by thermodynamic considerations.^1–6 The microscopic violation of detailed balance in such systems can be used to design a wide range of functional materials with enhanced thermomechanical, optoelectronic, or drug-delivery properties.^7–9 The predictive inverse design to drive the assembly of target dissipative structures requires a dynamical description of the system.^10–14 We have developed a variational algorithm to automate the discovery of inverse design principles for colloidal self-assembly in a nonequilibrium steady-state in molecular dynamics simulations. The algorithm uses a variational principle arising from rare dynamical fluctuations of the system in a trajectory ensemble and optimizes the yield of target clusters with statistically estimated explicit gradients in the design parameter space. We demonstrate the performance of this algorithm by obtaining optimal design principles for the self-assembly of DNA-labeled colloids¹⁵ driven out of equilibrium by a shear flow. We expect that the ability to uncover general optimal inverse design principles away from equilibrium will enable bottom-up synthesis of new materials and elucidate the processes encoding structure in biological contexts.^16–18 

The self-assembly of nanoscale building blocks is increasingly used to engineer functional materials with novel properties arising from their complex nanostructures. Colloidal systems offer a versatile paradigm for inverse design toward a desired target structure due to the independent tunability of the shape, valency, and assembly environment.¹⁹ In thermodynamic equilibrium, stabilizing a target structure amounts to lowering its free energy, which is typically achieved by increasing both the strength and specificity of interactions. This principle has been exploited to achieve the self-assembly of a variety of clusters and superlattices from colloids and nanocrystals with crystal facets decorated with organic ligands or DNA.^20–22 However, employing these principles, in practice, requires mitigating dynamical effects such as slow coarsening and kinetic trapping.^11,12 Optimal forces for self-assembly must achieve a tradeoff between slow relaxation at high interaction strengths and slow growth at high interaction specificity.^23,24 The self-assembly of colloids and biomolecules in nonequilibrium steady-states provides a route to decouple kinetics from stability and mitigate this tradeoff. Directive self-assembly has been achieved by driving the system with a constant supply of chemical fuel or by applying external fields.^25–28 However, the design of such systems must confront the continuous supply of energy necessary to prevent the system from relaxing to equilibrium. Existing computational methods to discover inverse design principles for nonequilibrium self-assembly are limited due to the configurational probability not following the Boltzmann distribution and the corresponding variational structure afforded by the free energy no longer being valid under such dissipative conditions.

Recent advances in the theoretical treatment of the stochastic thermodynamics of nonequilibrium steady-states have made possible a trajectory ensemble description of self-assembly, treating structure, and dynamics on an equal statistical footing.^29,30 This has enabled basic principles governing assembly away from equilibrium to be formulated.^31,32 In this work, we develop a perspective and accompanying numerical technique based on these insights. Rather than considering the probability of observing a state and tuning its associated free energy, we consider the likelihood that a trajectory forms a specific structure as quantified by a stochastic action, and how that action is changed by modifying the intermolecular and applied forces. We show that fluctuations around a nonequilibrium steady-state encode the susceptibility of a system to assemble in a manner analogous to a fluctuation–dissipation relationship. Furthermore, optimal forces that assemble a target structure while minimizing the change to the stochastic action satisfy a variational principle.^33,34 We extend and apply an optimization algorithm that solves this variational expression and computes the optimal control force to sample rare dynamical phases.³⁵ We show that this algorithm can be used to solve the inverse design problem, deciphering how rare fluctuations encode stability away from equilibrium.

We outline below an inverse design algorithm for the self-assembly of sheared DNA-coated colloids into different target nanoclusters. The algorithm is based on a variational principle relating rare fluctuations in an ensemble of trajectories conditioned on evolving a target structure to effective forces achieving the target as the typical dynamical state. Working with a trajectory ensemble, where the probability distribution is known, circumvents the difficulty of not knowing the distribution of configurations within a nonequilibrium steady-state. To solve the variational problem, we have used generalized response relations for the gradients of the steady-state trajectory probability to a change in the inter-particle and externally driven forces.

For concreteness, we consider a model of N colloidal particles in a cubic box of length L, evolving with an overdamped Langevin equation of the form

where $ṙi$ are time derivatives of the coordinates of the ith particle and u_i are the forces acting on it. The friction coefficient of the colloids with the thermal bath is denoted γ, and η_i are Gaussian white noise that satisfy

where I₃ is the 3 × 3 identity matrix and k_BT is Boltzmann’s constant times the temperature. The angular brackets denote an averaging operation over the random noise distribution. As we consider dynamics in the presence of a shear flow, we use Lees–Edwards boundary conditions.³⁶ 

We use an ansatz of DNA-labeled spherical isotropic colloids as programmable building blocks for self-assembly. The interaction between these colloids, mediated by the DNA molecules attached to their surface, consists of a volume-exclusion repulsion and a short-range attraction.³⁷ The effective interaction strengths and the pairwise specificity can be independently tuned by varying the sequences of the grafted DNA molecules. During self-assembly, the short-range forces generate a competition between local and global order that leads to frustration and unique phase behavior and dynamical effects.^38,39 This system has been computationally and experimentally demonstrated to form finite nanoclusters with specific target structures.^40,41 The high-dimensional design space has the possibility to offer multiple pathways to stabilize any cluster out of the many nearly degenerate states formed without the specificity of the DNA-mediated attraction. To illustrate the performance of the variational algorithm, we consider the nonequilibrium self-assembly of 21 such rigid and nonrigid clusters, some examples of which are demonstrated in Fig. 1(a).

We examine the self-assembly of these colloidal particles under a constant linear shear flow. Shear flows are known phenomenologically to alter the stability of compact and extended colloidal structures.^42,43 A recent paradigm of colloidal assembly being increasingly explored is that in a microfluidic device, where the confining walls generate a strong shear on the assembling clusters.^44,45 This system offers a canonical nonequilibrium setting to explore inverse design principles. Taken together, the forces acting on the ith colloid are

where V(r) = V_WCA(r) + V_Morse(r) and V_WCA(r_ij) is a WCA pair potential representing the volume exclusion interactions and V_Morse(r_ij) denotes the DNA-mediated short-range pairwise attraction. The force due to a shear flow, $fiS(ri)$⁠, has the form

which has magnitude f and generates a constant gradient of the x component of the velocity along the z direction. The WCA pair potential has the functional form

with the particle diameter σ and energy scale ϵ. The attractive Morse potential has the functional form

where D_ij is the magnitude of the bond energy and α determines its width.

We work in units of k_BT = 1, γ = 1, and σ = 1. The natural time scale with these units is t₀ = γσ²/k_BT, and times are expressed in these units throughout. We set ϵ = 10k_BT and α⁻¹ = σ/10. The attractive energy scale D_ij and the shear flow rate f are tuned as variational parameters to induce self-assembly. They have been restricted to vary within the range 0 ≤ D_ij ≤ 10k_BT and 0 ≤ f ≤ 50k_BT/σ² to avoid large relaxation times and to stay within the overdamped regime. Figure 1(b) shows the potentials for the inter-particle interactions. The Morse potential and its force have both been truncated and shifted, using the Shifted Forces approximation,⁴⁶ to decay smoothly to zero at r_c = 2.12σ.

In order to avoid finite size effects in the formation of small clusters, we study a low packing fraction of ϕ = 0.01. We use a first order Euler discretization for the equation of motion in Eq. (1). Since the potentials in Eq. (3) are narrow and short-range, we have to use a small time step of 5 × 10⁻⁵t₀ in order to sample the potentials accurately. We have used trajectories of length ranging from τ/t₀ = 2.5 × 10³ to 10⁴.

In order to uncover design principles for self-assembly, we consider the task of finding the set of forces that fulfill the condition of assembling a target structure as the typical state of the system in the long time limit. Such tasks in stochastic dynamics are generalizations of Brownian bridges and known to have unique solutions.³⁴ They have played an important role recently in the application of large deviation theory to physical systems driven far from equilibrium.^35,47–50

We start by defining an observable A_τ as a time averaged indicator function for a target cluster,

where 1[r^N(t)] = 1 for a configuration satisfying a geometric criterion consistent with a target cluster and 1[r^N(t)] = 0 otherwise for each time t along a trajectory r^N(t) of total duration τ. The average value of the observable quantifies the yield of the target cluster. For all colloidal clusters considered, 1 is computed by constructing a bond-connectivity matrix. A cutoff of r_b = 1.35σ has been used to define a bond between two particles. Indicator functions for rigid target clusters are then uniquely determined by permutation-invariant measures of this connectivity matrix.⁵¹ For nonrigid target clusters, along with the bond-connectivity matrix, we additionally consider measures of the geometry of the cluster for defining the indicator function.

Rather than considering trajectories conditioned on a particular value of A_τ directly, which is numerically cumbersome, we work within an ensemble equivalent representation.⁵² Using a counting parameter λ, we can statistically bias a system toward a particular value of A_τ within a nonequilibrium steady-state. The cumulant generating function ψ(λ) is the partition function associated with the trajectory ensemble under the statistical action of λ,

where the angular brackets denote a path average over trajectory probability P₀[r^N(t)] as

The subscript 0 refers to the average being computed in a reference ensemble where the particles do not typically show the desired self-assembly behavior. For this reference system, we have chosen an equilibrium ensemble of colloids interacting only with the WCA repulsive forces, i.e., D_ij = f = 0, such that u_i = −$∇i$∑_j≠iV_WCA(r_ij), which is denoted as $FiWCA(rN)$⁠.

When the optimizable parameters are tuned to vary u, the trajectory probability P₀[r^N(t)] changes to P_u[r^N(t)]. The cumulant generating function can be estimated in the modified ensemble as

where the functional form of the relative action ΔS[u] can be derived from Onsager–Machlup theory,⁵³ 

with the integral being computed in the Ito sense. This change of measure analogous to a Girsanov transform³³ relates the original likelihood of self-assembly in the reference ensemble to the ensemble under the control force.

Since the exponential is a convex function, we apply Jensen’s inequality to Eq. (10),

to obtain a variational expression for the cumulant generating function. In the long time limit for τ → ∞, we can replace trajectory averages with static averages and simplify the relative action using the equation of motion for $ṙi$⁠. Hence, we arrive at our final variational expression

where ⟨Ω[u]⟩_u is the target function to optimize.

For a bounded observable like the indicator function, a large value of λ enforces the desired conditioning. The problem of saturating the variational inequality is known to have a unique solution when u(r^N) can take all possible functional forms, the optimal force being a generalization of Doob’s h-transform.³⁴ Optimizing ⟨Ω[u]⟩_u will lead to a set of many-body forces for assembling target clusters in high yield. In practice, the use of only one-body and two-body forces in Eqs. (3)–(6) need not saturate the inequality. The second term in the variational expression is associated with the Kullback–Leibler divergence between the reference and conditioned trajectory ensembles. This term enforces the smallest excess force out of all possible control forces and thus acts as a regularizer in the optimization process. While the solution to Eq. (13) uniquely selects the force that shows fluctuations closest to rare fluctuations in the original ensemble, it is not a unique inverse design criterion and alternatives can, in principle, be constructed.⁵⁴ However, the optimization scheme that we construct in Sec. II C can be generally extended to other functional forms of regularizers.

To numerically optimize Eq. (13), we derive explicit gradients of the variational estimator ⟨Ω[u]⟩_u using an algorithm that we have previously employed to estimate large deviation functions in nonequilibrium steady-states.³⁵ The general form of the gradient with respect to any variational parameter c ∈ {D_ij, f} is

where $Ṡ[u]$ is the time derivative of the action S[u]. Equation (14) is a generalized fluctuation–dissipation relation for a nonequilibrium response in the design parameter space. For computational purposes, we approximate the gradient expression by integrating the correlation function in the second term up to a fixed large time interval Δt = 5t₀. Due to the small density, we have to use a low variance estimate for the explicit gradient in Eq. (14) for the specific case of optimizing the shear flow rate f (see  Appendix A).

The corresponding variational algorithm consists of a stochastic gradient descent optimization⁵⁵ for a large value of λ in Eq. (13). Starting from an initial point in parameter space {D_ij, f}, we simulate the dynamics of the system using Eq. (1) and, after relaxation into a steady-state, statistically estimate the explicit gradient of the variational estimator [Eq. (14)]. We perform stochastic gradient descent updating all variational parameters c at every step of the optimization, with the update rule at the nth step being

where the stochastic gradients are evaluated within a steady-state with the current value of the parameters and α_n is the learning rate for any of the c parameters in the nth optimization step. The level of noisy fluctuations in each parameter during the optimization process can be tuned independently through the corresponding learning rates. If the variational surface changes sharply, the learning rate has to be decreased with increasing n to anneal to the optimal solution basin. The learning rates have also been chosen individually for each example such that in each optimization step, the rate of change in D_ij/k_BT is in the range [0.1, 0.5] and that of fσ²/k_BT is in the range [1, 5].

To illustrate the performance of the optimization algorithm, we study the assembly of six particles into an octahedral (O_h) target cluster. An octahedron is not the highest yield cluster formed in a system of six hard sphere colloids with infinitely short-range attractions⁴⁰ and is formed in only 6% yield with strong, nonspecific interactions. Figure 2(a) shows the yield as a function of optimization steps with different learning rates and different trajectory noise histories. For all these examples, the yield is optimized over multiple orders of magnitude with the final converged value being close to 100%. This change in the order parameter over several orders of magnitude arises from the observable being defined as the probability of forming the target cluster, and in cases where the change is more drastic, it would necessitate the use of two different learning rates in Eq. (15). At a constant learning rate, the learning curves show two distinct regions such that a gradual rise in yield is followed by a rapid convergence to the saturation value.

Figure 2(b) shows the convergence of D_ij for one of the optimization runs. For six particles, there are 15 distinct interactions, all of which are optimized. The starting point is a nonspecific attraction D_ij = 4k_BT for all ij pairs. The optimization curve shows two regions, an initial spreading of the D_ij values followed by a rapid permutation symmetry breaking and a clear segregation of the 15 interactions into 12 attractive and 3 repulsive parameters. The 12 attractive interactions are all statistically equal, as are the 3 repulsive interactions. The attractive interactions correspond to the 12 bonds in the connectivity matrix for the octahedron. The symmetry breaking is spontaneous and is aided by the initial noisy fluctuations during optimization. Different noise histories in the trajectory lead to a symmetry breaking for which different sets of D_ij parameters become attractive or repulsive. For the finite clusters considered, this symmetry breaking is general. We refer to the specific D_ij solutions for the optimal yield of a target cluster as an alphabet and the particular D_ij in which there is a pairwise attractive interaction for every contact in the target structure as a Maximal Alphabet (MA). This strategy has been previously shown to be effective in the self-assembly of short-range interacting colloids into small clusters.^15,41

For the octahedral cluster, we have studied the stability of the MA solution for varying values of λ. We find a bimodal structure of the variational surface, with the algorithm converging to either a MA solution or a trivial solution u = F^WCA, depending on the value of λ. Figure 3(a) demonstrates the convergence of the octahedral yield with different values of λt₀ varying in the range [10³, 4 × 10⁴], starting from an MA solution with 10k_BT bond energies and a yield of 100%. For moderate but decreasing values of λ, the algorithm remains stable in the MA solution but with monotonically decreasing yields and bond energies. Below a critical value of λ_c = 5 × 10³/t₀, the MA solution becomes unstable and the algorithm finds the u = F^WCA solution. Rather than optimizing the yield in Eq. (13), at small values of λ, the second term is optimized. For some moderate λ values, the MA basin is only a local optimum and the crossover behavior shows λ-dependent hysteresis. We find this bistability of the variational surface to be generic.

Figure 3(b) shows the dependence of the average optimized bond energies of the converged MA solutions from Fig. 3(a), as a function of varying λ. This illustrates the depth of the MA basin in the parameter space. In the limit that each bond is formed independently, $Dλ*$ is expected to asymptotically vary as ∼2k_BT ln(λt₀), as demonstrated in  Appendix B. Over the range of λ considered, we find a logarithmic dependence but with a different coefficient, ∼1.3 ln(λt₀). This suggests that the free energy is approximately pairwise additive.

In the limit of large λ, which, in practice, is chosen such that the estimate of the first term in Eq. (13) is at least an order of magnitude larger than the negative second term, the variational algorithm can be used to automate the discovery of optimal forces for the self-assembly of clusters of arbitrary shapes and sizes in a nonequilibrium steady-state. The optimal forces stabilize the target clusters in an arbitrary ensemble without accounting for the dynamics of transient relaxation toward its steady-state. For a fixed number of tunable parameters, the computational cost scales linearly with the system size since the only bottleneck is propagating a steady-state trajectory long enough to compute statistically converged gradients. The algorithm also scales linearly with the number of variational parameters but with a small proportionality constant as all the gradients are estimated from the same trajectory. The use of the statistically estimated gradients significantly lowers the computational cost in contrast to numerically estimating the gradients from finite difference techniques by propagating multiple trajectories at different points in the parameter space. We next use our variational algorithm to study and rationalize the optimal design principles for a collection of rigid and nonrigid clusters.

We have investigated the formation of small low-energy rigid and nonrigid clusters with six, seven, or eight particles. We discover distinct design principles of these clusters and rationalize our findings by analyzing the response function of yield to the shear flow rate. We also demonstrate that the variational algorithm can obtain high yield optimal solutions even with constraints imposed on the total number of experimentally realizable design parameters. The design principles we obtain are expected to be general for the nonequilibrium self-assembly of short-range interacting colloids in a sheared steady-state.

We study the formation of a family of rigid clusters that are all known to be the lowest energy structures for systems of hard sphere colloids with identical infinitely short-range attractions. These clusters have previously been systematically enumerated and tabulated,^56,57 and their free-energy landscapes have been theoretically and experimentally studied.⁴⁰ These finite clusters are the colloidal analogs to small molecules and have been shown to be involved in the controlled seeding and growth of polycrystalline phases and kinetically arrested gels.^58,59

For each of these clusters, we optimize {D_ij, f} to extremize the yield within this force ansatz. In the limit that the attractive interactions between the particles were infinitely short-range, there would be no internal low-energy distortion modes, and the bond-connectivity matrix would correspond to a unique geometry. For our optimization, the indicator function refers to the corresponding bond-connectivity matrix conditions being satisfied. Figure 4(a) summarizes the design principles discovered for these clusters. The point groups for the symmetries of each of these clusters have been indicated along with the highest yields obtained. For chiral C_n clusters, the yields are the racemic yield.

For each of these clusters, the corresponding optimal alphabet discovered by the variational algorithm has also been indicated. We find that for clusters (i–iv), the optimal solution for D_ij is the MA. For the chiral C₂ cluster (v), the optimal alphabet is closely related to the MA but has a higher symmetry and is equivalent to a smaller 3 × 3 alphabet while having the same yield. For clusters (vi–ix), all of which contain a radial five-fold motif, the optimal yields are much less than 100% and the optimal alphabets are not MA. The reason is the competition with structures with a higher number of bonds. These lower energy structures would be geometrically unfeasible for infinitely short-range attractions; however, in our model, the short-range bonds have nonzero vibrations, which is sufficient to lead to the formation of the extra bonds. Unlike MA, the optimal alphabet discovered by the variational algorithm penalizes the formation of these lower energy competing structures.

Figures 4(b) and 4(c) show two slices through the optimization landscape in the parameter space of {D_ij, f}. Figure 4(b) is a diagonal slice through D_ij such that all D_ij = D while fixing f = 0. We find that for the two six particle clusters (i) and (ii), there is a monotonic increase in yield with increasing D. This suggests that even when the attractive forces are not infinitely short-range, both of these clusters are energetically the most stable configurations, and the only competing structures are higher in energy. When D = 10k_BT, the C_2v cluster (i) is formed with a yield of 93% compared to the 6% yield of the O_h cluster (ii), which is consistent with the stabilization due to the rotational entropy in the former.⁴⁰ All the other clusters (iii–ix) in Fig. 4(b) show a turnover in yield with increasing D. This is due to competing low-energy structures that are formed at large enough nonspecific D. We expect this design principle of a turnover in yield with increasing nonspecific attractions to be general for larger clusters since most larger clusters built from short-range interacting particles will contain the radial fivefold motif.⁵⁷ The value of the attraction D at the yield turnover is determined by a competition between the energetic stabilization of the lower energy structure and the destabilization of structures with missing bonds.

Figure 4(c) shows a slice through the optimization landscape with varying the shear rate f while fixing D_ij at their optimal value found by the variational algorithm. We see that for this class of rigid clusters, the yield monotonically decreases with the increasing shear rate. We expect this feature to be general for rigid clusters since rigid clusters need no additional geometric stabilization that can be provided by shear, which only energetically destabilizes the bonds in the cluster. This perspective is confirmed in Sec. III C using a linear response theory.

For our ansatz of short-range interacting colloids, clusters with N particles but fewer than 3N − 6 bonds in total, and fewer than three bonds for every particle, are not minimally rigid in that they have zero energy deformation modes.^56,57 These clusters are not formed in high yield as stable ground state structures in equilibrium. We have used the variational algorithm to uncover optimal nonequilibrium design principles for a family of such nonrigid clusters. The clusters we have considered belong to a family of planar two-dimensional structures known as polyiamonds. They have been shown to self-assemble from colloids in the presence of a spatial heterogeneity like in hydrodynamically driven assembly of sedimenting colloids in the presence of a substrate.^60–62 Within our control force ansatz, we investigate whether the shear flow planes are sufficient to stabilize these clusters.

For the optimization process, the indicator function for the cluster yield has been defined using both the bond-connectivity matrix and the flatness of the clusters, as discussed in  Appendix C. Figure 5(a) shows the optimal design principles obtained for clusters (x–xiv). The variational algorithm converges on the MA solution for the D_ij parameters for all of these clusters. The yields, however, are not 100% due to contribution from competing buckled configurations where the polyiamonds fold over the triangular faces to form tetrahedral motifs. Moreover, with the MA fixed, the optimal yield is at a non-zero shear flow. Figure 5(b) shows the yield as a function of the shear rate for fixed optimal alphabets. The location of the turnover in yield depends on the competition between geometric stabilization of the planar structure from the shear flow lines and energetic destabilization of the bonds. This design principle of planar two-dimensional clusters being stabilized by a shear flow appears to be general and stands in contrast to the rigid clusters that are strictly destabilized by shear.

The origin of the response of yield to shear flow is related to the relaxation dynamics of order parameter fluctuations in the unperturbed system. The response coefficients for rigid and nonrigid clusters can be understood using a generalized linear response theory.^63,64 Keeping the D_ij parameters fixed, the linear response of the yield to a change in shear flow rate can be decomposed into two terms,

where we have used the low variance estimator for the correlation function in  Appendix A. Here, Γ is related to the dynamical torque acting on the cluster,

and M refers to the virial stress on the cluster due to internal forces,

where $Fi=−∇i∑j≠iVWCA(rij)+VMorse(rij)$ is the conservative force acting on the ith particle. Γ and M are time-reversal asymmetric and symmetric parts, respectively,⁶⁵ of the full stochastic action gradient, R = Γ + M. Decomposing in this form is necessary to preserve the deconvolution of the center-of-mass motion from the internal coordinates of the cluster.⁶⁶ 

This linear response function is 0 at equilibrium due to the spatial parity symmetry of the system. Hence, we have characterized the different components of this correlation function at a small value of shear f = 5k_BT/σ² for the rigid cluster (ii) and the nonrigid cluster (x), fixing D_ij to their corresponding MA interactions. The results are shown in Fig. 6. We find that the component coming from the virial stress has opposite signs at small times for the rigid and the nonrigid cluster, which accounts for the opposite signs of the gradient of the yield. The tumbling motion of the clusters in a shear flow couples positively with the internal order parameter in the case of a nonrigid cluster and leads to an increase in yield with the increasing shear flow rate at small values of shear. The shear flow planes function as a spatial heterogeneity that is generally a precondition for stabilizing these planar clusters during self-assembly. At large shear, the yields of both rigid and nonrigid clusters are decreased with increasing shear due to the larger anti-correlation between the torque and the indicator function.

These design principles in- and out-of-equilibrium are general in their scope of applicability for small clusters formed by DNA-coated colloids. Nevertheless, a key limitation of this approach is the linear system size scaling of the number of different kinds of DNA-labeled colloids required in order to assemble a cluster, evident in the corresponding quadratic scaling of the number of variational parameters. We have addressed this limitation in Sec. III D.

Engineering a system with an extensive number of specific interactions is difficult, even with DNA-coated colloids. It is advantageous in this regard to uncover alphabets that code for the minimal sufficient interactions to stabilize a target structure in such a way that does not increase with the increasing system size. For example, polymers and crystals are macroscopic structures that can be assembled with a finite number of specific interactions, as both only require a repeating microscopic number of components to be stabilized, either a sequence of monomers or a unit cell. For clusters that do not have a clear repeating unit, discovering optimal design principles with smaller alphabets is nontrivial.

We have studied this problem by considering the six different low-symmetry eight-particle clusters (xvi–xxi) shown in Fig. 7. For each of these clusters, there is no direct way to partition the interactions into two or four classes based on their bonding environment or symmetry. We have used the variational algorithm to optimize the yield of each of these clusters, with a two-particle, four-particle, and eight-particle alphabet, in which D_ij has 3, 10, and 28 independent variational parameters, respectively. Clusters (xvi–xix) have near 100% yields for the full sized alphabet. Clusters (xx) and (xxi) compete with higher bonded clusters containing the radial fivefold motif and so have an optimal yield of lower than 100% even with an eight-particle alphabet. Nevertheless, for all the clusters, a four-particle alphabet can give quite large yields in comparison to the maximum possible yield. The variational algorithm identifies the optimal way to partition the groups of interactions of these clusters despite the lack of clear symmetry.

For the C₁ cluster (xix), even a two-particle alphabet has a high yield, despite not having any exact twofold symmetry. In this case, the algorithm has recognized a near-symmetry in the cluster and has partitioned it into two groups. The symmetry of these letters is close to the symmetrical alphabet that was discovered by the algorithm in a related C₂ cluster (v) in Fig. 4. We expect this potential to discover optimal design principles for large clusters with a small number of groups to be promising toward the self-assembly of experimentally realizable systems with practical constraints on the limits of the bottom-up design.

We have developed an inverse design algorithm for the self-assembly of colloidal clusters in a nonequilibrium steady-state. The formalism exploits a variational structure originating from large deviation techniques for importance sampling in trajectory ensembles. The algorithm optimizes the yield of clusters of arbitrary shapes, sizes, and geometry by tuning control forces within an arbitrarily chosen ansatz, with statistically estimated explicit gradients. We have demonstrated the performance of the algorithm using an ansatz of DNA-labeled colloidal clusters self-assembled in a shear flow and have obtained design rules for different families of rigid and nonrigid clusters. This algorithm scales linearly both in the system size and in the number of optimizable parameters in the force ansatz, but its performance is independent of the specific order parameter chosen for the self-assembly process. For example, the choice of a locally defined structural order parameter such as the density or the degree of crystallinity as the optimized observable can produce design principles for the assembly of extended dispersed or periodic structures out of equilibrium. Similarly, dynamical order parameters such as the instantaneous flux between two stable states can also be optimized using the same variational procedure in a suitable trajectory ensemble. Hence, this algorithm can be used to tune both structural and dynamical order parameters of clusters in a nonequilibrium steady-state to produce dynamical phases having no equilibrium analogs.

This variational algorithm differs from other available inverse design algorithms for soft matter in and out of equilibrium. The equivalent variational structure in configuration space for systems in thermal equilibrium, where the potential energy function is optimizable and explicit gradients can be statistically estimated by autodifferentiation, has been used extensively as the basis of both importance sampling and inverse design algorithms.^67–69 This configuration space approach with explicit gradients is not feasible in nonequilibrium systems due to the probability measure being non-Boltzmann. Out of equilibrium, there have been theoretical approaches to rationalizing design principles in one or two component systems.^31,70,71

In the absence of a closed form expression for the configuration space measure, machine learning algorithms have been previously used to identity optimal design principles by tracking an order parameter during or at the end of finite-duration trajectories.^72–75 Machine learning or neuroevolution based approaches are equivalent to estimating numerical gradients in the design space using finite-difference methods and have similar convergence properties as explicit gradient based methods in the limit of small optimization steps.⁷⁶ However, since our multidimensional statistical gradient estimates are obtained using information from the same trajectory, our explicit gradient based method is expected to be advantageous in a high-dimensional design space as typically encountered in the bottom-up design of soft materials.

Finally, a class of design algorithms employ a trajectory ensemble based approach to statistically estimate explicit gradients of the dynamical response in colloidal systems^77,78 and are formally closest to our approach. These algorithms probe the transient dynamical response of self-assembly trajectories, which, however, do not predict the long-time properties of the self-assembled cluster in a nonequilibrium steady-state. In our algorithm, we have explicitly evaluated the long-time steady-state limit and arrived at the novel fluctuation–dissipation relation in design space in Eq. (14). This will enable our explicit-gradient based method to be directly used to optimize both structural and dynamic properties of driven phases of soft matter and automate the discovery of new functional materials.

A.D. and D.T.L. were supported by the U.S. Department of Energy, Office of Basic Energy Sciences through Award No. DE-SC0019375.

The data that support the findings of this study are openly available in Zenodo at https://doi.org/10.5281/zenodo.4289235 (Ref. 79).

Using Eq. (14), the second term in the explicit gradient with respect to the shear flow rate f is

where we have simplified the stochastic action using the equation of motion. Since z_i appears in the expression independently for each particle and the density of the particles is vanishingly small, the z-diffusion timescale of the cluster diverges, and the correlation function takes a long time to converge. Thus, any gradient estimate we obtain by integrating the correlation function to a finite time Δt will contain a systematic error. In order to obtain an unbiased gradient, we recognize that in the large λ limit we are working in, the major part of Ω(t) is from 1(t), which by our definition depends only on the internal coordinates of the particle, and due to the spatial translation symmetry in our system, it is decoupled from the center-of-mass diffusion. This decoupling is directly expressed by a regrouping of terms in the sum over particles,

where in the first term the z coordinate of the center of mass has been explicitly factored out. We identify that the indicator function does not correlate with the center-of-mass motion and so the first term is 0. We use only the second term to approximately evaluate the gradients of Eq. (13) with respect to the shear flow rate f.

We analytically solve for the asymptotic dependence of the bond energy D as a function of the bias λ, for the formation of one bond, independent from the dynamics of the other bonds. We can represent this simplified system as an equilibrium, two-state Markov model with the rates of transition between the bonded and unbonded states k_b and k_u, respectively. These rates are determined by their mean ν = (k_u + k_b)/2 and their ratio k_u/k_b = exp(ΔF/k_BT), where ΔF is the free energy difference of the two states.

The cumulant generating function associated with fluctuations in the bonded state satisfies a eigenvalue equation of the form,

where r_λ is a right eigenvector and operator W_λ is given by

which is equal to the adjoint of the transition rate matrix when λ = 0. For the optimal rates that generate the statistics equivalent to this rate matrix, we need Doob’s transform of the matrix.³³ For this purpose, we diagonalize the matrix to find the eigenvector corresponding to the dominant eigenvalue as (r₁, 1), where

and k_u and k_b have been rewritten with ΔF and ν.

The modified rates that generate the optimal dynamics are given by $k̃u=ku/r1t0$ and $k̃b=kbr1t0$⁠. The modified free energy difference corresponding to these rates is $ΔF̃=kBT⁡ln(k̃u/k̃b)$⁠. Using Eq. (B3), in the λ → ∞ limit, the optimal free energy goes as

In this limit, the free energy is dominated by the negative of the bond-formation energy D, and hence, the latter is asymptotically D ∼ 2k_BT ln(λt₀).

Shear flow enhances the yield of nonrigid clusters (x–xv) by stabilizing a planar geometry and suppressing buckling modes. Here, we have computed the probability distribution $P(cos⁡θ¯)$ of the average flatness, defined as $cos⁡θ¯=[cos(θABC)+cos(θDEF)]/2$⁠, conditioned on the correct bond-connectivity matrix for the cluster being satisfied. The angles, θ_ABC and θ_DEF, are defined in Fig. 8. The flatness is −1 for a perfectly planar geometry but increases due to buckling and bending of the nonrigid cluster. For defining the indicator function for the nonrigid cluster, we used a flatness cutoff of $cos⁡θ¯≤−0.8$⁠. We have looked at a population where the bond-connectivity matrix condition is satisfied, but the flatness is unconditioned. This is illustrated in Fig. 8 for the C_2h cluster (x). We have kept the D_ij forces fixed at the optimized alphabet and plotted the distribution of flatness at two values of shear, at equilibrium with f = 0 and also at f = 20k_BT/σ², which is close to the optimal value for the highest yield. We find that the planar geometry is a transient state at equilibrium, with the most probable states corresponding to the buckling of one or both of the angles θ_ABC and θ_DEF. The shear flow destabilizes the buckled conformations and stabilizes the planar state instead so that at f = 20k_BT/σ², the most probable conformation is the correct planar geometry.