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.
I. INTRODUCTION
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.
II. COMPUTATIONAL METHODS
A. Relative entropy optimization in other ensembles
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 (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
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 , 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.
B. Cluster target phase
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.
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.
C. Functional form for the relative entropy optimization
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,
with σ being the nominal particle diameter. The potential also had an auxiliary interaction that furnished an attractive well and a repulsive barrier,
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).
D. Simulation details
All simulations were performed in gromacs, version 5.1.2.76,77 The time step was , 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
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.
III. RESULTS AND DISCUSSION
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).
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.)
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).
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.
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.
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.
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.
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 |
IV. CONCLUSIONS
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.
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.
ACKNOWLEDGMENTS
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.
APPENDIX: RELATIVE ENTROPY OPTIMIZATION IN THE NPT ENSEMBLE
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:
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),
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
which can be written as
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
where η is the learning rate that is empirically tuned to maintain a stable optimization. Application of the gradient operator to Eq. (A4) yields
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 . Substituting this into Eq. (A4) and integrating over all of the coordinates except those in the pair potential yield
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:
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.