# designing_mechanical_metamaterials_by_learning_equivariant_flows__0a519f39.pdf Published as a conference paper at ICLR 2025 DESIGNING MECHANICAL META-MATERIALS BY LEARNING EQUIVARIANT FLOWS Mehran Mirramezani Department of Computer Science Princeton University m.mirramezani@berkeley.edu Anne S. Meeussen & Katia Bertoldi School of Engineering and Applied Sciences Harvard University {meeussen,bertoldi}@seas.harvard.edu Peter Orbanz Gatsby Computational Neuroscience Unit University College London p.orbanz@ucl.ac.uk Ryan P. Adams Department of Computer Science Princeton University rpa@princeton.edu Mechanical meta-materials are porous solids whose geometric structure results in exotic nonlinear mechanical behaviors that are not typically achievable via homogeneous materials. We show how to drastically expand the design space of a class of mechanical meta-materials known as cellular solids, by generalizing beyond translational symmetry of a unit pore cell. This is made possible by transforming a reference geometry according to a divergence free flow that is parameterized by a neural network and equivariant under the relevant symmetry group. We show how to construct flows equivariant to the space groups, despite the fact that these groups are not compact. Coupling this flow with a differentiable nonlinear mechanics simulator allows us to represent a much richer set of cellular solids than was previously possible. These materials can be optimized to exhibit desirable mechanical properties such as negative Poisson s ratios or to match target stress-strain curves. We validate simulated mechanical behaviors of these new designs against fabricated real-world prototypes. We find that designs with higher-order symmetries can exhibit a wider range of behaviors. 1 INTRODUCTION Mechanical meta-materials are engineered porous structures in which geometric features of the pores lead to uncommon mechanical behaviors not achievable from natural materials (Bertoldi et al., 2017) with various applications in soft robotics (Khajehtourian & Kochmann, 2021), biomedical devices (Wang et al., 2023), etc. One promising class of mechanical meta-materials are cellular solids, characterized by an array of pores, that exhibit exotic nonlinear properties such as auxeticity (expanding under tension) (Bertoldi et al., 2010) and reversible shape morphing (Wenz et al., 2021). The nonlinear functionalities of cellular solids can be programmed by: i) optimizing the geometry of the pores, and/or ii) manipulating the arrangement of the pores. Existing approaches to cellular solids design have been focused on the former, optimizing the shape (Overvelde & Bertoldi, 2014; Xue & Mao, 2022; Medina et al., 2023) or topology (Wang et al., 2014; Gao et al., 2019) of the pores of a cellular solid constructed by translational symmetric arrangement of a single pore unit cell. The design of such cellular solids has also begun to attract substantial attention from the machine learning community, e.g., Beatson et al. (2020); Ma et al. (2020); Kollmann et al. (2020); Ha et al. (2023). Figure 1(a) shows examples of such cellular solids designs from Overvelde et al. (2012). However, the space of potential symmetric arrangement of the pores is not limited only to translation, and exploiting other spatial patterns of pores can potentially unlock new cellular solids with unexpected mechanical responses. Representing and optimizing this larger space of cellular solids is challenging, but recent developments in machine learning of crystallographically-invariant functions (Adams & Orbanz, 2023) make it possible to construct unconstrained parameterizations of the symmetric arrangement of the pores that can be used for Published as a conference paper at ICLR 2025 (b) (c) (d) Figure 1: (a) Example of three shape functions s : Ω R2 {0, 1} constructed by translational symmetric arrangement of different unit pore cells, from Overvelde et al. (2012). (b) Real-world experiment demonstrating different nonlinear behaviors of fabricated cellular solids specified by the shape functions in (a), from Overvelde et al. (2012). (c) A vector field defined by a flow equivariant under the space group p4, constructed as in Section 3. (d) A design based on a level set of the p4-invariant function in (c). Note that it is disconnected. shape optimization. In this work, we leverage these models and introduce a novel framework for learning richer classes of cellular solids, in which we not only optimize pore shapes but also explore all possible arrangements of the pores from two-dimensional crystallographic symmetry groups. The approach we take to the problem is to define the geometry implicitly via a neural network flow from a reference configuration. By carefully constructing this flow to preserve group invariance, that enables all 17 crystallographically symmetric arrangement of the pores including translation, simple symmetric cellular structures can be transformed parametrically into complex cellular structures, as illustrated in Figure 3. Another advantage of this approach is that it ensures the solution avoids disconnected regions and associated numerical challenges observed in existing techniques such as level set method (Wang et al., 2014) shown in Figure 1(d) or moving morphable voids method (Du et al., 2024). We additionally show that such neural flows can be made divergence-free, which makes it possible to avoid modifying the total volume of the structure, allowing the optimization to focus purely on geometry as the means of solving the task. The paper is structured as follows. We first introduce cellular solids, framed in terms of a design space that can exhibit symmetries. Section 3 then develops a constrained optimization problem for learning shapes which are symmetric with respect to a space group while achieving a target volume fraction. Our approach uses a novel formalism based on equivariant neural network flows. Next we provide an overview of our implementation, including details of the neural network, the differentiable mechanics solver, and the mechanical material model used. We then demonstrate the application of our proposed approach in designing cellular solids with nontrivial functionalities in both simulation and real experiments. We conclude with related work and a discussion on limitations and future steps. 2 CELLULAR SOLIDS AND SYMMETRIC SHAPES An example of a cellular solid is illustrated in Figure 1(a). It consists of a homogeneous solid material, interrupted by empty regions that are called pores. Once the material is chosen, the solid is entirely characterized by its geometry. This geometry can be specified by a function s : Ω {0, 1}, where Ωis the entire volume occupied by the solid, including pores. This function is called a shape. A value s(x) = 1 specifies there is material at location x, a value s(x) = 0 that x is in a pore. We are particularly interested in cellular solids that are like the example in Figure 1(a) completely determined by a two-dimensional cross-section. We can then reduce the shape to a two-dimensional function s : Ω R2 {0, 1}, where Ωis now the area of the cross-section. The methods developed here can be used to construct shape functions on Ω Rn for both n = 2 and n = 3. Our simulations and real-world experiments assume n = 2. Periodic designs as tilings. Clearly visible in Figure 1(a) is the periodic structure. This periodicity can be formalized as a tiling: If Π Rn is a convex polytope, and T is a group of translations of Rn, then T is said to tile the space with Π if the images ϕΠ of the polytope under elements ϕ T cover the entire space and only their boundaries overlap formally, if ϕ T ϕΠ = Rn and ϕΠ ψΠ = empty set or face of ϕΠ (1) Published as a conference paper at ICLR 2025 for all ϕ, ψ T. If Π is, for example, an axis-aligned square in R2, the group T would consist of all horizontal and vertical shifts by multiples of the edge length of Π. A shape s is periodic if there is a function σ : Π {0, 1} such that s(x) = σ(ϕx) if x ϕΠ for some ϕ T . (2) In words, σ specifies a pattern of pores on a small area Π, and this pattern is then replicated over Ωby shifts. To avoid dealing with the boundaries of Ωexplicitly, we can define s as a function s : R2 {0, 1} on the entire plane, and restrict it to Ωwhere necessary (see Section 4 for details); observe that (2) indeed specifies s on all of R2, if T tiles R2 with Π. We also note that s satisfies (2) if and only if it is invariant under T, s(ϕx) = s(x) for all ϕ T and x R2 , in which case σ is the restriction of s to Π. Expanding the design space. Our point of departure from existing designs is the observation that tilings cannot only be generated by translations. In general, one may substitute T by a group G of isometries of Rn. Such a group is said to tile Rn with a convex polytope Π if (1) holds with G substituted for T. Any group of isometries that tiles Rn with a convex polytope Π is called a space group or crystallographic group on Rn. We are particularly interested in the case n = 2 that are also known as wallpaper groups. In addition to translations, space groups may also contain rotations, reflections, and other isometries. Our approach is to construct cellular solids with shape functions that satisfy the generalized periodicity s(ϕx) = s(x) for all ϕ G and x R2 , (3) for a space group G. We call such a shape symmetric with respect to G. Properties of space groups. The two-dimensional space groups are completely classified: Two space groups G1 and G2 are isomorphic if one can be transformed into the other by an affine deformation of the underlying space, i.e., if there is an affine, order-preserving bijection α : Rn Rn such that ϕ G1 if and only if αϕα 1 G2. Up to isomorphism, there are only finitely many space groups for each dimension n. For n = 2, there are 17 such groups. These are illustrated in Appendix A.1 Some properties of space groups relevant in the following are: (i) Each space group G is countably infinite. (ii) Since each element ϕ of G is an isometry, it is of the form ϕx = Aϕx + bϕ, where Aϕ is an orthogonal matrix and bϕ Rn. (iii) G contains a countably infinite subgroup of pure translations, namely TG := {ϕ G | Aϕ = identity} = {ϕ G | ϕ(x) = x + bϕ for some bϕ Rn} . This group is always infinite, and generated by n linearly independent shifts: There exist n linearly independent vectors b1, . . . , bn Rn such that ϕ TG ϕ(x) = x + Xn i=1 ci(ϕ)bi for some c1(ϕ), . . . , cn(ϕ) Z . We refer to Vinberg & Shvartsman (1993); Ratcliffe (2006) for more on the geometry of such groups, and Hahn et al. (1983) for their crystallographic properties. 3 SYMMETRIC SHAPES FROM SYMMETRY-PRESERVING FLOWS We select a class of symmetries by fixing a space group G. Our goal is to solve two problems: I. Constructively define a set S = {sθ | θ Θ} of shapes that satisfy (3), indexed by parameters in some parameter space Θ. II. Given a loss function ℓ: shape R, solve θ := arg minθ Θ ℓ(sθ). The loss represents how well the shape achieves some physical property, such as the Poisson s ratio illustrated in Figure 5. To solve (I), we proceed as follows: We start with a symmetric reference shape s0. This is a G-invariant function s0 : Rn {0, 1} that satisfies (3). Published as a conference paper at ICLR 2025 Figure 2: Examples of the reference shape s0 for different space groups. The labels p1, pg, etc. refer to the crystallographic naming standard for such groups see Appendix A.1. We then construct a class {Fθ|θ Θ} of flows that are symmetryand volume-preserving. Each shape sθ is defined as the deformation of s0 by Fθ. That is, we fix a terminal time tmax > 0, and define sθ(x) := s0(Fθ(x, tmax)) . These functions are shapes on Rn that satisfy (3), and can be restricted to shapes on Ω. The remainder of this section explains this construction in detail. This solution of (I) can then be combined with a standard differentiable mechanics simulator to solve (II), as described in Section 4. 3.1 SYMMETRYAND VOLUME-PRESERVING FLOWS A flow is a function F : Rn R+ Rn that satisfies F(x, 0) = x and F(x, s + t) = F(F(x, s), t) for x Rn and s, t R+ . For our purposes, we can restrict the second argument of F to an interval I = [0, tmax]. We call a flow symmetry-preserving for a space group G if it is G-equivariant pointwise in time, ϕ F(x, t) = F(ϕx, t) for all ϕ G and (x, t) Rn I . If that is the case, then sθ is G-invariant, since s0 is G-invariant and so sθ(ϕx) = s0(Fθ(ϕx, tmax)) = s0(ϕFθ(x, tmax)) = s0(Fθ(x, tmax)) = sθ(x) . We also require F to be volume-preserving. A continuous function g : Rn Rn is volumepreserving if vol(g(R)) = vol(R) for every (measurable) set R Rn . For a flow F, this means vol(F(R, t)) = vol(R) for every (measurable) set R Rn and t I . This ensures that the ratio of volume occupied by pores does not change if a shape is deformed using F. It is well known that a flow F can be constructed as the solution of a differential equation. If H : Rn I Rn is a (sufficiently smooth) function, the constrained equation d dt F(x0, t) = H(F(x0, t), t) subject to F(x0, 0) = x0 for t I, x0 Rn (4) has a unique solution F, and this solution is a flow. Each smooth function H hence determines a flow F, and we can specify a parameterized class of flows Fθ by specifying a parameterized class of smooth functions Hθ. Theorem 2 below shows that the flow F is volumeand symmetry-preserving if H is volume-preserving and satisfies H(ϕx, t) = AϕH(x, t) for ϕ G, x Rn, t I . (5) To construct such a function H, we define a symmetrization operator Γ and an operator Λ that turns functions into volume-preserving functions. 3.2 VOLUME PRESERVATION To define Λ, we can draw on existing work: A continuously differentiable function g is volumepreserving if and only if it is divergence-free, (div g)(x) = ( g)(x) = 0 for all x Rn . Published as a conference paper at ICLR 2025 It was noted in Richter-Powell et al. (2022) that, if g : Rn Rn(n 1)/2 is a function that is twice continuously differentiable then 0 2g1(x) 3g2(x) ngn 1(x) 1g1(x) 0 3gn(x) ng2n 3(x) 1g2(x) 2gn(x) 0 ng3n 6(x) ... ... ... ... ... 1gn 1(x) 2g2n 3(x) 3g3n 6(x) 0 1 1 1 ... 1 is a divergence-free function Λg : Rn Rn. 3.3 SYMMETRIZATION For finite groups, invariant or equivariant functions are typically constructed by the summation trick , which starts with a function f and sums over all compositions f ϕ with group elements. This is not possible for space groups, since these groups are countably infinite. We can, however, decompose the symmetrization into two steps. This decomposition is based on the following property. Theorem 1. Let G be a space group on Rn. For each ϕ G, denote by ci(ϕ) the linear coefficient of bϕ with respect to the shift basis vector bi, that is, ϕ(x) = Aϕx + P i n ci(ϕ)bi. Then b G := {ϕ G | c1(ϕ), . . . , cn(ϕ) [0, 1)} is a finite subset of G. For each ϕ G, there are unique elements ˆϕ b G and τϕ TG such that ϕ = ˆϕ + τϕ. If f : Rn Rn is invariant under the shift subgroup TG of G, the function (Γf)(x) := 1 ϕ b G A 1 ϕ f(ϕx) satisfies Aϕ(Γf) = (Γf) ϕ for each ϕ G. Proof. See Appendix A.2.1. To construct the TG-invariant function required by the result, we use maximal invariants: Define an n n-matrix B as B := matrix whose columns are the lattice basis vectors b1, . . . , bn . This matrix describes a change of basis from the cartesian basis to b1, . . . , bn. The set {B 1ϕ | ϕ TG} therefore consists of all shifts by vectors with integer entries, or in other words, B 1TG = Zn. Define a function ρ : Rn R2n as ρ(x) := ρ(B 1x) where ρ(u) := (cos(2πu1), sin(2πu1), . . . , cos(2πun), sin(2πun))T . This function maps Rn onto the 2n-dimensional torus. It is well known that ρ is a maximal invariant for the shift group Zn. This means that a function g on Rn is Zn-invariant if and only if it can be represented as g = g ρ for some function g on R2n. It follows that ρ is a maximal invariant for TG. 3.4 FLOW CONSTRUCTION We can now combine Γ, Λ and ρ to construct a divergence-free function H that satisfies (5) from a smooth function h, as the next result shows. We can therefore use a neural network with parameter vector θ to represent a class {hθ|θ Θ} of such functions. For each θ, we then obtain a unique flow Fθ that preserves volume and symmetry. Theorem 2. Let G be a space group on Rn. Choose a function h : R2n I Rn(n 1)/2 that is Lipschitz-continuous, and (k + 1)-times continuously differentiable on I, for some k N. Then H(x, t) := Γ(Λ(h(ρ(x), t))) for (x, t) Rn I is a divergence-free function Rn I Rn that satisfies (5). For this function H, the differential equation (4) has a unique solution F. The function F is a flow, is volume-preserving, and is symmetrypreserving for G. It is continuous in its first argument, and k times continuously differentiable in the second. Proof. See Appendix A.2.2. Published as a conference paper at ICLR 2025 differentiable mechanics solver symmetry/volume-preserving initial model cellular meta-material initial/target shape flowed model w/ sym. final design neural network flow( ) Figure 3: An overview of our modeling framework implemented in JAX. An equivariant neural network parameterized by θ works in tandem with a differentiable mechanics solver to flow an initial geometry to learn cellular solids with desired functionalities from a rich set of volumeand symmetry-preserving geometries. The upper left shape is s0 and the upper right shape is sθ. 4 IMPLEMENTATION To design cellular solids, we combine the shape representation defined above with a mechanics simulator. We proceed as follows: 1) Start with a reference shape s0. This is generated by selecting a space group G and a tiling convex polytope Π for this group. We then define a set of pores on Π and tile with the images ϕΠ under group elements ϕ G such that the entire region Ωis covered. 2) The region {x : s0(x) = 1} is represented by a triangular mesh using pygmsh, an open-source python package for geometry and mesh processing. This yields a discretization of the shape in terms of N points in R2. 3) The image of the mesh points under a flow Fθ is constructed as in Section 3 to define a symmetric shape covering the same volume as s0 for each parameter value θ. 4) Using a mechanics simulator, we simulate the behavior of sθ under physical effects, such as lateral deformation. We define a cost function ℓthat measures the response to a given physical effect, in such a way that small values of ℓencode desired behavior. 5) We minimize ℓwith respect to θ using gradient descent, as described in more detail below. The overall setup is illustrated in Figure 3. Differentiable mechanics simulator. To solve the optimization problem, we must backpropagate gradients through both the mechanics solver and the flow representation. This involves solving a static nonlinear elasticity problem, described in Appendix A.3. To do so, we use an end-to-end differentiable continuum mechanics solver (Oktay, 2024). This solver is based on a variant of the finite element method called isogeometric analysis (IGA) (Hughes et al., 2005; Hughes, 2012). IGA represents both geometry and solutions in the same smooth B-spline basis functions. This allows us to directly construct differentiable maps from (i) geometry parameters to the solution, and (ii) to a specified loss function defined by the solution. The optimizer is implemented in JAX (Bradbury et al., 2018), and uses automatic differentiation and adjoint methods. Neural network ansatz. The function hθ in Section 3 is represented by a fully-connected neural network with two hidden layers of size 10, with tanh nonlinearity. To optimize these parameters, we use ADAM, with a learning rate of 0.001. Each optimization is performed several times with different initialization of the neural network parameters. Published as a conference paper at ICLR 2025 Boundary conditions. In Section 3, shapes are constructed as functions on Rn, here for n = 2. To represent cellular solids, the function must be restricted to the finite volume Ω, and we must ensure this restriction respects boundary conditions. To ensure preservation of the boundaries regardless of the parameters θ, we use an envelope function which takes H(x, t) 0 at the edges of Ω. 5 APPLICATIONS Our experiments design novel cellular solids with two types of nonlinear behaviors: 1) nontrivial force-displacement responses, and 2) effective Poisson s ratios1 under tension and compression. These types of behaviors are important in engineering applications such as soft robotics and programmable actuations. Here we explore examples of such behavior that would be challenging or impossible to achieve with homogeneous materials. In uniaxial loading we apply a displacement uniformly to the top edge of the material while fixing the bottom edge; the left and right boundaries are free to move. The force-displacement response is usually described by a nominal stress-strain curve, in which nominal stresses (S) in our experiments are computed as the resulting vertical reaction forces per thickness (in the bottom boundary) normalized by the width of the cellular solids. Strain (ϵ) values are the applied displacements per meta-material height (H) that in this study we use a displacement of equal to 10% of the height of the structure resulting in the final applied strain of 0.1. The effective Poisson s ratios νef are approximated using horizontal displacement in a region of size H/2 in the middle of left and right boundaries (Ωl and Ωr, see Figure 5 top left) as in Medina et al. (2023): Ωl uxds 1 |Ωr| Ωr uxds . (7) When targeting force-displacement responses we use curves that are obtained at equally spaced displacement loading increments applied at the top edge to minimize a loss L = P i[S(ϵi) St(ϵi)]2 with respect to a target St response. We also design cellular solids with target Poisson s ratios at ϵ = 0.1. The mechanical behavior of cellular solids in this study is captured with a nearly incompressible neo-Hookean model with strain energy density function given in (13) that has material properties of Young s modulus E = 50 MPa and Poisson s ratio ν = 0.46. In our experiments, cellular solids are of size 7 7 with initial material volume fraction of 50% that is preserved during designs. In all force-displacement experiments, neural network parameters are optimized until a loss value less than 0.0001 is achieved, where for Poisson s ratio designs a loss value less than 0.01 was used for stopping simulations. Using a single NVIDIA Ge Force RTX 2080 Ti GPU, each optimization step (which requires forward simulation and gradient computation) requires approximately 60 seconds; there are at least 25,000 degrees of freedom in all mechanics models that found satisfactory through a mesh independence analysis described in Appendix A.4, such that simulation results were not affected by insufficient discretization resolution. 5.1 DESIGNING UNDER UNIAXIAL TENSION Cellular solid structures built from hyperelastic materials usually exhibit nonlinear force-displacement responses. Here, we aim at proposing cellular solids with linear responses under uniaxially pulling experiments i.e., St(ϵ) = Cϵ such that St(ϵ = 0.1) = βS0 where β changes from 0.1 to 1.5 with an interval of 0.1 for a given normalized nominal stress S0. We find that for all considered 17 symmetry groups, p2gg group indicated the best performance by being able to learn cellular solids with linear force-displacement responses for all β values except β = 1.5. The learned meta-materials with their linear responses are depicted in Figure 4 for β = 0.1, 0.5 and 1.4. While we are able to achieve a wide range of linear responses with p2gg group, with p1 translational symmetry we can design meta-materials that have mechanical responses with β up to 1.1 (see Appendix A.5). 1Informally, Poisson s ratio is how much a material expands laterally when squeezed from the top, or contracts when pulled. Almost all materials have a positive Poisson s ratio. Published as a conference paper at ICLR 2025 Figure 4: Undeformed configurations (i.e., ϵ = 0) of three cellular solids designs with 50% volume fraction and pore shapes respecting p2gg symmetry group, and their linear stress-strain responses with the corresponding target curves during a uniaxial tension of the top edge up to ϵ = 0.1. Next we consider designing auxetic meta-materials with negative Poisson s ratios. We explored cellular solids that can exhibit a large νef = 0.5 during tension. With p1 symmetry, the best design can only achieve νef = 0.05, whereas pg symmetry can improve this to νef = 0.12 (Figure 5). Our framework was able to learn a cellular solid with p4 symmetry that has νef = 0.45 that is very close to the target negative Poisson s ratio. The eight best-performing cellular solids designs among the 17 group-constrained designs in achieving νef = 0.5 are shown in Table 1. This indicates the benefit of exploring a larger space of cellular solids than those exhibiting simple translational symmetries. In this case, six symmetry groups have better performances than p1 symmetry. All final eight designs are depicted in Appendix A.6. Real-world experiments: To demonstrate transfer to the real world, we manufactured the final designs with pg and p4 symmetry groups and quantitatively verified their behavior under a uniaxial tension test (detailed in Appendix A.7). As predicted by our simulations in Figure 5 (second column) and confirmed with experimental observations (third column), the two proposed designs have negative Poisson s ratios. The experimental measurements νex = 0.14 and νex = 0.49 under pulling for pg and p4 models, respectively, are in strong agreement with simulation results. The video footage of the pg cellular solid sample during pulling experiment with the corresponding simulation result are included in the supplementary materials. 5.2 DESIGNING UNDER UNIAXIAL COMPRESSION The design of cellular solids under uniaxial compression is a more challenging task because of nonlinear buckling instabilities observed in such structures (Overvelde & Bertoldi, 2014). In such cases, the nominal stress-strain responses are characterized by a critical strain value (ϵcr) at which mechanical responses vary. Here, we seek to design cellular solids that behave linearly before an ϵcr while having a constant nominal stress as the applied strain increases to the final value of 0.1: S(ϵ) = Cϵ if ϵ ϵcr, Cϵcr if ϵcr ϵ 0.1. (8) Figure 6 illustrates mechanical responses of three cellular solid designs under uniaxial compression. We were able to design materials with various compression responses by controlling the plateau forces for a fixed ϵcr (designs 3 and 4) and tuning both the plateau force and the location of ϵcr (design 2) using p1 symmetry, which found to be a challenging design factor in previous works (Medina et al., 2023). Interestingly, the only symmetry group that enabled such designs was p1. Examples of designs with other symmetries that force-displacement responses do not match target values are depicted in Appendix A.8. Next, we focus on tailoring cellular solids with negative Poisson s ratios under compression. Here we aim to identify structures with νef = 0.2, i.e., requiring the material to shrink. To ensure this Table 1: Poisson s ratios of cellular solids from eight symmetry groups with top performances when aiming for νef = 0.5. groups p1 p2 pg p2mg p2gg p4 p3 p6 νef 0.05 0.04 0.12 0.40 0.10 0.45 0.11 0.14 Published as a conference paper at ICLR 2025 Figure 5: Undeformed configurations (i.e., ϵ = 0) of two cellular solids designs that have negative Poisson s ratios with 50% volume fraction and pore shapes respecting pg and p4 symmetry groups (first column), and their deformed configurations from simulations during a uniaxial pulling of the top edge up to ϵ = 0.1 (second column). Experimental realization of the same cellular solids under pulling confirms negative Poisson s ratios for both designs. The experimental measurements are also in strong numerical agreement with simulation results (third and fourth columns). Scale bars: 5cm. design criterion, we customized the loss function to enforce points on Ωl and Ωr edges to have mean positive and negative displacements, respectively. We were able to optimize a p6 and pg symmetric cellular solids shown in Figure 7 with νef = 0.15 and νef = 0.2 under compression, respectively, that they also shrink. Finally, we examine the ability of our framework in designing cellular meta-materials with a more complicated mechanical response: achieving zero effective Poisson s ratios under both uniaxial compression and tension via having zero displacements on Ωl and Ωr edges. Figure 8 shows the resulting successful structure that leverages p6 symmetry group, indicating zero displacement at both left and right edges for either pushing or pulling of 10% of the cellular solid height. 6 RELATED WORK Existing cellular solids design are mainly focused on structures constructed from translational symmetry of an optimized pore unit cell. For example, Wang et al. (2014) utilized a level set method Figure 6: Undeformed configurations (i.e., ϵ = 0) of three cellular solids designs with 50% volume fraction and pore shapes respecting p1 symmetry group, and their stress-strain responses with the corresponding target curves during a uniaxial compression of the top edge up to ϵ = 0.1. Published as a conference paper at ICLR 2025 Figure 7: Undeformed and deformed configurations of a p6 and pg symmetric cellular solids designs that have negative Poisson s ratios of 0.15 and 0.2, respectively, under uniaxial compression. Figure 8: Undeformed and deformed configurations of a cellular solid design with 50% volume fraction and pore shapes respecting p6 symmetry group indicating zero effective Poisson s ratios under both 10% uniaxial pushing and pulling, with loss values during optimization. to design cellular meta-materials via unit cells topology optimization. Xue & Mao (2022) proposed an FEM-based mapped shape optimization technique to design periodic meta-materials. Medina et al. (2023) employed a higher order moving-mesh method to optimize the pore shapes of cellular solids. Du et al. (2024) optimized pore structures through describing them by a moving morphable voids method. Machine learning methods are becoming an emerging tool for both accelerating the design process and discovering new meta-materials beyond intuitions. For example, Xue et al. (2020b) introduced a data-driven neural network homogenization approach for cellular meta-materials design. Ma et al. (2020) applied machine learning tools for characterization of non-uniform arrangement of cells in meta-materials. Ha et al. (2023) leveraged generative machine learning techniques for inverse design of meta-materials with target stress-strain responses. In this study, we propose an equivariant neural network flow to learn volume-preserving cellular solids via performing shape optimization of the pores while exploring all 2D crystallographic symmetry groups to find the optimal arrangement of the pores. 7 LIMITATIONS AND FUTURE WORK We introduce a machine learning framework that enables an efficient approach for inverse design of two-dimensional cellular solids with various nontrivial mechanical behaviors beyond the reach of existing methods. Scaling up to three-dimensional cellular solids design while considering crystallographic symmetric arrangements could potentially unlock uncommon mechanical functionalities not observed in two dimensional meta-materials. But this requires efficient and faster solvers that our framework suffers from at large number of degrees of freedom in the mechanics modeling. A future work could focus on speeding up the simulations leveraging learned surrogate models (Sun et al., 2021), neural network-based order reduction techniques (Beatson et al., 2020), and also amortized optimization methods (Xue et al., 2020a). Since all crystallographic symmetry groups include translational symmetry, further computational cost reduction can also be achieved by simulating a single unit cell with periodic boundary conditions instead of modeling the entire cellular solid (Mizzi et al., 2021). Published as a conference paper at ICLR 2025 ACKNOWLEDGMENTS We would like to thank Deniz Oktay for his computational mechanics insights. This work was partially supported by NSF grants IIS-2007278 and OAC-2118201, NSF under grant number 2118201, the NSF under grant number 2127309 to the Computing Research Association for the CIFellows 2021 Project. PO is supported by the Gatsby Charitable Foundation. Ryan P. Adams and Peter Orbanz. Representing and learning functions invariant under crystallographic groups. ar Xiv preprint ar Xiv:2306.05261, 2023. Ludwig Arnold. Random Dynamical Systems. Springer, 1998. Alex Beatson, Jordan T. Ash, Geoffrey Roeder, Tianju Xue, and Ryan P. Adams. Learning composable energy surrogates for PDE order reduction. In Proceedings of the 34th International Conference on Neural Information Processing Systems, pp. 338 348, 2020. Katia Bertoldi, Pedro M. Reis, Stephen Willshaw, and Tom Mullin. Negative Poisson s ratio behavior induced by an elastic instability. Advanced Materials, 22(3):361 366, 2010. Katia Bertoldi, Vincenzo Vitelli, Johan Christensen, and Martin van Hecke. Flexible mechanical metamaterials. Nature Reviews Materials, 2(11):17066, 2017. James Bradbury, Roy Frostig, Peter Hawkins, Matthew J. Johnson, Chris Leary, Dougal Maclaurin, and Skye Wanderman-Milne. JAX: composable transformations of python+numpy programs. 2018. URL https://github.com/google/jax. Peter Deuflhard and Folkmar Bornemann. Scientific Computing with Ordinary Differential Equations. Springer, 2002. Zongliang Du, Tanghuai Bian, Xiaoqiang Ren, Yibo Jia, Shan Tang, Tianchen Cui, and Xu Guo. Inverse design of mechanical metamaterial achieving a prescribed constitutive curve. Theoretical and Applied Mechanics Letters, 14(1):100486, 2024. Jie Gao, Huipeng Xue, Liang Gao, and Zhen Luo. Topology optimization for auxetic metamaterials based on isogeometric analysis. Computer Methods in Applied Mechanics and Engineering, 352: 211 236, 2019. Chan Soo Ha, Desheng Yao, Zhenpeng Xu, Chenang Liu, Han Liu, Daniel Elkins, Matthew Kile, Vikram Deshpande, Zhenyu Kong, Mathieu Bauchy, and Xiaoyu (Rayne) Zheng. Rapid inverse design of metamaterials based on prescribed mechanical behavior through machine learning. Nature Communications, 14(1):5765, 2023. Theo Hahn, Uri Shmueli, and James C. W. Arthur. International tables for crystallography, volume 1. Reidel, 1983. Thomas J.R. Hughes. The Finite Element Method: Linear Static and Dynamic Finite Element Analysis. Courier Corporation, 2012. Thomas J.R. Hughes, Austin Cottrell, and Yuri Bazilevs. Isogeometric analysis: Cad, finite elements, nurbs, exact geometry and mesh refinement. Computer Methods in Applied Mechanics and Engineering, 194(39):4135 4195, 2005. Romik Khajehtourian and Dennis M. Kochmann. Soft adaptive mechanical metamaterials. Frontiers in Robotics and AI, 8, 2021. Hunter T. Kollmann, Diab W. Abueidda, Seid Koric, Erman Guleryuz, and Nahil A. Sobh. Deep learning for topology optimization of 2D metamaterials. Materials Design, 196:109098, 2020. Chunping Ma, Zhiwei Zhang, Benjamin Luce, Simon Pusateri, Binglin Xie, Mohammad H. Rafiei, and Nan Hu. Accelerated design and characterization of non-uniform cellular materials via a machine-learning based framework. NPJ Computational Materials, 6(1):40, 2020. Published as a conference paper at ICLR 2025 Eder Medina, Chris Rycroft, and Katia Bertoldi. Nonlinear shape optimization of flexible mechanical metamaterials. Extreme Mechanics Letters, 61:102015, 2023. Luke Mizzi, Daphne Attard, Ruben Gatt, Krzysztof K. Dudek, Brian Ellul, and Joseph N. Grima. Implementation of periodic boundary conditions for loading of mechanical metamaterials and other complex geometric microstructures using finite element analysis. Engineering with Computers, 37 (3):1765 1779, 2021. Deniz Oktay. Translating between Scientific Computing and Machine Learning with Automatic Differentiation. Ph D thesis, Princeton University, 2024. Johannes T. B. Overvelde and Katia Bertoldi. Relating pore shape to the non-linear response of periodic elastomeric structures. Journal of the Mechanics and Physics of Solids, 64:351 366, 2014. Johannes T. B. Overvelde, Sicong Shan, and Katia Bertoldi. Compaction through buckling in 2D periodic, soft and porous structures: effect of pore shape. Advanced Materials, 2012. John G. Ratcliffe. Foundations of Hyperbolic Manifolds. Springer, 2nd edition, 2006. Jack Richter-Powell, Yaron Lipman, and Ricky TQ Chen. Neural conservation laws: A divergencefree perspective. Advances in Neural Information Processing Systems, 35:38075 38088, 2022. Xingyuan Sun, Tianju Xue, Szymon Rusinkiewicz, and Ryan P. Adams. Amortized synthesis of constrained configurations using a differentiable surrogate. Advances in Neural Information Processing Systems, 34:18891 18906, 2021. Efim B. Vinberg and Oleg V. Shvartsman. Discrete groups of motions of spaces of constant curvature. In E. B. Vinberg (ed.), Geometry II. Spaces of Constant Curvature. Springer, 1993. Hao Wang, Yongtao Lyu, Sergei Bosiakov, Hanxing Zhu, and Yuanfei Ren. A review on the mechanical metamaterials and their applications in the field of biomedical engineering. Frontiers in Materials, 10, 2023. Yiqiang Wang, Zhen Luo, Nong Zhang, and Zhan Kang. Topological shape optimization of microstructural metamaterials using a level set method. Computational Materials Science, 87: 178 186, 2014. Franziska Wenz, Ingo Schmidt, Alexander Leichner, Tobias Lichti, Sascha Baumann, Heiko Andrae, and Christoph Eberl. Designing shape morphing behavior through local programming of mechanical metamaterials. Advanced Materials, 33(37):2008617, 2021. Tianju Xue and Sheng Mao. Mapped shape optimization method for the rational design of cellular mechanical metamaterials under large deformation. International Journal for Numerical Methods in Engineering, 123(10):2357 2380, 2022. Tianju Xue, Alex Beatson, Sigrid Adriaenssens, and Ryan P. Adams. Amortized finite element analysis for fast PDE-constrained optimization. In Proceedings of the 37th International Conference on Machine Learning, volume 119 of Proceedings of Machine Learning Research, pp. 10638 10647, 13 18 Jul 2020a. Tianju Xue, Alex Beatson, Maurizio Chiaramonte, Geoffrey Roeder, Jordan T. Ash, Yigit Menguc, Sigrid Adriaenssens, Ryan P. Adams, and Sheng Mao. A data-driven computational scheme for the nonlinear mechanical properties of cellular mechanical metamaterials under large deformation. Soft Matter, 16(32):7524 7534, 2020b. Published as a conference paper at ICLR 2025 A.1 THE 17 WALLPAPER GROUPS The following figure illustrates the 17 distinct (up to isomorphy) space groups on R2, also known as the wallpaper groups. In each figure, the region marked red corresponds to the polytope Π in (1). The letter F is inscribed to mark orientation. Figures from Adams & Orbanz (2023). p1 p2 pm pg cm p2mm p2mg p2gg c2mm p4 p4mm p4gm p3 p3m1 p31m p6 p6mm A.2.1 PROOF OF THEOREM 1 This proof draws on one of the fundamental properties of space groups, namely that, if G is a space group and TG its subgroup of pure translations, the quotient group G/TG is finite. See for example (Vinberg & Shvartsman, 1993). Step 1. We first show b G is finite. Consider the relation on G defined as ϕ ψ : ϕ = ψ + τ for some τ TG . Since TG is an equivalence relation, and by definition of quotients, the equivalence classes of correspond to the elements of the quotient group G/TG. Since G is a space groups, this implies has a finite number of equivalence classes. Recall that each element ϕ of G is of the form ϕ(x) = Aϕx + bϕ = Aϕx + X i n ci(ϕ)bi , where c(ϕ) = (c1(ϕ), . . . , cn(ϕ)) is a vector in Rn. It follows that there are unique vectors ˆc(ϕ) [0, 1)n and c(ϕ) Zn such that c(ϕ) = ˆc(ϕ) + c(ϕ). The transformations ˆϕ(x) := Aϕx + X i n ˆci(ϕ)bi and τϕ(x) := x + X i n ci(ϕ)bi are hence indeed elements ˆϕ b G and τϕ TG that satisfy ϕ = ˆϕ + τϕ, and are the only such elements. It follows that each equivalence class of contains exactly one element of b G, which shows the set b G is indeed finite. Step 2. To verify the properties of Γ, we must first establish two simple properties of b G, namely (i) d Gψ = b G and (ii) ˆϕψ(x) = c ϕψ(x) + bτ for some τ T0 . (9) Here, Gψ is the set obtained by composing each element of G with ψ from the right, and d Gψ denotes application of the definition of b G to Gψ. Since G is a group and ψ one of its elements, we have Published as a conference paper at ICLR 2025 Gψ = G, so property (i) holds. To verify (ii), observe that the hat operation ˆϕ preserves the equivalence class of ϕ: Since ϕ = ˆϕ + τϕ and τϕ TG, we have ˆϕ ϕ. If ψ is a further group element, we also have ˆϕψ(x) = ϕ(ψx) bτϕ and therefore ˆϕψ ϕψ. In summary, we have ˆϕψ ϕψ c ϕψ which shows (ii) Step 3. Now consider a TG-invariant function f. For any ψ G, we have ˆϕ b G A 1 ϕ f(ˆϕψx) = X π b Gψ A 1 πψ 1f(πx) substitute π := ˆϕψ π b Gψ A 1 πψ 1f(bπx + bτ) by (9i), for some τ T0 π b Gψ A 1 πψ 1f(bπx) since h is T0-invariant π b Gψ A 1 π f(bπx) since A 1 πψ 1 = A 1 ψ 1A 1 π = AψA 1 π bπ b G A 1 π f(bπx) by (9ii) . It follows that (Γf)(ψx) = X ˆϕ b G A 1 ϕ f(ˆϕψx) = Aψ X bϕ b G A 1 ϕ f(bϕx) = Aψ(Γf)(x) for any ψ G. A.2.2 PROOF OF THEOREM 2 The existence and uniqueness of solutions of (4) is given by the Picard-Lindelöf theorem (e.g. Deuflhard & Bornemann, 2002). Here is a statement that suffices for our purposes: Fact 1 (Version of the Picard-Lindelöf theorem). Fix n, k N and a bounded interval I := [0, tmax]. Let H : Rn I Rd be a Lipschitz function that is k 1 times continuously differentiable in its second argument. Then F(x0, t) := x0 + Z t 0 H(F(x0, s), s)ds defines a function F : Rn I Rd that is continuous in its first argument, k + 1 times continuously differentiable in the second, and is the unique solution of (4). We break the main ingredients of the proof of Theorem 2 up into three lemmas. First, we show that a function H that satisfies (5) yields flows with the symmetry properties we require. Lemma 1. Let G be a group of isometries of Rn, and fix some tmax > 0. Let H : Rn I Rn be a function is Lipschitz, and k times continuously differentiable in its second argument, for k 1. Then there is a unique function F : Rn [0, tmax] Rn that satisfies (4) for all (x, t) in its domain, and this function satisfies F(ϕx, t) = ϕF(x, t) for all ϕ G, all x Rn and 0 t tmax (10) if and only if H satisfies (5). The solution F is continuous in its first argument, and (k + 1) times continuously differentiable in the second. Proof. We first construct the solution F : Rn I Rn on the larger interval I = [ 1, tmax]. (This is to ensure that 0 is in the interior of I, which we need for step 3 below.) The restriction of F to Rn [0, tmax] is then solution in the statement of the result. Step 1. Existence, uniqueness, and smoothness of F on Rn I then follow immediately from Fact 1, where we set tmin = 1 and t0 = 0. Fact 1 also shows that F solves d dt F(x, t) = H(F(x, t), t) with F(x, 0) = x . (11) What we have to show is that F is equivariant if and only if H satisfies (5). Published as a conference paper at ICLR 2025 Step 2. Suppose H satisfies (5). Write F ϕ(x, t) := ϕ 1F(ϕx, t). Thus, F is equivariant iff F ϕ = F for all ϕ. d dt F ϕ(x, t) = d dt(Aϕ 1F(ϕx, t) + bϕ 1) since ϕ 1x = Aϕ 1x + bϕ 1 dt F(ϕx, t) Aϕ 1 and bϕ 1 do not depend on t = Aϕ 1H(F(ϕx, t), t) by (11) = H(ϕ 1F(ϕx, t), t) by (5) = H(F ϕ(x, t), t) definition of F ϕ. Substituting the initial condition into the definition of F ϕ shows F ϕ(x, 0) = ϕ 1(F(ϕx, 0)) = ϕ 1ϕx = x Thus, both F and F ϕ solve the differential equation (11) for the initial condition x. Since the solution is unique by Fact 1, it follows that F(x, t) = F ϕ(x, t) for all t. In summary, (5) implies equivariance of F. Step 3. Conversely, assume F is equivariant, so F ϕ = F. Since H is continuously differentiable, Fact 1 shows t 7 F(x, t) is twice continuously differentiable. Its Taylor expansion around any t in the interior ( 1, tmax) is hence F(x, t + ϵ) = F(x, t) + d dt F(x, t) ϵ + o(ϵ2) = F(x, t) + H(F(x, t), t) ϵ + o(ϵ2) for each x Rn and ϵ > 0 . (12) If we set t = 0 and use the initial condition F(x, 0) = x, we obtain ϕF(x, ϵ) = Aϕ(x + H(x, 0) ϵ + o(ϵ2)) + bϕ = ϕx + AϕH(x, 0) ϵ + o(ϵ2) , where we note that Aϕo(ϵ2) = o(ϵ2) since Aϕ is an isometry. Substituting ϕx into (12) shows F(ϕx, ϵ) = ϕx + H(ϕx, 0) + o(ϵ2) Since F ϕ = F, we hence have ϕx + H(ϕx, 0) ϵ + o(ϵ2) = ϕx + AϕH(x, 0) ϵ + o(ϵ2) and therefore H(ϕx, 0) = AϕH(x, 0) + o(ϵ) Since that is true for all ϵ, it follows that H(ϕx, 0) = AϕH(x, 0). In summary, equivariance of F implies that (5) holds at t = 0. Since F is a flow, the flow F t(x, s) := F(x, t + s) satisfies (4) for the function Ht(x, s) := H(x, t + s) and the initial value xt 0 := F(x, t). If F is equivariant, so is F t. We can hence apply the same argument at s = 0, which shows H(ϕx, t) = AϕH(x, t). Given a differentiable function f : Rn Rm, we denote by Df its differential. If m = n, this is the Jacobian matrix, and if m = 1, the transpose (Df)T is the gradient of f. Lemma 2. For any differentiable function g : R2n Rm and any m N, the differential D(g ρ) is TG-invariant. In particular, Dρ is TG-invariant. Proof. Since the Jacobian matrix of ρ has entries (Dρ(x))ij = 2πρ2i 1(x) if j = 2i 1 2πρ2i(x) if j = 2i 0 otherwise it can be written as a function Dρ = M ρ of ρ, where M : R2n Rn 2n. That shows D(g ρ)(x) = (Dg)(ρ(x)) (Dρ)(x) = (Dg M)(ρ(x)) , and hence D(g ρ) ϕ = (Dg M) ρ ϕ = (Dg M) ϕ. Published as a conference paper at ICLR 2025 Our final auxiliary result shows that the symmetrization operator Γ does not increase divergence. For Theorem 2, this means that if a function f is divergence-free, then so is Γf. Lemma 3 (Γ does not increase divergence). If G is a space group on Rn, and f : Rn Rn is a continuously differentiable function, then div(Γf) sup divf sup . Proof. Since the divergence is the trace of the differential, and both trace and differential are linear, div(Γf)(x) = tr D 1 ϕ ˆG A 1 ϕ f(ϕx) = 1 ϕ ˆG tr D A 1 ϕ f(ϕx) Each summand has differential D(A 1 ϕ f ϕ)(x) = A 1 ϕ (Df)(ϕx)(Dϕ)(x) = A 1 ϕ (Df)(ϕx)Aϕ Since the trace is linear and orthogonally invariant, it follows that | tr D(A 1 ϕ f ϕ)(x)| = | tr((Df)(ϕx))| = |div(f)(ϕx)| div(f) sup The result then follows via the triangle inequality, div(Γf)(x) sup = 1 ϕ ˆG divf sup = divf sup , which is all we had to show. Proof of Theorem 2. Since ρ is TG-invariant and infinitely often differentiable, the function h ρ is TG-invariant and (k + 1) times continuously differentiable. The function Λ(h ρ), by definition of Λ, a function of the gradient (h ρ). Since this gradient is TG-invariant by Lemma 2, Λ(h ρ) is TGinvariant. By Theorem 1, ΓΛ(h ρ) satisfies (5). It is also divergence-free, since Λ(h ρ) is divergencefree and Γ does not increase divergence, by Lemma 3. We can therefore choose H := ΓΛ(h ρ) in Lemma 1, and the lemma shows that there is a unique solution F that is an equivariant flow. Since H is also divergence-free, as suitable version of Liouville s theorem (e.g. Arnold (1998), Theorem B.3.1) shows F is volume-preserving. A.3 MECHANICAL MODEL We use a nearly incompressible neo-Hookean hyperelastic model, which effectively describes the mechanical behavior of the elastomeric materials commonly employed in manufacturing of flexible meta-materials. The mechanical behavior of such hyperelastic material can be described by introducing a strain energy density function Ψ that is independent of the path of deformation and is a function of the deformation gradient tensor F = u+I. Here, u is the displacement field, which is a function u : Ω R2. A representation of this function is computed by the differentiable mechanics solver described in Section 4. Since the mesh representation of the shape sθ uses N mesh points, and the vector field has two degrees of freedom at each point, the representation of u is an element of R2N. With these ingredients, the function Ψ is given by tr(FT F) 2 2 ln det(F) + λ ln det(F) 2 , (13) where µ = E/2(1+ν) and λ = νE/(1+ν)(1 2ν) are Lamé parameters of a material with Young s modulus E, and ν is the Poisson s ratio. The displacement field of the structure at the equilibrium u is computed from the space of all admissible displacement fields H that satisfy Dirichlet boundary conditions (that is, u takes a suitable prescribed value on the boundary), by minimizing the total stored potential energy from material deformation under mechanical loading u = arg min u H Ω Ψ(F)d X. (14) The minimizer is computed using a second order Newton method for sparse Hessians. Published as a conference paper at ICLR 2025 A.4 MESH INDEPENDENCE ANALYSIS We performed mesh independence analysis for all the cases to ensure simulation results from mechanics solver were not affected by insufficient mesh resolution. The details of an example of these analyses with associated time for each optimization step is shown in Table 2 for designing under uniaxial tension of p1 group. Mesh independence is studied in an initial reference model. Table 2: Mesh independence analysis based on the deformation of the middle right edge of a cellular solid with p1 group symmetry under uniaxial tension. number of Do F 9,000 14,000 20,000 25,000 29,000 middle deformation (mm) 3.05 3.31 3.49 3.65 3.65 optimization time (s) 15 34 49 61 95 A.5 META-MATERIAL DESIGNS WITH P1 SYMMETRY UNDER UNIAXIAL TENSION Three cellular solids designs while pores respecting p1 symmetry group under uniaxial tension when aiming for linear responses with β equals to 0.1, 0.5, and 1.5. The design 4 is with p1 symmetry group when targeting linear response with β = 1.5. A.6 META-MATERIAL DESIGNS WITH NEGATIVE POISSON S RATIOS All eight cellular solids designs with best performances in achieving νef = 0.5 under uniaxial tension while pore shapes respecting different symmetry groups. Published as a conference paper at ICLR 2025 A.7 FABRICATION AND EXPERIMENTAL DETAILS 140mm-by-140mm samples with pull tabs were laser cut from butyl-based laser-engravable rubber sheets of 2.38 mm thickness. Laser cutting was performed using a Universal Laser Systems PLS6.150D machine with two 75-W CO2 lasers emitting at 10.6 µm. Cut patterns were programmed using the Universal Laser Systems software suite with the following settings: 3 % travel speed, 20 % power. Cutting was repeated three times. Experimental tensile tests were performed using an Instron 5969 UTM outfitted with an Instron 2530-series 50N static load cell. Samples were mounted in the UTM s clamps by their pull tabs; a controlled displacement of the top clamp between 5mm and 20mm was applied over the course of three cycles. Video footage of the sample was used to extract the vertical and lateral strain of each sample midway between the clamps, allowing calculation of the average experimental Poisson s ratio νex at 0.1 vertical strain. A.8 META-MATERIAL DESIGNS UNDER UNIAXIAL COMPRESSION Examples of cellular solids designs with target force-displacement responses under uniaxial compression respecting symmetry groups pm, p2mg, and p4 that do not match target values. These simulations ran with several different initialization of the neural network parameters.