Design of two-dimensional particle assemblies using isotropic pair interactions with an attractive well

The manufacture of functional materials with specific nanoscale structural features presents considerable scientific and engineering challenges. While top-down fabrication approaches (e.g., lithography) have been advanced to address such challenges, they are often too expensive or time consuming for adoption in industrial manufacturing applications.^1–3 Bottom-up approaches such as self-assembly, on the other hand, stand as promising alternatives in which material building blocks (nanoparticles, block copolymers, etc.) might be designed–through modification of their effective mutual interactions^4–7–to spontaneously self-organize into a state that exhibits the desired microstructural features.^8–11 To determine which interactions stabilize a desired self-assembled structure, a design strategy is needed. Traditional ‘forward’ approaches discover promising material combinations through screens (i.e., combinatorial searches) or more limited testing of candidate systems judiciously chosen based on physical intuition. Alternatively, ‘inverse’ design strategies provide a more direct means for discovering interactions suitable for stabilizing the target structure, typically via solution of a constrained optimization problem.^12–15 

A classic example of an inverse design problem for materials is the discovery of an isotropic pair potential ϕ(r; {α}) between particles that stabilizes a specified crystalline lattice arrangement as the ground state (GS). Here, r is the distance between particle centers and {α} is the set of optimizable parameters that, along with the specified functional form, define the pair potential. Researchers have developed robust inverse design strategies that have been applied to discover isotropic interactions that stabilize a variety of open structures in two and three dimensions under various constraints (including honeycomb,^16–19 kagome,^20–22 simple cubic,^19,23 and diamond^23–25 lattices, to mention a few). Recently, a reformulation of this type of GS optimization problem was introduced²⁶ which significantly improved performance and allowed for (1) exploration of how design goals affect trade offs for the target phase (e.g., thermal versus volumetric stability²⁶) and (2) discovery of interactions that stabilize very challenging target structures (e.g., snub square,²⁷ truncated square and truncated hexagonal²⁸ lattices). An advantage of GS-focused optimization is that, because it requires the a priori identification of structures that most closely compete with the target lattice, it can offer insights into how those competitions influence the functional form of the optimized interactions. On the other hand, the identification of such competitors can be highly nontrivial, especially for target lattices with open structures. Furthermore, stabilization of the target as the GS does not guarantee that the structure will be kinetically accessible via a self-assembly process at higher temperature.

An alternative, simulation-based design method focused on relative-entropy (RE) maximization was recently shown^29,30 capable of discovering interactions that promote spontaneous assembly of a target lattice from the disordered fluid upon cooling. Notably, the RE optimization computes interactions “on-the-fly” based on structures accessed in a simulation and thus avoids the need to explicitly identify possible competing structures in advance. As a result, the RE approach has been able to help design isotropic interactions that favor unusually open ordered phases (e.g., truncated hexagonal²⁹ and truncated tri-hexagonal³⁰ lattices, which naturally compete with a variety of stripe-phase structures that are difficult to identify a priori) as well as disordered hierarchical structures (e.g., porous mesophases^30,31 and cluster fluids,^30,32 which are not easily designable with GS-based optimization strategies). However, because RE maximization does not require identification of competing structures, it cannot provide direct insights into how the optimized interactions stabilize the target relative to its competitors.

Here, we use both GS and RE optimization strategies to investigate a thus far unanswered question: which open two-dimensional crystal structures can be stabilized by an isotropic pair interaction comprising a repulsive core and a single attractive well? Interactions of this form are ubiquitous in nature, e.g., those present in noble gases and liquid metals,³³ the former of which have thermodynamic properties that are well captured by the familiar Lennard-Jones model. Single-well (effective) interactions can also arise in colloidal systems where, e.g., depletion interactions, van der Waals forces, and screened Coulomb interactions are present and can be tuned via the material selections made for the particle core, surface-passivating polymers or ligands, co-solutes, and the solvent.^5,34–36 Given the diverse contexts in which single-well interactions naturally emerge and can even be systematically modified, it is of interest to theoretically explore how the combination of a repulsive core and an attractive well might be chosen to favor various targeted crystal structures.

Our specific approach is to study the inverse design of isotropic, single-well pair potentials to stabilize 2D crystal structures (square, honeycomb, kagome, truncated square, and truncated hexagonal) using GS and RE methodologies described in section II. In section III, we present our results in the context of ‘design rules’ that can be inferred from our study of square, honeycomb, and kagome lattices. We then apply these design rules to discover single-well interactions that stabilize the more challenging truncated square and truncated hexagonal lattices. We conclude with section IV, where we summarize the rules and provide some further thoughts about their application to the inverse crystal design process.

In this method, we seek an isotropic, single-well interparticle pair potential ϕ(r; {α}) that stabilizes a target lattice l_t as the ground state over (at least) a narrow density range. Broadly speaking, we formulate the design as an analytical nonlinear program and determine the pair potential parameters {α} via its numerical solution using GAMS (General Algebraic Modeling System).^37–39 

For the model pair interaction, we choose the functional form ϕ(r; {α}),

where r = r/σ, A, n, σ₀, B_i, r_i and σ_i are design parameters (i.e. {α}), r_c is the cut-off radius and $fshift$(r) = Pr² + Qr + R is a quadratic shift function added to enforce ϕ(r_c) = ϕ^′(r_c) = ϕ^″(r_c) = 0. Location of the potential minimum r_min is fixed as

and an appropriate well profile is achieved by means of the following constraint equation

where r_i denotes an appropriately distributed set of radial points left(i = l) or right(i = r) of r_min. The magnitude of the minimum was fixed such that ϕ(r_min) = −ϵ. Finally, we constrain B₁ and σ₀, σ₂ as

To simplify notation, quantities are implicitly nondimensionalized from this point forward in terms of the usual (dimensionally appropriate) combination of parameters, including the energy scale ϵ, the length scale σ, or the Boltzmann constant k_B.

We compute quantities of interest over a narrow density range using the canonical ensemble. In particular, we consider a grid of densities {ρ_j} that spans a range Δρ = 1.21 − 1.23 to compute the potential energy U of the target lattice and equidensity competing structures,

where the sum is over all lattice-dependent coordination shells n_i at distances r_i up to the pair-potential cut-off. Similarly, the potential energy for phase-separated competitors U_PS was computed as⁴⁰ 

where x is the molar fraction for the crystal phase l₁, U_PS; (ρ_j) is the total system energy at net density ρ_j and U_i is the potential energy for crystal l_i at density ρ_i (see supplementary material for further discussion of this equation). Due to the narrow density range, energies of phase separated structures at only the midpoint density of the range ρ_m were sufficient to consider for optimization purposes. Stability of the target lattice relative to these competitors for the remainder of density points within the range was verified by means of forward ground-state phase diagram calculations, as discussed below.

Finally, we define an objective function F that maximizes the sum of the energy differences of the target and competitors at the midpoint density ρ_m = 1.22 as

where U_t(ρ_m) denotes the target’s energy and $Uli(ρm)$ the energy of a competitor lattice l_i. Stability of the target lattice for the remaining grid points {ρ_j} was ensured by constraining the energy differences for every competitor l_i as

Below, we discuss how structures that closely compete with the target are identified.

The pool of closely competing structures for the GS optimization was determined in the following way. First, an initial competing pool comprising lattices commonly stabilized by isotropic potentials was chosen, and the optimization was performed based on comparisons of the target lattice with that pool. The resulting optimized pair potential was then used to a carry out a ‘forward calculation’ in which the stability of the target lattice was determined relative to a comprehensive list of possible competing lattices (including those with variable parameters chosen via optimization using GAMS as well as lattices co-existing at different densities). Any structures found to be more energetically stable than the target were then added to the competing pool for a subsequent optimization. This process was repeated until the forward calculation of the optimized pair potential yielded no new competitors and established the target as the energy minimum, a point at which the optimization procedure was considered converged and completed. The detailed calculation procedure for single phases has been explained in detail previously.²⁷ A discussion of how to consider phase-separated crystals is provided in the supplementary material.

Following this approach, the final competing pools for square, honeycomb and kagome lattices were determined to be as follows:

• Square target: hexagonal, kagome, snub hexagonal, snub square, elongated triangular, and rectangular (b/a = 1.17) lattices, as well as phase-separated lattice structures comprising square and hexagonal lattices which we denote by (square lattice density)/(hexagonal lattice density) and the corresponding square lattice molar fraction: (1) 1.1479/1.4344, 0.7040; (2) 1.2675/0.7008, 0.9518; (3) 1.2512/0.7095, 0.9665; (4) 1.1833/1.44961, 0.83636; (5) 1.2279/0.7050, 0.9912.

• Honeycomb target: hexagonal, square, kagome, snub hexagonal, elongated triangular, rectangular (b/a = 1.59) lattices, and three specialized competitors, two of which resembled a distorted honeycomb and a staggered elongated triangular (see figure 1S in the supplementary material). There was a single phase-separated structure composed of honeycomb and hexagonal crystals which we denote by (honeycomb density)/(hexagonal density) and honeycomb molar fraction: 0.49180/1.86031 0.18861.

• Kagome target: hexagonal, square, snub hexagonal, elongated triangular, honeycomb, kagome (b/a = 1.02), and twisted kagome lattices, as well as a specialized competitor resembling rows of staggered rectangular boxes (see Fig. 2S in the supplementary material). No phase-separated competitors were found for this target.

Relative entropy course graining relies on maximizing the probability of achieving a desired equilibrium configurational ensemble by optimizing a set of parameters θ = {θ₁, θ₂, …, θ_m} of the interaction pair potential u(r|θ). This can be achieved by direct (“on-the-fly”) optimization in a simulation,²⁹ where potential parameters are periodically updated as

where i indicates iteration step, g(r|θⁱ) is the current radial distribution function, g_tgt(r) is the target lattice radial distribution function and γ is an adjustable scalar parameter to ensure stability and convergence of the optimization. Detailed derivation of this update procedure is available in previous work.^29,30 In this paper, we use u(r|θ) consisting of Akima splines whose knots are the parameters. We constrain the splines to display the features of a single attractive well with a minimum located at r_min. If we consider a set of scalars r = {r₁, r₂, …, r_m} as denoting the position of the Akima knots, the well constraint is easier to implement by considering the parameters θ ≡ θ(r) as differences of the pair potential

where the sign has been inverted to resemble a force term and for simulation convenience. The desired interaction form can then be achieved by constraining positions to the left (r_l) and right (r_r) of the minimum as

Lastly, the minimum position itself is free to vary and is updated as follows

where |…| indicates absolute value, the i + 1 superscript indicates the next iteration round, and k subscripts indicates the knot position immediately to the left (k − 1) or right (k + 1) of the r_min position at the ith iteration. Derivation of this result is provided in the supplementary material.

The Relative Entropy (RE) optimization is implemented with the GROMACS 4.6.5 molecular dynamics package.^41,42 Briefly, a system of particles (N = 200 − 1000) interacting via a single-well spline pair potential is simulated in the canonical ensemble using a periodically replicated rectangular simulation cell with aspect ratio chosen to accommodate the target lattice. The system is initiated in the target lattice at T = 1, melted via a run of 2 × 10⁶ time steps at T = 1.5, and then slowly cooled back to T = 1 (over 5 × 10⁶ time steps). Radial distribution function statistics g(r) are collected over the last 10⁶ time steps and compared to the target’s structure g_tgt(r). The spline potential is then updated per Eq. 9 and the process is iterated in this manner. In practice, γ = 0.005 − 0.01 is sufficient for all crystal target structures considered here. Convergence is typically achieved in about 100-200 iterations, and optimization is considered complete once the self-assembled crystal remains stable up to T = 1.5.

Monte Carlo simulations were used to test for stability of target lattices of particles interacting via the potentials obtained from the ground-state optimizations. In particular, systems of N particles (N = 64, 40, 108 for square, honeycomb, and kagome targets, respectively) interacting via Eq. 1 with their respective ground-state optimized parameters were simulated in the canonical ensemble using a periodically replicated simulation cell with densities and aspect ratios chosen to best accommodate the target lattices. In each case, the target was melted to form a liquid at a high temperature and then either quenched (square lattice) or slowly cooled (over ∼ 5 × 10⁶ Monte Carlo steps for the honeycomb and kagome lattices) to low temperature to check for assembly of the target crystal. For the square, honeycomb, and kagome lattices, the aforementioned high temperatures were T = 0.5091, 1.2091, and 2 and the low temperatures were T = 0.1627, 0.5, and 0.5. For the 24 Monte Carlo simulations we performed for the three targets, we observed perfect assembly in 50% (square), 17% (honeycomb), and 100% (kagome) of the runs.

As demonstrated below (see Fig. 1–5), we were able to discover, using the inverse design methodologies of RE and GS optimization and the Monte Carlo simulations described in Sec. II, isotropic pair potentials with a single well that promote self-assembly of all targeted lattices. Below, we examine the resulting optimized pair potentials and briefly discuss them within the context of the ‘design rules’ that they suggest for this class of systems.

This rule of thumb states the intuitive idea that the range of an isotropic interaction potential must be large enough to distinguish (and hence provide energetic advantage to) the target structure compared to its competitors. While only strictly true at zero temperature, we find that it provides a helpful guide for designing isotropic pair potentials to stabilize and promote assembly of our target lattices. Evidence of this rule’s utility was recently demonstrated in an RE optimization study where the shortest range isotropic, repulsive pair potential for stabilizing each of a wide variety of two-dimensional lattices was determined.²⁹ In that work, it was shown that the square and honeycomb lattices could be stabilized by an interaction spanning only two coordination shells, while the kagome lattice required a potential spanning three shells. These results are consistent with the single-well pair potentials we obtain here from optimizations for the corresponding target lattices (see black bars denoting coordination-shell locations in figures 1–3 for square, honeycomb and kagome lattices respectively). As discussed below, the single-well potentials designed to stabilize truncated square and truncated hexagonal lattices also span the minimum number of shells previously reported for these structures.²⁹ Consistent with this rule, our attempts to design single-well pair potentials that stabilize any of these structures using shorter range interactions were unsuccessful, and simulations employing the resulting optimized potentials led to either amorphous (glassy) or competing crystalline structures. Consistency with this rule may also be found retrospectively in the other published works featuring the same design targets (e.g., square^17,19,26 and honeycomb^17–19) and was briefly hinted at in another recent inverse design study.²⁰ Thus, while the results of the present study are not the first to suggest this design rule, they solidify the notion that a minimum interaction range must typically be spanned by an isotropic potential to favor a target structure compared to its competitors, an observation which acts as the foundation for the next two rules.

Given the observation that an isotropic pair potential requires a minimal range to stabilize a given target structure relative to its competitors, one might ask whether one can expect that a simple (e.g., Lennard-Jones-like) pair potential with a wider range of negative energy (i.e., a wider well) that spans the first and more distant coordination shells in the target lattice can stabilize any of the open structures considered here. Based on the results of our inverse design study, the answer is no, save for the square lattice. The optimized potentials for the square lattice target (see figure 1) have their minimum in the first-neighbor shell, which may at first glance appear counterintuitive since the first shell of the closely competing triangular lattice has more neighbors (six versus four) and hence a more favorable first-shell contribution to the energy. However an energetic advantage of the square lattice is still accomplished with this interaction by having a well that is wide enough to reach the second neighbors of the square lattice but not the more distant second neighbors of the triangular lattice.

For the more complex honeycomb and kagome structures, we found that the attractive well feature had to be constrained to the last one or two (closely spaced) shells of the interaction range (see plots in figure 2 and 3). If the pair potential wells were allowed to incorporate nearer coordination shells as well, the optimizations and subsequent forward calculations would inevitably find that phase-separated structures involving a more highly coordinated and hence more energetically stable (typically triangular) lattice were more favorable than the single-phase target structure. This was true for either GS or RE optimization strategies, and the qualitative similarity of the designed pair potentials (and corresponding forces shown in the figure insets) points toward the generality of this observation.

Given the constraints on the location and width of the attractive well discussed in the rules above, there is still a question of how simple the repulsive profile (r < r_min) can be in a single-well pair potential and simultaneously stabilize the target structure. As can be seen from the pair potentials and forces shown in figures 2 and 3 for the honeycomb and kagome lattices, nontrivial features in the repulsive interactions naturally emerge from the GS and RE optimizations for low-coordinated target structures. As emphasized in recent inverse design studies using repulsive isotropic potentials,^26–28 these features are tied to the location of coordination shells in competing lattices and necessary to stabilize target structures relative to more highly coordinated lattices, including the triangular crystal, which are otherwise naturally favored by smooth, short-range repulsive interactions. Instead of shoulders, the presence of multiple positive wells in the potential (if allowed in the optimization) could also have the same effect. Importantly, the inclusion of a single well in the pair interaction does not obviate the need for shoulders or other nontrivial repulsive features in the potential for r < r_min.

Having identified the three qualitative rules above, we conclude by presenting their application to the design of the more challenging, highly open structures of the truncated square lattice and the truncated hexagonal lattice. Pair potentials stabilizing these structures with purely repulsive interactions have very recently been reported using both the RE^29,30 and the GS²⁸ optimization strategies. However, due to the known complexity and tediousness involved in targeting these open structures with GS calculations (specifically, in identifying the large set of relevant competing structures which include a multiplicity of stripe phases), here we limit the pair potential discovery to the RE method. As such, following rule 1, we know in advance from the published studies the expected minimum interaction range to stabilize these targets. Next, we anticipate that the single well feature must lie towards the end of the interaction well and not be too broad as per rule 2. Finally, we expect the stabilization to require the introduction of multiple repulsive shoulders in the pair potential before the well, as per rule 3.

Our designed single-well interactions for these target structures are displayed in figure 4. Indeed, as can be seen, both potentials display the expected features as per the applied ‘design rules’ and, save for the single well, closely resemble those published for the same targets in previous RE work (figure [3]²⁹). The resulting self-assembled structures using the optimized interactions are shown in figure 5. Note that when we constrained the attractive well to be centered in a shell closer to the origin but still spanning the minimum interaction range (i.e. a broad well encompassing more than the outermost coordination shell(s)), the design strategy failed to stabilize any of these target structures, which further confirms the validity of design rule 2.

Isotropic pair interactions comprising a repulsive core and a single attractive well are models for effective interactions that are ubiquitous in colloidal fluids and interesting for material design applications due to their possible tunability. Using GS and RE inverse design methods, we have inferred a set of ‘design rules’ that help to understand the properties of single-well pair potentials that can stabilize three distinct two-dimensional target lattices (square, honeycomb and kagome). These rules can be summarized as (1) the interaction range must span a minimum number of coordination shells to differentiate the target structure from its competitors, (2) the well must be relatively narrow and located toward the end of the minimally required interaction range, and (3) nontrivial repulsive features (e.g., shoulders) in the potential spanning specific coordination shells in the target and competing lattices are required to impose energetic advantages to the target structure.

We further examined two challenging, low-coordinated target structures using the RE inverse design approach: truncated square and truncated hexagonal lattices. The optimized pair potentials, which displayed features consistent with the aforementioned design rules, readily assembled into the target structures upon cooling in computer simulations. Given their applicability across a wide variety of target structures, potential functional forms, and design methods, we believe these rules will be equally applicable to other crystal design targets in 3D, where transferability of 2D results to 3D targets are known to apply in repulsive systems¹⁹ and greatly simplify the design process.

See supplementary material for an illustration of specialized honeycomb and kagome competitors, parameters of reported GS potentials as well as derivations of the binary phase competitor calculations and the well minimum position update scheme used in the RE method.

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