Breadth versus depth: Interactions that stabilize particle assemblies to changes in density or temperature Free

The synthesis of matter with precise and well-characterized structures at nanometer length scales is critical to the discovery and manufacture of new material systems with desirable optical,¹ mechanical,² and other physical³ properties. Unfortunately, despite recent progress,^4–6 direct fabrication of such materials using top-down approaches can be technologically challenging, expensive, and time-consuming. An alternative bottom-up strategy is to synthesize systems of nanoscale particles with effective interactions^7,8 that favor spontaneous or directed self-assembly into the targeted structure (e.g., a periodic crystal or superlattice configuration).^9–12 

Historically, design in self-assembly has focused on forward approaches, whereby the interactions between particles are altered in some rationally guided or Edisonian fashion, and the resulting equilibrium structures that such interactions produce are identified and catalogued according to their ability to meet various design goals concerning their material properties. However, attention has begun to shift toward inverse design strategies, whereby interactions that favor the targeted structure or properties in the thermodynamically stable state are formally discovered through a constrained, statistical mechanical optimization.^13,14

Several recently introduced computational methods for inverse design focus on optimization of interparticle interactions to stabilize a targeted ground-state configuration with the assumption of an isotropic pair potential ϕ(r, {α_i}) with variable parameter set {α_i}. Such approaches have found various interactions that stabilize two-dimensional square, honeycomb and kagome lattices^15–18 as well as the three-dimensional diamond crystal structure.^19,20 Furthermore, it was demonstrated that particles with the optimized interactions indeed assembled into the targeted lattice phases at a higher temperature using molecular simulations.^15,16,20

In this same vein, we have used inverse methods to design convex-repulsive, isotropic pair potentials of the following form:^21,22

that stabilize open crystal structures. The final optimized potentials qualitatively resembled effective pair interactions observed in various soft colloidal systems (e.g., star polymers, ligand-passivated nanocrystals, and microgels).⁷ Here, σ and ϵ represent characteristic length and energy scales, respectively; H is the Heaviside function; {A, n, λ_i, k_i, δ_i} are variable parameters (i.e., α_i), one of which is fixed to ensure ϕ(1) = ϵ; r_cut is a cut-off radius; and f_shift is a quadratic function f_shift(r/σ) = P(r/σ)² + Qr/σ + R added to enforce ϕ(r_cut/σ) = ϕ′(r_cut/σ) = ϕ″(r_cut/σ) = 0. Using a simulated annealing optimization approach, parameters for this potential that stabilize, over a very wide range of densities, square and honeycomb lattices in two dimensions²² as well as simple cubic and diamond ground-state structures in three dimensions were determined.²¹ Complete phase diagrams were also calculated for the three-dimensional systems,²³ which illustrated rich and complex phase equilibria with the targeted assemblies exhibiting good thermal stability relative to other competing crystal lattices.

One basic question that has not yet been addressed is which features of a pair potential would tend to maximize the melting temperature of a given target structure? Moreover, how would encoding high thermal stability into the interactions affect the corresponding range of densities for which the target structure is favored? In other words, is there a natural compromise between designing for robustness to changes in temperature versus volume? Such questions are challenging to answer directly via inverse methods because they would require incorporating full molecular simulations (for a wide range of model parameters and thermodynamic conditions) into the optimization problem, which is computationally unfeasible at present. A pragmatic alternative is to search for features of the ground-state behavior that, while easier to compute than higher temperature properties, correlate with thermal stability of the target phase. In the present work, we find that placing constraints on the minimum chemical potential advantage that the target structure would exhibit over selected equi-pressure competing lattices at zero temperature helps determine optimized interactions with higher target-phase melting temperatures.

The specific structure that we target via ground-state inverse optimizations in this work is the two-dimensional square lattice ground state, which has attracted considerable theoretical interest in recent years,^22,24–26 and the class of pair potentials we consider are those described by Eq. (1). The stable equilibrium ground-state structure can be established by determining the global minimum of the potential energy U at a fixed density and zero temperature or the minimum of the chemical potential μ at a fixed pressure and zero temperature (amongst other possibilities that follow from classical thermodynamics²⁷). Following other works,^21,22,28,29 we adopt the latter fixed pressure framework for our optimizations for convenience because any coexistence between the target structure and another lattice also requires equality of pressure between phases. Through our optimizations, we test how the maximum achievable range of densities for stability of the square-lattice ground state is affected by constraining the differences between its chemical potential at a prescribed state point and those of selected competing lattice structures at the same pressure. To do this systematically requires the solution of a series of optimization problems, each utilizing different constraints. Given the considerable computational expense of using stochastic optimizers (e.g., simulated annealing and genetic algorithms) for even a single optimization, we instead formulate each optimization problem of interest as a constrained mixed-integer nonlinear problem, and we solve it numerically using the General Algebraic Modeling System (GAMS).³⁰ We then explore the consequences of our imposed chemical potential advantage of the target ground state for the resulting thermal stability (i.e., melting temperature) of the resulting lattices.

The paper is organized as follows. We first introduce the computational methods used for carrying out the pair-potential optimization and the melting point estimation. Next, we show the relationship between the density range of ground-state stability for the target lattice and the minimum chemical potential difference between the target and selected competitors. We conclude by discussing how these trends manifest in the resulting pair potentials and lattice melting temperatures.

In this section, we describe how we formulate and solve the inverse design problem of interest in this work: finding isotropic pair interactions that maximize the density range over which the targeted square lattice is the ground-state configuration given a constraint on its chemical potential advantage over selected competing lattices. We further detail the implementation of molecular dynamics (MD) simulation methods for estimating the melting point to characterize the corresponding thermal stability of the designed lattice structures.

We formulate the design optimization problem following the general paradigm

where f(x) is an objective function and g_i(x) are constraint equality or inequality equations with variables x. The mathematical forms of f(x) and g_i(x) define the type of problem to be solved (e.g., if integer variables or non-linear functions are necessary). For the inverse design calculations of interest here, the set of equations g_i incorporate any desired constraints to be imposed on the interparticle pair potential and f is formulated to ensure optimization of the desired thermodynamic property. To optimize for specific ground states, one needs to consider not only the target lattice l_t of the design but also other lattices {l} that naturally compete with it for thermodynamic stability (the procedure to determine the pool of competing lattices is discussed separately in Sec. II A 3). Using Eq. (1) as the model pair potential, we ultimately seek potential parameters {A, n, l_i, k_i, d_i} (i.e., the decision variables) that maximize the density range for which the target lattice l_t has a chemical potential lower than that of lattices in {l} at the same pressure such that a minimum specified chemical potential advantage of the target structure over select competitors is obtained at a given state point.

Specifically, to incorporate the pair potential of Eq. (1), we introduce constraint equations that ensure that the potential is appropriately normalized, repulsive, convex, and continuous (we implicitly nondimensionalize energies by ϵ, lengths by σ, and omit parameter notation below for brevity). The normality condition is given by

and the other constraints are given by

and

We set r_cut = 2.271 83 as motivated by previous work considering square lattices designed via this potential form.²² As a practical matter, r is discretized over a finite set of points in (0, r_cut); we use ten uniform points in ϕ(r) and 60 points distributed in a 1:6:5 ratio from ranges [0.2,0.8), [0.8,1.2], (1.2,r_cut) for ϕ″(r), which we find sufficient to enforce the constraints. It is not necessary to include the constraint on ϕ′(r) so long as the other constraints are fulfilled.

Next, we specify the equations describing the physics of the ground state. The first is for the internal energy per particle, which can be expressed as

Here, r_i,l(ρ_l) are the density-dependent coordination distances for each lattice l and n_i,l are the number of neighbors at those distances.³¹ The pressure P_l of lattice l is related to its density ρ_l by the virial expression,

For our purposes, the relevant density of a competing lattice l, ρ_l, is that which leads to equality of pressure with the target lattice l_t of density ρ_t. In other words, ρ_l(ρ_t) can be determined from knowledge of ρ_t via the relation P_l(ρ_l) = P_t(ρ_t), and thus from Eq. (6), we have

The chemical potential of a ground-state lattice l is, in turn, given by

Finally, an auxiliary equation is used,

where r₀(ρ_l) represents the nearest neighbor distance for competing lattices at density ρ_l. This helps tighten the optimization formulation by keeping density within a reasonable range.

An objective function f(x) that fulfills our optimization goals must also be specified. We choose such an objective function to evaluate to a finite scalar value f(x) → f and to be conducive to maximizing the range of densities Δρ_t = ρ_t,f − ρ_t,i for which the target lattice exhibits a chemical potential lower than that of the competing lattices. This is then defined as

where the sum is over discretized target lattice densities (each spaced a distance δ apart), ρ_l(ρ_t) is computed from Eq. (7), and H is again the Heaviside step function.

An additional constraint equation,

is introduced to specify the minimum acceptable chemical potential difference Δ between the target lattice at an intermediate density point ρ_t,o and selected competing lattices at the same pressure. Here, we use ρ_t,o = 1.39, which was found in an earlier study²² to be in the middle of the density range of stability for a square lattice designed for Δ = 0 and the same pair potential form. In practice, we have found that the post-optimization chemical potential difference between the target lattice and its closest selected competitor Δμ ≡ min{μ_l(ρ_l(ρ_t,o))} − μ_t(ρ_t,o) is approximately equal to the constraint Δ in all cases.

We implemented the optimization problem described above in GAMS,³⁰ using the Basic Open-source Nonlinear Mixed INteger (BONMIN)^32,33 solver with Interior Point OPTimizer (IPOPT)³⁴ as the non-linear sub-solver. This choice of solver permits us to use integer valued functions such as in Eq. (10) (i.e., the Heaviside function) as well as the remaining non-linear functions present in the potential and system physics formulation.

In practice, each optimization begins by inputting an initial guess for the pair potential parameter set that does not violate the constraints of Eqs. (2)-(4) and specifying a narrow target lattice density range [ρ_t,i, ρ_t,f] containing ρ_t,o to consider.

If the maximum attainable value of f is realized in the optimization (i.e., the maximum number of feasible density points Δρ_t/δ = (ρ_t,f − ρ_t,i)/δ is achieved), the boundaries of the density range are widened and the previously attained potential is used as the initial guess for a new optimization. This procedure is repeated until the optimization returns f < Δρ_t/δ, indicating that the maximum density range of stability for a given chemical potential constraint Δ was attained in the previous optimization. We carry out the optimizations described above for different values of Δ to explore how an imposed chemical potential advantage of the target lattice affects the maximum attainable Δρ_t. As discussed in Sec. III, there is a maximum value of Δ above which a feasible solution does not seem to exist for any density range. While found values are not verifiably global due to the local nature of the optimizer, they are optimal to the best of our efforts.

In addition to the explicit constraints described above, only pair potentials that result in mechanically stable target ground-state structures (as determined from phonon spectra analysis) were considered. Spontaneous assembly of particles interacting via the optimized potentials from the fluid state into the target structure upon temperature quenching was also verified at ρ_t,o using Monte Carlo simulations (see the supplementary material³⁵).

In our previous work on inverse design of targeted lattices,^21,22 we sought pair potentials that simply maximized the density range of stability of the desired structure (Δ = 0). For that type of optimization, it was necessary to choose a finite pool of competitive structures to compare with the target lattice, ideally those with the lowest values of chemical potential at the pressures of interest (which are not generally known in advance). We determined the composition of this competitive pool from an iterative procedure. An initial set of structures was selected (e.g., Bravais lattices plus a small number of non-Bravais lattices or tilings) based on intuition and knowledge obtained from earlier simulation studies on similar pair potentials. An optimization was then performed using the chosen competitive pool, followed by a forward calculation of the ground-state phase diagram with the optimized pair potential for densities in the targeted range. Any new structures that appeared were subsequently added to the previous competitive pool, and a second optimization with the updated list of competitive structures was completed. This process–updating the competitive pool and optimizing the pair potential considering the expanded list of possible lattices found in forward calculations–was repeated until no new competing lattice structures emerged.

In the present study, we repeat similar optimizations but with the added requirement of a minimum chemical potential difference between the target and selected competitors. The hypothesis is that such a constraint will find potentials displaying enhanced thermal stability of the targeted phase. Note that one cannot enforce a fixed chemical potential difference between the target lattice and all possible competitors. To understand why, consider a representation of lattice structure defined by a set of primitive and basis vectors {v}. If {v} can be modified continuously in some way (without adding or removing particles), e.g., by a set of suitable parameters {Θ}, then {v({Θ})} will then define a hyperspace of continuously connected lattices with the target structure l_t representing a specific point in this space. For ground-state systems of a given pair potential at a specific pressure, state quantities such as μ depend on the lattice structure (i.e., μ({v})) such that μ itself can be represented as a hypersurface of continuously connected lattices v({Θ}). Thus, one can always find structures in the neighborhood of the target lattice on the hypersurface with chemical potentials arbitrarily close to that of the target.

Considering this, it is clear that one cannot enforce a nonzero minimum chemical potential difference Δ between the target and all possible competitive structures. However, one can meaningfully constrain the μ hypersurface in the optimization by enforcing a minimum chemical potential difference Δ between the target and a chosen set of lattices {l_g} that define “flag points” on the landscape. This helps to achieve a standardized and well defined constraint depth that is feasible for the optimization. Indeed, a similar approach was introduced by Zhang et al. for a related Δμ optimization and justified under similar premises.²⁹ 

We provide an example for concreteness. The chosen target square lattice can be represented as a point in a larger Bravais subspace spanned by oblique (OBL) primitive vectors {v_B({Θ})} with {Θ} consisting of an aspect ratio b/a and primitive vector angle θ. Thus, the square lattice is represented by v(1, π/2), while other Bravais lattices like the triangular lattice are given by {v(1, π/3)} and so on. The corresponding μ({v}) landscape for this Bravais subset is then a function of (b/a, θ) (i.e., μ({v(b/a, θ)}). As such, the flag-point lattices we choose for enforcing the depth constraints here are the triangular lattice and a rectangular (REC) lattice which capture independent variations along the θ and b/a directions in the neighborhood of the target (for an illustration of the resulting Bravais chemical potential landscape from one of the optimized potentials, see Figure 1S of the supplementary material).³⁵ Similar subspace arguments can be made to account for elongated triangular (ET) and snub-square (SS) non-Bravais lattices, where the former can be transformed into square by a row-shift and the latter by a rotation of a single tile around its next neighboring square tiles. Other tilings or non-Bravais lattices can, in principle, be important for target lattice stability, but we did not find others that were relevant in the present square-lattice design problem.

Given the above considerations, our final competing pool consisted of lattices determined from the iterative forward procedure, some of which (those which naturally belonged to a subspace that continuously deformed into the square lattice) were also chosen as flag-point lattices. From the Bravais subspace, the final competing pool consisted of triangular, REC b/a = 1.17, and OBL b/a = 1.1, θ = 1.09 lattices, with triangular and REC also serving as flag points for the chemical potential constraint. Similarly, for the other relevant non-Bravais subspaces, the competitive pool included one SS lattice (b/a = 1.0) and three ET lattices (b/a = 1.07, b/a = 1.20, and b/a = 1.23), with all but the last ET lattice serving as flag points for the chemical potential constraint.

The Z-method is a microcanonical MD simulation strategy for estimating the melting point of a crystal that does not require free energy calculations.³⁶ The approach is based on the idea that a crystal remains metastable upon raising the temperature until it reaches its superheating limit, where it is hypothesized to have the same internal energy as the liquid at the equilibrium melting temperature T_m. The name is due to the fact that the estimate comes from the Z-shaped (zig-zag) graph that one observes for the system in the temperature T vs pressure P plane as it transforms from the solid phase upon raising the energy at a constant volume. It has been applied for a variety of model systems and has been repeatedly tested for both accuracy and variability.^37–39 For our purposes here, where we seek only estimates of melting temperatures to compare the widely varying thermal stabilities of targeted assemblies designed under various constraints, the Z-method provides an adequate guide.

Energy sweeps for the Z-method are carried out as follows. Initial particle positions are set in a perfect square lattice, and initial velocities are chosen from a random distribution and scaled to achieve a desired initial kinetic energy. For a series of progressively increasing energies, microcanonical MD simulation trajectories of N = 1024 particles (and a periodic square cell of length V^1/2 chosen to set ρ = N/V = 1.39) are initiated with a time step value of 0.001. After an initial pre-equilibration period at each energy, averages of static quantities like temperature T and pressure P are taken every 1000 time steps for at least 10⁶ steps. Near the transition region, averages for liquid and solid properties are taken separately with the phase being determined by the translational order parameter τ,

Here, r_i denotes particle positions’ vectors and k is a reciprocal lattice vector. We chose $k=2πl(1,1)$ for this purpose, where l denotes the lattice constant value at density ρ. We use τ ≥ 0.5 to indicate solid configurations and τ ≤ 0.1 to denote liquid configurations. These assignments were additionally supported by monitoring the mean square displacement of the particles as a function of time.

Reported estimates of T_m are averages of the temperature of the superheating limit of the solid obtained from twelve independent energy sweeps.

As further corroboration of the estimates obtained from the Z-method described above, we also carry out melting point estimations by the hysteresis method. This method is based on analysis of superheating and supercooling processes in the framework of nucleation theory and validated through molecular dynamics simulations.^40,41 The basic approach is to carry out a simple heating and cooling sweep of the system near the melting point to determine the temperature of superheating T₊ and supercooling T₋. The melting point T_m is then estimated from

As such, we carry out Monte Carlo simulations in the canonical ensemble for N = 400 particles in a periodic square cell of length V^1/2 adjusted to fix density at ρ = N/V = 1.39. Simulations are started from the crystal phase and heated until melting is achieved. The system is then cooled from the liquid back into the crystal. The T₊ and T₋ points are obtained from the resulting hysteresis loop in an energy vs temperature diagram.

Using our described ground state optimization procedure, we were able to obtain pair potential parameters for Eq. (1) that satisfied all of our objective goals. That is, we found potentials that (a) were convex repulsive, (b) maximized the density range Δρ_t for which the square lattice is the stable structure, and (c) were such that the target at density ρ_t,o displayed a specified minimum chemical potential advantage Δμ over the flag-point competitors (as elaborated in Sec. II). The resulting relationship between Δρ_t subject to increasing values of Δμ for the optimized potentials is plotted in Figure 1.

As seen, there is a clear negative correlation between Δμ and Δρ_t. While the exact values of Δρ_t may change based on the choice of non-linear subsolver (also given local nature of the solutions), test runs using a different subsolver showed values that yielded a very similar trend (not shown). In other words, there appears to be a clear compromise with this potential form between designing for high stability at a given density and designing for stability with respect to changes in density. There also appears to be a limit with this potential form to how stable one can make the square lattice ground state at ρ_t,o relative to the flag-point lattices (Δμ ≈ 0.23). For instance, we were only able to find solutions consistent with larger Δμ than those shown in Figure 1 if we allowed the pair potential to violate the convexity constraint.

In terms of judging the overall quality of the optimizations, we can compare to one result from a previous study,²² where a simulated annealing algorithm was used to find parameters for the potential of Eq. (1) that maximized the range of densities for which the square lattice was the stable ground-state structure (with no chemical potential constraint). In that paper, Δρ_t = 0.39 was found for the optimized potential, which displayed a minimum chemical potential advantage of Δμ ≈ 0.01 over the flag-point lattices considered here. This can be compared to that of the potential obtained in this study with a Δμ = 0.01 constraint, which exhibits a 50% wider density range, Δρ_t = 0.58. While reported solutions are not verifiably global, the fact that such a large improvement in the objective function was obtained points to one of the advantages that the present rigorous framework has over heuristic optimization approaches like simulated annealing (a global optimizer in principle).

We now explore how features of the optimized interparticle potentials help to explain the observed trade-off associated with designing for a large chemical potential advantage of the target ground-state structure at a given density versus designing for target stability over a wide range of densities. In Figure 2, the pair potentials corresponding to Δμ = [0.01-0.23] are shown (for the full list of potential parameters values see Tables S1 and S2 of the supplementary material).³⁵ While no pronounced features can be expected for strictly convex-repulsive interactions, two important aspects of the potential do manifest. As Δμ increases, so does the steepness of the core repulsion (for r ≲ 0.8) as well as the rate of radial decay towards the cut-off point (for r ≳ 1.2). The latter part can be seen more clearly in the log-log inset where intermediate features of core repulsion and radial decay can be seen to lie approximately between the two extrema potentials corresponding to Δμ = 0.01 and Δμ = 0.23. As we discuss next, this sharpening of radial-dependent features with increasing Δμ is what provides the chemical potential advantage of the target over its competitors, but at the cost of target lattice stability at other densities.

To look closer into the relation between pair potential form and target stability, it is helpful to recall that the chemical potential expression for a ground state system is given as μ_l = U_l + P_l/ρ_l. Using the energy and pressure expressions in (5) and (6), it is possible to recast this expression as

where ψ(r) has been defined as

As such, we see chemical potential depends not only on the pair potential but also on its gradient. Analyzing the radial dependence of ψ(r) will thus help to understand how the various lattice coordination shells at their respective radial separations contribute to the chemical potential and how they bias the functional form of the optimized potentials leading to the observed negative correlation between Δμ and Δρ_t.

We illustrate these points by plotting ψ(r) for optimized interactions corresponding to the limiting cases of strongly (Δμ = 0.23) and weakly (Δμ = 0.01) constrained chemical potential advantages of the square lattice ground state over the flag point structures. The plot in Figure 3(a) compares both ψ(r) and the fractional coordination-shell contributions to the chemical potential of the square lattice for the two potentials. As can be seen, interactions obtained with the larger Δμ constraint impart greater emphasis on first-shell contributions that translate into potentials with harder cores and faster decays at these distances. These ψ(r) features help the square lattice realize a lower chemical potential than that of the triangular lattice whose more densely packed first-coordination shell lies at a separation similar to that of the square lattice. Equally important is the shoulder-like region that decays between the square lattice’s first and second coordination shells. The role of this shoulder is to destabilize the closely competitive rectangular and elongated triangular lattices that have second coordination shells at separations within the shoulder region and thus contribute to their higher values of chemical potential compared to that of the square lattice (See Table S3 of the supplementary material for a list of μ_l,i/μ values of selected lattice competitors shells up to the third coordination).³⁵ 

The potential shape trends obtained from optimizations with the high Δμ constraint described above can be contrasted to the muted features that manifest when a smaller Δμ constraint is applied (which leads to considerably larger Δρ_t). Shown in Figure 3(b) is also ψ(r) for the two cases, but now presented along with shaded areas to indicate the range of the first- and second-coordination shell distances of the corresponding stable square lattices. The key point is that small changes in coordination distances (due to changes in density) would have very different consequences for the chemical potential of the Δμ = 0.23 system as compared to the Δμ = 0.01 system due to their different forms for ψ(r). For the Δμ = 0.23 system, small changes in density and coordination distances will produce pronounced changes in ψ(r) and hence the chemical potential. As a result, the specific shape that provided a great chemical potential advantage for the square lattice at ρ_t,o is no longer able to favor the structure at even modestly lower or higher densities. In contrast, the slower varying form of ψ(r) for the Δμ = 0.01 system, while providing a reduced chemical potential advantage at ρ_t,o, is able to keep the square lattice stable over a wider density range. The inverse relationship in Figure 1 emerges as a natural consequence of this trade-off.

Moving on to understand how designing potentials for large Δμ for the square lattice ground state at ρ_t,o = 1.39 affects the thermal stability of the target structure, we use the Z-method and the hysteresis method to estimate the corresponding melting temperatures T_m at that density. As can be seen in Figure 4, the potentials optimized with larger Δμ constraints also show higher T_m irrespective of the estimation method. For instance, while the Δμ = 0.01 system has a melting point at around T_m = 0.02, the melting point for the Δμ = 0.23 potential is approximately T_m = 0.2, an order of magnitude greater. A slightly more pronounced but largely similar result is found from the hysteresis method (T_m = 0.02 to T_m = 0.22 at Δμ = 0.01 to Δμ = 0.23, respectively). Considering the ψ(r) analysis presented above, this trend makes intuitive sense. Potentials designed with larger Δμ constraints impose greater penalties to target lattice deformation, and hence a higher average energy is required to move particles from their perfect lattice positions. This apparently translates directly to a higher melting point for the structure. A similar argument can be made based on the discussion of the μ hypersurface in Sec. II. Since Δμ captures an effective “well depth” for the target structure, imposing higher Δμ has the effect of creating greater “restoring forces” on the target (i.e., higher eigenvalues of the μ Hessian).²⁹ This results in increased mechanical stability at the ground state and a correspondingly higher melting point as shown here.

Finally, an important question arises when comparing to Figure 1. Since we probed the melting points along a path where both “range” (Δρ_t) and “depth” (Δμ) change simultaneously, how does T_m change if we hold a particular depth constant and vary the range or vice versa? From our discussion so far, we expect that depth alone will determine the thermal trend while the range will be largely inconsequential. Indeed, test runs where we probed systems with the same depth but different ranges yielded scatter around a mean value, whereas holding range constant and varying depth produced melting points consistent with Fig. 4 (not shown). Thus, for our inverse optimized pair potential, Δμ of the target in the ground state appears to strongly correlate with the thermal stability of the assembly, while the corresponding density range of stability has no such clear connection to the melting temperature.

We have used inverse methods of statistical mechanics to gain new insights into the trade-off between designing interactions for stability of a target structure and changes in temperature versus density. Specifically, we have explored the consequences of constraining the minimum chemical potential advantage of a target square lattice ground state at a prescribed density ρ_t,o over select competitors (Δμ or “depth” on the μ landscape) while designing potentials that maximize the range of density where the target ground state is stable (Δρ_t). The resulting constrained nonlinear optimization problem was solved numerically. For the isotropic, convex-repulsive pair interactions considered here, pair potentials designed with a larger Δμ constraint exhibited a narrower range of density stability Δρ_t. The reasons for this compromise are apparent when examining the radially dependent forms of the optimized pair potentials and their gradients. To enable high stability at a given density, features in the potential and its derivative must align with specific coordination shells to help produce the desired differences in chemical potential. When such features are present, however, the resulting target structures can lose stability with even modest changes in density.

We have also verified, via MD and Monte Carlo simulations, that potentials exhibiting ground states designed with larger Δμ constraints have higher melting temperatures at the target density. Preliminary tests further suggest that it is Δμ alone, and not Δρ_t, that correlates with T_m. Both results are in accord with the idea that Δμ constraints ensure restoring forces on the μ hypersurface that resist deformation (and ultimately melting) of the target structure.

T.M.T. acknowledges support of the Welch Foundation (F-1696) and the National Science Foundation (No. CBET-1403768). We also gratefully acknowledge many insightful and productive conversations with Dr. Avni Jain about the ground-state calculations in this work, as well as her help with comparing our results to those obtained from earlier simulated annealing approaches to inverse design.