Isotropic pairwise interactions that promote the self-assembly of complex particle morphologies have been discovered by inverse design strategies derived from the molecular coarse-graining literature. While such approaches provide an avenue to reproduce structural correlations, thermodynamic quantities such as the pressure have typically not been considered in self-assembly applications. In this work, we demonstrate that relative entropy optimization can be used to discover potentials that self-assemble into targeted cluster morphologies with a prescribed pressure when the iterative simulations are performed in the isothermal-isobaric ensemble. The benefits of this approach are twofold. First, the structure and the thermodynamics associated with the optimized interaction can be controlled simultaneously. Second, by varying the pressure in the optimization, a family of interparticle potentials that all self-assemble the same structure can be systematically discovered, allowing for a deeper understanding of self-assembly of a given target structure and providing multiple assembly routes for its realization. Selecting an appropriate simulation ensemble to control the thermodynamic properties of interest is a general design strategy that could also be used to discover interaction potentials that self-assemble structures having, for example, a specified chemical potential.

Inverse design strategies are powerful tools for discovering simple interactions that promote spontaneous assembly of particles into states exhibiting specific structural motifs.1–27 Recently introduced approaches build upon molecular coarse-graining methods such as iterative Boltzmann inversion (IBI)28–30 or relative entropy (RE) optimization,28,31,32 which prescribe how to map a detailed model onto one with reduced degrees of freedom, e.g., replacing a multibody potential with a pairwise interaction or a group of atoms with a single bead. In design for assembly, this allows one to discover simple pair potentials that drive particles to sample equilibrium configurations that closely match those exhibited by a model with considerably more complex interactions that favor a target structure. Such inverse strategies have enabled the design of isotropic pairwise interactions that self-assemble particles into a rich variety of phases including equilibrium cluster fluids,6,7,10 porous mesophases,7,8,10 colloidal crystals,9–13 and fluids and gels comprising particle strings (“colloidomers”).14 

The simplest and most naïve coarse-graining strategy for determining such an interaction is direct Boltzmann inversion, which sets the optimized pair interaction between particles equal to the potential of mean force obtained from the target configurations.28 Unsurprisingly, such optimized potentials are unable to adequately reproduce the structure of most targets of interest at equilibrium,28 and more sophisticated, iterative coarse-graining strategies such as IBI or RE optimization are required for successful design. Moreover, since the coarse-graining process changes the interparticle interactions, there is an inherent trade-off between reproducing the structural characteristics of the target and matching its corresponding thermodynamic properties.33–38 For example, standard IBI seeks a pair potential that matches the radial distribution function g(r) of the target but not its pressure.35 It is possible to largely preserve agreement between the structure of the optimized and target ensembles while also bringing their pressures into accord by incorporating a correction term into the interaction.33–38 For instance, the pressure can be corrected in IBI by adding a linear ramp to the pair potential,34,35 while force-matching (another coarse-graining strategy) has a different, volume-dependent pressure correction term.36–38 Often, the specific form of the optimized potential is not particularly important for molecular coarse-graining applications, and modifying the pair potential in these ways does not generally pose problems.

In the context of inverse design for self-assembly, however, it may be necessary to constrain the functional form of the pair interaction in order to make contact with a particular experimental system or to learn about the fundamental requirements for a self-assembly process. The RE optimization framework is particularly well suited for this purpose because it optimizes the parameters of a specified potential. The form of the resulting potential, however, would not generally be preserved after application of a numerical pressure correction, motivating the need for other strategies to control the pressure for self-assembly design applications. Fortunately, the RE framework is not specific to any particular statistical mechanical ensemble,31,32 and while RE optimization is often performed in the canonical ensemble, it can also be carried out in the isothermal-isobaric ensemble to yield an optimized interaction that possesses the desired pressure by construction. For colloidal self-assembly, RE-optimized potentials are often thought of as effective interactions between colloids mediated by an implicit solvent; in this context, the corresponding osmotic pressure39 is controlled by the optimization. Perhaps more significantly, in applications of inverse design for assembly, the pressure can serve as a control parameter that can be tuned to bias the RE optimization toward a more net attractive (lower pressure) or repulsive (higher pressure) effective interaction, yielding a family of pair potentials that favor assembly of the same target morphology.

In this work, we demonstrate how the RE optimization strategy for inverse design is modified when simulated in the isothermal-isobaric ensemble. We validate the approach by revisiting a prior study where we optimized pairwise interparticle interactions in the canonical ensemble to assemble a rich variety of equilibrium mesophases (porous dispersions, lamellae, clusters, etc.).7 For this article, we use the equilibrium cluster phase that emerged from this earlier study as our target ensemble. Clusters, characterized by self-limited growth and preferred finite aggregate size, are of fundamental interest in the context of self-assembly40–65 but also are of practical relevance for understanding and controlling the viscosity of concentrated solutions of therapeutic proteins for subcutaneous injection.66–68 

The balance of the article is organized as follows. In Sec. II, we give the RE update equation applicable when the coarse-grained simulations are performed in the isothermal-isobaric ensemble instead of the canonical ensemble, in addition to describing other relevant computational details. In Sec. III, we first demonstrate that, for an appropriately chosen pressure, the protocol described in Sec. II recovers the same pair interaction as when the coarse-grained simulations are performed in the canonical ensemble. We then use the isothermal-isobaric optimizations to discover a family of cluster-forming potentials at different pressures and characterize the morphologies associated with the resulting interactions before concluding in Sec. IV.

RE optimization is equivalent to maximizing the likelihood that a given simulation protocol will sample the configurations characteristic of a target ensemble.9,10,31,32,69 The protocol is specified by both the total interaction potential UR|θ (where θ denotes a set of adjustable parameters and R indicates a particle configuration) and the thermodynamic ensemble in which the optimization simulations are carried out. When the target ensemble has a constant number of particles N and volume V, the canonical (“NVT”) ensemble is a natural choice for the RE simulation protocol, but other ensembles are also permitted. In the present work, we perform the optimization simulations in the isothermal-isobaric (“NPT”) ensemble, where the pressure P is constant, and the volume fluctuates. Our target ensemble is at a constant volume, so its volume distribution is effectively a delta function. The RE optimized model will approach a similar volume distribution to avoid poor overlap of configurations between the two ensembles; the pressure will match the desired value by construction.

We previously provided a derivation of the RE update scheme when the optimization simulations are performed in the canonical ensemble.9 The derivation for other ensembles proceeds analogously. The primary differences are (1) the probability of observing a given configuration R depends on the ensemble, (2) the averages that are taken over the simulation data correspond to the chosen ensemble, and (3) prefactors that depend on N or V in the update scheme may need to be explicitly accounted for instead of absorbed into the optimization learning rate. When maximum likelihood fitting is performed in the NPT ensemble and the form of the potential is taken to be an isotropic pairwise interaction denoted u(r|θ), the parameters at step i in the optimization, θ(i), are updated by

(1)

where η is the learning rate, g and gtgt are the radial distribution functions of the trial system in the current step of the optimization at density ρ and that of the target ensemble at density ρtgt, respectively, Vtgt is the volume of the target ensemble, and β=(kbT)1, where kb is Boltzmann’s constant and T is the temperature. A complete derivation for Eq. (1) can be found in the  Appendix.

In prior work, we have either used a static value for the learning rate η or infrequently adjusted η manually in response to the behavior of the optimization. However, for optimizations in the NPT ensemble, we have found that automated adjustment of the learning rate is essential. The magnitude of the update to the potential is determined both by η and by the difference between ρ and ρtgt. As a result, a value of η that might be appropriate when ρ and ρtgt are very close may yield an update that is so large that the optimization becomes unstable when ρ deviates more strongly from ρtgt. Since ρ can fluctuate frequently and rapidly over the course of an optimization, automated control over η is needed to maintain stability. In the  Appendix, we describe an empirical, automated procedure to update the learning rate in a way that leads to stable optimization.

To test the RE framework in the isothermal-isobaric ensemble, we optimized pair interactions to promote self-assembly of clusters. By working at the relatively low densities where clusters readily form, we can use a large simulation box with a modest number of particles. Clusters have the additional benefit of being straightforward to characterize via the cluster size distribution (CSD).64,65

To compute the CSD, two particles were defined as neighbors if their interparticle separation was less than a specified cutoff distance, rcut. Two particles are members of the same cluster if they are neighbors or if there is a contiguous pathway of neighboring pairs that connects them. The CSD is defined as the probability, p(n), that an aggregate contained a specified total number of particles, n. The pairwise nature of this analysis can lead to artifacts, particularly for larger values of rcut. For instance, two aggregates that are largely spatially separated could be classified as a single cluster if a single particle bridges them. To better reflect the concept of a cluster as a compact and relatively spherical object, we resolved groups of clusters that are defined as a single aggregate in the CSD on the basis of k-means clustering70–72 on the unwrapped particle coordinates of each aggregate identified by the CSD. The number of individual clusters in a given group (k) was taken to be the nearest integer to the quotient of the number of particles in the group by the most probable cluster size (primary peak in the CSD). The clustering was performed using the k-means++ algorithm for selecting initial conditions.73 

From the identities of the clusters, we computed their centers of mass in order to determine the cluster-cluster radial distribution function, gcl-cl(r). We also determined which lattice (if any) was formed by the clusters using Polyhedral Template Matching (PTM) with a cutoff of 0.15 for the root-mean-square deviation relative to an ideal lattice.74 PTM can distinguish between face-centered cubic (FCC), hexagonal close-packed (HCP), body-centered cubic (BCC), simple cubic, and icosahedral crystals and works well even at higher temperatures. The PTM calculations and visualization of the simulation configurations were performed using the Open Visualization Tool (OVITO), version 2.9.0.75 

The specific cluster target ensemble was inspired by prior work where we used IBI to discover an unconstrained pair interaction that prompted self-assembly of a variety of microphase-segregated states as a function of ρ.7 At ρσ3 values between 0.076 and 0.172 (where σ is the nominal particle diameter), the IBI interaction [Fig. 1(a)] prompted self-assembly of clusters. Over most of that range, the average cluster size increased as the system was densified. However, between ρσ3 values of 0.144 and 0.172, the average cluster size was relatively insensitive to density,7 perhaps reflecting the preferred cluster size encoded by the pairwise interaction in the absence of entropic effects. We therefore defined ρtgtσ3 as 0.144, using the low end of preceding density range for computational convenience.

FIG. 1.

(a) The dimensionless isotropic pairwise interaction, βu(r), discovered via inverse design using IBI in Ref. 7. This interaction drives self-assembly of particles into various microphase-separated states as a function of number density (ρ) in the canonical ensemble. [(b) and (c)] The corresponding radial distribution function and cluster size distribution, respectively, at ρσ3 = 0.144. (d) A snapshot of clusters that assemble onto a lattice using the interaction shown in panel (a) at ρσ3 = 0.144. See Ref. 7 for further simulation details.

FIG. 1.

(a) The dimensionless isotropic pairwise interaction, βu(r), discovered via inverse design using IBI in Ref. 7. This interaction drives self-assembly of particles into various microphase-separated states as a function of number density (ρ) in the canonical ensemble. [(b) and (c)] The corresponding radial distribution function and cluster size distribution, respectively, at ρσ3 = 0.144. (d) A snapshot of clusters that assemble onto a lattice using the interaction shown in panel (a) at ρσ3 = 0.144. See Ref. 7 for further simulation details.

Close modal

We generated the target ensemble by simulating this IBI-derived interaction at ρtgtσ3 = 0.144 and computed gtgt(r) [Fig. 1(b)]. The intracluster region extended out from r = σ to approximately r = 3σ, followed by a depletion zone between distinct clusters that led into cluster-cluster correlation peaks, the first of which occurred near r = 8σ. The minimum at r = 1.71σ signaled the termination of the first coordination shell on the particle level, and therefore, we used this value for rcut in the CSD calculations throughout. The target-ensemble clusters were relatively large and had good size-specificity [Fig. 1(c)]; the median cluster size was 55 particles. Additionally, the clusters themselves were arranged onto a lattice [Fig. 1(d)]. PTM indicated that the structure had large co-existing domains of BCC and HCP lattices that account for 43% and 40% of the clusters, respectively.

In our prior work,10 the target ensembles were usually obtained from simulations with physically unrealistic or multibody interactions, and RE optimization was used to replace these interactions with simpler pair potentials. Here, the target ensemble is obtained from an isotropic pair potential, but the form of the interaction is relatively complex and its features reflect many length scales. We can use RE optimization to replace this complex pair interaction with a more constrained functional form having simpler features.

To this end, we employed the following potential form in the RE optimization that was inspired by prior work7,10 but has been further modified to reproduce some of the qualitative features of the IBI potential more closely. The potential had a hard-core-like contribution given by the Weeks–Chandler–Anderson (WCA) potential,

(2)

with σ being the nominal particle diameter. The potential also had an auxiliary interaction that furnished an attractive well and a repulsive barrier,

(3)

where ϵ1 and ϵ2 are energy scales implicitly nondimensionalized by the thermal energy and α1 and α2 are also implicitly nondimensionalized by the nominal particle diameter σ. The first term provides an asymmetric attractive well that is deeper at smaller values of r and smoothly connects to the peak of a Gaussian repulsive barrier given by the second term; together these terms model the main qualitative features of the interaction shown in Fig. 1(a). The variables ϵ1, ϵ2, α1, and α2 are the adjustable parameters that comprise θ.

To facilitate a comparison between RE optimizations in different ensembles, we took special care to monitor the evolution of θ to determine when the optimization was complete, as opposed to simply stopping the optimization when the desired structure was obtained. Some fluctuations in the parameters always persisted as a consequence of the stochastic nature of the simulations, but we manually identified the point at which the parameters no longer systematically evolved. We then selected the step from this region that had the best g(r) matching, as quantified by the lowest value of γ(i) given by Eq. (A8).

All simulations were performed in gromacs, version 5.1.2.76,77 The time step was Δt=0.001βmσ2, where m is the particle mass, and periodic boundary conditions were enforced in all three dimensions. A canonical velocity-rescaling thermostat78 was employed to control the temperature with a time constant of τT = 100Δt. The NPT simulations used a Berendsen barostat with a time constant of τP = 500Δt and compressibility 7.2 × 10−4βσ3 to maintain a constant pressure. The Berendsen barostat was chosen because it is particularly stable even when the starting box size is not consistent with the pressure of the barostat, which was often the case as the potential evolved over the course of the optimization. However, unlike the employed thermostat,78 the Berendsen barostat does not produce fluctuations consistent with the appropriate ensemble.79 Therefore, we verified all NPT simulation results with NVT simulations after the optimization was complete, fixing the volume based on the average over the NPT simulation.

Since the simulations within the iterative optimization framework must be performed repeatedly (for each of the potentials presented here, hundreds of steps were generally required), we carried out relatively short and modestly sized simulations during the optimization. We then verified the results with larger and longer simulations outside of the optimization framework. During the optimization, the simulations contained 2000 particles and were run for 106 time steps each. The first 7.5 × 105 steps are taken to be the equilibration period, and any required quantities were measured over the remainder of the simulation. The NVT simulations were performed at ρtgtσ3 = 0.144; the NPT simulations were also started from ρσ3 = 0.144 but ρ evolved over the simulation. In order to further expedite the optimization simulations, we also used a somewhat aggressive value of ≈0.0025kbT/σ as the force threshold to determine the distance to truncate the pairwise interactions.

In order to verify the results of the optimization, we characterized the optimized potentials using NVT simulations with 10 000 particles and 2 × 107 time steps, where the first 107 time steps were defined as the equilibration period. Configurations were saved every 104 steps for subsequent analysis. We also used a more conservative value of 0.001kbT/σ in the force to determine the truncation distance, effectively including more of the repulsive barrier in the simulation from Eq. (3). As a result, the pressure measured from these NVT simulations was systematically slightly higher than the pressure input into the NPT optimization. To confirm that the pressure in the validation simulation was consistent with the NPT optimization, we computed a cutoff correction to the pressure,80 

(4)

where θopt indicates the optimal parameters from the optimization, rNVT is the cutoff in the NVT validation simulation, and rNPT is the cutoff from the NPT optimization. When the above term was included, the pressure from the validation simulation matched the targeted pressure to within 1.4% in all cases, except for the lowest pressure optimization, which we will discuss in Sec. III.

Finally, preliminary calculations indicated that using the final configuration from the previous step as the initial configuration in the iterative simulations caused kinetic issues with respect to the organization of the clusters. Clusters tended to become trapped in kinetically arrested and disordered configurations over the course of the optimization, whereas the clusters would facilely crystallize when initialized from a nonaggregated fluid state. Since the pressures associated with these two different arrangements of the clusters were different, the measured pressure from the NVT validation simulations did not match the pressure imposed by the barostat in the NPT simulations. In past work,9,12 we have used temperature annealing to avoid such issues, but large changes to the temperature at constant pressure change the box size significantly, which can cause practical issues in the simulation. Instead, we reset the initial configuration for every optimization simulation to a snapshot taken from an equilibrated WCA fluid simulation at ρtgt.

To provide a basis for comparison for the NPT RE optimizations, we first performed a RE optimization in the NVT ensemble. Since the paths through parameter space may be different between the two ensembles, we monitored the convergence of the parameters to determine when the optimization was complete. Both the parameters [Fig. 2(a)] and the pressure [Fig. 2(b)] stopped evolving meaningfully after about 210 steps. We also show the convergence criterion at each step, γ, in Fig. 2(c). Smaller values of γ indicate better overlap between the target ensemble and the simulation performed with the optimized potential [see Eq. (A8)]. We took the simulation having the lowest value of γ as the final result of the optimization, denoted by the square in Fig. 2(c).

FIG. 2.

The evolution of (a) Δθ(i) = θ(i)θopt, where θopt are the optimal parameters, (b) the nondimensionalized pressure P(i)βσ3, and (c) the convergence criterion γ(i) defined in Eq. (A8) as a function of step i in the NVT optimization at ρσ3 = 0.144. The step corresponding to the optimal interaction is denoted by the open square in panel (c). The (d) interaction potentials and (e) corresponding CSDs are shown for the self-assembled clusters for the steps denoted by the open symbols in (c).

FIG. 2.

The evolution of (a) Δθ(i) = θ(i)θopt, where θopt are the optimal parameters, (b) the nondimensionalized pressure P(i)βσ3, and (c) the convergence criterion γ(i) defined in Eq. (A8) as a function of step i in the NVT optimization at ρσ3 = 0.144. The step corresponding to the optimal interaction is denoted by the open square in panel (c). The (d) interaction potentials and (e) corresponding CSDs are shown for the self-assembled clusters for the steps denoted by the open symbols in (c).

Close modal

Prior to the convergence of the parameters, we found that many values of γ were comparable to that of the optimized potential, providing many interactions that appear consistent with the targeted cluster morphology. In Figs. 2(d) and 2(e), we show a subset of these interactions and their corresponding CSDs from the steps demarcated by all the open symbols in Fig. 2(c). The CSDs are all mutually consistent; therefore, we see that there is some flexibility in the interaction potential in terms of promoting self-assembly of clusters with the correct size. However, the pressure associated with the potentials [Fig. 2(b)] varies from P = 0.76–1.01kbT/σ3.

The final optimized interaction and its corresponding CSD and g(r) are shown in Figs. 3(a)–3(c), where they are compared to their respective measurements from the target simulation. The restricted functional form of the interaction prevented perfect matching between the two, with the optimized interaction producing slightly smaller clusters than the target (the median values are at 52 and 55, respectively). This is reminiscent of prior work on porous mesophases, where a similarly restricted functional form produced slightly smaller pores than the unrestricted functional form.7 The optimized g(r) was slightly less structured overall, particularly near r = σ, which is probably a consequence of the deeper and sharper attractive well present in the target interaction. Such a sharp feature in the potential is not realizable for the constrained functional form. The optimized interaction prompted self-assembly of clusters that crystallized onto a BCC lattice [Fig. 3(d)]. On average, 70% of clusters were classified as BCC. The coordination structure of most of the remaining clusters (25% of all clusters) was unknown in the PTM analysis, which is likely at least partially attributed to defects in the lattice that occur due to the finite size of the simulation box. (There were slightly fewer than 200 clusters in the larger simulations.)

FIG. 3.

At ρσ3 = 0.144, comparison of (a) the IBI-derived interaction potential used to generate the target ensemble (gray) and the RE-optimized potential in the NVT ensemble with the constrained functional form given by Eqs. (2) and (3) (dark orange), (b) CSD, and (c) g(r) from the target ensemble and the optimized simulation. (d) A snapshot of clusters that are arranged onto a BCC lattice using the optimized RE interaction shown in panel (a).

FIG. 3.

At ρσ3 = 0.144, comparison of (a) the IBI-derived interaction potential used to generate the target ensemble (gray) and the RE-optimized potential in the NVT ensemble with the constrained functional form given by Eqs. (2) and (3) (dark orange), (b) CSD, and (c) g(r) from the target ensemble and the optimized simulation. (d) A snapshot of clusters that are arranged onto a BCC lattice using the optimized RE interaction shown in panel (a).

Close modal

In order to compare between RE in the NVT and NPT ensembles, we measured the pressure associated with the optimized interaction in Fig. 3(a). We extended the simulation with the optimized interaction for another 106 time steps, from which we computed a pressure of P* = 0.788kbT/σ3. Using this value to control the pressure and setting ρtgtσ3 equal to 0.144, we carried out a RE optimization in the NPT ensemble, the results of which are shown in Figs. 4(a)–4(c). While the parameters in Fig. 4(a) stopped evolving more quickly in the NPT ensemble than the NVT ensemble as a function of optimization step, we note that the learning rate evolved with step differently in the two ensembles. In the case of the NPT calculation, the optimization seemed to heavily favor approximate matching of the density ρ to ρtgt [Fig. 4(b)]. The step associated with the final optimized interaction is shown by an open square in Fig. 4(c).

FIG. 4.

(a) The evolution of (a) Δθ(i) = θ(i)θopt, where θopt are the optimal parameters, (b) Δρ(i)σ3 = ρ(i)σ3ρtgtσ3, and (c) the convergence criterion γ(i) defined in Eq. (A8) as a function of step i in the NPT optimization, where P = P* = 0.788kbT/σ3.

FIG. 4.

(a) The evolution of (a) Δθ(i) = θ(i)θopt, where θopt are the optimal parameters, (b) Δρ(i)σ3 = ρ(i)σ3ρtgtσ3, and (c) the convergence criterion γ(i) defined in Eq. (A8) as a function of step i in the NPT optimization, where P = P* = 0.788kbT/σ3.

Close modal

Despite the different paths through parameter space, the NVT and NPT RE optimizations arrived at effectively identical solutions. Each final optimized potential and the corresponding CSD and g(r) are directly compared in Figs. 5(a)–5(c), where the results are nearly indistinguishable. The median cluster size was 52 particles in both cases. Both interactions promoted self-assembly of the same morphology, with 72% of the clusters from the NPT optimization assigned via PTM to a BCC structure [Fig. 5(d)] and 23% of the clusters unclassified. As should be the case if the simulation box is sufficiently large, performing the optimization in the distinct ensembles but with compatible choices for the thermodynamic variables yielded equivalent results.

FIG. 5.

Comparison of (a) the optimized potential from the NVT (ρσ3 = 0.144, dark orange) and NPT (P = P* = 0.788kbT/σ3, teal) optimizations and the corresponding (b) CSD and (c) g(r). (d) A snapshot of clusters that are arranged onto a BCC lattice using the optimized NPT interaction shown in panel (a).

FIG. 5.

Comparison of (a) the optimized potential from the NVT (ρσ3 = 0.144, dark orange) and NPT (P = P* = 0.788kbT/σ3, teal) optimizations and the corresponding (b) CSD and (c) g(r). (d) A snapshot of clusters that are arranged onto a BCC lattice using the optimized NPT interaction shown in panel (a).

Close modal

After validating the NPT optimization strategy, we subsequently varied the pressure to discover different interaction potentials that also resulted in a clustered microphase. The results for P from 0.5P* to 4.0P* are shown in Figs. 6(a)–6(c). As expected, potentials optimized at P < P* had a deeper attractive well and a more muted repulsive barrier; the inverse is true for P > P*. On the whole, there is significant flexibility with respect to the parameters if either the volume or the pressure of the system can be tuned although it is reasonable to assume that the parameters must be optimized in a coordinated fashion to produce the desired morphology. The CSDs resulting from the potentials in Fig. 6(a) are shown in Fig. 6(b), where we see that the fidelity with respect to cluster size is consistent with the NVT optimization (median cluster size of 52 or 53 particles in all cases), with polydispersity decreasing as the pressure increases.

FIG. 6.

(a) The optimized interaction potentials, βu(r), (b) the corresponding CSDs, and (c) the radial distribution functions between the center of masses of the clusters, gcl-cl(r), as a function of the pressure of the barostat in the NPT optimization for P = 0.5P*, 1.0P*, 2.0P*, 3.0P*, and 4.0P*, where the lighter colors indicate higher pressure and the arrow is drawn in the direction of increasing pressure.

FIG. 6.

(a) The optimized interaction potentials, βu(r), (b) the corresponding CSDs, and (c) the radial distribution functions between the center of masses of the clusters, gcl-cl(r), as a function of the pressure of the barostat in the NPT optimization for P = 0.5P*, 1.0P*, 2.0P*, 3.0P*, and 4.0P*, where the lighter colors indicate higher pressure and the arrow is drawn in the direction of increasing pressure.

Close modal

All of the above interactions resulted in clusters that primarily crystallized onto a BCC lattice. PTM indicated that between 55% and 72% of the cluster centers were in a BCC environment, and most of the remaining clusters were not classified as crystalline by PTM. The pair distribution functions of the cluster center-of-masses, gcl-cl(r) [Fig. 6(c)], were also consistent with self-assembly into a BCC lattice.

From the P = 4.0P* result, we note that attractions are not required to form a clustered phase, consistent with prior work.55–63 Analogously to previous work on pores and crystals comprised of single particles,10,12 the attractive well can be completely replaced with a shoulder while retaining the targeted morphology. While cluster-forming potentials possessing competing interactions and those characterized by purely repulsive interactions are sometimes discussed as two distinct classes of potentials, by tuning the pressure, we see that there is a continuous family of interactions that spans this space. With the exception of the P = P* result, the potentials shown in Fig. 6(a) do not correspond to stationary points for RE optimization in the NVT ensemble for the given gtgt(r). Therefore, performing the RE optimization in the NPT ensemble is critical to systematically obtain a family of interactions that may in general yield insights into the fundamental requirements for self-assembly of a given target structure.

We did not perform calculations at a higher pressure than P = 4.0P* because we found that the optimization did not converge (i.e., the parameters did not stop evolving) even after 1000 iterations. This may be due to some frustration in the optimization: the functional form qualitatively changes as ϵ1 becomes negative, i.e., the part of the function that generates the attractive well inverts to form another repulsive barrier and in doing so generates a rather abrupt attractive well between this barrier and the WCA component. Negative values for ϵ1 are disfavored in the optimization so that at sufficiently high pressures, ϵ1 stays close to zero and other parameters must compensate to maintain the appropriate box size for the given pressure. Other functional forms tailored to generate purely repulsive potentials may be a better choice for higher pressures than those investigated here.

We also performed NPT optimizations at even lower pressures than shown in Fig. 6. The results for P = 0.125P* and 0.25P* are in Fig. 7, where the P = 1.0P* result is also plotted for context. We see a similar deepening of the attractive well as the pressure decreases in Fig. 7(a). Conversely, while the repulsive barrier moves into lower r with lowering pressure, the magnitude of the repulsive barrier no longer significantly decreases from P = 0.25P* to P = 0.125P*. It is understood that for interactions possessing competitive attractions and repulsions that a minimum degree of repulsion is needed to form clusters in order to thwart macroscopic phase separation,52–54 and the form of the P = 0.125P* potential appears to be influenced by this limitation.

FIG. 7.

(a) The dimensionless optimized interaction potentials, βu(r), (b) the corresponding CSDs, and (c) the radial distribution functions between the center of masses of the clusters, gcl-cl(r), as a function of the pressure of the barostat in the NPT optimization for P = 0.125P* and 0.25P*, where P = 1.0P* is shown for comparison. The lighter colors indicate higher pressure, and the arrow is drawn in the direction of increasing pressure.

FIG. 7.

(a) The dimensionless optimized interaction potentials, βu(r), (b) the corresponding CSDs, and (c) the radial distribution functions between the center of masses of the clusters, gcl-cl(r), as a function of the pressure of the barostat in the NPT optimization for P = 0.125P* and 0.25P*, where P = 1.0P* is shown for comparison. The lighter colors indicate higher pressure, and the arrow is drawn in the direction of increasing pressure.

Close modal

Perhaps as a consequence of the above considerations, the self-assembled structures corresponding to the lower pressure potentials in Fig. 7(a) had some notable qualitative deviations from the optimizations at higher pressures. First, the primary peak in the CSD was shifted to slightly larger clusters, and the clusters were significantly more polydisperse [Fig. 7(b)]. For the case of P = 0.125P*, the clusters never crystallized during the validation simulation, continuing to exist as a fluid. (The structure of over 99% of the clusters was not identifiable in the PTM analysis.) At P = 0.25P*, the assembly did order somewhat (although a plurality of the clusters, 46%, were unclassified by PTM) into a polymorphic structure with significant degrees of BCC, HCP, and FCC crystal structures (34%, 12%, and 8% of the clusters, respectively). gcl-cl(r) for both low pressure optimizations are shown in Fig. 7(c), where the morphological differences are apparent.

We also noticed that at the lower pressures, the optimizations took noticeably longer to converge, with the P = 0.25P* optimization requiring about 500 iterations (200–250 iterations were typical at the higher pressures) and the parameters in the P = 0.125P* optimization failing to converge after 1000 iterations (although the convergence criterion did not improve over the last ≈800 steps). Additionally, while the optimal ρ value measured from the higher pressure results generally matched ρtgt closely (less than 1% difference), the final ρ value measured from the optimal P = 0.125P* simulation deviated more strongly (over 11%). The discrepancy between the optimal ρ value and ρtgt posed a challenge to our optimization strategy since we began each simulation at ρtgt; this issue manifests as the pressure measured from the NVT validation simulation varying more strongly from the corresponding barostat pressure than in any other optimization (≈7% difference). By contrast, when we begin the NPT simulation from a disordered state at the correct density for the interaction, thereby allowing the clusters to self-assemble while ρ is not evolving, we restore close agreement between the measured pressure of the NVT simulation and the barostat of the NPT simulation (within 1.5%). Given that we observed kinetic issues when the iterative simulations were not initialized from a disordered state (see Sec. II), it seems reasonable that equilibration may also be inhibited if self-assembly occurs while ρ is also evolving, which could explain the above discrepancy with respect to pressure.

On the whole, the low pressure regime investigated in this work underscores some of the limitations of tuning the pressure and still achieving self-assembly of a targeted structure. Based on the above, we suggest that the convergence of the parameters and the closeness of ρ for the optimized interaction to ρtgt may serve as indicators for what pressures are compatible with a given targeted structure. On the other hand, there may be applications where tuning either the pressure or the interaction potential to modulate structural morphology may be desirable, and locating regions of phase-space where the morphology of the optimal self-assembled structures varies in response to pressure may be useful toward this end. Furthermore, a systematic study of NPT optimizations as a function of pressure provides an avenue to determine the necessary range in pressure required in order to self-assemble a target structure via any specified form for the interaction potential—a practically useful guideline for realization of such self-assembled structures.

The primary aim of this article has been validation of an approach to optimize for interactions that possess desired thermodynamic properties as well as structural correlations. In the course of doing so, we have also discovered a series of interaction potentials that self-assemble into clusters that may be useful in future work in the field of clustering. In Table I of the  Appendix, we report the parameters that define all of the optimized potentials from this work and the corresponding ρ values at which the validation simulations were performed.

TABLE I.

Optimal parameters for the potential form defined by Eqs. (2) and (3) and the corresponding density as a function of pressure.

Pressureϵ1α1ϵ2α2ρσ3
0.125P* 1.345 1.201 0.095 2.111 0.128 
0.25P* 1.189 1.293 0.098 2.260 0.144 
0.5P* 0.980 1.441 0.121 2.461 0.144 
0.75P* 0.844 1.550 0.138 2.591 0.145 
1.0P* 0.727 1.621 0.159 2.671 0.144 
1.25P* 0.646 1.659 0.177 2.745 0.144 
1.5P* 0.567 1.682 0.194 2.799 0.145 
2.0P* 0.427 1.728 0.223 2.904 0.145 
3.0P* 0.179 1.742 0.294 3.035 0.145 
4.0P* −0.009 1.478 0.373 3.352 0.145 
Pressureϵ1α1ϵ2α2ρσ3
0.125P* 1.345 1.201 0.095 2.111 0.128 
0.25P* 1.189 1.293 0.098 2.260 0.144 
0.5P* 0.980 1.441 0.121 2.461 0.144 
0.75P* 0.844 1.550 0.138 2.591 0.145 
1.0P* 0.727 1.621 0.159 2.671 0.144 
1.25P* 0.646 1.659 0.177 2.745 0.144 
1.5P* 0.567 1.682 0.194 2.799 0.145 
2.0P* 0.427 1.728 0.223 2.904 0.145 
3.0P* 0.179 1.742 0.294 3.035 0.145 
4.0P* −0.009 1.478 0.373 3.352 0.145 

In this article, we have demonstrated how, by manipulating the simulation protocol used in RE optimization, we can directly optimize interaction potentials that possess a desired pressure while maintaining a constrained functional form. We showed that the RE optimization converges to the same solution in both the NVT and NPT ensembles when the pressure for the latter is chosen to be consistent with the output of the NVT optimization. Moreover, we have tuned the pressure in the NPT optimizations to generate a family of potentials that all form clusters of the correct size. As noted in prior work, specifying the pressure inherently degrades the matching between gtgt(r) and g(r) when the potential is infinitely flexible;35 however, this is not necessarily the case when the functional form of the interaction is restricted. We found comparable g(r) matching at all state points except for the two lowest pressure results that did not form a BCC lattice of clusters—these matched the target less well. In Fig. 8, we compare the g(r) computed from the P = 0.5P* optimization to that from the P = 4.0P* optimization. Despite the obvious differences in the optimized potentials [see Fig. 6(a)], the radial distribution functions are well-matched, with the higher pressure result being only slightly more structured.

FIG. 8.

Comparison of the radial distribution functions g(r) using the optimized interaction potentials from the P = 0.5P* and P = 4.0P* optimizations.

FIG. 8.

Comparison of the radial distribution functions g(r) using the optimized interaction potentials from the P = 0.5P* and P = 4.0P* optimizations.

Close modal

Using a barostatted simulation protocol within the RE framework provides a straightforward and computationally convenient means to control the pressure associated with the interaction that is outputted by the RE optimization. This protocol allows us to discover the interaction form that is most likely to reproduce the target ensemble subject to the pressure constraint as opposed to tuning an interaction potential that matches the pressure in an ad hoc fashion. Finally, control of thermodynamic quantities via RE optimization is not restricted to the pressure. For instance, the chemical potential of a coarse-grained interaction could be tuned by performing the iterative simulations in the grand canonical ensemble.

This research was primarily supported by the National Science Foundation through the Center for Dynamics and Control of Materials: an NSF MRSEC under Cooperative Agreement No. DMR-1720595. This work was also partially supported by the Welch Foundation (F-1696). We acknowledge the Texas Advanced Computing Center (TACC) at The University of Texas at Austin for providing HPC resources. B.A.L. acknowledges support from the Darleane Christian Hoffman Distinguished Postdoctoral Fellowship at Los Alamos National Laboratory.

1. Relative entropy update equation in the NPT ensemble

The update scheme for RE optimization within the NPT ensemble follows analogous mathematical steps to the procedure outlined in prior work for the NVT ensemble.9 As before, we assume a target ensemble of M statistically independent and identically distributed (IID) configurations, R1:M, all of which have fixed volume Vtgt. The target ensemble should reflect the desired structural correlations to be realized via the optimized interaction. Adopting the NPT ensemble as the simulation protocol, the probability of observing a configuration (R) with a volume V is given by the following Boltzmann factor:

(A1)

where U(R|θ) is the potential energy for configuration R, θ is a vector of the tunable values that parameterize the potential, β is the inverse thermal energy, P is the pressure, and Z(θ, P) is the partition function. Similarly, the probability of observing the set of IID target configurations, R1:M, with volume Vtgt, is given by the product of probabilities with the form of Eq. (A1),

(A2)

This quantity measures how likely a simulation in the NPT ensemble with parameters θ is to sample the configurations R1:M at the target volume. Maximization of this quantity is called maximum likelihood fitting.69 

As in the NVT derivation, it is easier to maximize the log-likelihood. Taking the natural log of Eq. (A2) and dividing by M yields

(A3)

which can be written as

(A4)

in the large configuration limit (i.e., M → ∞), where ptgt(R|Vtgt) is the probability distribution of the target simulation. This expression only differs explicitly from the NVT version through the additional pressure-volume term.

Within the gradient ascent optimization algorithm, updates to θ follow from

(A5)

where η is the learning rate that is empirically tuned to maintain a stable optimization. Application of the gradient operator to Eq. (A4) yields

(A6)

which is formally equivalent to the NVT result, as the pressure-volume contribution in Eq. (A4) is a constant with respect to θ.

Differences due to the RE optimization ensemble are only implicitly contained within the average that is performed over the simulation data [the right-most term in Eq. (A6)]. As implied by Ref. 31, the general relative entropy update in any other ensemble is given by the following replacements: ptgt(R|Vtgt) → ptgt(R|Etgt) and p(R, V |θ, P) → p(R, E|θ, I), where E are the fluctuating extensive variables coupled to the fixed intensive variables I and Etgt are the fixed values of the extensive variables in the target ensemble.

Equation (A6) can be greatly simplified in the case of isotropic pair interactions. For a one-component system with particles interacting via the pair potential u(r|θ), the total potential energy is U(R|θ)12ijNu(ri,j|θ). Substituting this into Eq. (A4) and integrating over all of the coordinates except those in the pair potential yield

(A7)

where ρ and ρtgt are the ensemble-averaged optimized and target densities, respectively, and g(r|θ, P) and gtgt(r|Vtgt) are the corresponding radial distribution functions. Importantly, in writing Eq. (A7), we have assumed that the pair interaction is finite in range and the system is macroscopic in size. This allows us to extend the integrals—which are technically over different volumes—arising from both terms in Eq. (A6) to infinity and add them together. This expression is similar to the NVT analog apart from the different values of density attached to each radial distribution function. The factors in front of the integral are constants and can be absorbed into the learning rate of Eq. (A5).

2. Additional details for relative entropy optimization

We use the evolution of the following convergence criterion:

(A8)

which is an integrated measure of the similarity of the current step i to the target ensemble, to govern the learning rate η. [See Eq. (A7) for definitions of the symbols.] In particular, we compute δ(i) = (γ(i−1)γ(i))/γ(i−1). If γ significantly decreases from one step to the next (i.e., δ(i) < −0.1), then the optimization is converging as desired. Therefore, η remains the same (η(i+1) = η(i)). If the convergence criterion increases, we have empirically observed that the optimization might be entering an unstable region in parameter space and ρ may begin to fluctuate. If the increase in γ is modest (i.e., 0.1 < δ(i) < 0.25), then we maintain the same value for η as above. However, if δ(i) is larger, then the learning rate is decreased (η(i+1) = 0.5η(i) if 0.25 ≤ δ(i) < 2.0 and η(i+1) = 0.25η(i) if δ(i) ≥ 2.0) to stabilize the optimization. Finally, if γ changes little (−0.1 ≤ δ(i) ≥ 0.1), then it is possible that the optimization is not exploring parameter space efficiently, and so η is increased accordingly (η(i+1) = 2η(i)). We also include minimum and maximum values for η (0.001 ≤ η ≤ 0.16 for most cases in this work), and we initialized the optimization using the minimum value for η.

For the initial guess for the NVT optimization and the subsequent NPT calculation where the pressure is set to the measured pressure from the preceding NVT result (P*), we use a WCA potential [practically done by setting ϵ1 and ϵ2 equal to zero in Eq. (3)]. As can be seen from Eq. (3), the length scales α1 and α2 are nonlinear parameters and therefore require reasonable initial guesses for the optimization to converge. Although the exact functional form differs, α1 and α2 have the same qualitative interpretation as in related prior work. Therefore, we use initial guesses for the length scales from prior RE calculations (α1 = 1.46 and α2 = 2.55).7 Preliminary calculations demonstrated that the optimization did not appear sensitive to the choice of initial guess so long as the α values were reasonable. One other choice of initial guess was arrived at by fitting the functional form of Eqs. (2) and (3) to the IBI result shown in Fig. 1(a). This IBI-fit initial guess gave nearly identical optimization results to the WCA initial guess. As we altered the pressure in the NPT optimizations, we used results from optimizations at other nearby pressures for the initial guess. We performed optimizations at P = 0.125P*, 0.25P*, 0.5P*, 0.75P*, 1.0P*, 1.25P*, 1.5P*, 2.0P*, 3.0P*, and 4.0P*. As the pressure was lowered below P*, the initial guess was given by the closest optimization result at higher pressure, and as the pressure was raised above P*, the initial guess corresponded to the nearest optimization result performed at lower pressure. The quality of the optimizations at P = 0.75P*, 1.25P*, and 1.5P* was consistent with the results presented in the main text and only omitted for clarity in the figures.

The optimized parameters at each pressure and the density at which the validation simulations were performed are given in Table I.

1.
H.
Chao
and
R. A.
Riggleman
, “
Inverse design of grafted nanoparticles for targeted self-assembly
,”
Mol. Syst. Des. Eng.
3
,
214
222
(
2018
).
2.
S.
Torquato
, “
Inverse optimization techniques for targeted self-assembly
,”
Soft Matter
5
,
1157
1173
(
2009
).
3.
A.
Jain
,
J. A.
Bollinger
, and
T. M.
Truskett
, “
Inverse methods for material design
,”
AIChE J.
60
,
2732
2740
(
2014
).
4.
E.
Marcotte
,
F. H.
Stillinger
, and
S.
Torquato
, “
Optimized monotonic convex pair potentials stabilize low-coordinated crystals
,”
Soft Matter
7
,
2332
2335
(
2011
).
5.
A.
Jain
,
J. R.
Errington
, and
T. M.
Truskett
, “
Inverse design of simple pairwise interactions with low coordinated 3D lattice ground states
,”
Soft Matter
9
,
3866
3870
(
2013
).
6.
R. B.
Jadrich
,
J. A.
Bollinger
,
B. A.
Lindquist
, and
T. M.
Truskett
, “
Equilibrium cluster fluids: Pair interactions via inverse design
,”
Soft Matter
11
,
9342
9354
(
2015
).
7.
B. A.
Lindquist
,
R. B.
Jadrich
, and
T. M.
Truskett
, “
Assembly of nothing: Equilibrium fluids with designed structured porosity
,”
Soft Matter
12
,
2663
2667
(
2016
).
8.
B. A.
Lindquist
,
S.
Dutta
,
R. B.
Jadrich
,
D. J.
Milliron
, and
T. M.
Truskett
, “
Interactions and design rules for assembly of porous colloidal mesophases
,”
Soft Matter
13
,
1335
1343
(
2017
).
9.
B. A.
Lindquist
,
R. B.
Jadrich
, and
T. M.
Truskett
, “
Communication: Inverse design for self-assembly via on-the-fly optimization
,”
J. Chem. Phys.
145
,
111101
(
2016
).
10.
R. B.
Jadrich
,
B. A.
Lindquist
, and
T. M.
Truskett
, “
Probabilistic inverse design for self-assembling materials
,”
J. Chem. Phys.
146
,
184103
(
2017
).
11.
W.
Piñeros
,
B. A.
Lindquist
,
R. B.
Jadrich
, and
T. M.
Truskett
, “
Inverse design of multicomponent assemblies
,”
J. Chem. Phys.
148
,
104509
(
2018
).
12.
B. A.
Lindquist
,
R. B.
Jadrich
,
W. D.
Piñeros
, and
T. M.
Truskett
, “
Inverse design of self-assembling Frank-Kasper phases and insights into emergent quasicrystals
,”
J. Phys. Chem. B
122
,
5547
5556
(
2018
).
13.
C. S.
Adorf
,
J.
Antonaglia
,
J.
Dshemuchadse
, and
S. C.
Glotzer
, “
Inverse design of simple pair potentials for the self-assembly of complex structures
,”
J. Chem. Phys.
149
,
204102
(
2018
).
14.
D.
Banerjee
,
B. A.
Lindquist
,
R. B.
Jadrich
, and
T. M.
Truskett
, “
Assembly of particle strings via isotropic potentials
,”
J. Chem. Phys.
150
,
124903
(
2019
).
15.
D.
Chen
,
G.
Zhang
, and
S.
Torquato
, “
Inverse design of colloidal crystals via optimized patchy interactions
,”
J. Phys. Chem. B
122
,
8462
8468
(
2018
).
16.
W. L.
Miller
and
A.
Cacciuto
, “
Exploiting classical nucleation theory for reverse self-assembly
,”
J. Chem. Phys.
133
,
234108
(
2010
).
17.
E.
Jankowski
and
S. C.
Glotzer
, “
Screening and designing patchy particles for optimized self-assembly propensity through assembly pathway engineering
,”
Soft Matter
8
,
2852
2859
(
2012
).
18.
A. W.
Long
and
A. L.
Ferguson
, “
Nonlinear machine learning of patchy colloid self-assembly pathways and mechanisms
,”
J. Phys. Chem. B
118
,
4228
4244
(
2014
).
19.
A. W.
Long
and
A. L.
Ferguson
, “
Rational design of patchy colloids via landscape engineering
,”
Mol. Syst. Des. Eng.
3
,
49
65
(
2018
).
20.
A. F.
Hannon
,
Y.
Ding
,
W.
Bai
,
C. A.
Ross
, and
A.
Alexander-Katz
, “
Optimizing topographical templates for directed self-assembly of block copolymers via inverse design simulations
,”
Nano Lett.
14
,
318
325
(
2014
).
21.
F. A.
Escobedo
, “
Optimizing the formation of colloidal compounds with components of different shapes
,”
J. Chem. Phys.
147
,
214501
(
2017
).
22.
J.
Madge
and
M. A.
Miller
, “
Optimising minimal building blocks for addressable self-assembly
,”
Soft Matter
13
,
7780
7792
(
2017
).
23.
M. R.
Khadilkar
,
S.
Paradiso
,
K. T.
Delaney
, and
G. H.
Fredrickson
, “
Inverse design of bulk morphologies in multiblock polymers using particle swarm optimization
,”
Macromolecules
50
,
6702
6709
(
2017
).
24.
M. Z.
Miskin
,
G.
Khaira
,
J. J.
de Pablo
, and
H. M.
Jaeger
, “
Turning statistical physics models into materials design engines
,”
Proc. Natl. Acad. Sci. U. S. A.
113
,
34
39
(
2016
).
25.
A. L.
Ferguson
, “
Machine learning and data science in soft materials engineering
,”
J. Phys.: Condens. Matter
30
,
043002
(
2017
).
26.
H. M.
Jaeger
and
J. J.
de Pablo
, “
Perspective: Evolutionary design of granular media and block copolymer patterns
,”
APL Mater.
4
,
053209
(
2016
).
27.
K. R.
Gadelrab
,
A. F.
Hannon
,
C. A.
Ross
, and
A.
Alexander-Katz
, “
Inverting the design path for self-assembled block copolymers
,”
Mol. Syst. Des. Eng.
2
,
539
548
(
2017
).
28.
W. G.
Noid
, “
Perspective: Coarse-grained models for biomolecular systems
,”
J. Chem. Phys.
139
,
090901
(
2013
).
29.
J. F.
Rudzinski
and
W. G.
Noid
, “
Coarse-graining entropy, forces, and structures
,”
J. Chem. Phys.
135
,
214101
(
2011
).
30.
F.
Müller-Plathe
, “
Coarse-graining in polymer simulation: From the atomistic to the mesoscopic scale and back
,”
ChemPhysChem
3
,
754
769
(
2002
).
31.
A.
Chaimovich
and
M. S.
Shell
, “
Coarse-graining errors and numerical optimization using a relative entropy framework
,”
J. Chem. Phys.
134
,
094112
(
2011
).
32.
M. S.
Shell
, “
The relative entropy is fundamental to multiscale and inverse thermodynamic problems
,”
J. Chem. Phys.
129
,
144108
(
2008
).
33.
C.-C.
Fu
,
P. M.
Kulkarni
,
M. S.
Shell
, and
L.
Gary Leal
, “
A test of systematic coarse-graining of molecular dynamics simulations: Thermodynamic properties
,”
J. Chem. Phys.
137
,
164106
(
2012
).
34.
D.
Reith
,
M.
Pütz
, and
F.
Müller-Plathe
, “
Deriving effective mesoscale potentials from atomistic simulations
,”
J. Comput. Chem.
24
,
1624
1636
(
2003
).
35.
H.
Wang
,
C.
Junghans
, and
K.
Kremer
, “
Comparative atomistic and coarse-grained study of water: What do we lose by coarse-graining?
Eur. Phys. J. E
28
,
221
229
(
2009
).
36.
A.
Das
and
H. C.
Andersen
, “
The multiscale coarse-graining method. V. Isothermal-isobaric ensemble
,”
J. Chem. Phys.
132
,
164106
(
2010
).
37.
N. J. H.
Dunn
and
W. G.
Noid
, “
Bottom-up coarse-grained models that accurately describe the structure, pressure, and compressibility of molecular liquids
,”
J. Chem. Phys.
143
,
243148
(
2015
).
38.
K. M.
Lebold
and
W. G.
Noid
, “
Systematic study of temperature and density variations in effective potentials for coarse-grained models of molecular liquids
,”
J. Chem. Phys.
150
,
014104
(
2019
).
39.
J.
Mewis
and
N. J.
Wagner
,
Colloidal Suspension Rheology
(
Cambridge University Press
,
2013
).
40.
A. J.
Archer
and
N. B.
Wilding
, “
Phase behavior of a fluid with competing attractive and repulsive interactions
,”
Phys. Rev. E
76
,
031501
(
2007
).
41.
Y.
Zhuang
and
P.
Charbonneau
, “
Equilibrium phase behavior of the square-well linear microphase-forming model
,”
J. Phys. Chem. B
120
,
6178
6188
(
2016
).
42.
Y.
Zhuang
and
P.
Charbonneau
, “
Recent advances in the theory and simulation of model colloidal microphase formers
,”
J. Phys. Chem. B
120
,
7775
7782
(
2016
).
43.
E.
Mani
,
W.
Lechner
,
W. K.
Kegel
, and
P. G.
Bolhuis
, “
Equilibrium and non-equilibrium cluster phases in colloids with competing interactions
,”
Soft Matter
10
,
4479
4486
(
2014
).
44.
M. B.
Sweatman
,
R.
Fartaria
, and
L.
Lue
, “
Cluster formation in fluids with competing short-range and long-range interactions
,”
J. Chem. Phys.
140
,
124508
(
2014
).
45.
S.
Mossa
,
F.
Sciortino
,
P.
Tartaglia
, and
E.
Zaccarelli
, “
Ground-state clusters for short-range attractive and long-range repulsive potentials
,”
Langmuir
20
,
10756
10763
(
2004
).
46.
F.
Sciortino
,
S.
Mossa
,
E.
Zaccarelli
, and
P.
Tartaglia
, “
Equilibrium cluster phases and low-density arrested disordered states: The role of short-range attraction and long-range repulsion
,”
Phys. Rev. Lett.
93
,
055701
(
2004
).
47.
J.
Groenewold
and
W. K.
Kegel
, “
Colloidal cluster phases, gelation and nuclear matter
,”
J. Phys.: Condens. Matter
16
,
S4877
S4886
(
2004
).
48.
G. F.
Kendrick
,
T. J.
Sluckin
, and
M. J.
Grimson
, “
Macrocrystal phases in charged colloidal suspensions
,”
Europhys. Lett.
6
,
567
572
(
1988
).
49.
J.-M.
Bomont
,
J.-L.
Bretonnet
,
D.
Costa
, and
J.-P.
Hansen
, “
Communication: Thermodynamic signatures of cluster formation in fluids with competing interactions
,”
J. Chem. Phys.
137
,
011101
(
2012
).
50.
A. J.
Archer
,
C.
Ionescu
,
D.
Pini
, and
L.
Reatto
, “
Theory for the phase behaviour of a colloidal fluid with competing interactions
,”
J. Phys.: Condens. Matter
20
,
415106
(
2008
).
51.
D.
Pini
,
G.
Jialin
,
A.
Parola
, and
L.
Reatto
, “
Enhanced density fluctuations in fluid systems with competing interactions
,”
Chem. Phys. Lett.
327
,
209
215
(
2000
).
52.
R. P.
Sear
and
W. M.
Gelbart
, “
Microphase separation versus the vapor-liquid transition in systems of spherical particles
,”
J. Chem. Phys.
110
,
4582
4588
(
1999
).
53.
R. B.
Jadrich
,
J. A.
Bollinger
,
K. P.
Johnston
, and
T. M.
Truskett
, “
Origin and detection of microstructural clustering in fluids with spatial-range competitive interactions
,”
Phys. Rev. E
91
,
042312
(
2015
).
54.
A. I.
Campbell
,
V. J.
Anderson
,
J. S.
van Duijneveldt
, and
P.
Bartlett
, “
Dynamical arrest in attractive colloids: The effect of long-range repulsion
,”
Phys. Rev. Lett.
94
,
208301
(
2005
).
55.
H.
Shin
,
G. M.
Grason
, and
C. D.
Santangelo
, “
Mesophases of soft-sphere aggregates
,”
Soft Matter
5
,
3629
3638
(
2009
).
56.
B.
Mladek
,
D.
Gottwald
,
G.
Kahl
,
M.
Neumann
, and
C.
Likos
, “
Formation of polymorphic cluster phases for a class of models of purely repulsive soft spheres
,”
Phys. Rev. Lett.
96
,
045701
(
2006
).
57.
P.
Camp
, “
Structure and phase behavior of a two-dimensional system with core-softened and long-range repulsive interactions
,”
Phys. Rev. E
68
,
061506
(
2003
).
58.
R.
Sear
and
D.
Frenkel
, “
Continuous freezing in three dimensions
,”
Phys. Rev. Lett.
90
,
195701
(
2003
).
59.
C.
Likos
,
N.
Hoffmann
,
H.
Lowen
, and
A.
Louis
, “
Exotic fluids and crystals of soft polymeric colloids
,”
J. Phys.: Condens. Matter
14
,
7681
7698
(
2002
).
60.
C. N.
Likos
,
B. M.
Mladek
,
D.
Gottwald
, and
G.
Kahl
, “
Why do ultrasoft repulsive particles cluster and crystallize? Analytical results from density-functional theory
,”
J. Chem. Phys.
126
,
224502
(
2007
).
61.
E.
Lascaris
,
G.
Malescio
,
S. V.
Buldyrev
, and
H. E.
Stanley
, “
Cluster formation, waterlike anomalies, and re-entrant melting for a family of bounded repulsive interaction potentials
,”
Phys. Rev. E
81
,
031201
(
2010
).
62.
K.
Zhang
,
P.
Charbonneau
, and
B. M.
Mladek
, “
Reentrant and isostructural transitions in a cluster-crystal former
,”
Phys. Rev. Lett.
105
,
245701
(
2010
).
63.
W.
Klein
,
H.
Gould
,
R. A.
Ramos
,
I.
Clejan
, and
A. I.
Mel’cuk
, “
Repulsive potentials, clumps and the metastable glass phase
,”
Physica A
205
,
738
746
(
1994
).
64.
P. D.
Godfrin
,
N. E.
Valadez-Pérez
,
R.
Castañeda-Priego
,
N. J.
Wagner
, and
Y.
Liu
, “
Generalized phase behavior of cluster formation in colloidal dispersions with competing interactions
,”
Soft Matter
10
,
5061
5071
(
2014
).
65.
J. C. F.
Toledano
,
F.
Sciortino
, and
E.
Zaccarelli
, “
Colloidal systems with competing interactions: From an arrested repulsive cluster phase to a gel
,”
Soft Matter
5
,
2390
2398
(
2009
).
66.
A.
Stradner
,
H.
Sedgwick
,
F.
Cardinaux
,
W. C. K.
Poon
,
S. U.
Egelhaaf
, and
P.
Schurtenberger
, “
Equilibrium cluster formation in concentrated protein solutions and colloids
,”
Nature
432
,
492
495
(
2004
).
67.
E. J.
Yearley
,
P. D.
Godfrin
,
T.
Perevozchikova
,
H.
Zhang
,
P.
Falus
,
L.
Porcar
,
M.
Nagao
,
J. E.
Curtis
,
P.
Gawande
,
R.
Taing
,
I. E.
Zarraga
,
N. J.
Wagner
, and
Y.
Liu
, “
Observation of small cluster formation in concentrated monoclonal antibody solutions and its implications to solution viscosity
,”
Biophys. J.
106
,
1763
1770
(
2014
).
68.
S.
Xu
,
C.
Sun
,
J.
Guo
,
K.
Xu
, and
C.
Wang
, “
Biopolymer-directed synthesis of high-surface-area magnetite colloidal nanocrystal clusters for dual drug delivery in prostate cancer
,”
J. Mater. Chem.
22
,
19067
19075
(
2012
).
69.
D.
Barber
,
Bayesian Reasoning and Machine Learning
(
Cambridge University Press
,
2012
).
70.
H.
Steinhaus
, “
Sur la division des corps matériels en parties
,”
Bull. Acad. Pol. Sci.
IV
,
801
804
(
1956
).
71.
S. P.
Lloyd
, “
Least squares quantization in PCM
,”
IEEE Trans. Inf. Theory
28
,
129
137
(
1982
).
72.
J.
MacQueen
, “
Some methods for classification and analysis of multivariate observations
,” in
Proceedings of the Fifth Berkeley Symposium on Mathematical Statistics and Probability, Volume 1: Statistics
(
University of California Press
,
Berkeley, CA
,
1967
), pp.
281
297
.
73.
D.
Arthur
and
S.
Vassilvitskii
, “
K-means++: The advantages of careful seeding
,” in
SODA’07: Proceedings of the Eighteenth Annual ACM-SIAM Symposium on Discrete Algorithms
(
Society for Industrial and Applied Mathematics
,
Philadelphia, PA, USA
,
2007
), pp.
1027
1035
.
74.
P. M.
Larsen
,
S.
Schmidt
, and
J.
Schiøtz
, “
Robust structural identification via polyhedral template matching
,”
Modell. Simul. Mater. Sci. Eng.
24
,
055007
(
2016
).
75.
A.
Stukowski
, “
Visualization and analysis of atomistic simulation data with OVITO–the open visualization tool
,”
Modell. Simul. Mater. Sci. Eng.
18
,
015012
(
2009
).
76.
D.
Van Der Spoel
,
E.
Lindahl
,
B.
Hess
,
G.
Groenhof
,
A. E.
Mark
, and
H. J. C.
Berendsen
, “
GROMACS: Fast, flexible, and free
,”
J. Comput. Chem.
26
,
1701
1718
(
2005
).
77.
S.
Plimpton
, “
Fast parallel algorithms for short-range molecular dynamics
,”
J. Comput. Phys.
117
,
1
19
(
1995
).
78.
G.
Bussi
,
D.
Donadio
, and
M.
Parrinello
, “
Canonical sampling through velocity rescaling
,”
J. Chem. Phys.
126
,
014101
(
2007
).
79.
S.
Rogge
,
L.
Vanduyfhuys
,
A.
Ghysels
,
M.
Waroquier
,
T.
Verstraelen
,
G.
Maurin
, and
V.
Van Speybroeck
, “
A comparison of barostats for the mechanical characterization of metalorganic frameworks
,”
J. Chem. Theory Comput.
11
,
5583
5597
(
2015
).
80.
M. P.
Allen
and
D. J.
Tildesley
,
Computer Simulation of Liquids
, 2nd ed. (
Oxford University Press, Inc.
,
New York, NY, USA
,
2017
).