Communication: Inverse design for self-assembly via on-the-fly optimization

There is a growing appreciation that a diverse array of structural motifs can be stabilized in systems of particles interacting via isotropic pair potentials, including various microphases,^1–3 open crystalline lattices,^4–7 and quasi-crystals.^8,9 Discovery of such potentials has been facilitated by inverse methods of statistical mechanics, most commonly optimization algorithms that iteratively refine the form of the interaction to attain increasingly favorable structural or thermodynamic properties.^10,11 For cases where such strategies have been employed to find isotropic pair potentials that assemble into specific two- and three-dimensional (2D and 3D) lattices, optimization typically builds on a ground-state calculation, wherein interactions are sought that stabilize the target structure at zero temperature relative to relevant competing structures. Despite successful application of such methods to design purely repulsive interactions stabilizing multiple open lattices (including honeycomb,^5–7 square,^5–7 rectangular,⁴ and kagome⁴ in 2D and diamond,^5,12,13 simple cubic,^5,12 and fluorite⁴ in 3D), there are some notable drawbacks to these approaches. First, they are encumbered by the requirement to specify the pool of relevant competing structures, a list that is not fully known in advance and thus must be modified as the potential is updated in the optimization. Moreover, the target structure must be checked explicitly for mechanical stability with the pair potential. Finally, interactions designed to stabilize the target structure in the ground state are not guaranteed to readily assemble the target from the fluid state upon cooling from higher temperature.

In this communication, we report a molecular dynamics (MD) simulation-based, inverse optimization strategy—carried out at nonzero temperature—that iteratively and systematically refines the pair interaction according to the fluid and crystalline pair structures dynamically encountered during the assembly process. The approach, while encoding practical aspects of assembly of the target from the fluid, is also technically simple and easy to implement. To illustrate the power of the methodology, we successfully employ it to design isotropic and purely repulsive pair potentials to assemble six distinct 2D lattices, including two structures which—to our knowledge—have never been stabilized via a pair potential before.

Machine learning-based techniques have been successfully applied to the self-assembly of colloids, from the rational design of building blocks and templates needed to fabricate nanomaterials^15–17 to the elucidation of pathways involved in the assembly process.¹⁸ In this work, the optimization scheme we present is general and based on maximum-likelihood machine learning¹⁹ (called “relative entropy coarse-graining” in the statistical mechanics community^20–22). Within it, particle-particle interactions are tuned in order to maximize the likelihood of reproducing desired configurations (in this specific case, 2D periodic lattices). We have previously used a related optimization strategy based upon relative entropy coarse-graining to discover isotropic pair potentials that favor microphases with porous architectures of a prescribed size.³ 

One advantage of relative entropy coarse-graining is the expression of the pair potential u(r|θ) in terms of a functional form that is parametrized by m arbitrary, scalar parameters, θ ≡ [θ₁, θ₂, …, θ_m]; this allows for a variety of constraints to be straightforwardly placed on u(r|θ). For flexibility of the potential, we optimize the amplitudes of the knots for an Akima spline spaced at an interval (δr) of δr/σ ≈ 0.0056, where σ is the nearest neighbor crystal distance. In order to avoid unphysical oscillatory potentials,²³ we constrain the knots to be monotonically increasing with decreasing r (i.e., purely repulsive interactions), though other more complex functional forms and constraints are possible. For the case of isotropic pair potentials, the updates to the parameters that characterize the potential are derived from the difference in the radial distribution functions associated with the present potential g(r|θ) and the target simulation, g_tgt(r),

where i indexes the iteration and α is the learning rate to be set empirically by observing the optimization stability.²⁴ A derivation of the update scheme and additional details pertaining to the optimization algorithm can be found in the supplementary material.

The target configurations are generated by an MD simulation where the particles are pinned to their lattice positions via a quadratic confining potential;²⁵ all simulations contained ≥1000 particles and were performed in Gromacs 4.6.5^26,27 in the NVT ensemble²⁸ with periodic boundary conditions in the x and y directions. For notational convenience, we define our optimization temperature as T^∗; however, the actual outcome of the optimization is βu(r), which uniquely defines the u(r) at any given temperature (β = (k_BT)⁻¹) that yields the desired pair structure. For the analysis below, the u(r) given by T^∗ is used, where k_BT^∗ assumes a value of unity. Further details regarding simulations can be found in the supplementary material.

In contrast to prior work based on finding a potential for which the desired lattice is the ground state, we target potentials that self-assemble into the targeted lattice at finite temperature. Therefore, a key step in our procedure is to initiate each simulation from a disordered fluid state, thereby allowing all relevant structural motifs to compete as needed prior to crystallization, including structures such as disordered microphases that might not be easily amenable to inclusion in an explicit competitor pool. Moreover, mechanical stability is directly incorporated into the optimization framework by the presence of finite temperature thermal motion. Practically, the configuration from the previous step in the optimization is melted via heating to an empirically determined temperature of T/T^∗ = 1.5, and then subsequently cooled to T/T^∗ = 1.0 prior to collecting statistics for the g(r). We consider the optimization complete when the crystal is sufficiently stable that it no longer melts at T/T^∗ = 1.5.

In Figs. 1(a) and 1(b), we compare the radial distribution functions of the optimized (g(r|θ)) and target (g_tgt(r)) structures for both the honeycomb (HC) and square (SQ) lattices, respectively. Both show excellent matching in the peak positions over many coordination shells, even though the ranges of the potentials (u(r), in black) only span the first two coordination shells. We also show g(r) from the last step before any crystallization occurred in the optimization to demonstrate that the disordered fluid locally contains muted structural signatures of the lattice. This correspondence between disordered fluid and crystalline lattice structure allows for the former to provide implicit “competitor pool” information for the optimization. As the potential is optimized prior to crystallization and therefore the fluid state evolves, fluids that have structural signatures corresponding to other lattices (and therefore do not match g_tgt(r)) are suppressed (any incorrectly assembled structures are also explicitly penalized). However, because the fluid is globally disordered and highly mobile, issues such as phase boundaries and metastability that become problematic upon lattice formation do not inhibit proper sampling of phase space in the fluid.

In addition to HC and SQ, this optimization procedure was carried out for the rectangular lattice with an aspect ratio of three (R3), kagome (KG), truncated square (TS), and truncated hexagonal (TH) lattices. The resulting potentials were then simulated shortly at a sufficiently high temperature to fully melt the crystal, and then the simulation was very slowly cooled through the empirically determined melting-freezing transition range and then further cooled to T = 0. The resulting structures are shown in Figs. 2(a)-2(f), where it is clear that the optimizations were successful. Only small expected defects are present due to finite size of the periodically replicated simulation cell and any misalignment of the nucleated crystal with the cell. The corresponding potentials are shown in Fig. 3. The potentials are essentially repulsive shoulders, with one to three such features beyond the core. Shoulders originate from the monotonicity constraint, which prevents the development of attractive wells at specific coordination shells. Instead, u(r) develops features with stiff repulsive forces (1) to penalize the surmounting of the shoulder, while (2) yielding a strong “thermodynamic” pressure (i.e., the ensemble averaged force due to the influence of every other particle) that pushes two particles into the shoulder. In essence, attractive wells are forgone for stiff repulsions and higher pressures to build the strong, specific coordination needed for targeting a crystal phase.

While the optimization framework presented here does not guarantee that the assembled crystal is also the ground state, we can confirm that the desired lattice, where the particles are in their ideal positions, is the lowest energy state at T = 0 of the lattices featured in this work. In Figs. 4(a) and 4(b), we show the total energy, U, as a function of ρ for the six lattices studied here in addition to the triangular (TR) lattice, using u_HC(r) and u_SQ(r), respectively. The circle denotes the optimization density for the lattice, and we see that the desired lattice is the lowest in energy for a reasonable range about ρ in the optimization, though the values of ρ outside of the optimization point are not guaranteed to be mechanically stable. The other potentials display similar behavior (data not shown).

The importance of heating the system into the fluid phase at each optimization step, thereby incorporating the assembly process into the optimization, can be demonstrated in the context of the TH lattice. If optimizations are performed beginning from the crystalline state, then mechanical stability is accounted for since the crystal can fall apart in the MD simulation, but the role of competitor phases to the self-assembly process is not included. Simulated annealing of a potential optimized for the TH lattice without melting the crystal at the outset of each iteration resulted in assembly of a crystalline stripe phase as shown in Fig. 5. The energies at T = 0 of the resulting (defective) striped phase and the perfect TH lattice are similar, with the latter being more stable by 5%. However, the striped phase forms first and is sufficiently kinetically stable that it never transitions into the TH lattice, even with very slow annealing schedules. Therefore, the TH phase appears to be kinetically inaccessible via simulation, though it might in fact be the ground state. Beginning with a disordered state in every optimization step circumvents this difficulty entirely because the g(r) is collected from the assembled phase, and the update scheme therefore adjusts the potential accordingly if an incorrect structure is encountered.

In summary, we have introduced a simple, relative entropy based inverse design approach, which we demonstrated can successfully discover purely repulsive, isotropic pair potentials that favor assembly of particles into a wide variety of open 2D lattices. Because this method is based on standard, MD simulation techniques, there is no need to construct or update a large competitor pool during the optimization. Instead, structural information is simply encountered in the simulation and utilized in the optimization on-the-fly, via self-assembly, to target a given crystal–including states that would otherwise be hard to predict a priori or to incorporate into a competitor pool (such as periodic microphases). In addition to the major simplification from the implicit competitor pool, mechanical stability, and likewise some finite degree of thermal stability, is naturally encoded in the method. However in this approach, it is not known whether the desired lattice is the ground state; rather, the optimization gives insight into what structure will result from self-assembly at nonzero temperature for a given potential. Also, it should be noted that crystals with large kinetic barriers to self-assembly will likely be difficult to treat with this strategy, though it is also reasonable to suspect that such lattices will encounter real-world complications that may preclude their assembly in general—particularly for large nanoparticle to micron-sized colloidal systems.

Intriguing avenues for future work include optimizing other constrained functional forms for u(r), and the use of different types of target structures, such as 3D lattices or binary lattices with multiple interactions. With respect to the former, one strategy to move towards realization of such assemblies would be to constrain u(r|θ) to experimentally motivated models that describe, for example, ligand-coated colloids.^29,30 Some existing interaction models, for micelles,⁸ for instance, bear resemblance to the simpler potentials presented in this work (HC and SQ in particular). However, simplification of the remaining potentials via greater constraint on the functional form might be necessary to realize the corresponding lattices.

See supplementary material for derivation of Eq. (1) and additional details regarding the simulation and optimization protocol.

T.M.T. acknowledges support of the Welch Foundation (F-1696) and the National Science Foundation (Grant No. 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.