Nathan A.
Mahynski
*^{a},
Runfang
Mao
^{b},
Evan
Pretti
^{b},
Vincent K.
Shen
^{a} and
Jeetain
Mittal
^{b}
^{a}Chemical Sciences Division, National Institute of Standards and Technology, Gaithersburg, Maryland 20899-8320, USA. E-mail: nathan.mahynski@nist.gov; Tel: +1 301-975-6836
^{b}Department of Chemical and Biomolecular Engineering, Lehigh University, 111 Research Dr., Bethlehem, Pennsylvania 18015-4791, USA

Received
9th December 2019
, Accepted 18th February 2020

First published on 19th February 2020

Inverse design methods are powerful computational approaches for creating colloidal systems which self-assemble into a target morphology by reverse engineering the Hamiltonian of the system. Despite this, these optimization procedures tend to yield Hamiltonians which are too complex to be experimentally realized. An alternative route to complex structures involves the use of several different components, however, conventional inverse design methods do not explicitly account for the possibility of phase separation into compositionally distinct structures. Here, we present an inverse design scheme for multicomponent colloidal systems by combining active learning with a method to directly compute their ground state phase diagrams. This explicitly accounts for phase separation and can locate stable regions of Hamiltonian parameter space which grid-based surveys are prone to miss. Using this we design low-density, binary structures with Lennard-Jones-like pairwise interactions that are simpler than in the single component case and potentially realizable in an experimental setting. This reinforces the concept that ground states of simple, multicomponent systems might be rich with previously unappreciated diversity, enabling the assembly of non-trivial structures with only few simple components instead of a single complex one.

The direct link between a material's macroscopic properties and its microscopic structure has driven extensive research into the design rules underlying rational bottom-up self-assembly.

An alternative and possibly simpler way to assemble complex colloidal structures involves the use of several “simple” components rather than a single “complex” one.^{24–27} By this we mean a set of several different colloidal species each with a relatively simple set of pairwise interactions as opposed to a system with a single species whose self-interaction is very complicated. While multicomponent mixtures can be tuned with conventional inverse design,^{28} this generally proceeds by optimizing a canonical system with a fixed number of particles at the same composition as the target; however, under experimental conditions this number can fluctuate and is more accurately represented by the grand canonical ensemble. Accurately representing the phase behavior becomes significantly more challenging than in single component systems. Gibbs' phase rule stipulates that, when temperature and pressure are fixed, coexistence may occur between as many different compositionally distinct phases as there are components. In spite of this, it is reasonable to expect that when the interparticle potentials are permitted to be arbitrarily complex, multicomponent systems can be found that yield the desired target and are stable with respect to phase separation. However, as constraints on interparticle potentials are imposed to make the designs more experimentally tractable, this becomes less plausible. As a result, there is a need for an inverse design scheme that can account for the possibility of phase separation in multicomponent mixtures.

Here, we report an inverse design scheme to enable the self-assembly of open, exotic lattices which derive their complexity from diversity rather than through chemically infeasible pair potentials. This scheme relies on a recently developed method to compute ground state phase diagrams for multicomponent colloidal systems,^{29} which allows the fitness of a desired target structure to be determined for a given Hamiltonian. In turn, this may be optimized using an active learning algorithm to design multicomponent systems with simple potentials that assemble into desirable porous structures. To illustrate this inverse design methodology, here we focus on two-dimensional (2D), binary systems with structures that have potential applications in, e.g., creating size-selective porous membranes; however, the underlying principles are extensible to three-dimensional and more multicomponent systems as well. We demonstrate this approach for several different lattices and show how it enables us to identify systems with simple, isotropic interactions that can be used to assemble non-trivial, open lattices.

U_{i,j}(r) = U^{r}_{i,j}(r) +λ_{i,j}U^{a}_{i,j}(r), | (1) |

(2) |

(3) |

Fig. 1 Construction of phase diagrams and a simple Hamiltonian. (a) Example of a ternary phase diagram constructed from the convex hull of specific free energy, A/N_{tot}, using the lowest energy candidates found at each mole fraction, or stoichiometry, searched. A k-component system yields a k-dimensional hull, and this construction is valid for all multicomponent systems. (b) Depiction of the pair potential (cf.eqn (1)) between species with a chosen λ value which varies smoothly from long-range attractive (λ = 1) to long-range repulsive (λ = −1). |

Pairwise interactions are varied by changing the value of λ_{i,j}, so that the Hamiltonian can be described by three parameters for a binary system, H = f() = f(λ_{1,1}, λ_{1,2}, λ_{2,2}). Otherwise, ε_{i,j} = σ_{i,j} = 1, which were used to non-dimensionalize all energy and length scales, respectively. Interactions were truncated at r_{cut} = 3σ.

Using these phase diagrams (cf.Fig. 2) we compute a target's fitness, F, to assess its degree of (meta)stability as a function of parameters in the Hamiltonian. We consider F to be a function of three contributions determined by the phase diagram: (i) the height above the convex hull of (free) energy which defines a structure's degree of (meta)stability, B; (ii) the composition range over which the target is stable, if it falls on the hull, Ω; and (iii) the convexity of the hull, C, at the target's composition. Thus, F = B × Ω × C.

The first and most important factor, B, we define akin to a Boltzmann factor of the (free) energy difference between the target and the convex hull at the given .

B = exp(−Δu), | (4) |

The second factor, Ω, weights B by accounting for the range of compositions over which a system will exhibit stable coexistence with the target structure. This is introduced to break “ties” between hulls which all contain the target, and would otherwise all have F = B = 1. A k-component system has a k-dimensional convex hull in (x_{a},x_{b},…,x_{k−1}, U/N_{tot})-space which provides a phase diagram by projecting into the first (k − 1) dimensions [cf.Fig. 1(a)]. We choose to define this weight as the fraction of the k-component system that should form the target integrated over mole fraction space, .

(5) |

This fraction, f_{target}() is simply given by the lever-rule. For example, the shaded blue region in Fig. 2(a) denotes the f_{target}() for the cyan target at x_{a} = 0.5; in panels (b) and (c), the shaded red regions denote the results for the red target at x_{a} = 0.2. In the binary case, S is given by the sum of half the differences in mole fraction between that of the target and each of its neighboring vertices on the hull. Ideally, only the target should exist on the hull and solutions with any composition in will exhibit coexistence with the target structure. In this case, S attains its highest value S_{max} = 1/k! (cf. ESI†):

(6) |

This function is bounded , where A ≥ 0. The lower bound occurs when S = 0 and corresponds to the case where the target structure is metastable so it does not belong to the convex hull of (free) energy.^{32–34} This mathematical form confers the important property that the fitness has a well-defined reference value,

(7) |

C = max[δu,0] + 1, | (8) |

The fitness and uncertainty may be calculated at any _{*} desired, after computing the covariance between all pairs of points in the prior, K, between the prior points and a point(s) of interest, K_{*}, and between the point(s) of interest, K_{**}. The predicted fitness of our target structure at _{*} is simply given by:

(9) |

(10) |

The hyperparameters in the kernel, , are determined by maximizing the marginal likelihood of the model^{36} to find the values that lead to the most likely of description of the observations. Active learning seeks to determine the global maximum fitness with a minimum number of new experiments (measured phase diagrams). We employed the regret-free upper confidence bound acquisition function given by:^{39,40}

(11) |

Thus, after fitting the current observations the next query is performed at _{next} = argmax(a). We set η = 2, thus using the 95% confidence interval to determine this. This cycle of fitting followed by prediction is repeated until a certain computational budget is exhausted or we have reached sufficient convergence, i.e., the optimal fitness value and location does not change significantly over time.

The very first point predicted [cf.Fig. 3(b)] was found to have F > 0.4 revealing that regression alone, with a reasonable prior, can be sufficient to locate regions of OHC stability. As the optimization proceeds, it explores regions of higher uncertainty, gradually finding marginally better coordinates until it converges to its optimum = (−1.0, 0.5017, −0.9325) ≈ (−1, 0.5, −1), which turns out to be the missing point from the prior grid. Fig. 3(a) shows representative hulls found during this process. Fig. 3(b) tracks the best fitness found over the course of the optimization procedure; in the inset, convex hulls are drawn around points with increasing F, revealing a structure to the data that shows a single peak around this point, reflecting reasonable behavior of the proposed fitness function. The hulls for the red, black, and blue points correspond to the same colors in Fig. 3(a). The red hull in Fig. 3(a) reveals that an alternating square lattice with x_{1} = 0.5, rather than the target OHC, is stable at that . Although the OHC does fall on the black hull, there are a number of competitors that manifest closer to x_{1} = 0.5 than in the optimal hull found (blue).

To test the accuracy of the ground state phase diagrams, MD simulations were performed, some of which are summarized in Fig. 3(c). More detailed results are available in the ESI.† The target OHC lattice is stable at the native stoichiometry of the lattice (x_{1} = 0.5) for the optimal , but as predicted by the phase diagram, loses its stability to a compositionally distinct competitor, namely a Kagome ring, at x_{1} = 0.4 (and also at x_{1} = 0.6). Other structures found along the hull are also reflected in the MD simulations, validating the accuracy of the ground state phase diagrams used in grand canonical inverse design.

However, we found that this single component system, containing multiple characteristic length scales in its optimal single-component pair potential, can be exchanged for a binary mixture whose compositional diversity replaces this complexity with a set of simpler interactions. In Fig. 4 we show the inverse design results for this lattice. For validation purposes, we have verified that this lattice can be successfully formed in MD simulations. We emphasize that it can be formed at the same or even higher temperatures used for the single-component case designed via relative entropy optimization (T* = k_{B}T/ε = 0.005).^{43} For 3 replicates, less than 5 optimization cycles were needed to initially locate stable regions of -space for this lattice, after which only marginal further improvement of the fitness was found. Notably, there are two stable regions found in Fig. 4(c), located symmetrically about the λ_{1,1} = λ_{2,2} plane. In this case, the fitness is multimodal, and the structure is stable only over a very narrow region of -space which would almost certainly be missed by a regular grid search.

The stabilization of this tile with simple, isotropic pair potentials is intriguing if one considers the smallest geometric element of the crystal. Another way to view this 4.8.8 Archimedean tile is as an assembly of square tetramers attached to their neighboring squares by their corners [cf.Fig. 4(a)]. Clearly, the 4-fold square symmetry is stable locally within each tetramer, and yet it is not stable globally; if it were, a square lattice, as depicted in red in Fig. 3(a), would result instead. We performed MD simulations using the optimal parameters to investigate further and found that the 4.8.8 tile could be reliably formed across a range of different temperatures; however, it was nearly always accompanied by the formation of some amount of OHC. Fig. 4(d) shows a representative result as identified using the fast neighbor graph analysis algorithm.^{44} While the final result varies depending on the starting conditions, in the image shown, roughly half of the system has crystallized into each competing polymorph.

To understand this polymorphism thermodynamically, we conducted free energy calculations for the periodic 4.8.8 Archimedean tile, the OHC, and the square lattice using the Einstein molecule approach^{45,46} as shown in Fig. 4(f). The free energy per particle, A/N_{tot}, is always lower for the 4.8.8 tile compared to both the OHC and square lattice over the entire temperature range. However, as the ground state is approached, the difference between the OHC and 4.8.8 tile specific free energies converges to approximately 10^{−3}ε/N_{tot}, in agreement with ground state optimization of their potential energies. This difference is similar to that separating the close-packed HCP and FCC polymorphs,^{47,48} which are known to crystallize into a polymorphic mixture due to this small free energy difference. Therefore, we attribute the apparent coexistence of the OHC and 4.8.8 tile to kinetic effects which trap some of the system in the less stable OHC phase. Regardless, the free energy calculations validate our inverse design, which accurately predicts the ground state energy rankings to relatively high precision.

As an alternative route to reducing polymorphism, it is possible to imagine an assembly procedure using prefabricated “metaparticles”^{24,49,50} composed of a set of different monomers irreversibly attached to one another. This multi-stage approach would first require the assembly of these building blocks from individual monomer units, then subsequent mixing of these metaparticles instead of a single stage approach in which the monomers are simply combined. We examined self-assembly using prefabricated rigid dimers and tetramers instead of simple monomers and found that nearly defect-free 4.8.8 Archimedean tiles can be formed by using rigid tetramers, as shown in Fig. 4(e) (cf. ESI†). This can be achieved either with rigid square tetramers, or with tetramers which have been minimized at the chosen Hamiltonian resulting in a slightly oblique rhombus. For comparison, we repeated the simulations with the parameters used for stabilizing square lattices (λ_{1,1} = −1, λ_{1,2} = 1, λ_{2,2} = −1) over the same temperature range. Only clusters of squares were observed regardless of the building block used (cf. ESI†). Thus, the 4.8.8 Archimedean tile can be formed not only in the ground state but also at finite temperature using the potentials found corresponding to the single-stage assembly directly from monomers, though polymorphism may be reduced by using a multi-stage approach.

- N. Vogel, M. Retsch, C.-A. Fustin, A. del Campo and U. Jonas, Chem. Rev., 2015, 115, 6265–6311 CrossRef CAS PubMed.
- M. N. O'Brien, M. R. Jones and C. A. Mirkin, Proc. Natl. Acad. Sci. U. S. A., 2016, 113, 11717–11725 CrossRef PubMed.
- M. A. Boles, M. Engel and D. V. Talapin, Chem. Rev., 2016, 116, 11220–11289 CrossRef CAS PubMed.
- M. E. Leunissen, C. G. Christova, A.-P. Hynninen, C. P. Royall, A. I. Campbell, A. Imhof, M. Dijkstra, R. van Roij and A. van Blaaderen, Nature, 2005, 437, 235–240 CrossRef CAS PubMed.
- S. C. Glotzer and M. J. Solomon, Nat. Mater., 2007, 6, 557–562 CrossRef PubMed.
- S. Sacanna, W. T. M. Irvine, P. M. Chaikin and D. J. Pine, Nature, 2010, 464, 575–578 CrossRef CAS PubMed.
- D. Klotsa and R. L. Jack, J. Chem. Phys., 2013, 138, 094502 CrossRef PubMed.
- V. N. Manoharan, Science, 2015, 349, 1253751 CrossRef PubMed.
- S. Whitelam, 2016, arXiv preprint arXiv:1606.00493.
- A. W. Long, C. L. Phillips, E. Jankowski and A. L. Ferguson, Soft Matter, 2016, 12, 7119 RSC.
- X. Zheng, Y. Wang, Y. Wang, D. J. Pine and M. Weck, Chem. Mater., 2016, 28, 3984–3989 CrossRef CAS.
- M. J. Solomon, Langmuir, 2018, 34, 11205–11219 CrossRef CAS PubMed.
- P. K. Bommineni, N. R. Varela-Rosales, M. Klement and M. Engel, Phys. Rev. Lett., 2019, 122, 128005 CrossRef CAS PubMed.
- E. Pretti, H. Zerze, M. Song, Y. Ding, R. Mao and J. Mittal, Sci. Adv., 2019, 5, eaaw5912 CrossRef PubMed.
- M. Rechtsman, F. Stillinger and S. Torquato, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2006, 73, 011406 CrossRef PubMed.
- M. Torikai, J. Chem. Phys., 2015, 142, 144102 CrossRef PubMed.
- B. A. Lindquist, R. B. Jadrich and T. M. Truskett, J. Chem. Phys., 2016, 145, 111101 CrossRef.
- R. B. Jadrich, B. A. Lindquist and T. M. Truskett, J. Chem. Phys., 2017, 146, 184103 CrossRef.
- A. Jain, J. R. Errington and T. M. Truskett, Soft Matter, 2013, 9, 3866–3870 RSC.
- A. Jain, J. R. Errington and T. M. Truskett, Phys. Rev. X, 2014, 4, 031049 Search PubMed.
- B. A. Lindquist, R. B. Jadrich, W. D. Piñeros and T. M. Truskett, J. Phys. Chem. B, 2018, 122, 5547–5556 CrossRef CAS PubMed.
- C. S. Adorf, J. Antonaglia, J. Dshemuchadse and S. C. Glotzer, J. Chem. Phys., 2018, 149, 204102 CrossRef PubMed.
- B. A. Lindquist, R. B. Jadrich, M. P. Howard and T. M. Truskett, J. Chem. Phys., 2019, 151, 104104 CrossRef PubMed.
- M. Gruenwald and P. L. Geissler, ACS Nano, 2014, 8, 5891–5897 CrossRef CAS PubMed.
- K. S. Khalil, A. Sagastegui, Y. Li, M. A. Tahrir, J. E. S. Socolar, B. J. Wiley and B. B. Yellen, Nat. Commun., 2012, 3, 794 CrossRef PubMed.
- W. M. Jacobs, A. Reinhardt and D. Frenkel, Proc. Natl. Acad. Sci. U. S. A., 2015, 112, 6313–6318 CrossRef CAS PubMed.
- M. Song, Y. Ding, H. Zerze, M. A. Snyder and J. Mittal, Langmuir, 2018, 34, 991–998 CrossRef CAS PubMed.
- W. D. Piñeros, B. A. Lindquist and T. M. Truskett, J. Chem. Phys., 2018, 148, 104509 CrossRef PubMed.
- N. A. Mahynski, E. Pretti, V. K. Shen and J. Mittal, Nat. Commun., 2019, 10, 2028 CrossRef PubMed.
- H. S. Ashbaugh and H. W. Hatch, J. Am. Chem. Soc., 2008, 130, 9536–9542 CrossRef CAS PubMed.
- N. A. Mahynski, H. Zerze, H. W. Hatch, V. K. Shen and J. Mittal, Soft Matter, 2017, 13, 5397–5408 RSC.
- A. Tkachenko, Proc. Natl. Acad. Sci. U. S. A., 2016, 113, 10269–10274 CrossRef CAS PubMed.
- P. G. Debenedetti, J. Chem. Phys., 1986, 84, 1778 CrossRef CAS.
- P. G. Debenedetti, Metastable liquids: concepts and principles, Princeton University Press, 1996 Search PubMed.
- S. Plimpton, J. Comput. Phys., 1995, 117, 1–19 CrossRef CAS.
- C. E. Rasmussen and K. I. Williams, Gaussian Processes for Machine Learning, MIT Press, 2006 Search PubMed.
- R. Garnett, M. A. Osborne and P. Hennig, 2013, arXiv preprint arXiv:1310.6740.
- D. R. Jones, J. Global Optimization, 2001, 21, 345–383 CrossRef.
- P. Auer, J. Mach. Learn. Res., 2002, 3, 397–422 Search PubMed.
- N. Srinivas, A. Krause, S. Kakade and M. Seeger, Proceedings of the 27th International Conference on Machine Learning, 2010.
- Q. Chen, S. C. Bae and S. Granick, Nature, 2011, 469, 381–384 CrossRef CAS PubMed.
- E. W. Cochran, C. J. Garcia-Cervera and G. H. Frederickson, Macromolecules, 2006, 39, 2449–2451 CrossRef CAS.
- W. D. Piñeros and T. M. Truskett, J. Chem. Phys., 2017, 146, 144501 CrossRef PubMed.
- W. F. Reinhart and A. Z. Panagiotopoulos, Soft Matter, 2018, 14, 6083–6089 RSC.
- C. Vega and E. G. Noya, J. Chem. Phys., 2007, 127, 154113 CrossRef PubMed.
- E. Pretti and J. Mittal, J. Chem. Phys., 2019, 151, 054105 CrossRef.
- L. V. Woodcock, Nature, 1997, 385, 141–143 CrossRef CAS.
- P. G. Bolhuis, D. Frenkel, S.-C. Mau and D. A. Huse, Nature, 1997, 388, 235–236 CrossRef CAS.
- W. M. Jacobs and D. Frenkel, J. Am. Chem. Soc., 2016, 138, 2457–2467 CrossRef CAS PubMed.
- M. B. Zanjani, J. C. Crocker and T. Sinno, Soft Matter, 2017, 13, 7098–7105 RSC.

## Footnote |

† Electronic supplementary information (ESI) available: More details on the fitness function, active learning algorithm, movies, and comparison of phase diagrams to MD simulations. See DOI: 10.1039/c9sm02426c |

This journal is © The Royal Society of Chemistry 2020 |