Turning Statistical Physics Models Into Materials Design Engines

Despite the success statistical physics has enjoyed at predicting the properties of materials for given parameters, the inverse problem, identifying which material parameters produce given, desired properties, is only beginning to be addressed. Recently, several methods have emerged across disciplines that draw upon optimization and simulation to create computer programs that tailor material responses to specified behaviors. However, so far the methods developed either involve black-box techniques, in which the optimizer operates without explicit knowledge of the material's configuration space, or they require carefully tuned algorithms with applicability limited to a narrow subclass of materials. Here we introduce a formalism that can generate optimizers automatically by extending statistical mechanics into the realm of design. The strength of this new approach lies in its capability to transform statistical models that describe materials into optimizers to tailor them. By comparing against standard black-box optimization methods, we demonstrate how optimizers generated by this formalism can be faster and more effective, while remaining straightforward to implement. The scope of our approach includes new possibilities for solving a variety of complex optimization and design problems concerning materials both in and out of equilibrium.

Despite the success statistical physics has enjoyed at predicting the properties of materials for given parameters, the inverse problem, identifying which material parameters produce given, desired properties, is only beginning to be addressed. Recently, several methods have emerged across disciplines that draw upon optimization and simulation to create computer programs that tailor material responses to specified behaviors. However, so far the methods developed either involve black-box techniques, in which the optimizer operates without explicit knowledge of the material's configuration space, or require carefully tuned algorithms with applicability limited to a narrow subclass of materials. Here we introduce a formalism that can generate optimizers automatically by extending statistical mechanics into the realm of design. The strength of this approach lies in its capability to transform statistical models that describe materials into optimizers to tailor them. By comparing against standard black-box optimization methods, we demonstrate how optimizers generated by this formalism can be faster and more effective, while remaining straightforward to implement. The scope of our approach includes possibilities for solving a variety of complex optimization and design problems concerning materials both in and out of equilibrium. materials design | directed self-assembly | optimization | inverse design C omputer programs that can design material properties have led to exciting, new directions for materials science (1)(2)(3). Computational methods have been used to predict crystal (4) and protein (5,6) structures, yielding the toughest crystals known to mankind (4) and de novo protein configurations unseen in nature (5). Applied to polymers, Monte Carlo methods (7)(8)(9) and evolutionary algorithms (10,11) have paved the way toward optimizing directed self-assembly. Similar methods have been used to identify the crystal structures of patchy, colloidal particles (12). For far-from-equilibrium systems like jammed, metastable aggregates of particles (3), simulation-based optimization has been successfully used to design bulk properties like stiffness (13) and packing density (14) by way of tuning complicated microscale features like particle shape.
However, despite these successes, most of the existing methods work only for narrowly defined classes of materials: Optimization techniques that prove successful at designing one class of materials may struggle or fail on other systems. Thus, designing new materials can require a large investment in trial and error at the level of the algorithm itself, even if, for given parameters, the material's behavior can be simulated easily.
In black-box approaches, the algorithm tunes the material by adjusting control parameters without considering the likelihood of finding the material in microscale configurations. Instead, the optimizer operates in some auxiliary space, defined outside the physical model, and remains ignorant of the statistics in the physical configuration space. On the other hand, for the overwhelming majority of materials, an accurate description of macroscale behavior comes about by explicitly considering the probability of finding the system in certain microscopic configurations. Several materials optimization approaches exist that take this statistical nature into account, examples include optimizers that design spin configurations (15,16), patchy colloidal particles (17), and self-assembly driven by short-ranged interactions (18). These approaches build heavily upon the specific model defining the material. As a consequence, it is difficult to extend them beyond the material class they are designed for. Evidently, microscale configurations present key statistical information about a material, which is completely ignored by black-box approaches, yet there is no formalism that generically incorporates this information into materials design.
The question we ask is whether microstate information not only can be used to enhance an optimizer's speed and range of applicability, but also can become the cornerstone of an approach that automatically transforms a design goal and a statistical model into an optimizer. One can then construct optimizers that can work on material classes that are as broad as those described by statistical mechanics, without the need for ad hoc modification.
Here we take some first steps toward such a framework. We present a formalism that can be used to transform the capacity to predict material behavior into an optimizer that tunes it. Furthermore we find that our formalism often solves optimization problems faster and more reliably than approaches built around black-box methods.

Deriving the Optimization Equations
Our approach starts by assuming a given model ρðxjλ i Þ that predicts the probability of finding a material in some configuration, x. The model depends on the adjustable parameters λ i and by tuning λ i , a user can impact the emergent, bulk properties, averaged over configuration space. Thus, design proceeds by tweaking λ i to promote a desired, user-specified response.
We work toward this goal, starting from the heuristic equation Significance A fundamental tenet of science is that the properties of a material are intimately linked to the nature of the constituent components. Although there are powerful methods to predict such properties for given components, a key challenge for materials design is the inverse process: identifying the required components and their structural configuration for given target properties. This paper presents a new approach to this challenge. A formalism is introduced that generates algorithms for materials design both under equilibrium and under nonequilibrium conditions and operates without the need for user input beyond a design goal. This formalism is broadly applicable, fast, and robust, and it provides a powerful tool for materials optimization as well as discovery.
where the angle brackets denote an average over configurations weighted by ρ, the overdot denotes a derivative with respect to an artificial time that indexes the optimization steps, and f ðxÞ is a function that quantifies how well configuration x represents the user-specified goal. In this way, Eq. 1 attempts to increase the probability of finding the system in states with better than average values of f ðxÞ: If a configuration x has a value of f ðxÞ greater than the average, then the probability of finding the system near x increases, for average f ðxÞ it does not change, and for worse than average it decreases. As written, Eq. 1 assumes that ρðxjλ i Þ can be independently set for every possible point in the configuration space, despite the fact that ρðxjλ i Þ is constrained by the physics behind the material of interest. In actuality, ρðxjλ i Þ can offer only a limited flexibility through the parameters λ i . Thus, for many problems it will not be possible to exactly satisfy Eq. 1, although it is possible to make a best approximation to Eq. 1, given a particular physical distribution. We achieve this by setting the changes to λ i such that they minimize the average error between the updates implied by Eq. 1 and the actual changes to ρðxjλ i Þ. Explicitly, we rewrite Eq. 1 as _ λ i ∂ λi log½ρ = f ðxÞ − hf i and select the _ λ i that minimize the average value of the squared error, e = hð _ λ i ∂ λi log½ρ − ½f ðxÞ − hf iÞ 2 i. The equations of motion for λ i that minimize « at a given instant in time are found by setting the partial derivatives with respect to _ λ i equal to zero (Supporting Information): Eq. 2 is now a closed expression for λ i that depends only on expectation values. Thus, it can function as an algorithm: one can draw samples from ρðxjλ i Þ, use the samples to evaluate the right-hand side of Eq. 2, and then integrate the equations of motion to generate new, improved parameter settings. Eq. 2, and its motivating Eq. 1, overlap with a surprising number of different fields. For example, the matrix elements in Eq. 2 resemble kinetic coefficients, suggesting the interpretation that f ðxÞ generates a thermodynamic force that pushes the system to solve the design goal (19). Alternatively, Eq. 2 appears in the optimization and mathematics literature as natural gradient descent: It bears the interpretation of a gradient descent method that takes the steepest step such that the change in entropy stays bounded (20)(21)(22). Indeed, the matrix in front of Eq. 2 is the Fisher information metric and is constraining the driving force to move in directions of small entropy change (21). This interpretation is also associated with state-of-the-art optimizers like the covariance matrix adaptation evolution strategy (CMA-ES) (22)(23)(24); however, here the design parameters λ i are treated as random variables drawn from a Gaussian distribution irrespective of the design problem, and a version of Eq. 2 is used to update the mean and covariance of this auxiliary distribution. This is in contrast to our proposition that λ i should be treated as deterministic variables that evolve according to Eq. 2 with randomness entering only at the level of material configurations. If the task considered is changed to finding the best-fit model parameters to a given set of data, then Eq. 2 represents the direction in parameter space that decreases the fit error most efficiently per unit of behavioral change in the model (25). In this scenario, one can consider regularizing Eq. 2 to produce the Levenburg-Marquart algorithm modified to account for the geometric aspects of the optimization. Finally, one can note that the motivating equation, Eq. 1, is the replicator equation from game theory and evolutionary biology (26)(27)(28). Thus, one could also interpret the dynamics as a process of reproduction and competition in a continuous parameter space (26,27), projected onto ρðxjλÞ.
Whatever the picture, Eq. 2 has a number of powerful properties. In particular, Eq. 2 is invariant to any invertible reparameterization of λ i , including rotations, dilations, and translations in parameter space. If the reward function, f ðxÞ, is chosen correctly, the velocity flow is also invariant to rank preserving changes in the design goal (22). Thus, there will be no difference in performance between two design problems that differ by coordinate choice over λ i and/or a rank preserving change in the design goal [e.g., gðxÞ vs. expðgðxÞÞ], provided the initial values of λ i are the same. These invariances also provide stability to the algorithm: by making the search algorithm invariant to both the goal function magnitude and the parameterization, the effect of sampling errors gets bounded in a parameterization invariant way. Thus, errors from sampling parameters in Eq. 2 will not cascade, even if the matrix in Eq. 2 becomes illconditioned. Altogether, these features greatly simplify the optimization task because, now, the designer is free from worrying about trivial choices surrounding λ i and the goal function. Details on these points can be found in Supporting Information.
For the task of optimizing materials, we stress one further property: By using Eq. 2 optimization takes place in configuration space, rather than in an auxiliary space introduced to define an optimizer. As we will show, this gives a unique advantage in applying Eq. 2 to materials design: More information is used by the optimizer when updating parameters, without incurring an increase in computational cost. Perhaps best of all, this optimizer is constructed by straightforwardly applying the formalism encoded in Eq. 2 to the relevant statistical model. For example, when ρðxjλ i Þ is given by the canonical ensemble, ρðxjλ i Þ ∝ exp½−λ i h i ðxÞ, the optimizer follows immediately from Eq. 2 as _ λ i = −Cov½h i ðxÞ, h j ðxÞ −1 Cov½h j ðxÞ, f ðxÞ. Ultimately, Eq. 2 makes the transition from describing a material to designing a material in just one step.

Comparison Against Black-Box Optimizers
To demonstrate these strengths, we test our method against standard approaches that feed simulation parameters into a model by way of a black-box optimizer. We use adaptive simulated annealing (ASA) (29) and the CMA-ES (24). In each test, we allow the optimizers a fixed budget of material simulations, and each simulation requires a fixed amount of computational power. Thus, in our comparisons, computational cost and number of simulations are equivalent and an efficient optimizer is one that solves a design problem simulating as few candidate materials as possible. Implementation details are in Supporting Information.
As a first example, we task these two black-box algorithms and our approach with designing a square-lattice Ising model to maximize the magnitude of its magnetization. To do so, each method is allowed to vary the coupling constants that define the energy of spin alignments in the horizontal (J x ) and vertical (J y ) directions. To implement the black-box methods, we allow each optimizer to guess a set of coupling constants and evaluate the quality of that guess by computing the average magnetization. We find that both ASA and the CMA-ES struggle when searching entirely in the zero magnetization phase. This is an obvious consequence of the fact that the optimizer sees no variation in the quality for each parameter setting. Consequently, it receives no guidance about how to update its parameter guesses and can at best walk randomly until finding the phase boundary ( Fig. 1 A and B).
By contrast, the updates encoded in Eq. 2 navigate a path that links one phase to the other: Fig. 1C shows the flow field generated by Eq. 2 upon taking ρðxjλÞ to be the canonical ensemble with the Ising Hamiltonian, where ½ij x denotes summing over nearest neighbors along the x direction, likewise for ½ij y , and s i denote the spin variables. Given this statistical model, the control parameters λ i for optimization become λ x = J x =kT and λ y = J y =kT. Finally, the quality function f ðsÞ is set to reward states with higher magnetizations (Supporting Information). If, for shorthand, we call the individual energy components h x = P ½ij x s i s j and h y = P ½ij y s i s j , then Eq. 2 gives the velocity A similar equation holds for _ λ y but with the variables x and y appropriately interchanged. In this form, it is clear that our method will optimize as long as there is covariation between the quality function, f, and the energy components, h i . Because magnetization and energy are correlated, even if the average magnetization is zero, Eq. 2 can purposefully optimize even when operating in regions of parameter space where the black-box methods fail. The difference between these approaches lies in the fact that black-box methods are trying to solve a problem defined over the space of λ i , whereas our approach is tasked to solve a problem defined on the space of configurations, x, via Eq. 1. The CMA-ES generates multiple guesses of parameters from a Gaussian distributed over all possible λ and ASA samples by assigning an energy value to each choice of λ. These methods associate only one piece of information, the quality function, to full ensembles of configurations defined by each choice of parameters. This is in contrast to our approach that tries to solve the problem of reproducing configurations, x, that are better than average. Consequently, Eq. 2 is able to use information about how fluctuations in configurations correspond to fluctuations in quality because it has been built to exploit the extra fact that the simulation data were generated from a known distribution, ρ.
As a second example, we consider a thermalized particle trapped on a 2D substrate, defined on the x 1 − x 2 plane and at thermal equilibrium. The substrate applies a potential to the particle, making some positions more likely than others. We task the optimizer with trapping the particle in a specific potential well. To do so, we give the optimizer the freedom to tune the interaction strength with the substrate, and we give it control over a linear electric field to drive the particle. To make the problem interesting, we use a rough substrate potential: With the field included, the total Hamiltonian becomes where v x 1 and v x2 represent the field strength in the two coordinate directions. To simplify the form of Eq. 2, we represent these effects to the optimizer by defining λ s = 1=kT and absorb a factor of kT into the field coupling constants so that ρ ∝ exp½−λ s h s − λ x 1 x 1 − λ x2 x 2 . In discussing the optimizer's performance, however, we will convert the results to the original, physical variables kT, v x1 , and v x2 . The solution to this problem requires the design engine to make the target well the global minimum, and cool the system to zero temperature to trap the particle. For definiteness, the target well is the point ð5,5Þ and we initialize the algorithm with the substrate at 1kT and the field parameters set to zero.
In Fig. 2A we plot the energy landscapes generated by the optimizer, as well as the points sampled during each iteration. Indeed, the optimizer quickly learns to tilt the landscape, correctly making the target well the global minimum, and then cools the system (Fig. 2B), trapping the particle in the well. By comparing the performance against black-box approaches (Fig. 2C), our approach is both faster and more reliable: It correctly tilts the well after only 35 iterations, whereas it takes ∼ 100 simulations for the CMA-ES and ASA. Further, neither of the blackbox methods learns to completely cool the system in the allotted 1,000 simulated ensembles.
We speculate that this shortcoming is again a consequence of indirect problem representation. We note that when kT ≈ 0.1, the particle is almost evenly distributed between both the central well and the four nearest neighbors that surround it. Thus, for black-box methods, noise in the average particle position can play an overpowering role during parameter updates. By contrast, Eq. 2 considers covariances at the finer scale of configuration space. As with our prior example, this extra information makes our method appreciably more robust to flat regions in the search landscape and in this case yields an essentially exponential convergence to the optimized state (Fig. 2B).

Designing a Polymer to Fold into an Octahedron
The success of the two simple examples in the previous section invites more complicated design problems. As an example, we consider a basic model for a polymer: a string of hard, colored balls interconnected by rigid rods. The balls are weakly attractive, interaction strengths determined by the colors. For example, red and blue may be attracted more strongly than blue and green.
In principle, by tuning the color interactions, it should be possible to fold the chain into specific, desired shapes. To make a concrete task, we take a chain of six particles and create an optimizer to fold them into an octahedron, defined by minimizing the sum of the distances to the center of mass. Note that the search space is appreciably Here we treat the problem of a thermalized particle, trapped on a sinusoidal energy landscape superimposed on a quadratic background with a minima at the origin, depicted by the energy contours in A. The optimizer is given control over the system temperature kT plus a linear field potential in the two coordinates v x1 and v x2 . Its task is to trap the particle in a specific well located off center at ð5,5Þ in the x 1 − x 2 plane, marked by a cross in A. The sampled particle locations, given each choice of parameters generated, after evaluating 0, 10, and 100 iterations of Eq. 2 are plotted as red points. The optimizer uses the field to first tilt the potential and make the target well the global minimum, and then it cools the system. In tasking the same problem to black-box optimizers, we find that both ASA and the CMA-ES are able to tilt the well, but never learn to cool the system (B). Comparing how the temperature and field parameters change at each iteration (C) against d ave in B shows that the exponential convergence occurs in concert with an exponential decrease in system temperature. We note that the field parameters v x1 and v x2 track one another, reflecting the fact that the optimizer is invariant to rotations in the configuration space: The optimizer moves the field along the direction associated with the greatest improvement in solution quality and is insensitive to the fact that the problem was parameterized in the arbitrary coordinates x 1 and x 2 .
larger than in the prior examples (dimension 10), and that simply setting all of the interaction strengths to large values will not produce the optimal solution. By choosing the interaction appropriately, the same energy can be given to the octahedral and polytetrahedral configurations for identical coupling constants. Entropic arguments imply that the polytetrahedron will dominate the chain configurations unless the optimizer carefully adjusts the coupling constants to take on unique values (30,31). Fig. 3A shows a typical chain configuration from each generation, and Fig. 3B shows the median sum of distances to the center of mass, normalized relative to that of a perfect octahedron. Initially, the coupling constants are set to 1kT, and random chain configurations are typical. However, as the optimizer drives the interaction energies to larger values, the shapes become compact and structured. Around 200 generations, virtually every shape generated is octahedral (Fig. 3A).
By plotting the values of the interaction strengths against iteration number, we find the optimizer's solution is simple, logical, and arguably optimal. Early on, the optimizer attempts to meet the design goal by simply increasing the coupling strengths to make more compact objects (Fig. 3A). However, as the coupling constants are undifferentiated, the results are predominately polytetrahedral geometries. To compensate, the optimizer deactivates three coupling constants around 100 generations and sends the remainder to infinity (Fig. 3C). The logic behind this maneuver becomes clear by plotting the interactions as a network: The active interactions plus the polymer backbone form the contact graph of an octahedron. This strategy, transforming the contact matrix to an interaction matrix, has been identified as an approach to programing, by hand, the optimal interaction parameters for selfassembly (18). In fact, the specific problem of creating a selfassembling octahedron has been solved using a virtually identical motif (32). Evidently our optimizer can reproduce a well thought-out approach to self-assembly, and it does so automatically, requiring only a model and a design goal.
Optimization of an Out-of-Equilibrium System Because Eq. 2 holds for any parameterized probability distribution function, it can be used to create optimization schemes beyond the canonical ensemble in the prior examples. The only essential ingredients are a model that predicts the probability of microstates, an engine that samples configurations from said model, and a design goal. As simple extensions, chemical potentials or constraints on pressure could be included as tunable parameters (33). A new theoretical concept, termed "digital alchemy," extends statistical mechanics to account for microscale geometric parameters, such as the particle shape in a colloidalnanoparticle assembly (34). Thus, by coupling this approach with our optimization formalism, particle geometry can be tuned to produce optimized bulk responses. One can also note that the range of parameters to design is at the user's discretion: Eq. 1 can be used to rederive Eq. 2, assuming that some of the model parameters are not controllable by taking them to be time independent. Indeed, for the particle in a well problem, the wavelength of corrugation was taken as a fixed parameter and the resulting optimizer was quite effective. One can also consider optimization for global quality functions that exist over multiple a range of parameters (35). For instance, our approach can optimize the density of a crystal lattice over a range of pressures and system volumes by defining ρ = ρ 0 ðxjλ, V , PÞUðV , PÞ, where UðV , PÞ is a uniform distribution for V and P over a range of consideration and ρ 0 is the appropriate distribution for microstates given a fixed volume and pressure. Eq. 2 can then be applied to optimize in this extended parameter space to find choices of λ that work well over a range of possible densities and pressures. Abstracting the concept, statistical models could be based on calculations like self-consistent field theory or experiments with the code directly measuring correlation functions in the laboratory and tuning physical parameters.
Moreover, nonequilibrium processes, provided they have a statistical description, are fair game. For example, if one simulates diffusion by adding white noise to a mean drift, then the paths are distributed by a product of Gaussian distributions conditioned on the prior steps. Clearly, the paths are statistical objects, with diffusion and drift as the distribution parameters, λ. Thus, one can build an optimizer that tunes these control parameters using Eq. 2, even if they are time-dependent functions.
As proof of this point, we return to the problem of a particle trapped on a substrate, but now simulate the particle dynamics. The applied field and the system temperature are treated as functions of time and the optimizer is tasked with moving the particle from one well to another in a given interval. Fig. 4A shows the median distance to the target well after executing the optimizer's processing protocol in each generation. In the first 60 generations, the optimizer learns to transport the particle from its starting location to the target well via a large, deterministic driving force. It then spends the remaining iterations monotonically decreasing the system temperature while developing a trapping protocol with the field. After 2,000 iterations, the optimizer seems to traps the particle by oscillating the driving force, changing its direction before the particle can transition to another well. In effect, the optimizer learns to drag the particle to the target and trap it in place, using both the temperature and the field. By 2,000 iterations, ∼ 90% of the points in the path fall within the target well. When left to run longer, the optimizer continues to improve the quality of solutions, but at the cost of becoming unphysical, the optimizer generates extremely large fields that move the particle to the well faster and faster. To optimize beyond the proof of concept demonstrated here, one may have to restrict the range of parameters allowed to the optimizer or account for arbitrary velocities by increasing the number of steps in the walk.
Optimization of Directed Self-Assembly As a final demonstration of our approach we consider the realworld problem of designing the directed self-assembly of block copolymers on a chemically patterned substrate. This represents a task at the forefront of both materials design and sublithographic Early on, the polymer configurations are dominated by thermal energy and random, yet as the interaction strengths increase and differentiate, more structured objects appear, and ultimately only the octahedron configuration exists (iterations 190-210). (B) Plotting the median percentage of deviation between the polymer's radius of gyration R g and the radius of gyration for an octahedron R oct at each iteration shows that the optimizer again produces a monotonic decrease at each step, with an effectively exponential convergence in the last 30 iterations. By the last 10 iterations, the median deviation from a perfect octahedron is roughly 1%. (C) In plotting the coupling constants against iteration number, we find that the optimizer adjusts coupling constants in groups. It is clear that the active coupling constants form the contact network for an octahedron (C, Inset).
fabrication (1, 7-11, 36, 37). The goal is to lithographically pattern a substrate with a small number of chemical features such that these features promote block copolymers to self-assemble into a desired target morphology. Here, we consider a task that has been identified as a promising candidate for the manufacture of nextgeneration semiconductor devices and high-density storage media: self-assembly of AB-diblock copolymers into an ordered striped or lamellar morphology (1,(9)(10)(11). We use a theoretically informed coarse-grain model for block copolymer simulations (36). Polymer chains are simulated as beads that are linked together. The system is considered at fixed temperature and volume, and thus the probability of finding a given microstate configuration is defined in terms of an energy given by three parts. The first is a linear spring bond energy between beads in each polymer chain. The second is a nonbonded interaction energy that characterizes repulsion from unlike species and the material compressibility. Details are in Supporting Information (36). Both the model and the parameters in it represent a polystyrene-block-poly(methyl methacrylate) (PS-b-PMMA) diblock copolymer with a number-averaged molecular weight of 22,000-b-22,000 and a stripe period of 28 nm.
The final contribution to the energy is the substrate interaction. The substrate consists of two regions: the patterned stripes of width w and the background. Both are defined to have shortranged effects on the polymer beads and assume the form H s = kT = ΛðαÞ=d s exp½−ðz=2d s Þ 2 , where d s defines the decay length of the interaction, z is the distance from the plane of the substrate, and ΛðαÞ is the interaction strength between the substrate and a bead of type α. Thus, if the particle is over the guiding stripe and of type A, ΛðAÞ = Λ s . If the particle is of type A and over the background region, ΛðAÞ = Λ b . Following ref. 11, we simplify our model by assuming the interactions are antisymmetric: Λ s ðAÞ = −Λ s ðBÞ and Λ b ðAÞ = −Λ b ðBÞ. The design problem is to adjust the width of the strips w and the two energy parameters Λ s and Λ b so the target stripe phase replicates itself m times between two guiding stripes spaced by the polymer period multiplied by m.
The results for m = 3 and m = 6 density multiplication are shown in Fig. 5. For the m = 3 problem, we ran the optimization four times, varying the time step used in integrating Eq. 2. In every instance, the optimizer not only brought the system to a state that successfully met the design goal, but also, within noise, converged to the same state. The optimized parameters suggest that directed self-assembly is best achieved by setting the stripe width equal to roughly half the polymer period, Λ s ≈ −1kT and Λ b = 0.05kT. All of these parameters agree with simulation results obtained by a brute-force solution to the problem and experimental verifications performed on the real polymer system (37) and are physically consistent with optimization results obtained for triblock copolymer pattern multiplication (11). These results can be explained by considering the interfacial energies in the system. The background interaction is required to be weak because the background region has roughly equal coverage between the A and B phases and is significantly larger in area than the size of the stripe. Moreover, the interaction strength for the stripe components is larger in the m = 6 problem (Fig. 5D), because the larger distance between patterned regions requires stronger anchoring to guide assembly.
When the optimizer was run at the most aggressive time step, we were able to achieve convergence for m = 3, (m = 6) in less than 10 (18) iterations, despite the fact that the material required the simulation of roughly 50,000 (100,000) polymer beads. We stress that the performance obtained here is not a consequence of initializing the system too close to an optimal state, but rather evidence of the power behind Eq. 2. Fig. 5 shows that indeed the  Fig. 4. Optimizing a nonequilibrium process. By using Eq. 2 we can tune processing protocols for a Brownian particle walking on a rough energy landscape controlled by a time-dependent temperature, kT, and linear mean drift components, μ x and μ y . The optimizer has been tasked to adjust to place and trap the random walker in a well located at the x-y coordinates ð10,10Þ. (A) Ensemble median distance to the objective well after executing a processing protocol at each iteration of the algorithm. Callouts show representative paths taken by the particle, and contours in the callouts show lines of constant energy over the substrate potential. The large image represents the protocol executed after 2,000 iterations. Every protocol attempted at 10-iteration intervals is illustrated in B-D, where the temperature, kT (D) as well as the applied fields in the x direction and y direction normalized by temperature, μ x =kT (B) and μ y =kT (C), are plotted against time. At t = 0 the particle is released from its initial position at ð−10,0Þ and allowed to wander and the processing protocol is executed until the simulation is stopped at t = 1. At each iteration, the optimizer works to monotonically decrease the temperature, while arriving at a field protocol that quickly drives the particle to the target well and then oscillates the fields to trap it there.

A B
C D E Fig. 5. Optimizing directed self-assembly. Here we design the width of the stripe, the strength of its attraction to the red polymer beads Λ s , and the attraction strength of the background substrate Λ b toward the blue polymer beads, to match the self-assembled phase as closely as possible to the target of alternating stripes. We quantify the success of our optimizer by comparing an order parameter ΨðxÞ = ðn a Þ=ðn a + n b Þ binned along the x axis of the box and averaged over y and z to the target stripe pattern. With Eq. 2, we are able to produce optimized parameters for 3× (A) and 6× (B) density multiplication after simulating 10 and 20 parameter choices. Two characteristic configurations are plotted in A and B, Insets, separated by just a handful of iterations, yet displaying markedly different phases of the polymer. Asymptotic configurations depicting Ψ (solid, marked line) and the target (dashed line) (A and B, Insets) show that the optimized parameters match the desired morphology, typically within 80% or better. In plotting the parameters generated by Eq. 2 (C-E), we find that the interaction with the background brush is the most relevant parameter in directed self-assembly. For both the 6× and 3× problems, the rapid convergence toward the optimized state takes place once the background strength is reduced to Λ b ≈ 0.05. After a weak background is established, the strip width and strength function as fine-tuning parameters that facilitate defect-free assembly.
initial parameters do not produce a solution to the design problem, let alone a stripe pattern of any kind. This particular problem also has been attempted in some variant using other materials design methods. In fact, our initial conditions were selected to match those given to an implementation of the CMA-ES, solving the same design problem but using triblock copolymers (11). For that problem, the CMA-ES took roughly 50 generations to converge, simulating 32 ensembles in parallel per iteration, each requiring 200,000 Monte Carlo steps. Our approach also used 200,000 Monte Carlo steps per iteration, but required only 10 iterations to converge. If algorithm performance is measured in terms of the number of microstates simulated, then solving directed self-assembly problems by way of the CMA-ES requires at least 5 times as much compute power as the approach proposed here. If it is not possible to run the ensemble simulations in parallel, our approach is roughly 130 times faster than the CMA-ES and completes a full optimization process before the CMA-ES has completed a single iteration. Additionally, inverse Monte Carlo methods have been used to solve directed selfassembly problems involving the placement of guiding posts instead of stripes (8). Although there are relevant physical differences, we note that the results presented for inverse Monte Carlo converge after simulating roughly 30 million microstates. Because this number of microstates simulated is roughly 15 times larger than what was used here, we can speculate that our proposed methods could also be faster for such applications.

Conclusions
To the extent that the goal of materials design is a unified framework that handles a wide range of complex inverse problems, we believe the formalism introduced here represents a significant step forward. By applying Eq. 2, we can solve problems with flat search landscapes (Fig. 1) and multipleinteraction types (Fig. 2), incorporate constraints (Fig. 3), tune processing conditions (Fig. 4), and address application scale design and optimization tasks (Fig. 5). Furthermore, in all of the examples presented, the end result is intuitive even though it was achieved in a complicated search landscape where other optimization schemes struggle or fail. Finally, the fact that processing conditions such as applied fields or temperature protocols and model parameters like internal interaction energies can be optimized with the very same framework presents an unexplored direction for materials design. Because these are the essential aspects that determine the properties of any material, the capacity to tune both simultaneously, one accounting for the other, opens the doors to a more coherent and conceptually complete design program.