De Novo Atomistic Discovery of Disordered Mechanical Metamaterials by Machine Learning

Abstract Architected materials design across orders of magnitude length scale intrigues exceptional mechanical responses nonexistent in their natural bulk state. However, the so‐termed mechanical metamaterials, when scaling bottom down to the atomistic or microparticle level, remain largely unexplored and conventionally fall out of their coarse‐resolution, ordered‐pattern design space. Here, combining high‐throughput molecular dynamics (MD) simulations and machine learning (ML) strategies, some intriguing atomistic families of disordered mechanical metamaterials are discovered, as fabricated by melt quenching and exemplified herein by lightweight‐yet‐stiff cellular materials featuring a theoretical limit of linear stiffness–density scaling, whose structural disorder—rather than order—is key to reduce the scaling exponent and is simply controlled by the bonding interactions and their directionality that enable flexible tunability experimentally. Importantly, a systematic navigation in the forcefield landscape reveals that, in‐between directional and non‐directional bonding such as covalent and ionic bonds, modest bond directionality is most likely to promotes disordered packing of polyhedral, stretching‐dominated structures responsible for the formation of metamaterials. This work pioneers a bottom‐down atomistic scheme to design mechanical metamaterials formatted disorderly, unlocking a largely untapped field in leveraging structural disorder in devising metamaterials atomistically and, potentially, generic to conventional upscaled designs.


S2. Forcefield effect on two-body excess entropy
Figure S2 shows the topography of two-body excess entropy S2 in selected planes of the forcefield parameter space.Indeed, S2 shows a complex dependance on the forcefield parameters.Upon closer inspection, the S2 value generally declines with respect to the increase of bond directionality, in accordance with the loose packing of directional covalent-type bonding (see Fig. 7A in the main text).Interestingly, we find that the S2 topography shows a diagonal symmetry at the slice plane of  = 1.2, where larger bond angle is more likely to yield high-S2 structure, suggesting the formation of some loose-packed open structures when large bond angle is energy-favored.

S3. Stress shift per atom in disordered versus crystalline porous network
Taking the "diamond" forcefield as an example, its crystalline and disordered network share the same forcefield but exhibit a pronounced difference in exponent n.According to the reasoning chain (see Fig. 4A in the main text), the different exponent can be ascribed to a significant difference of local stress state in disordered versus crystalline network, thus leading to distinct local dynamics (see Fig. 4B in the main text).
Figure S3 provides the normal and shear stress shift per atom as a function of normal strain in disordered versus crystalline network upon tensile deformation.Notably, the disordered network exhibits higher normal stress shift but much lower shear stress per atom than its crystalline counterpart.This local stress state influences exponent n in two ways: Firstly, the normal stress shift per atom directly promotes the normal stress and stiffness thereof, while the shear stress shift per atom has very limited contribution.Second, the shear stress shift per atom promotes local deviatoric deformation that causes the onset of fluctuating particles and the stiffness loss thereof.From this perspective, structure disorder can linearize stiffness-density scaling through (i) more bond-stretching response that promote normal stress shift per atom and (ii) less bond-bending response that harmfully facilitate shear stress shift per atom.

S4. Fraction of fluctuating particle detected by TimeSOAP approach
Here, we conduct a TimeSOAP analysis to unveil the connection between structural disorder and the observed linearization in stiffness-density scaling (see Fig. 3 in the main text).The TimeSOAP approach excels at detecting the highly fluctuating particles from a time series of configuration evolution [ [1,2] ].
Following this approach, we first compute the vector of SOAP many-body descriptors qsoap(i) for the local neighborhood environment of particle i, that is, the normalized form of partial power spectrum vector pi [ [3,4] ], namely, where the vector pi consists of a set of elements { "" !# ()} defined in the following [ [3,4] ]: where n and n' are indices for the different radial basis function gn(r) up to nmax (herein, nmax=8 and gn(r) adopts spherical gaussian type orbitals), and l is the angular degree of the spherical harmonics Ylm(, ) up to lmax (herein, lmax = 6).
The coefficient  "#( is defined as the following inner product [ [3,4] ]: where R 3 denotes the integral sphere centered at particle i with a cutoff radius Rcut (herein, Rcut = 2.2), and (i, r) is the gaussian smoothed particle density at a distance r for all particles j within the cutoff sphere [ [3,4] ]: where rij is the distance between the center particle i and its neighbor particle j.Here, we implement the SOAP calculation using the Dscribe package [5] , and the output SOAP many-body descriptor vector qsoap(i) contains 252 normalized elements {pnn'l} for each particle i.The environment similarity between local packing of particle i and j is measured by the SOAP distance dsoap(i, j) [ [3,4] ]: where dsoap(i, j) ranges from 0 to 1 that represent, respectively, identical or zero-overlapping local environment between atom i and j.
Based on the SOAP distance, we calculate the TimeSOAP value  !5&65 of particle i at the normal strain  by [ [1,2] ]: where the SOAP distance  1234 is calculated between two normal strain states  and  + Δ for the same particle i, and the normal strain increment Δ is set as 0.1%.This TimeSOAP quantity captures the degree of dissimilarity between two local packing structures at normal strain  and  + Δ, respectively.Figure S4 shows the evolution of TimeSOAP signals with respect to normal strain for all particles in the silicon-type disordered and crystalline porous network upon uniaxial tensile deformation, wherein the TimeSOAP values are normalized by multiplying 0.1%.At each normal strain, we find that some particles exhibit higher TimeSOAP than others, and crystalline network shows very low TimeSOAP value in the initial deformation but quickly increases afterwards.By decreasing the packing density, we find that the patterns of TimeSOAP signals remain similar for disordered networks, while the patterns of crystalline networks show distinct difference, with lower packing density resulting in high-TimeSOAP particles at smaller normal strain.
To further characterize the particle dynamics, we compute the first derivative of TimeSOAP signal  ̇! 5&65 for each particle i at the normal strain  by [ [1,2] ]: Here, we quantify the fraction of fluctuating particles at each strain by counting the number of particles that exhibit signal derivative values larger than 0.01, that is, a criterion determined as the position between two peaks in the distribution density of signal derivative values (see Fig. S5).The main text has discussed the

S5. Stability of atomistic configurations at finite temperature
We examine herein the structural stability of different typical structures present in this work.Figure S6 shows the fictive temperature of different structures and their particles' mean square displacement (MSD) with time at various temperatures below Tf.The simulation settings are described in the main text.We find that Tf is generally above 0.03 in reduced unit, which is equivalent to 755 K in the silicon energy scale of SW potential [ [6] ].When the system temperature T = 0.03, we find that the particle MSD in different structures enter into their diffusion regime (or subdiffusive regime) but remain a very confined level within ~1 particle diameter displacement after a long simulation time of 1000 in reduced unit, equivalent to ~100 picoseconds in the silicon time scale of SW potential.When the system temperature T < 0.03, the particles   The right panel shows the particle mean square displacement as a function of MD time in logarithmic scale at system temperature T = 0.01, 0.02, and 0.03, respectively.

S6. Simulation reliability of SW potential
To validate the simulation reliability of SW potential, we could examine two extreme disordered networks, i.e., a network with no angular constraint and a strong covalent silicon network.The former network is solely described by the two-body interaction term of SW potential [ [6] ], that is, a generalized format of Lennard-Jones potential suitable for describing short-range two-body interaction.By adding strong angular interaction term of SW potential, the disordered silicon exhibits tetrahedral packing structure similar to that prepared by an accurate EDIP potential as a reference-i.e., a well-established forcefield to describe disordered silicon phase [ [7] ], and their stiffness-density scaling comparison is provided in Fig. S7.Note that both the SW and EIDP potential adopts real units of silicon attributes for fair comparison.The timestep is fixed to 1 fs.The melt quenching simulation is the same as the reduced-unit simulation described in the main text, except that the melt-state temperature is set as 5000 K.We find that the two forcefields offer very close scaling results, suggesting that the SW potential is capable of simulating disordered silicon and, generically, other similar covalent systems consisting of short-range radial and angular interactions.by SW [ [6] ] and EDIP [ [7] ] potential.
Based on the simulation reliability of the two extreme scenarios, we then intentionally select the key parameters that tune bond directionality, i.e., angle value and angle strength, instead of blindly unlocking the variability of all parameters in SW potential.This allows us to filter out and keep away from any unrealistic regions of SW potential, considering that both the angular value and strength parameters show no dependance or impose no influence on any other parameters of SW potential.In other words, the SW

Fig. S1 :
Fig. S1: Initial dataset of forcefield parameters.The red points represent the forcefields

Fig. S2 :
Fig. S2: Topography of two-body excess entropy S2 as a function of forcefield parameters

Fig. S3 :
Fig. S3: Comparison of normal (left) and shear (right) stress shift per atom as a function rate of local environment change for each particle i during tensile deformation.FigureS5shows the first derivative of TimeSOAP signals for silicon-type disordered versus crystalline porous network at different packing densities, wherein the derivative values are normalized by multiplying 0.1%.Indeed, the disordered network has a fraction of highly fluctuating particles at each normal strain, while particles in the crystalline network remains affine deformation with negligible structural fluctuation at initial strains.By decreasing the packing density, crystalline network shows earlier onset of fluctuating particles, while disordered network has more or less similar signal derivative patterns.

Fig. S5 :
Fig. S5: First derivative of TimeSOAP signals for (A) silicon-type disordered versus (B) generally show very limited mobility through confined vibration motions around their initial positions.These results demonstrate the extensive stability of different structures presents in this work at modest finite temperature.

Fig. S7 :
Fig. S7: Comparison of stiffness-density scaling for disordered silicon network generated