Designing athermal disordered solids with automatic differentiation
Mengjie Zu Carl Goodrich


APS/123-QED





 Institute of Science and Technology Austria
 mengjie.zu@ist.ac.at

 carl.goodrich@ist.ac.at

 Institute of Science and Technology Austria






















             


The ability to control forces between sub-micron-scale building blocks offers considerable potential for designing new materials through self-assembly. A typical paradigm is to first identify a particular (crystal) structure that has some desired property, and then design building-block interactions so that this structure assembles spontaneously. While significant theoretical and experimental progress has been made in assembling complicated structures in a variety of systems, this two-step paradigm fundamentally fails for structurally disordered solids, which lack a well-defined structure to use as a target. Here we show that disordered solids can still be treated from an inverse self-assembly perspective by targeting material properties directly. Using the Poisson's ratio, ν, as a primary example, we show how differentiable programming connects experimentally relevant interaction parameters with emergent behavior, allowing us to iteratively “train" the system until we find the set of interactions that leads to the Poisson's ratio we desire. Beyond the Poisson's ratio, we also tune the pressure and a measure of local 8-fold structural order, as well as multiple properties simultaneously, demonstrating the potential for nontrivial design in disordered solids. This approach is highly robust, transferable, and scalable, can handle a wide variety of model systems, properties of interest, and preparation dynamics, and can optimize over 100s or even 1000s of parameters. This result connects the fields of disordered solids and inverse self-assembly, indicating that many of the tools and ideas that have been developed to understand the assembly of crystals can also be used to control the properties of disordered solids. 


  keywords
self-assembly | disordered solids | auxetic materials | simulation | automatic differentiation



Designing athermal disordered solids with automatic differentiation
    Carl GOODRICH
    April 29, 2024
===================================================================







§  INTRODUCTION
























To what extent can the material properties of disordered solids be controlled? Recent results show that the properties of random spring networks close to the isostatic point can be radically and precisely tuned through slight adjustments to the network topology <cit.>. More specifically, the Poisson's ratio, ν, of a randomly generated network can be tuned to either the upper or lower bounds by removing only ∼ 1% of the springs – the choice of which springs to remove determines the final value of ν <cit.>. This is possible because each spring's contribution to the bulk modulus is independent of its contribution to the shear modulus, meaning that their ratio, which determines ν, can be tuned by the choice of removed springs.
However, in a real material, it is usually not possible to make precise, targeted alterations to structure. In fact, the notion of removing springs from a network simply does not translate to materials made up of particulate building blocks, ranging from atoms or molecules to larger objects like proteins or colloids. Unlike a spring network, one cannot simply remove the interaction between two particular particles, let alone do so in a scalable way. 

Nevertheless, this paper proposes a strategy for tuning the properties of particulate-based disordered solids. This strategy is akin to an inverse self-assembly approach, except that interaction parameters are adjusted to tune properties rather than structure. The necessary connection between interaction parameters and emergent properties is made by exploiting a class of numerical techniques called Automatic Differentiation (AD), enabling the exploration of high-dimensional and complicated design spaces. Once this connection is made, the system is “trained" similar to how one trains a neural network. 



We demonstrate this strategy by considering the simple case of athermal sticky spheres. Specifically, we consider a two-dimensional system of N particles divided evenly into  species, where the diameter of each species, σ_α, and the binding energy between each pair of species, B_αβ=B_βα, can be continuously varied. 
Motivated by DNA-coated colloids, the particles interact via a Morse potential with the short-ranged repulsive part replaced with a finite soft repulsion (see Methods)  <cit.>.

The system is prepared at zero temperature following the protocol developed in the study of the jamming of soft spheres <cit.>, where particles are placed randomly (corresponding to infinite temperature) and then quenched to the nearest local energy minimum. 





The final material properties, such as the Poisson's ratio, clearly depend on the + (+1)/2 values of σ_α and B_αβ, but we are left with two critically important questions. First, can these parameters be adjusted in order to accurately and precisely control material properties? For example, how close can ν be tuned to a particular target ? Second, can multiple properties be controlled simultaneously and independently, enabling highly nontrivial design? We will show that, even for our exceedingly simple model, the answers to both questions are a definitive yes. 

This result is obtained by directly connecting the + (+1)/2 model parameters to changes in material properties. More precisely, we construct an objective function , e.g. = (ν - )^2, that indicates how far a property (e.g. ν) is from a target (e.g. ), and use Automatic Differentiation (AD) <cit.> to calculate the gradient ∇_θ of  with respect to the parameters θ = {σ_α}∪{B_αβ}. ∇_θ indicates how changes to the parameters affects the objective, and there are numerous gradient-descent-based algorithms <cit.> for using this information to minimize , thus tuning the property of interest. See Methods for more details.

At zero temperature, the N particles arrange themselves into one of many possible “configurations," or local minima of the potential energy landscape. In this paper, we will manipulate material properties first on the level of individual configurations, meaning we will pick a random configuration b and tune the parameters until,  e.g., ν_b =. We find that this works surprisingly well all the way to the perfectly auxetic limit → -1. Interestingly, the Poisson's ratio of multiple configurations can also be tuned precisely and simultaneously, although this becomes more challenging as more configurations are optimized at the same time. We then consider the ensemble level, and again find a surprising ability to manipulate the average Poisson's ratio  within the range 0.2 <  < 0.7. Importantly, this approach is highly general, and we use it to tune other properties and even multiple properties simultaneously.














§ RESULTS



 §.§ Training individual configurations



To begin, we select a configuration b by placing N particles randomly in a 2d periodic box and minimizing the energy to the nearest minimum. We then define the objective = (ν_b - )^2, where ν_b is the measured Poisson's ratio and  is the desired target chosen from between -1 and 1, which are the theoretical bounds for isotropic systems in two dimensions.
Figure <ref>b confirms that our AD-based calculation of ∇_θν_b accurately predicts the change in ν_b over finite changes in the parameters θ. This is shown for a particular representative example where N=368 particles are evenly divided into =2 species with initial diameters of 0.8 and 1.0, a number density of ρ=1.6, and with constant binding energies B_αβ=0.1. Before training, we measure ν_b≈ 0.453. We train the system using standard gradient-descent-based algorithms (see Methods) to iteritively adjust θ before recalculating ν_b and ∇_θν_b. At each step, we reminimize the energy with respect to the particle positions, allowing us to track the configuration as parameters change. Figure <ref>c shows how  and ν change during 1600 iterations. 

Despite only having 5 parameters, the Poisson's ratio is tuned exactly to the target, meaning that the Poisson's ratio at this configuration is successfully and accurately controlled.


Figure <ref>d shows how the parameters evolve during training for 35 successfully trained configurations. Clearly, the final parameters depend strongly on the particular configuration; while the systems remain bidisperse and thus disordered, there is no obvious trend in the relative particle diameters nor in the three binding energies. 

Interestingly, for  less than approximately 0.2 or greater than approximately 0.7, the final Poisson's ratio seems to be dominated by a small localized region. To see this, we calculate K_ij and G_ij, which are the contributions of the interaction between particles i and j to the bulk and shear moduli, respectively, so that K=∑_ijK_ij and G=∑_ijG_ij  <cit.>. We observe a typical probability distribution for P(K_ij) and P(G_ij) when 0.2 < ν < 0.7 (orange data in Fig. <ref>e), with no noticeable spatial correlations in the largest positive and negative values. However, for ν < 0.2, the positive and negative tails of P(K_ij) widen considerably (blue data in Fig. <ref>e) but equally so they largely offset each other, and we find that the largest values are spatially localized (inset). We hypothesize that the broad tails in P(K_ij) enable more dramatic tuning of K, allowing K/G to become small and thus ν to become negative. Conversely, for ν > 0.7, the tails of P(G_ij) broaden, enabling K/G to become large. See SI for more comprehensive data. While we do not fully understand the origin of these extreme and localized regions, we note that the range of ν where they occur coincides with our ability to train the ensemble average Poisson's ratio, see below. 













Figure <ref> shows the probability P_success that we are able to successfully train a configuration of N=368 particles for various target Poisson's ratios and different numbers of species. Training is considered successful if the objective  decreases below _thresh = 10^-6 within 10^4 optimization steps. Not surprisingly, P_success increases with the number of species since that gives access to additional parameters, and with more conservative targets that are closer to the initial Poisson's ratio of approximately 0.5. We also find that the range of the attractive interaction matters, with shorter interactions making it slightly harder to train. 

These results show that we can accurately and precisely control ν_b over a wide range. Interestingly, Fig. <ref>a-b show that the parameters learned through this process influence the entire ensemble of configurations even though they were never considered in the training process. For example, the parameters obtained from training a single system to =-1 results in = 0.33± 0.02, well below the initial untrained value of =0.49. This is a first indication that we can manipulate ensemble averages.











This effect can be amplified by training on multiple configurations simultaneously. Specifically, we next choose a set of  configurations, with  ranging from 2 to 7, and define the objective = ∑_b^' ( ν_b^' - )^2, where here the sum is over the  chosen configurations b^'. Figure <ref>c shows that training becomes harder (less likely to succeed) as  increases, especially for aggressive targets (<0). Despite this, however, increasing  does increase the effect on , as shown by Fig. <ref>d. 



 §.§ Training ensemble-average quantities

We now directly consider the ensemble average Poisson's ratio, , by defining the objective = ( - )^2. Note that this is not the same as the →∞ limit of the objective from the previous section because here we are only concerned with the mean of the Poisson's ratio, not its value at every individual configuration. The other conceptual difference is that we can only ever estimate  by averaging over a finite number of (randomly chosen) configurations, meaning that all calculations of  and ∇_θ are necessarily stochastic. This is analogous to training in many machine learning models, where data is “batched" and gradients are highly noisy. Importantly, as demonstrated in Fig.S1, average gradients of  are predictive in the same way as in Fig. <ref>b, allowing the use of stochastic gradient descent optimization. 























Figure <ref>a shows  during optimization with a target of =0.35 using =2 species. At each step, we average over only 8 systems of size N=368, leading to large noise in . This leads to a systematic overestimation of  whenever the noise is larger than -. Therefore, in order to ascertain how well the training has done, we also perform “validation runs," where we use the final parameters obtained after 5000 optimization steps and calculate  over a fresh set of 2000 configurations. As indicated by the green circle, _val≈ 10^-6, meaning that we have trained  to  with an accuracy of 0.001. The purple star in Fig. <ref>a shows an alternative validation where  is calculated from 100 configurations of N=2048 particles each, demonstrating that even though our results are trained using small systems, the solution nevertheless applies to much larger systems. 



Figure <ref>b shows the 5 parameters at each optimization step (thick lines). The decrease in parameter fluctuations is caused by our variable learning rate (see SI), but the lack of systematic trends over the final 500 optimization steps suggests that the optimization has converged. 

The thin lines show the same data for 29 other training runs that also target =0.35 and all achieve similarly small _val. Clearly the final parameters are not unique, meaning that the solution is degenerate, but there are still clear trends in the parameters (e.g. small particles bind more strongly to large particles than to other small particles). 


Figure <ref>c shows the final validated _val as a function of the target , showing that we are able to successfully tune the Poisson's ratio between roughly 0.2 and 0.7. 
The existence of such bounds in our ability to tune  is expected, especially for low  as creating low- and negative-Poisson's ratio materials is notoriously challenging <cit.>. We expect that more complicated systems, for example with non-spherical particles, could expand these bounds.

Notably, this range of roughly 0.2 to 0.7 coincides exactly with the onset of broad and localized bond-level response observed in Fig. <ref>e. While training the Poisson's ratio of an individual configuration outside of this range is possible, it appears to rely on configuration-specific structural motifs that do not generalize to the full ensemble. This is also consistent with the observation in Fig. <ref>d that the final parameters vary dramatically from one configuration to another. 





Interestingly, we find that considering more species does not increase the range over which we can tune the system (blue triangles in Fig. <ref>c). In principle, increasing the number of species should increase tunability, but there are no guarantees of this. It is unclear if improved training methodology could lead to better results with more species or if these bounds are imposed by physical constraints. Finally, Fig. <ref>d shows the final parameters after 5000 optimization steps for different , suggesting a general design strategy for tuning  over the range 0.2 <  < 0.7. 




















 §.§ Training multiple ensemble-average quantities simultaneously

So far, we have focused on tuning the Poisson's ratio, but Fig. <ref> shows that we can tune other quantities in exactly the same way, specifically the pressure and a structural order parameter , which measures small amounts of local 8-fold symmetry in neighbor orientations (see SI). Now we ask whether we can tune multiple average properties simultaneously. Figure <ref>a shows the simultaneous optimization of the Poisson's ratio (=0.5) and pressure (=0.4), which results in a final validated loss of _val = 6× 10^-4, demonstrating success. To be precise, we use a combined objective function = ( - )^2 + ( - )^2, where the two components are plotted separately. 

Not all values of (, ) can be obtained, and Fig. <ref>b shows the final (_val, _val) obtained from a series of systematic trials with different (,) (indicated by light grey dots). This data shows a well-defined two-dimensional region of tunability, within which simultaneous optimization can be obtained. 



§ DISCUSSION

We have shown that the properties of 2d athermal disordered solids can be inverse designed using Automatic Differentiation to connect properties of interest to particle-particle interactions. We have focused primarily on the Poisson's ratio, ν, because as a unitless measure of elasticity it cannot be adjusted by trivially changing energy scales. Using a simple model of sticky spheres with tunable species-level interactions that is motivated by DNA-coated colloids <cit.>, we have shown that the Poisson's ratio of a single configuration can be tuned with a high success rate anywhere in the range -1 ≤ν≤ 1 (the theoretical bounds for two-dimensional isotropic systems), and the ensemble average  can be tuned in the range 0.2 <  < 0.7. Furthermore, we can tune the pressure and local structural order in a similar way, and even design for multiple properties simultaneously. 

This can be compared with typical self-assembly paradigms, where one manipulates interactions so that, for example, a desired crystal structure becomes the thermodynamic ground state <cit.>. Instead, our approach targets the statistics of the metastable states. This has the distinct advantage that the assembly protocol (athermal relaxation in our case) is baked into the calculation of these statistics, and so we do not need to worry about kinetic traps or other barriers to assembly. 

While other systems exhibit a variable Poisson's ratio – for example, jammed packings of soft spheres display → 1 near the unjamming transition [Note that while our model reduces to the soft-sphere model by setting B_αβ=0 and D_1 = 1.4 D_2, the → 1 result is only obtained for finite systems in a fixed-pressure ensemble, and is therefore inaccessible in our fixed-volume ensemble.] <cit.>, our work demonstrates a new level of systematic and robust control. 

In addition, we have shown that directly connecting objectives to parameters through the gradient ∇_θ enables the targeted design of multiple properties, which is otherwise challenging. 
Together, this presents a scalable approach to inverse design and reveals a direct strategy for the targeted manipulation of disordered solids. 

Unlike for ensemble averages, we are able to precisely tune the Poisson's ratio of individual configurations with high success even in the auxetic or fully incompressible regimes. This is analogous to the bond-level tuning in Refs <cit.> except it is at the level of species-species interactions and maintains the constraint of force-balance on every particle. Figure <ref>e also hints at an unexpected localized mechanism for creating auxetic materials. While these individual configurations are likely difficult to obtain experimentally due to the overwhelming number of competing structures, it is nevertheless telling how designable individual configurations are because it implies that the lower bound of  is inherently a collective effect. At ≈ 0.2, any particular ν_b can still be lowered, but doing so necessarily increase the Poisson's ratio of other configurations so that  is unchanged. 

One of the more surprising results is that increasing the number of species, and thus the number of parameters, does not increase the range over which we can tune . This implies a fundamental constraint that is perhaps related to the general rarity of low and negative Poisson's ratio materials and leaves open the challenge of using our approach in more sophisticated models, with other types of design parameters, to obtain < 0.




This paper has only scratched the surface of how disordered solids can be designed with Automatic Differentiation. Exciting extensions of this include 1) considering more complicated building blocks, for example with nontrivial shapes or anisotropic potentials, 
2) considering different and more complicated preparation protocols, and 3) focusing on a wider range of material properties (e.g. nonlinear stress-strain behavior <cit.>, density of states, and allosteric responses <cit.>). Recent advances in differentiable programming ecosystems, including the Python packages JAX and JAX-MD used in this work, make it possible to incorporate Automatic Differentiation (AD) into many existing calculations. Important considerations, which are discussed in the context of the present work in the Methods section, include the predictiveness of gradients, strategies for avoiding large memory loads, and the importance of initial parameter guesses.

































§ MODEL AND METHODS


Model.
We consider a system composed of N particles divided evenly into  species. 
Particles interact via a “harmonic-Morse" pairwise potential given by

    V(r)=
         k/2(r-σ̅)^2 - B,     r < σ̅
    
        B(e^-2a(r-σ̅)-2e^-a(r-σ̅)),     r ≥σ̅

where r is the center-center distance between two particles, σ̅ is the mean of their diameters, k characterizes the short-ranged repulsions, and B determines the strength of the medium-range attractions, whose extent is proportional to 1/a. Unless otherwise stated, we use k=5.0 and a=5 for all pairs of particles. However, B and σ̅ depend on the species type of the particles in question. Specifically, we independently vary the attractive strength B_αβ for every pair of species α and β (provided B_αβ=B_βα), as well as the particle diameter σ_α (so that σ̅_αβ = (σ_α + σ_β)/2). We then use the XPLOR smoothing function to truncate V(r) at a distance r_cut = σ̅+ 9.9/a. This model is commonly used to describe, for example, DNA-coated colloids where diameters and binding affinities can be manipulated at the species level <cit.>. Note that we use a harmonic repulsive force to ensure numerical stability at a random set of initial position, and that the strength k of this repulsion is decoupled from the attractive forces to allow B_αβ→ 0. 














For a given set of parameters θ = {σ_α}∪{B_αβ}, we obtain stable athermal structures by using the FIRE algorithm <cit.> to minimize the total potential energy E = ∑_< ij > V(r_ij) starting from a random set of initial positions. We always consider an even number of species, and set the initial diameters for half the species to be 0.8 times that of the other half, with an overall number density of 1.6. The initial binding strengths are all set to B_αβ=0.1. 



Optimization for individual configurations.
The original configurations are generated by quenching random packings to local energy minima using FIRE algorithm. Once an initial energy-minimized state with original setting are reached, the Poisson's ratio ν is calculated explicitly using linear response <cit.>, and combined with the target Poisson's ratio  to construct the objective function =(ν-)^2.  

Our model, the energy minimization, and the calculation of ν and  are performed entirely using the library JAXMD <cit.> and the Automatic Differentiation library JAX <cit.> on which it is built. As a result, the entire calculation is end-to-end differentiable, allowing us to accurately and efficiently calculate ∇_θ. Implicit differentiation, implemented using the JAXOPT library <cit.>, is used to propagate derivatives through the energy minimization as well as a key step in the linear response calculation, which ensures gradient accuracy and mitigates any extensive memory overhead. We then employ the RMSProp algorithm (implemented as part of JAX) to iteratively update the parameters θ based on the gradients. To improve the convergence speed, we dynamically adjust the learning rate using a meta-learning program that adaptively generates per-step hyperparameters. At each iteration, we use the final configuration of the previous step for the initial positions of the energy minimization process, ensuring that we consistently track a given energy minimum. Occasionally, the energy minimum transforms into a saddle point and the system undergoes a rearrangement. This results in a small spike in , e.g. see Fig. <ref>c at around 1300 steps.




Optimization for ensemble-averaged quantities.
The primary difference in optimization between individual and ensemble systems lies in the input configurations. In the optimization process for ensemble-averaged quantities, we generate 8 new configurations (by minimizing from random initial positions) at every step. Given our relatively small system size, this is well below the number of configurations needed to obtain accurate ensemble averages, meaning our estimates of  and ∇_θ are highly noisy. However, stochastic optimization algorithms like RMSProp are well-suited to very noisy gradients and this does not prevent convergence. 






