Energy-conserving equivariant GNN for elasticity of lattice architected metamaterials
Authors: Ivan Grega
Ilyes Batatia
Gabor Csanyi
Sri Karlapati
Vikram Deshpande
Venue: ICLR 2024 poster (2024)
Keywords: mechanical metamaterials
lattices
elasticity
GNN
equivariant
positive definite
energy conservation
Primary area: applications to physical sciences (physics, chemistry, biology, etc.)
TL;DR: We present equivariant GNN model for prediction of fourth-order stiffness tensor of mechanical metamaterials which is consistent with the law of energy conservation.

Abstract: Lattices are architected metamaterials whose properties strongly depend on their geometrical design. The analogy between lattices and graphs enables the use of graph neural networks (GNNs) as a faster surrogate model compared to traditional methods such as finite element modelling. In this work, we generate a big dataset of structure-property relationships for strut-based lattices. The dataset is made available to the community which can fuel the development of methods anchored in physical principles for the fitting of fourth-order tensors. In addition, we present a higher-order GNN model trained on this dataset. The key features of the model are (i) SE(3) equivariance, and (ii) consistency with the thermodynamic law of conservation of energy. We compare the model to non-equivariant models based on a number of error metrics and demonstrate its benefits in terms of predictive performance and reduced training requirements. Finally, we demonstrate an example application of the model to an architected material design task. The methods which we developed are applicable to fourth-order tensors beyond elasticity such as piezo-optical tensor etc.

Published as a conference paper at ICLR 2024
ENERGY-CONSERVING EQUIVARIANT GNN FOR ELAS- TICITY OF LATTICE ARCHITECTED METAMATERIALS
Ivan Grega' Ilyes Batatia', Gábor Csányi', Sri Karlapati .2., Vikram S. Deshpande" Department of Engineering, University of Cambridge; 2 Amazon Science
ABSTRACT
Lattices are architected metamaterials whose properties strongly depend on their goometrical design. The analogy between lattices and graphs enables the use graph neural networks (GNNs) tional methods such as finite element modelling. In this work, we generate a big dataset of pruuctue-property relationships for strut-based lattices. The dataset is made available to the community which can fuel the development of methods an- chored in physical principles for the fitting of fourth-order tensors. In addition, we present a higher-order GNN model trained on this dataset. The key features of the model are (i) SE(3) equivariance, and (ii) consistency with the thermody- namic law of conservation of energy. We compare the model to -quivvrivat models based on number of eттоr metrics and demonstrate its benefits in terms of predictive performance and reduced training requirements. Finally, we demon- strate an example application of the model to an architected material design task. The methods which we developed are applicable to fourth-order tensors beyond elasticity such as piezo-optical tensor etc.
faster surrogate model compared to tradi-
INTRODUCTION
a
Edge erdeccing
Vode enbedding
Syne
Laken
Nemootss
UAGE
Iraenedion
Era proke
No'l neer mediot
Clach poeling
zt
PSD stckk
Vandel rop representatio
FSD lundlce
1mm
Ien pwwppea
Figure I: (a) X-ray CT scan of 3d-printed lattice. A computer model of the unit cell is shown as an inset. (b) Model schematic. The dimensionality of intermediary quantities is noted between layers
using e 3nn convention. We omit simple linear layers from the diagram for clarity.
A relatively new class of materials, architected Imaterila,,,,,ls, emerged in the last century. (Fleck jet al. 2010) These materials draw inspiration from nature, where many materials are light, yet strong. because of their porosity and microscopic architecture. As subclass of architected ma- terials, lattices are a collection of struts (edges) which are connected at nodes. See Figure T and Figure5inn in the Appendix. Lattices are especially mechanically efficient, offering very high specific stiffness (stiffness divided by density). For instance, is possible to make materials with the density of water and the strength of steel.
"ig3 4 8.8.ammm....mmmmmm....mmmm...... uk
Work done outside of Amazon Science through informal collaboration.

---

Published as a conference paper at ICLR 2024
The established tool for computational analysis of lattices is the finite element (FE) method, which is also the industry standard for other structures from buildings to cars and airplanes. There are a num- ber of principles which the FE solution satisfies (subject to a suitable PDE and constitutive model). First, force equilibrium satisfied at all nodes and the computed displacements are compatible. Second, the strain energy under any deformation is nonnegative as required by energy conservation. Third, results are equivariant to rigid body transformations: rotating the lattice does not change its fundamental properties; they rotate accordingly.
Although the FE method is robust and physically grounded, its high computational cost can be prohibitive: for example, if each unit cell of size lcm is discretized into 100 elements, a wing structure of ~ 20m length would require ~ 10 elements.
Machine learning methods have been used overcome the computational cost of FE methods. In- durkar et al. (2022) employed message-passing GNN to classify lattices based on their mechanical response. Karapiperis Kochmann (2023) used GNN to predict the crack path in disordered lat- tices. Xue et al. (2023) build a GNN to learn the non-linear dynamics of mechanical metamaterials. Maurizi et al. (2022) use GNN to predict the mechanical response composites and lattices. |Myer et al. (2022) have presented GNN framework to predict the stiffness tensor of shell-based lattices. Machine learning methods have also been used to do inverse design of materials. Kumar et al.
2020 couple inverse and forward models to design spinodoid materials with orthotropic symmeetry Bastek et al. (2022) use a similar models for strut-based lattices with fully tailorable 3D anisotropic stiffness. /Zheng al. (2023) build a VAE model for generation of lattices with up to 27 nodes and
cubic symmetry.
While these machine learning models offer a much higher speed than FE, they lack grounding in physical principles. This might not be an issue when the application is restricted to a particular lattice symmetry class and when the model is deployed for data that are close to the training distri- bution. However, models without encoded equivariance and energy conservation principles could fail dramatically if deployed lattice at different orientations might not be ceff-consistent, and negative deformation energy could be predicted, implying the ability to extract energy from passive material.
cuttof-distribution lattice topologies the predictions for the same
In this work, we rely on the equivariant methods which have been introduced by the computational chemistry community. (Thomas et al.) 2018 Batzner et al.) 2022 Batatia et al. 2022) As key contributions:
We introduce a new task into the ML community and provide real-world dataset which can used by researchers in the future to improve higher order physics-focused models.
We present one such model the first equivariant model trained for prediction of the fourth- order elasticity tensor whose predictions are always energy conserving (consistent with the laws of physics).
We benchmark the model against non-equivariant models and show the benefits of key model components and training strategies.
2 BACKGROUND
2.1 EQUIVARIANCE
In the domain of physical sciences, invariance or equivariance under some physical transformations an important property. For example, in chemistry, the energy of a molecule needs to be the same regardless of the coordinate system chosen to represent the coordinates of the atoms. In our modeling of lattices, model predictions need to satisfy similar rules. Let represent a lattice that has attributes of the following types: scalars (e.g. edge lengths L), vectors (e.g. edge directions w = Uz), tensors (e.g. 4th order stiffness tensor C Cijka).
In the discussion of equivariance,
focus two main actions: rigid body rotation and translation. Let Q(: R) represent the rotation given by the rotation matrix R applied on the object, which is the first argument of the function. The following analytical transformation rules apply for scalars T.,
2

---

Published as a conference paper at ICLR 2024
vectors v, and higher order tensors K:
Q(1;
0 = Q(U: R) Rayus
Kigk... Q(K; R)- Ria RiaRBs Rkc. Kabc...
All these objects are invariant to rigid body translation.
The notation (C; R) represents the rotation of the lattice L, which implies the rotation of all attributes of the lattice according to the aforementioned transformation rules. Translation of the lattice T(C; t) simply means displacing the nodal positions by the vector t: f- t.
Our task is to predict the stiffness tensor C for the lattice L. The prediction of model M is M(C). The equivariance requirement is then
M(T(C:t) = M(C) Vt Q(M(C); R) M(Q(E; R)) VR
2.2 EUCLIDEAN EQUIVARIANT MESSAGE PASSING NEURAL NETWORKS
Message Passing Neural Networks Euclidean Equivariant Message Passing Neural Networks (MPNNs) (Liao & Smidt 2023 Thomas al. 2018, Weiler et al. 2018 Kondor al.) [0188] Batzner et al. 2022 Brandstetter et al.) 2022 Batatia et al., 2022, Satorras et al., 2021) are graph neural networks that are equivariant to rotations and translations. MPNNs map a graph with labels
called states on each node i, to a target y. At each layer t, MPNNs operate in four successive steps, the edge embedding, the pooling, the updute and the readout,
= Mc(a;") (t) (.)),, h ((+1) U,(o,") mm)), y= )
(1)
JEN(()
where M, the edge embedding function, + (EN(() the pooling operation (usually just sum) over the neighborhood of the node (i). U2 is the update function. These steps are repeated T times. Finally, the readout R, maps the states to the target quantity.
Most Euclidean MPNNs expand their internal features in a spherical basis. Node features carry an index im specifying the order of the basis expansion.
Equivariant MPNNs
mm (R (ri...., = D' m (1)) (TI,...,TN) TN)
(2)
mid
with D m R the Wigner-D matrices corresponding to the action of the rotation group on the spherical basis. Therefore, this im index is carried over to all internal features of the model.
Higher order MPNNS In full generality, the message can be simultaneous function of all neigh- boring atoms of the central atoms i. Therefore, one can expand the message in many-body expan- sion of the states,
m Σ un ( (0)
Σ S1.ja
(1) u2 (1
Σ S1r..l
(0 :0
(8) J
The number of simultaneous dependency is called the body-order. MPNN potentials were shown to increase the body order of messages by stacking layers. An alternative route is include higher- order terms in the message construction. The MACE (Batatia et al. 2022) architecture, on which we will be building in this work, introduced a systematic way to efficiently approximate equivariant messages of an ordered arbitrary body.
2.3 SOLID MECHANICS
We consider lattices as infinite periodic tessellations of a unit cell. The resulting metamaterial can be characterized by macroscopic (homogenized) properties. The key variables in solid mechanics
3

---

Published as a conference paper at ICLR 2024
under the assumption of small deformations stress, dzj. which is measure of force, and strain e = Ekt. which is a measure of deformation. Both stress and strain are symmetric 3 3 second order tensors (0ig jis Ext E.,)
A third key component in solving a solid mechanics problem is the constitutive law, which relates stress and strain. Linear elasticity postulates that ojj CijkiCki where C Cista is the fourth-order stiffness tensor.
Deformation energy vo under deformation is given by the following tensor contraction. Impor- tantly, thermodynamic laws prescribe non-negative deformation energy & 2 for any admissible deformation e:
1 2 C. 2 CiSCEEENEE
1
0
Ve
Since this has to be true for all strains e, all eigenvalues of the stiffness tensor must be non- negative. The stiffness tensor must be positive semi- semi-definite.
The 4th order stiffness tensor CAju is 3 x 3 x 3 3 tensor. While a tensor with such dimensionality could have up to 81 components, it can be easily shown that the tensor has only 21 independent components, because possesses both minor and major symmetries (Section A.1
Cajui = Cjika = Cijtk = CkRs
3 RELATED WORK
Finite element (FE) modelling The gold standard in computational methods in mechanics has been finite element modelling. In the FE framework, the properties of constituent material (e.g- steel) are known, and FE is used to calculate the structural response. The structure has degrees of freedom új. and loaded by external forces fi. The first step to solving FE problem is to assemble stiffness matrix Kij. which relates displacements and forces: fi Kiyuj- The next step to solve this matrix equation for u (usually by LU factorization). If properties of the material are known, FE provides very accurate predictions of the overall structural response. However, the computational complexity of the matrix inversion (or LU factorization) can very high.
Dataset Lumpe & Stankovic (202]) explored the property space of large dataset of mechanical lattices. The dataset which they used and made available comes from two crystallographic databases (Ramsden et al., 2009; O Keeffe et al.) 2008). The assembled dataset includes nodal positions, edge connectivity, crystal constants, and some elastic properties (Young's 's moduli, shear moduli and Poisson's 's ratios in the three principal directions). We use the crystallographic structures from this dataset basis for our dataset here.
Crystal Graph Convolutions (CGC) and modified CGC (mCGC) To our knowledge, the only existing GNN model used to predict the stiffness tensor of architected materials due to Meyer al. (2022). Instead of ream- lattices, the authors fitted the stiffiness tensor of shell-based lattices. Their model is not equivariant. They use a form of data augmentation whereby each lattice rotated 90º around the z= y- and 2- axis, and mirrored about the x 1--20 y. y 2, and x: planes. This increased the size the training dataset 7-fold. The loss used in training is component-wisg smooth_L. L1 on the 21 independent components of the stiffness tensor.
NNConv for 3d lattices Ross & Hambleton (0220) use GNN to model cubic lattices with 48 rotational and reflectional symmetries. Their model is based on the NNCony layer, which was introduced by Gilmer et al. (2017). In the NNCony model, messages between nodes are formed as a matrix-vector product, where the entries of the matrix are not constant but rather depend on the features of the edge connecting the two nodes.
E(3)- ((3-)aqqivar Message Passing Neural Networks Equivariant Message Passing Neural Net- works (MPNNs) Anderson et al. (019); Thomas et al. (2018); Brandstetter et al. (202);; Batzner et al. )0 Batatia et al. (2022) are class of GNNs that respect Euclidean symmetries (rotations, reflections, and translations). Messages are usually expanded spherical basis, and depending on the order of expansion, not only vectors but also higher-order features such as tensors can be
4

---

Published as a conference paper at ICLR 2024
passed between layers. Equivariant MPNNs have emerged powerful architecture for learning on geometric point clouds.
Methods for enforcing positive (semi-)definiteness Jekel et al. (2022) review a number of meth- ods that can be used ensure that the output of model is positive semi-definite. These include methods based on Cholesky factorization by (Xu et al. 2021; Van 't Sant et al. (003) and meth- ods based on eigenvalue decomposition. Note that these methods cannot be used in our framework because the eigenvalue decomposition often has unstable gradients and assembling the matrix by Cholesky factorization not SO(3) equivariant.
4 METHODS
4.1 DATASET
We created dataset on the basis of the dataset from Lumpe & Stankovic (2021). We process the dataset fix or avoid problematic lattices which reduces the dataset size to 8954 hase lattices. We augment the dataset by introducing nodal perturbations: c.g. for perturbation level 0. each node of a lattice is displaced from the original position by distance 0.1 a random direction. After the new perturbed lattice is obtained, its elastic properties have be computed using FE analysis. Nodal perturbations are applied at levels 0.01, 0.02, 0.03, 0.04, 0.05, 0.07, 0.1. At each level, we formed 10 distinct realizations of nodal perturbations. (Perturbations could only be applied lattices which have at least 2 fundamental nodes.) This enlarged the entire database to 635 454 distinct geometries. For each geometry, FE analysis
relative densities (strut thicknesses).
For the machine learning tasks this paper, we selected from base lattices: (i) 7000 training base lattices, and (ii) 1296 validation/test base lattices. This split ensures that we do not have similar perturbations of the same lattice in both training and test sets.
See the Appendix for the various compositions of the training dataset. Validation and sets fixed for all training runs. Validation set consists of the 1296 lattices without any perturbations, Test set consists of 3 realizations of nodal perturbations level 0.1 for the 1296 lattices. Thus, testing done on OOD data.
4.2 ARCHITECTURE
The diagram in Figure depicts the architecture of the model. Further details are explained in the Appendix. We highlight the main components here. The model relies on the MACE architecture for message passing which was adapted with minor changes. In particular, we used Gaussian embedding of edge scalars and all node features were initialised as ones and expanded using linear layer. We modified the nonlinear readout to enable the processing of higher order tensors.
The significant contribution of this work is the positive semi-definite (PSD) stack. As first step, the fourth-order tensor is transformed to Cartesian basis and then represented in Mandel notation as matrix. Subsequently, suitable PSD) function is applied to the matrix, which enforces its positive -emi-definiteness. This ensures energy conservation which was the key requirement of this work. The following section provides further details about the PSD stack.
4.3 MANDEL REPRESENTATION AND PSD LAYER
fourth-orde Cartesian tensor with major and minor symmetries CiyRı Cigik = Cjisi Cuisj has 21 independent components. Suppose this fourth-order tensor is map between second-order stress and strain:
Oij = CijkiCki
In Mandel notation, the second-order tensors can be written as 6-component vectors, and the fourth- order tensor C can be represented as a 6 x symmetric matrix:
T
(M) 011, 022, F33, V20za, V2 22012
(M) = E11, 622, eaa, V2eza, V2ela, / 22
5

---

Published as a conference paper at ICLR 2024
Cun
C112
C1133
/2C1133
v 2C
V2C1112
Cum
Czzza
Czzsa
22C221
22C2m
V2C221
C
Cas
Csa22
Cas33
/2ss23
V2Cs11
V2C'311
22C2n1
222222
22C2aa
2C2a23
2C2113
2C2312
22Cmm
2C122
v /21333
2C1323
2C1312
2C1312
22C1211
12C1222
V2C1233
2C1223
2C1213
2C1211
The energy conservation requirement can be rewritten as
2 (M)TT (M)
1 - (M)). ((M.T(MM)Ce (M)
Ye (MO)
2
Therefore, positive-defintte fourth-order C is equivalent to positive-definite matrix (MM),
We can apply various methods to enforce the positive definiteness of matrix C(M), These include taking even powers of the matrix, A and A matrix exponential, and its truncated versions, e (I + A/2)2, (T + A/4)a.
prove in the Appendix that the PSD stack maintains equivariance of the framework.
4.4 TRAINING AND EVALUATION DETAILS
Base CGC and mCGC models are trained according to the procedure described in ref Meyer et al. 2022). Model MACE is plain version of MACE that is trained without data augmentation. Where "tr" is added to the model name, it denotes that the model was trained using training method as outlined below. The suffix "+ve" denotes a model which includes the positive semi-definite layer and was trained using our training method. In training method, augmentation, whereby every time a lattice is accessed from the dataset, is returned at a different random orientation (and target stiffness transformed accordingly). Further details about training including the equations for the various types of loss (Leomps Ldit, Ldir,rels Lequiv, Ax)) are in the Appendix.
dynamic data
RESULTS
In this section compare the performance of our model with other models and identify the key components of both the model and training procedures. We also show an example of a downstream application of the GNN model in design task.
5.1 EQUIVARIANT MODELS OUTPERFORM NON-FQUIVARIANT MODELS
In Table I we show the performance of three main classes of models: CGC, NNConv, and MACE for dataset 1 1imp (find dataset details in the Appendix) Crystal graph convolution (CGC) is based on works by Xie & Grossman
(018)) and Meyer et al. (2022). In CGC, a constant learnt matrix multiplies node and edge features to create messages. Models NNConv and MACE are different in the nature of their message passing: the matrix which multiplies node features obtain messages is a function of edge features. Details of all the models are explained the appendix. Comparing errors, is evident that the equivariant MACE model class achieves lowest errors by all metrics.
We observe the following. (i) By adding data augmentation to CGC and NNConv models, models CGC+tr and NNConv+rr achieve substantially reduced stiffiness-bas based errors (Loomp- Ldirs Ldirret), as well as equivariance loss, Lequivv However, these models are more prone to predict negative eigen- values (Xry) ). (ii) Adding data augmentation to equivariant MACE model leads to a much smaller improvement. The percentage of predicted negative eigenvalues, A is not significantly affected. (iii) Models CGC+ve, NNConv+ve, and MACE+ve with encoded positive semi-definite output suf- fer in terms of increased stiffness-based emrors. (iv) CGC-based models outperform NNConv-based models. NNConv models will therefore not be considered in the following sections.
'Te choice of lattice representation. Meyer et al. 2022) combined message passing on primal and dual graph in model mCGC, We do not observe any performance gain from using dual graph representation on data, therefore this model is omitted from the main discussion. We report more details in the Appendix (Table

---

Published as a conference paper at ICLR 2024
Table 1: Performance of the different models and training strategies
CGC
CGC+t
CGC+ve
NNConv
NNConv+tr
NNConvtve
MACE
MACE*t
MACE+ve
Teamp
8.37
4.47
5.47
8.99
5.57
7.07
3.88
3.47
3.61
Ldir
8.77
9.65
6.90
4.17
4.11
4.21
Laiaad
0.44
0.38
0.35
0.21
0.20
0.21
Lequiv
0.33
0.15
0.17
0.39
0.25
0.22
0
Ag
30
1.3
1000
CGC MACE
MACF MACE
0.8
CGC
0.5
60-
Validation
COC-tr
MACE
40-
MACC-+b
0.2
20-
0.0
10P
11F-
10
10a
10F
10a
10"
10"
10*
Training step
Training step
Training step
The evolution of (a) component loss, Trcomp (b) equivariance loss, Teequiv and (c) per- centage of negative eigenvalues, X during training.
Figure 2:
5.2 INDUCTIVE BIASES SUPERIOR TO OBSERVATION AND LEARNING BIASES
As outlined by Kamiaaaaii et al. 2021 there are three conceptual pathways to embedding physics knowledge into machine learning models: observation bias, learning bias and inductive bias. Here we evaluate the efficacy of these biases from the viewpoint equivariance and energy conservation.
Table 2: Performance of various bias types for the learning of equivariance and energy conservation
Equivariance observation CGC+t
Energy conservation
none CGC
inductive MACE
none MACE
learning MACF+Ih
inductive
MACF+ve
Lequiee
0.64
0.13
Ag
30
29
Leomp
9.31
4.47
3.88
Leoop
3,88
3,63
3.61
Equivariance learning We achieve observation bias for equivariance by rotating the data that the model is trained on. As explained in section 4 models which end in "+tr" suffix were trained using data augmentation by rotation. Table shows that the equivariance error reduces dramatically when incorporate rotation augmentation of data. During training of model CGC in Table 2 all lattices were processed at a single orientation. Note that this is different from model CGC in Table I which was trained with 7-fold data augmentation presented Meyer et al. (022)). Figure 2p shows that the equivariance error reduces during training for both CGC and CGC+tr models. Not only is the final equivariance loss for model CGC+tr lower, but also the rate which Lequiv reduces faster. is instructive to note further that while validation loss keeps reducing during training (Figure 2h). the equivariance loss is not a monotonously decreasing function. Model MACE equivariant by design, therefore equivariance loss, Lequiv is zero both in Table and Figure Furthermore, the cquivariant MACE model also has lower component loss, Leump
Energy conservation learning A second physical principle that our model should be aligned with is the positive semi-definitents of stiffness. We evaluate this for the MACE model in Table 2 and Figure Pb. The base MACE model trained on the data predicts negative eigenvalues for 30% of lattices. We attempt to introduce learning bias in model MACE+Ib as follows. Loss is modified include penalty which is calculated from directional projections oy the of predicted stiffness tensor, Ć, into 250 random directions d,; e d: The penalty is then computed as l x relu(-c,) where & is suitably chosen multiplier. This penalty is added loss during training.

---

Published as a conference paper at ICLR 2024
6-
Lown
Loue
3
3-
2
CDC
2-
MACE ve
10'
10
10
2
2
Training dataset
Lescs
Figure 3: (a) Convergence with the amount of data for various model classes. (b) Sensitivity to the maximum 'frequency' Lman and (c) correlation order for MACE model.
Table 3: Various methods for ensuring positive semi- comi-definiteness for equivariant model
eA
A²
A
(I A/2)²
(I A/4)"
Leoop
4.58
3.61
4.41
3.87
3.62
Ldir
5.12
4.21
5.17
4.63
4.41
Tdbraael
0.35
0.21
0.28
0.24
0.26
Table shows that this learning bias is not very effective in guiding the model towards positive semi- semi-definite predictions. Figure 2 shows that the learning bias can have a positive effect during the dynamics of learning, but the final values Ang are similar whether or not learning bias is used.
Scaling with dataset size Using more training data is effectively an observation bias which should lead to better results for all models. In Figure 3a we plot component loss, Leomp- for base CGC model, CGC with data augmentation and positive semi-definite layer (CGC+ve) and MACE model with positive semi-definite layer (MACE+ve). The composition of training datasets explained in the Appendix. At any dataset size, the equivariant MACE+ve model outperforms the CGC-based models. The CGC model with data augmentation outperforms the base CGC model. The scaling slope was calculated as linear fit log-log axes and is displayed on the graph. It is evident that the MACE+ve model has the most favourable scaling.
5.3 CHOICE OF POSITIVE (SENI-)DEFNITEE LAYER AND MACE-SPECIFI PARAMETERS
We evaluate number of methods for making the output positive (semi-)definice for MACE model class. The results are displayed in Table 3 We empirically observe that the matrix square method, A achieves most favourable results. Moreover, this method also has the lowest associated compu- tational cost. For these reasons, we use the matrix square method throughout the paper whenever "4ve" suffix is used, unless stated otherwise.
Spherical frequency Liaz and degeneracy One of the most important hyperparameters of the MACE model is the maximum frequency of expansion in spherical basis, I-max- In Figure 3 we show the sensitivity of model accuracy to the degree of expansion Tamax- Empirically, we observe that degree Liaa = 4 is required to achieve good accuracy. This is in line with the spherical form of the fourth-order stiffness tensor, which contains L = 44 components. Moreover, it has been remarked by Joshi et al. (2023) that certain types of highly symmetric graphs require high order of tensors L. More specifically, to identify the orientation neighbourhood with L -fold symmetry, at least L- order tensors required. In the Appendix, show how model which is internally truncated to Immaz 2 or 3 is unable to predict anisotropic behaviour of simple cubic lattice (Section A.14 Figure 7
Posiitive definite vs semi-definite The physical principle of energy conservation requires nog-negative deformation energy き = ci C// 2 Ove. The case of zero eigenvalue C does not violate energy conservation. A structure whose stiffness tensor has zero eigenvalue is mechanism it is possible to deform it without exerting any work. However, it is a feature of our dataset that all eigenvalues are positive. Therefore, have the freedom try both positive definite and positive semi-definite layers.
8

---

Published as a conference paper at ICLR 2024
xy
z
orgnal
larget
optimized
Figure 4: (a) Unit cell of the lattice used as the starting point for optimization. (b) Original stiffness surface and (c) optimized stiffness surface. Inset shows projections into x y plane of the starting point, optimization target, and the result of ML-based optimization.
Correlation order v In Figure 3 we show the sensitivity of test error to the body-order of mes- sages. Model with == = 1 does not contain the equivariant product layer and equivalent to Tensor Field Network (TFN) with 2-body messages. We observe that increasing the order messages to three- and four-body (v 2 and 3) significantly reduces error over the test dataset.
5.4 SPEEDUP USING MACHINE LEARNING METHODS
Table shows a comparison of inference time for the three classes of investigated GNN models well for finite-element calculation. Time is reported for computation of stiffness tensor for 5000 lattices. The tests were run on desktop computer with Intel i7-1 CPU, 96GB RAM and Nvidia RTX3070 GPU. While the equivariant MACE-based architecture is slower than the more standard CGC- and NNConv-base models, all these models are 3 orders of magnitude faster than FE calculation.
Table 4: Comparison inference speed for the different ML models and the FE baseline
CGC
NNConv
tsuce
5
11x 10
5.5 EXAMPLE APPLICATION: DESIGN OF AN ARCHITECTED MATERIAL
An important application of architected solids is to achieve complex anisotropic stiffness tensor that cannot be found in existing materials. Figure we show an example application of our GNN model in gradient- based optimization scheme for the design of specific stiffness tensor. The starting unit cell, as correctly predicted by the GNN model, has the same stiffness in x- and y-directions. The task is to perturb nodal positions to break the y symmetry and reduce stiffness in the direcctoo.. We use gradients returned from backpropagation and execute 50 steps of gradient descent algorithm. The optimization scheme produced the desired result with great accuracy as verified post-optimization using FE baseline. Error between the desired output and FE-verified ground truth is Loomp 1.84. Based the speedup of the GNN model compared to FE, estimate the GNN- optimisation to also be orders magnitude faster than the FE baseline.
CONCLUSION
In this work, we present the application of Euclidean equivariant GNNs to the prediction of the 4th order stiffness tensor of architected lattice metamaterials, In addition to the intrinsic equiv- ariance to rigid body rotations and translations, we designed the model also preserve positive semi-definiteness of the predicted stiffness, in line with energy conservation. We benchmark the model against other architectures that were previously used for property prediction of lattice mate- rials and demonstrate superior performance by all the metrics studied. Finally, we demonstrate possible downstream use of the model in ML-based structural optimization. Fast and accurate prop- erty prediction models, such the the high computational cost of traditional FE methods and they are applicable to tensors beyond the stiffness tensor such piezo-optical, elasto-optical and the flexoelectric tensors.
presenting, achieve significant improvement over

---

To ensure reproducibility and completeness, we include detailed descriptions of the mod- els used, hyperparameters, and data sources in the Appendix. The code available pub- haat---eti--lttieggmee gnn.git Datasets are available publicly at Published as a conference paper at ICLR 2024
REPRODUCIBILITY STATEMENT
licly at Boo.org/10.7763//AMMMM0
ACKNOWLEDGMENTS
The authors acknowledge funding from the UKRI Frontier Research grant "Graph-based Learn- ing and design of Advanced Mechanical Metamaterials" with award number EP/X02394X/IL This work was performed using resources provided by the Cambridge Service for Data Driven Discovery (CSD3) operated by the University of Cambridge Research Computing Service www.csd..cam.a..uut provided by Dell EMC and Intel using Tier-2 funding from the Engineering and Physical Sciences Research Council (capital grant EP/T022159/1), and DiRAC funding from the Science and Technology Facilities Council www.:w.da..a.... Padmeya Prashant Indurkar for the initial exploration of graph neural networks for lattice metama- terials and are grateful to Angkur Shaikeea for the fruitful discussions throughout the project.
The authors further acknowledge
REFERENCES
Brandon Anderson, Truong Son Hy, and Risi Kondor. Cormorant: Covariant molecular neural networks. In H. Wallach, H. Larochelle, A. Beygelzimer, F. AlcheBuc, E. Fox, and R. Gar- nett (eds.), Advances in Neural Information Processing Systems, volume 32. Curran Asso-
ciates, Inc., 2019. URL 03573b3262 (6666/8888999911122220000---0000 http://a...rrrr.o/ppp/////////1
pdff
Peter J. Basser and Sinisa Pajevic. Spectral decomposition of 4th-order covariance tensor: Appli-
cations to diffusion tensor MRI. Signal Processing, 87(2):220 2007. ISSN 0165-1684. doi: http::/doi.org/100 101 oiw.::/o.iiipro..00...022222.
httpttss:: www. www.sciieessc. rect. com/
Jan-Hendrik Bastek, Siddhant Kumar, Bastian Telgen, Raphaël N. Glaesener, and Dennis M.
Kochmann. Inverting the structure- duructurepropertyy map of truss metamaterials by deep learning. Pro- ceedings of the National Academy of Sciences, (19())e2111505110) January 2022. ISSN 0027- 8424, 1091-6490. doi: 12.1177//paa.2211555111 URL https:
10 1073/ pnas 2111105111
Ilyes Batatia, David Kovacs, Gregor Simm, Christoph Ortner, and Gabor Csanyi. Mace: Higher order equivariant message passing neural networks for fast and accurate force fields. In S. Koyejo, S. Mohamed, A. Agarwal, D. Belgrave, K. Cho, and A. Oh (eds.), Advances in Neural Information Processing Systems, volume 35, pp. 11423-11433 1136.. Curran Associates, Inc., 2022. URL
ccpppe Files/paper/2027/ Saw.appper-oooefenec..mm pdf
file
Simon Batzner, Albert Musaelian, Lixin Sun, Mario Geiger, Jonathan P. Mailoa, Mordechai Ko- mbluth, Nicola Molinari, Tess E. Smidt, and Boris Kozinsky. E(3)-equivariant graph neural networks for data-efficient and accurate interatomic potentials. Nature Communications, 13 (1):2453, May 2022. ISSN 2041-1723. doi: 10. 10.1038/s41467-10 4.9-02299..... 7/doi. 10w.11038/1661. 022 29939-5
URL https:
Johannes Brandstetter, Rob Hesselink, Elise van der Pol, Erik J Bekkers, and Max Welling. Geomet- ric and physical quantities improve e(3) equivariant message passing, In International Confer-
ence on Leurning Representations, 2022. URL htttps://openreviest net/ forum?id
xwr8gOBeVI
N. A. Fleck, V. S. Deshpande, and M. F. Ashby. Micro-architectured materials: past, present and future. Proceedings of the Royal Society A: Mathematical,, Physical and Engineering Sciences, 466(2121):2495 2516, June 2010. doi: 10. 10.10/ssp..2000.0225 URL https:
10. 10.1008/rae rww.///spa... 2010 201000001110 0215
10

---

Published as a conference paper at ICLR 2024
Mario Geiger, Tess Smidt, Alby M., Benjamin Kurt Miller, Wouter Boomsma, Bradley Dice, Kostiantyn Lapchevskyi. Maurice Weiler, Michat Tyszkiewicz, Martin Uhrin, Simon Batzner, Dylan Madisetti, Jes Frellsen, Nuri Jung, Sophia Sanbom, jkh, Mingjaan Wen, Josh Rackers, Marcel Rød, and Michael Bailey. e3nn/e3nn: 2022-12-12, December 2022. URL https:
//doi. 10. /221/zenooo. 7430260
Justin Gilmer, Samuel S. Schoenholz, Patrick F. Riley, Oriol Vinyals, and George E. Dahl. Neural message passing for quantum chemistry. In Proceedings of the 34th International Conference on Machine Learning- Volume 70, ICML' 17, pp. 1263 1272. JMLR.org. 2017.
Peter Helnwein. Some remarks on the compressed matrix representation of symmet- ric second-order and fourth-order tensors. Computer Methods Applied Mechanics
and Engineering, 190122):7553277700 2001. ISSN 0045-7825. doi: https://doi.org//
1016/800045-7255500000662-22 URL https: atticlle/p11//00457882505
httppp::/www...iiet.t... /ttp://wwwcinncieeeetttco wwwpeeedddrreeeeeeeeeeeeeeeeeeee scinncci
Padmeya Prashant Indurkar, Sri Karlapati, Angkur Jyoti Dipanka Shaikeea, and Vikram S. Desh-
pande. Predicting deformation mechanisms architected metamaterials using gnn, 2022. URL
http::::/ddo..org//// 0ttp//100.85550/1 48550/ arxiv 2202. 09427
Charles F. Jekel, Kenneth E. Swartz, Daniel A. White, Daniel A. Tortorelli, and Seth E. Watts. Neural network layers for prediction of positive definite elastic stiffness tensors, 2022. URL
fittss//doo wr..orgg/oo //rgg/1000 48550/48555668 wiw.i.vv.... 2203 13938
Chaitanya K. Joshi, Cristian Bodnar, Simon V. Mathis, Taco Cohen, and Pietro Liò. On the expres- sive power of geometric graph neural networks, 2023.
Konstantinos Karapiperis and Dennis M. Kochmann. Prediction and control of fracture paths in disordered architected materials using graph neural networks. Communications Engineering, (1):32, June 2023. ISSN 2731-3395. doi: (0.1038//44122-023-00000
URLhttpps://ooi.
George Em Karniadakis, Ioannis G. Kevrekidis, Lu Lu, Paris Perdikaris, Sifan Wang, and Liu Yang. Physics-informed machine learning. Nature Reviews Physics, 36))422-440, June 2021. ISSN 2522-5820. doi: 10. 0011038//2225--01-0000111 942254-021-0031-
URL http:::tps//soo /ttp://22c..r11111133333
Risi Kondor, Zhen Lin, and Shubhendu Trivedi. Clebsch-gordan nets: fully fourier space spherical convolutional neural network. In S. Bengio, H. Wallach, H. Larochelle, K. Grauman, N. Cesa- Bianchi, and R. Garnett (eds.), Advances in Neural Information Processing Systems, volume 31.
Curran Associates, Inc., 2018 URL ht httt:::://ppccceediinig www.rieersssp..p cwwpuiip..ccppppp....ppppp
Siddhant Kumar, Stephanie Tan, Li Zheng, and Dennis M. Kochmann. Inverse-designed spinodoid metamaterials. upj Computational Materials, 6(1):73, June 2020. ISSN 2057-3960. doi: 10. 0008/441524-000004----- URL
Yi-Lun Liao and Tess Smidt. Equiformer: Equivariant graph attention transformer for 3d atomistic graphs. In The Eleventh International Conference on Learning Representations, 2023. URL
https: http///ppeeeeeeeeeeeeeeee net www.moiddorroooo f ARgOTD
Thomas S. Lumpe and Tino Stankovic. Exploring the property space of periodic cellular structures based on crystal networks. Proceedings of the National Academy of Sciences, 118(7), Febru- ary 2021. doi:
URL https: htttp:://ddi..og/1000
2003504118
Marco Maurizi, Chao Gao, and Filippo Berto. Predicting stress, strain and deformation fields in materials and structures with graph neural networks. Scientific Reports, 12(1):21834, December 2022. ISSN 2045-222 2322. doi: 11033/441198-02222644401 URL ht //doi. org/10.

---

Published as a conference paper at ICLR 2024
Paul P. Meyer, Colin Bonatti, Thomas Sancogne-Dejenn and Dirk Mohr. Graph-based meta-
materials: Deep learning of structure-property relations. Materials & Design, 223:111175,
2022. ISSN 0264-1275. doi: http://dcioogp//...06/mmmteeeee.
URL hetps:
// //www. sc encedi rect. cww
Tomáš Mánik. A natural vector/matrix notation applied in an efficient and robust return- mapping
algorithm for advanced yield functions. European Journal of Mechamics A/Solids, 90: 90:104357, 2021. ISSN 0997-7538. doi: https://dbi.org/100/ 101 1010000111...onn.. htt.:///gw.......ooooooooooooooa
URL (tLppsss:
// /wwww www.sciencedirect..om htt.a//ccw....aaao.oaaaaaaaaaata
Michael O'Keeffe, Maxim A. Peskov, Stuart J. Ramsden, and Omar M. Yaghi. The reticu- lar chemistry structure resource (RCSR) database of, and symbols for, crystal nets. Accounts
of Chemical Research, 41(12):1782 1789, October 2008. doi: 10.1021/ar801124u URL
Stttps://ol1 wtt.:/i......................... 11210//0111.111///00002111
S. J. Ramsden, V. Robins, and S. T. Hyde. Three-dimensional euclidean nets from two-dimensional hyperbolic tilings: kaleidoscopic examples. Acta Crystallographica Section Foundations of
Crystallography, 65(2):81 January 2009, doi: 10. 10.11077601188767380000065333331 7/doi.orq/ 10.
Elissa Ross and Daniel Hambleton. Using graph neural networks approximate mechanical response on 3d lattice structures. AAG2020, 2020. URL Stttps://chinkshelll. fr/ AAA2022 24 Ross pdf
wp htgp-conteerttaplllaas///oo
Víctor Garcia Satorras, Emiel Hoogeboom, and Max Welling. E(n) equivariant graph neural net- works. In Marina Mcila and Tong Zhang (eds.), Proceedings of the 38th International Con- ference on Machine Learning, volume 139 Proceedings of Machine Learning Research, pp. 9323 9332. PMLR, 18 24 Jul 2021. URL htttps://proceeeee aatorras2la.html
mwr. www.mllr.rr.rr pre.s////s///3///
Nathaniel Thomas, Tess Smidt, Steven Kearnes, Lusann Yang, Li Li, Kai Kohlhoff, and Patrick Riley. Tensor field networks: Rotation- and translation-equiivrraan neural networks for 3d point clouds, 2018.
William Thomson. Xxi. elements a mathematical theory of elasticity. Philosophical Transactions of the Royal Society of London, (146):481-498) 1856.
Sikko Van 't Sant, Prakash Thakolkaran, Jonàs Martínez, and Siddhant Kumar. Inverse-designed growth-based cellular metamaterials. Mechanics of Materials, 182: 104668, 2023. ISSN 0167-6636. doi: http:///doi.or//0..0000 mec.m///.mccmma222233..3 20..20223.1666. URL /https: //op:///www.. /www. sciencedirect. http.///eeeeeeeeaeeeeeeeeaaaaa/a
Maurice Weiler, Mario Geiger. Max Welling, Wouter Boomsma, and Taco S Cohen. 3d steerable cnns: Learning rotationally equivariant features in volumetric data. In S. Ben- gio, H. Wallach, H. Larochelle, K. Grauman, N. Cesa-Bianchi, and R. Gamett (eds.), Advances in Neural Information Processing Systems, volume 31. Curran Associates, Inc., 2018. URL httpp:://proceeeeeeseste igg wwt.:rt....reeccceoo.eeeee ccpppper..pp...peo//pp/e [ir ffiles/paperr/211 file/ TTle/A88641055000006222255321123
pdf
Tian Xie and Jeffrey C. Grossman. Crystal graph convolutional neural networks for an accurate and interpretable prediction of material properties. Phys. Rev. Leut, 120: 14531453,,, Apr 2018. doi: 10. (0.11103/hyssev..comm 120.144301. URL. httphht//////itttt aps wttp:/pp.o.o/o/....1111
PhysRevLett 120. 120.145301
Kailai Xu, Daniel Z. Huang, and Eric Darve. ric positive definite neural networks. Journal of Computational Physics, 428: 10072, 2021. ISSN 0021-9991.
Learning constitutive relations using symmet-
doi: http:///to.oo////....00gggg22200 wtt.i//iiiittcccceeeooooooaaaaas
URL httpp:://att //www.
Tianju Xue, Sigrid Adriaenssens, and Sheng Mao. Learning the nonlinear dynamics of mechanical
metamaterials with graph networks. International Journal Mechanical Sciences. 238: 107835, 2023. ISSN 0020- 7403. doi: hut.s//noooooono/aooooo/oooooonn
107835. URL https:
7/ 7/www. scwenceieeeeoo http//aawee..eee....c///eeeeeeaa htt.aaccc......c..occ..aaaaaaaaa
12

---

Published as a conference paper at ICLR 2024
Li Zheng, Siddhant Kumar, and Dennis M. Kochmann. Unifying the design space of truss metama- terials by generative modeling, 2023.
A APPENDIX
A.1
SYMMETRIES OF THE STIFFNESS TENSOR
From equation stiffness tensor can be expressed as the derivative of strain energy with respect to strain:
8h
C Cijkt dei
The order of 6ij and 66t is interchangeable, which results in the major symmetry for the stiffness tensor: Cejkı Clai--
Furthermore, strain is defined as the symmetric gradient of displacement, u:
1 Ou:
Ouj
Eij
Ox:
Therefore, Eji, which gives rise to the minor symmetry of the stiffness tensor. All in all, the stiffness tensor C has both minor and major symmetries:
Cijkt = Cjúat = Cijik = Cklij
As result, 21 of the 3 3 x x 3 81 components of the stiffness tensor are independent.
A.2 METHODS FOR ENFORCING POSITIVE (SEMI-) DEMII)DEFHNITENESS
Here we outline methods which can be used to enforce positive (semi- maniidddffinteness for n x n matrices R" + R Section A.4 explains how the thh order stiffness tensor can be efficiently represented in a matrix form which justifies the use of these methods.
A.2.1 CHOLESKY-BASED METHOD
Cholesky decomposition is defined for Hermitian positive-definite matrix R - R. as:
= -
where L is a lower diagonal matrix and L' is its conjugate transpose. The diagonal entries of positive.
Machine learning methods Xu al. 2021 Jekel et al. 2022 Van 't Sant al. 2023)
Cholesky factorization as follows. Suppose we require n n positive definite matrix A : R" + R neural network outputs & n(n + 1 entries: ao. ... ag. They are arranged into lower diagonal matrix L with the diagonal elements passed through suitable function : IR + Rso (such as p(x) exp(x)):
T((ao) a1
p(az)
a3
p(as)
Matrix product LL" then guarantees a symmetric (Hermitian) positive-definite matrix. Note that zero eigenvalues admissible, different function p: R R20 be used (e.g. relu).
Such construction, while simple, will not produce equivariant output because components of the matrix 6g, .. are treated independent scalars.
13

---

Published as a conference paper at ICLR 2024
A.2.2 EIGENVALUE- METHOD
Symmetric matrix is positive definite iff all its eigenvalues positive. The eigenvalue-based methods operate on this premise (Jekel et al. 2022).
Suppose we require
positive definite matrix A Rd + R" A neural network again outputs k = n(n + 1 entries: 40, ao..... ag. They are arranged into symmetric matrix M:
Tag
01
03
Q2
M = 03 04
05
Eigenvalue decomposition is performed on this matrix: M UAUT. Next, suitable function p: IR + Rso is applied to the eigenvalue matrix A:
(()1) 0
p(22)
A
and positive definite matrix A is assembled A UA+UT. Similarly, for positive semi- definiteness, function p R - R20 should be used.
The advantage of this method, opposed to the Cholesky- based method, is that the geometric rep- resentation of eigenvectors is maintained in other words, if the overall model had been equivariant with respect to vectors in U, it will remain equivariant after eigenvalues significant disadvantage of this method is that eigenvalue decomposition is not a stable operation with respect to gradients, which is also noted the official PyTorch documentation.
made positive. The
A.2.3 MATRIX POWER AND MATRIX EXPONENTIAL
To avoid the computational complexity and gradient instability of eigenvalue decomposition, can look for methods which will provide the same result - matrix with positive eigenvalues without ex- plicitly computing the eigenvalue decomposition. We have experimented with number of methods which are based on taking even powers of matrix and calculating matrix exponential.
Matrix exponential
The action of matrix exponential square symmetric matrix M is
A matrix exp( MM AT
iec eigenvalues of A are exponentiated eigenvalues of M.
The method is usually implemented as an iterative algorithm in which the explicit calcula- tion of eigenvectors not required. While it is stable with respect gradients, its exe- cution takes 1.5 times longer than computing eigenvalue decomposition (comparing PyTorch linalg. matrrix (xtrixxexp(M) and l sinalg. eig)(M)). The key difference between matrix ex- ponential and matrix powers is that it produces strictly positive eigenvalues (and hence positive definite matrix).
Matrix power
Even powers of symmetric x matrix ensure non-negative eigenvalues:
A = M = U
Therefore, has the same effect carrying out the eigenvalue decomposition and raising the eigen- values to power n. However, it has a lower complexity and could be up to 80 times faster (comparing PyTorch 1 linalg mat x powe M, 2) and linalg g. eigh(M)).
We evaluate the performance of ndd and th power in Section 5.3
Truncated matrix exponential
One of the ways to write matrix exponential
M
lim k-100
I+
14

---

Published as a conference paper at ICLR 2024
We evaluate the performance of positive semi-definite layer for and 4 in Section 5.3
An important advantage these methods as opposed to Cholesky-basee methods that they can be made equivariant. However, that is predicated on using the Mandel notation as opposed to the more traditional Voigt notation, as discussed in the following section.
A.3 SPECTRUM OF THE TtH ORDER TENSOR
As outlined by Lord Kelvin Thomson (856), there are principal strains EE() such that stress is parallel to strain under that deformation:
(<== = C: : E())
where ( is a scalar eigenvalue of the stiffness tensor and E 1,... 6 are the six 2nd order eigentensors. See also more recent text by Basser Pajevic (007)
The equivalence between tensor notation using C, fij. dij and vector/matrix notation using C(M), clAT), ((M) provides way to calculate the eigenvalues and eigentensors of the 4th der tensor C. The eigenvalues (+) for tensor C are the eigenvalues of matrix C(M), and the cigentensors ( are obtained by rearranging the cigenvectors of C(M), Suppose = ei 622, E33, /2 V2E13, 2212] is an eigenvector of ((M)), The corresponding eigentensor for C
[1
612
11
E
612
(22
C23
C13
23
333.
A.4 MANDEL/KELVIN NOTATION VOIGT NOTATION
Stress oij and strain Gy are 2nd order symmetric tensors:
0i
013
E11
612
613
012
022
023
612
622
623
161
623
633
They have 6 independent components: three direct components (indexed by 11,22,33), and three shear components (indexed by 12,13,23). They are often arranged in vector form using Voigt nota- tion.
011
611
0
C22
(v)
033
e(V)
C33
d23
2023
013
2313
012
12
The factor of 2 in front of shear components of strain is to preserve the dot product equivalence for strain energy: in 2nd order notation, strain energy can be written as the contraction of stress and strain:
2 : - 0 (V)
Following this notation, the 4th order stiffness tensor can be represented as 6 x 6 matrix (V) such that (V) = (V);
Cuu
C1122
C1133
C1123
C1113
Cun22
C2211
C2222
C2233
C2223
C2213
C2212
CV
C3311
C3322
C3333
C3323
C3313
C3312
C2311
C2322
C2333
C2523
C2313
C2312
Cran
C1322
C1333
C1323
C1312
C1312
C1211
C1222
C1233
C1223
Cia1a
C1212
The disadvantage of this approach is that stress and strain expressed in contravariant and co-
variant bases, respectively, which do not coincide (Helnwein) 2001) Mānik (02T). This makes the Voigt notation unsuitable for GNN model. In particular, if the equivariant GNN model outputs
15

---

Published as a conference paper at ICLR 2024
4th order tensor C which arranged into 6 x 6 matrix C(V using the Voigt notation, and we then apply a positive definite layer (e.g. matrix square), the output matrix loses equivariance property (as further explained in the following section).
This issue can be resolved using the Mandel/Kelvin notation. In Mandel notation, the second order stress and strain tensors also written 6-dimensional vectors, but they take the following form:
T1
E11
222
e22
((M)=
333
c(M)
E33
22
/2623
v 2013
V211
12
2212
Strain energy can still be expressed as contraction
(M (M)
Moreover, the norm of stress and strain is preserved under this notation
|||| = e: == (M);
|||| = o o ( (M)
The corresponding stiffness tensor (MM) can be written such that ((M) ((M)
Cun
1122
CHa3
V2C123
V2C113
22
Czu
C2222
C22a3
V2C2223
V2C2213
V2Cr22
((M)
Can
C3322
C3333
V2C3323
V2C3318
V2C3a12
V2C2s1
V2C2322
V2C2333
2C2323
2C2313
2C2312
22Cis11
V2C1322
V V2C1333
2C1323
2C1312
2C1312
/2C1211
V2C1222
12C1233
2C1223
2C1213
2C1212
Using this notation, both stress and strain are expressed in the same orthonormal basis. Contrary to using Voigt notation, GNN model outputs 4th order tensor C, we can arrange the components into x 6 stiffness matrix, and apply positive definite layer to this matrix. Importantly, this pipeline will satisfy equivariance as shown below.
the Mandel notation in equivariant framework. If equivariant
A.5 PROOF OF EQUIVARIANCE OF PSD LAYER IN MANDEL NOTATION
Let M be the output (arranged in Mandel notation) of equivariant MACE backbone M for lattice L
M(C) = M
and let Q( R) represent the rotation given by the rotation matrix R applied the object, which is the first argument of the function. We postulate that the representation of the rotation group in Mandel basis can be written terms of matrix R(M) such that the rotated output M given by
M = M(Q(C))) R)) = M) M M (M)T
After pass through the PSD layer (for instance using matrix square), the output and its rotated version are given by
A = 2
(5)
Á M (M)) M ((M) (
Therefore, the PSD layer constructed in this way is equivariant iff matrix R(M) is orthonormal: RMM) R(M) I. this section, we prove both the postulate that the rotation of stiffness tensor in Mandel basis can be expressed using matrix R(M) and the statement that this matrix is orthonormal.
16

---

Published as a conference paper at ICLR 2024
We first consider the basis of representation of stress in Mandel notation. For stress ((M) with components [u,,..., Ug]. the basis is formed by the following second order tensors
"ie(1) & e
u1
tuze e(22+
12
+uze(8) e
03
(e (()) ((3) e (2)
U5
(e(1)
60 ((8)
e())
u6
2g V2 (e(1)
2 e(2)
x e1))
We now proceed to derive the representation of SO(3) rotation group in Mandel notation. With- out loss of generality, we assume that the basis vectors e(1), e(2), e(3) are originally aligned with the Cartesian axes. Therefore, the i-th component of vector e is equivalent Kronecker delta: e(:) = ōg. The effect of rotation on the basis vectors e(1), e(2), (() can be expressed by matrix multiplication with the conventional rotation matrix Rij as
R
where ((P) the basis vector e()) expressed in the rotated basis and we use standard Einstein summation convention for repeated indices. The elements of the rotation matrix are therefore
Rij e(4) e
We now express stress the rotated frame ē in terms of the components of Mandel representation U1 ... Ug:
ōi = Rap Rip
uio c/ +
+
U1
(3)
+ +ep 1000000)
V2
= Rip undipóig +...
U4 (82p02g + d2g) + +...
v2
This can be written matrix- product in Mandel representation I (M)
(M) (M)
(M))
Ri
Ri2
Ri3
/ 2212R
22 R1a
/ /2RnRi
u1
R
Riu
R
2R
22 R23
v 12
Rsi
R3s
/2R Ray
/2Rn Ras
/2Rn Ryz
22 Rai
/2Rga Raz
22 Rea
R22 Raa + Reg Raz
R2i Raa + R2g
R21 Ra2 + Rez R8i
U4
/2Rn Rai
/2R12Rs2
/2R1
R12R33 + R1a
Rn Rss + Rig
Ry Rz2 + Ri2
Up
/2Rn Rai
V2R R22
/2R1aR2
R12 + Rig
Ru Ras + Ria
Ru Raz + R12 R21
Matrix R(M is the representation of SO(3) rotation in Mandel notation. We can now proceed to derive the corresponding rotation rule for the stiffness matrix C(M), In the original frame,
(MM) = (an)
(8)
while the rotated frame:
(M) = ((a))
(M)) ( (CM) (M)
g (M) (M R (M)
comparing equations equation g and equation 9 we obtain the rotation rule for stiffness matrix C(M);
(M)) R (M) (M.)-1
(10)
17

---

Published as a conference paper at ICLR 2024
We proceed to show that matrix R is orthonormal. This can be done by expanding the product R(M), R(M) and showing that R(A) R(M) diy, but we choose an alternative route: by considering the double contraction of stress dot product in Mandel basis.
pj
From equation equation] 7 is straightforward to show that
(M) (M)
Vai
We then consider this contraction in rotated basis:
(M) (M)
0 R a v
(M) Rig o (M) (M)) R.ATRR(M (M) g
Therefore, to show that matrix R(M) is orthonormal, it suffices show that of a( (M) (M)
o; 0;
o. (M) g
BaglgBaggi unddpprng +...
u4
(S2yday dapō2e) + ...) (undhnous + u4 (d2. S3ad2u) + ... )..)
V2
V2
6p
uuddppdg +... ... u4
v2 (dzpôro + )epd2a) +...) [uirradi +..
us
(d2adys + dsad2s) + t..)
v2
+
u2 ug
un
+ us up o. (M) g(1))
where we used the fact that the conventional x 3 rotation matrix R is orthonormal Rip Ria = 6pa- Therefore
R R (M)
We reformulate the rotation rule for stiffness in Mandel notation from equation equation 10
((M) R (M)) (M))-1
which validates the postulate of this proof. This combined with equation equation 6 proves that our PSD layer is equivariant.
A.6 PROOF OF NON-FQUTVARIANCE OF PSD LAYER IN VOIGT NOTATION
Analogous analysis can be performed for the stress/strain in Voigt basis. However, in Voigt notation, the bases for strain and stress are different:
(1e(1) @ e(1
t2) e2)
U2
(3ge e
(v)
U3
+ua (e(2)
(a) +e(*) & e
u4
Us
+us (e(1) (3) +e(s) (1)
16
+ug
(e(1)
(2) e(1)
(ee())
U11
(2)e(2)
U2
tuge(8) (e(+)
(V)
u3
(e(²) @ ( +e(() ((()) +
u4
us
(e(1) ((() +el) (e())) +
46
the
(e(1) e( +e e(1)
18

---

Published as a conference paper at ICLR 2024
a
Metamolerial
Micruercillcturee
Strul-based unit oel
Shellbaased unit cell
10
1mm
10 pm
Figure 5: (a) The assembly of millions of unit cells effectively behaves as continuum material, hence the name metamaterial. In this work assume all unit cells have cylindrical struts with radius r. (b) Different types of unit cells could used, such as triply periodic minimal surfaces (TPMS) which lead to shell-based based lattices.
It can be shown through analogous process that the corresponding representations of the rotation group are matrices R (o)) and R(V)):
g(V)
(V))))
Bi
Ria
Rīa
213
2Rn
2121
U1
Ri
Ry
Rig
2 2Rz Ria
2 2 Res
2 2221 Rez
u2
Rgi
Rga
Ras
2 Rs3
2 Rss
2Rs1 Rs2
ug
Ri Rai
Rez
Raa
RezRas + RaaR
Rn + RaaRa
Ru Ruz + RezRan
Rn Rai
Ryz
RigRas
Ry2Ras + R1aR
Rn + RIyRa
Rn Ryg + Rı
us
Rn Re
Ry2
Raa
R12R23 + RıaRza
R1R2 + RiaRa
R1 Rr2 + R12
UG
e(V) R (V)
Ria
Ria
R12R13
R11
Rn
M1
R2i
Ria
Re2R2s
Ru R2a
Ru Ru
R2y
R5z
R2g
Rsz Reg
Rsi R33
Rai Ryz
u3
24 Rai
2Rz
2R22
Re2Ra
+ Rag
Rn Ras + Reg
Ru Raz + Rez Ran
11
Rai
2R22Ra2
Ry2R33
R12
Ri Ras + Ria
Rn + Ri2Rai
28 R2i
2R12H22
2212
R12f23 RisR22
R1123 Ris
R1, R12
and can be shown that
R(Vo). (V.).T
We can now derive the rotation rule for stiffness in Voigt notation as
((V) (V)
8(V) = (()))(
R(V.0) (M)) -CV) (V)
( (VV))
O RR(V (Vo)-1 C(V p(
Importantly, matrices (Voo and RVV,e) are not orthonormal (R( R(oo + 7) which implies that using Voigt representation in PSD layer breaks equivariance.
A.7 LATTICE METAMATERIALS. GRAPH REPRESENTATION OF LATTICE UNIT CELLS AND FE
Figure SeeureBJoutinn oottlines the concept of metamaterials. We consider periodic lattices which constructed by repeating a predefined building block in three dimensions. This building block is called unit cell. When the unit cell is repeated over distance much longer than its size, there will be millions of unit cells in the sample of interest and its overall response can be characterized by effective material properties.
In this work, we model strut-based lattices. There is a clean analogy between such lattice and the mathematical concept of graph. The struts in lattice be thought of edges, while their
19

---

Published as a conference paper at ICLR 2024
intersections nodes. We wish model large samples of metamaterials, which comprise millions of unit cells. Rather than converting such sample into a graph with billions nodes and edges, we use the concept of periodicity: the lattice can be fully defined by its unit cell and wish to predict the stiffness of the material from the geometry of the unit cell.
All unit cells be represented in reduced coordinate system, where the unit cell is unit cube (0 S z; S 1). The real positions of nodes, transformed coordinates ε, can be expressed as an affine transformation of the reduced coordinates: æ Ax. where A is suitable matrix. Unit cells that we study in this work range from very simple geometries with just a few nodes to very complex unit cells with hundreds of nodes. We define four node types based on the following conditions:
1. inner nodes, where 0 < 24 < 141,
2. face nodes, where z; E {0, }} for only one index i,
3. edge nodes, where æ; {0, } for two indices i,
4. corner nodes, where z; {0. 1 for all three indices i.
4 node types. Further note that face, edge and corner nodes are shared by 2, 4 and 8 neighboring unit cells, respectively.
The unit cells are representation the infinite periodic lattice. Therefore, it is possible to shift the window of observation to perceive different unit cell of the same lattice. This is illustrated in Figure 66--c. The simple cubic lattice is typically represented as a square with four edges on the boundaries and four edge nodes Displacing the unit cell window, obtain view with 4 face nodes and 1 inner node.
The fundamental representation of a lattice is such where only the inner nodes are kept and edges are connected across periodic boundaries. In Figure we show the fundamental representation of the simple eubic lattice. The lattice has inner node (NO) and 2 fundamental edges (EO, B1). The edges are defined by edge adjacency, and edge shifts:
adjacency shift
FO:
NO > NO; [1,0]
E1:
NO - NO; [0.1]
edge has adjacency N₂ > Nj and shift (6)) then the edge vector will be v(()) x()) x ((
Graph representation is obtained when graph is connected according to the fundamental edge adja- cency and edge shifts are stored with the graph (Fig. (6).
Note that finite element simulations makes the handling of periodicity much simpler. In particular, when macroscopic strain 6j is applied to the material, the following equations are prescribed in FE setup:
the windowed representation of lattices, because
B u
u
Σ
Ei
(1
x11)
0 0.1 0
A.8 DATASET
The full dataset contains 8954 base lattices. We selected from this dataset: (i) 7000 training base names, and (ii) 1296 validation/test base names.
Training sets of various sizes formed follows:
note that intuitively, might call these nodes corner nodes, but to adhere to definitions above, in 2d 3% E {0, 1} for two indices
20

---

Published as a conference paper at ICLR 2024
Typical oall representation nodes, edges)
Fundamental representatin (1 node, 2 edges)
Infinite periodic latlice
NO
edge node
Windowed representation nodes, edges)
conner node
Graph representation node, edges)
face node
inner node
E0
Figure 6: The representation of lattice unit cell in our framework. For simplicity, we show a two- dimensional simple cubic lattice as an example (a). The typical representation of the unit cell is square (b). However, to simplify the periodicity conditions in FE calculations, we transform this unit cell to windowed representation (c). The fundamental representation connects edges across periodic boundaries and only inner node types remain (d). The graph used for GNN comes directly from the fundamental representation (e). For completeness, we illustrate the remaining edge node type in part (f).
graphs
Description
0imp quarter
1750
1750 lattices with no perturbations
0imp half
3500
3500 lattices with no perturbations
0imp
7000
7000 lattices with no perturbations
1imp
27847
7000 lattices
realization at 0.0.0.02.0.04.0000 levels
2imp
48681
7000 lattices with realization at and realizations 0.02.0.04.0.07 levels
4imp
90336
7000 lattices with realization at and realizations at 0.02.0.04.0.07 levels
Note that three relative densities (strut radii) of each distinct geometry are used.
A.9 GRAPH ATTRIBUTES
Vanilla CGC In the base CGC model, we use the same node, edge and graph features as Meyer ct al. (022) with the addition of strut radius. Node features of the graph are nodal positions:
h [21, 23)
Edge features are unit vector, length, and radius:
e = [1,, U2, U3, L.
Augmented GCG-based models We wish to have model which is invariant to rigid body trans- lation, therefore we drop nodal positions from node features. The following input features are used.
h = [1]
6 = [u1, U2 U3, L. r.
A. ARCHITECTURE
We consider lattices geometric graphs, with node positions, æ; E R33 edge adjacency, ((i,))) edge shift: El E 123 and edge thickness, Tt E Rt, Note that edge adjacency is multiset as there can multiple edges between the same set of nodes. The role of edge shifts account for periodic connections. Further detail can be found in the Appendix. The model we develop acts in three steps, first the embedding, then layers of MACE, and finally the readout.
As defined in Section A.7 of the Appendix
21

---

Published as a conference paper at ICLR 2024
Embeddings The length and thickness of the edges encoded using Gaussian embeddings with bases: GL. G, R + R They are then concatenated and used as edge attributes
Zij -e----llll
@ ff te-
(11)
where denotes concatenation and (Te. pe) are a collection of fixed parameters of the Gaussians and they depend on the number of bases. The edge vectors are expanded in a spherical basis up to Timaz Unlike atoms of various elements in chemistry, all our nodes are of the same type. There- fore, the node features are initialized as ones and are expanded to the desired number of dimensions
(0) using a linear layer: h Wi.
MACE layer At each layer s of MACE, edge embedding (x) Øif are formed by taking the tensor product between the node features h2 and the edge vectors expanded in spherical basis. This tensor product is weighted with non-linear function of the edge attributes Zij.
Clama milimggg laata YI (Lig)) j.kamm
R
(12)
Where & indexes feature channels, 1, m index angular momenta, Clymam Jim lamo are the Clebsch- Gordan coefficients that enforce the equivariancd and indexes combinations of Im which preserve equivariance. The Atomic Basis A of the node at layer is constructed by summing edge embeddings over the edges i:
Ew(el
Σ
(+) 25.kmplams
(13)
klepdaa la
JEN(()
tensor product applied to the Atomic Basis v times to increase the body order of the feature, and the resulting features are symmetrized using set generalized Clebsch-Gordae coefficients 7, Um
B
( i.nkIM
Σ
ECLM
II
(a) i,kigme
(14)
E=1
where B( MM are called sketched product basis. The B B-features then linearly mixed to form many-body message,
(a) mI TELLEM = ΣΣ (1 ()), senck K..LLM
(15)
Finally the message is used to update the next node features using an update function.
equivariant non-linear readout followed by global graph pooling. Invariance to tessellation is maintained by mean pooling operation which ensures that the predicted stiffness will not grow if nodes with identical neighborhoods graph. [6] Finally, another linear layer outputs two scalars, two = vectors and one / = 4 vector, corresponding the spherical component of a 4th order tensor with the correct permutation symme- try. This is converted to Cartesian basis and assembled into Mandel notation to form the final matrix output A (see the Appendix for details).
Readout After S layers of MACE,
added to the
Positive Semi-Definite Layer The positive semi-definite layer is the key step to ensure that the final output is positive semi-definite, line wih the law of energy conservation. We evaluate number of methods to make the stiffness tensor positive semi-definite. These include taking even powers of the matrix, 2 and A4, matrix exponential, and truncated versions, e (I + A/2)², (I + A/4)*.
$ For further details, see the original MACE paper by Batatia et al. 20022 TThis is another fundamental inductive hias of our model: the outputs for unit cell and its X tessellation are identical.
22

---

Published as a conference paper at ICLR 2024
A.11 TRAINING AND EVALUATION DETAILS
Definitions of loss metries If C, - Cpij and Cp Cpsj are the predicted and target stiffnesses for lattice (in Mandel notation), respectively, and Cp Cpijki and C Cpijkl the predicted and target stiffnesses for lattice p (in 4th order notation), respectively, the loss used during training calculated as:
1 7pp 36
Σ 3/
CpeyC
Leomp.p = =E (Cris -Cpts) )² Lurain
Σ P
Lempp.P
(16)
B
7p
where denotes the total number of lattices p. 7pp the mean stiffness, Leomp,p the component loss and Luain the overall loss. For testing purposes, in addition to loss Leamp. we define directional loss, Liis and relative directional loss, Lirree, which are calculated using N 250 random directions on the unit sphere (dg. 9 = I...N):
pprr,p
1 N
Σ
Σ iik!
(Crijkt - Epishaa durdesdaaedd
Lditzel.p Ldir.p Tp
(17)
with Lpir,p the directional loss and Lairpp the relative directional loss. We further report the percent- age of lattices with negative eigenvalues, Xx. and equivariance loss, Lequiv which is calculated as
follows. We choose 5 random orientations (here S 10) parameterized by corresponding rotation
matrices ((*) R (s = S). The predicted stiffiess tensor the original orientation is C (p) For each lattice in the test dataset C(P), calculate the predictions for each of the 10 rotations: C(p.a) Equivariance loss is defined as "jkll"
@
(((CT)
Lequiv
Σ
Σ
R((s) Gj
(p.s) Gjkl
SNB
pqs RgRI
Note that the equivariance loss is caleulated purely from predictions, disregarding the mismatch from the ground truth.
All models were trained single NVIDIA A100 GPU with 80GB of memory. The training routines were handled by Pytorch Lightning. Specifics vary between models and are outlined below.
A.11.1 CGC AND MCGC
Hyperparameters were searched on a grid (Table Every experiment was run with constant Iearn- ing rate for up to 100 000 steps. Optimizer AdamW used with settings (31, (22 (0.9, 0.999), == x 10-3, weight decay=1 x 10-8 Validation loss was checked every 100 steps and training was stopped by early stopping callback with patience 50 if validation loss has not reduced.
The setup of the vamilla CGC and mCGC followed as closely as possible the methods used by Meyer al. (2022). Key differences between training routines for CGC-based models and the vanilla versions are as follows.
Static 7-fold data augmentation was used for training set. Each lattice was rotated by 90 deg around the r-, y- and z axes, and mirrored about the y, x and y planes. This added six more versions of the lattice into the dataset.
The components of stiffness tensor were independently normalized to lie between 0 and
Optimizer RAdam was used
Loss smooth_11 was used
A.12 NNCONV
In model NNConv, linear learnable layers were used as node and edge feature embeddings to in- crease latent dimensionality. The message passing NNCony layer composed of 3-yeyer neural
23

---

Published as a conference paper at ICLR 2024
Table 5: Hyperparameter search for CGC and mCGC models
ndim
16, 32, 64, 128, 256, 512
message passes
2, 3,4,4 6 min, max, 256
aggregation
batch size
Ir
0.0003, 0.001, 0.003, 0.01
# parameters
16K - 10M
Table Hyperparameter search for NNConv models
ndim
16, 32, 64
message passes
3,4,
aggregation
min, max, mean
batch size
256
Ir
0.0003, 0.001, 0.003, 0.01
# parameters
22K 1.6M
network with ReLU nonlinearity. "Sum" aggregation was used to aggregate messages from neigh- bors, after which ReLU layer was applied. The layers of message passing were not shared, but inde- pendent, A residual connection between layers was used. Hyperparameters were searched (Table ( Every experiment was run with constant learning rate for up to 100 000 steps. Optimizer AdamW W was used with settings (91, 32)
grid
(0.9, 0.999), Validation loss was checked every 100 steps and training was stopped by early stopping callback with patience 50 if validation loss has not reduced.
10 weight decay=1 x 10
A.12.1 MACE
Hyperparameters for MACE-based models were searched grid in Table 7 Every experiment was run with constant learning rate for to 000 steps. Optimizer AdamW was used with settings (31, (2) (0.9, 0.999),
10-8, 8 weight decay=1 10-8. smaller batch size of 64 was used because of the higher memory requirements of the MACE model. To maintain consistency with the batch size of 256 from CGC models, gradient accumulation over batches was used. Validation loss was checked every 100 steps and training was stopped by early stopping callback with patience if validation loss has not reduced. Value-based gradient clipping was used with cutoff 10.0.
A. PRIMAL, DUAL AND COMBINED GRAPHS
The choice of graph over which message passing is run is important. For instance, some studies in the mechanics community use the dual graph where the centres of lattice cells are converted to nodes, and the neighbouring cells connected by graph edges. (Karapiprasi & Kochmann 2023)
Meyer et al. 2022) combine the primal graph, with a line graph. The original (primal) graph can be converted to line graph follows. The nodes of line graph are the edges of primal graph. Two
Table 7: Hyperparameter search for MACE models
hidden dim readout dim
8, 16, 32, 64
8, 16. 32
message passes aggregation batch size
0.0003, 0.001, 0.003
parameters
50K - 600K
24

---

Published as a conference paper at ICLR 2024
Table Performance of CGC model for primal, dual, and combined graph representation lattice unit cells
Static augmentation
Dynamic augmentation primal CGC
primal CGC
dual
combined mCGC
dual
combined
mCGC
Loamp
7.81
11.90
8.49
4.63
9.10
5.36
Ldie
8.31
15.46
8.71
5.29
12.03
5.84
Laiegrl
0.38
1.14
0.42
0.25
0.76
0.27
Ag
10
1
26
28
b
Z
y
Leus
Loue
Hwai:: Ground truth
Figure 7: (a) Unit cell and (b) the true stiffness surface of simple cubic lattice. (c) Projection into x y plane and comparison of models with various Lmaa-
nodes in the line graph are connected if the two corresponding edges of the primal graph meet at node.
Here investigate the performance of GNN based on the choice of graph passing is done. Table & shows the results for various models and training strategies. Primal and combined correspond to models CGCNN and mCGCNN from Meyer al. (2022). Dual message passing done purely on the line graph. Static augmentation corresponds to the training routine from Meyer et al. (2022) as described Section A.III.] Dynamic augmentation corresponds our training routine whereby each time a lattice is retrieved from the dataset, it is obtained at a different orientation.
which message
In summary, we empirically do not see any benefit of incorporating the line graph into our model. Therefore, we do consider these models in the main text.
A.14 TRAINING WITH Lmaz AND DEGENERACY OF HIGHLY-SYMMETRIC LATTICES
In Figure 7 we show the unit cell of the simple cubic lattice. The lattice has high degree of sym- metry which has profound consequences for message passing. If the maximum degree spherical expansion, Immazs inside the model is lower than 4, the model is restricted to fitting isotropic stiffness tensor for this lattice. When the Lmaz 2 4, the anisotropy can be captured.
Note that when the model is trained with Lanax 4, it still needs output fourth-order tensor whose spherical form includes L, = 4 component: x Oe + 2 x 2e - x 4e. [ We enable this by incorporating a tensor product expansion layer. Suppose the message passing done up to Lmaz = 2. After message passing and graph pooling, each graph has features of the form N Oe + x + N 2c, where is the number of channels (N chosen even). We split the channels into two sets of size m = N/2 and do tensor product between them:
[m mm+mm 1o+ 2e] & (m m x 1o+ x 2e] + [33m² x 0e (4m²) 2e (m²) x 4e]
The output is passed through a linear layer with learnable weights which reduces it to the correct dimensionality 2 x 0e 2 2e 4e.
Thhe notation used here is specific the software implementation of e3nn Geiger
2022
25


Reviews:
[{"review_id": "ANQUYzXtvl", "reviewer": "Reviewer_azcZ", "summary": "This paper proposes an application of equivariant GNNs for predicting the stiffness tensor of architected lattice metamaterials. To ensure the validity of physical significance, a layer is proposed to preserve the positive semi-definiteness of the predicted stiffness.", "questions": "Please refer to the weakness.", "limitations": "", "rating": 5, "confidence": 4, "soundness": 2, "presentation": 2, "contribution": 2, "strengths": "1.\tThere is less work focused on studying higher-order tensors beyond first-order tensors, such as coordinates, velocity, and force. Research on a 4th-order tensor represents a new situation for application.\n2.\tThis paper try to design a new module to ensure the validity of physical significance, which is inspiring. \n3.\tThe figures are pretty.", "weaknesses": "1.\tThis work chooses to predict a 4-th ordered tensor – “stiffness tensor”, whose symmetry has been maintain well with existing models like TFN, PAINN. This paper does not propose theoretical innovations in maintaining equivariance, so this methodology and the emphasis on “fourth-order” in abstract are not directly related (the Positive Semi-Definite Layer is based on physical meaning, regardless of whether it is a fourth-order tensor).\n2.\tSince the whole backbone is built with existing MACE layer, the main novelty is to propose the  “Positive Semi-Definite Layer”(let us call it “PSD-layer”). About this, here are some questions:\n     - a) To maintain the equivariance of the total model, we must ensure each layer in such a model is equivariant. But for the PSD-layer, the PSD matrix $A$ is based on $M$, which is created by arrange $n(n+1)/2$ entries of the output of an equivariant model.\n        - **i**. Is it an equivariant operator to arrange $n(n+1)/2$ entries to get $M$?  This requires a formal mathematical proof.\n        - **ii**. Is it an equivalent operator of the function $\\rho$ to turn $M$ into $A$? In the end of Section A.2.2, “If the overall model had been equivariant with respect to vectors in $U$, it will remain equivariant after eigenvalues are made positive” is not trivail, it may also require a formal mathematical proof.\n         - **iii**.\tFor now, we assume that the previous question (function $\\rho$ is equivariant) has a good proof. It is necessary to discuss whether this design will reduce the representational capacity of the entire neural network. Assume that this neural network is a bijection (we consider an equivalence class divided based on different perspectives as an input. But the function $\\rho$ may not a bijection, e.g. $(\\pm \\Lambda) ^2$ will get a same result. The Wigner-D based network build a faith representation of the group, but the induction of function $rho$ will lead to an unfaith representation, it requires more analysis. For example,  an analysis for the increased stiffness-based errors of CGC+ve, NNConv+ve, and MACE+ve in Table 1 may be a manifestation of the decline in the ability of network representation.\n    - b)\tLet assume the PSD-layer is equivariant, there are some questions about the experiments:\n        - **i**.\tIn Section 3 (Related Work) and Section 5(Conclusion), Finite element (FE) modelling is mentioned, but the experiment could not find the data comparison (e.g. time cost or accuracy) between your model and FE.\n        - **ii**.\tMore powerful baselines should be discussed, especially the geometric GNNs like TFN[1], NequIP [2], SCN[3], GMN[4], SE(3)-Transformer [5]. \n3.\tHere are some possible typos and recommended symbol modifications:\n    - a)\tAlmost all the inline formulas in the article lack punctuation at the end.\n    - b)\tIn Section 2.1, the rotation of the stiffness tensor $R_{ia}R_{jb}C_{abcd}R_{kc}R_{ld}$ may be better to generalize to $n$-th order tensors with the format of changing bases of a tensor($ R_{ia}R_{jb} R_{kc}R_{ld} C_{abcd}$).\n    - c)\tIn Eq. (2),  the Wigner-D matrix generally written in the form $D_{m’,m}^{(l)}(R)$.\n    - d)\tIn Section 4, the third line in the paragraph “Positive Semi-Definite Layer”, “in line wih”, may be “in line with”.\n\nReference:\n\n- [1] Tensor field networks: Rotation-and translation-equivariant neural networks for 3d point clouds. \n- [2] E(3)-equivariant graph neural networks for data-efficient and accurate interatomic potentials\n- [3] Spherical Channels for Modeling Atomic Interactions\n- [4] Equivariant Graph Mechanics Networks with Constraints\n- [5] SE(3)-Transformers: 3D Roto-Translation Equivariant Attention Networks", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "GnQvdg7ax7", "reviewer": "Reviewer_kbo4", "summary": "This paper introduces higher-order energy-conserving SE(3) equivariant GNNs which build upon the MACE architecture. These GNNs are applied to lattices, i.e. architected metamaterials. The new features of the model are conservation of energy and SE(3) equivariant predictions of the a 4th order stiffness tensor.  For the stiffness tensor positive-definiteness is ensured. Experiments compare against existing deep learning models of Crystal Graph Convolution (CGC) and NNConv.", "questions": "- When exactly is the positive semi-definite layer applied? after the readout?\n- In Figure 2 validation loss of CGC+tr seems to be still dropping sharply, is it possible that training was stopped too early?", "limitations": "", "rating": 8, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "- This paper introduces a novel and very interesting application of SE(3) equivariant GNNs.\n- Positive semi-definiteness of stiffness tensor is a interesting new tool for physics-based machine learning.\n- The paper contains all important information to follow (general physics background, no material-design expert).\n- The evaluation scheme is quite solid and the different training methods are sound.", "weaknesses": "- It seems that most of the important information which is novel to the deep learning community is put into the appendix. For example, definitions of stress and strain tensors, or their relation to the stiffness tensor. Even more importantly, only in the appendix one can read why the stiffness tensor can be represented as a matrix, and how the positive definiteness of a matrix can be ensured. \n- On the other hand, background and related work are a bit repetitive\n- Code / or pseudocode would be pretty helpful to understand output layers and how the stiffness tensor is ensembled.\n- There are no runtime comparisons in the paper, especially since e3nn based models are known to be slow this would be interesting to know.\n- Related to runtime comparisons, comparisons to FE methods are needed to get a better understanding of the presented performances. I understand that FE models are used as ground truth and that this might be tricky to obtain, but for example one could estimate runtime comparisons (especially since it is stated the FE methods are very slow) and report performance differences for nodal perturbations to get a perspective for the reported loss values? The latter should give a feeling to what nodal perturbations the presented losses are comparable. \n- As far as I can see there is no OOD experiments although this is claimed to be one of the main reasons why such models built on physical principles are built.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "7a7P2dOxd8", "reviewer": "Reviewer_TDG7", "summary": "This paper introduces a specially design graph neural network that well tackles the physical constraints in lattice architected metamaterials. To be specific, the model guarantees by design SE(3)-equivariance and energy conservation. The latter is fulfilled by the proposed positive semi-definite layer that overcomes the limits of previously proposed approaches using Cholesky factorization and eigenvalue decomposition. The model is examined on a dataset ``limp'' showcasing that MACE is superior when equipped with the proposed PSD layer.\n\nUpdate: contribution raised from fair to good, overall score raised from 3 to 5 due to additional experiments on material design.", "questions": "1. How does the method perform compared with other solvers, deep learning-based or even not?\n\n2. Are there other tasks/datasets that can benefit from the development of such GNN-based method that is specifically designed to regress on the stiffness?\n\n3. Even if the paper discusses the limitations of previous Cholesky/eigenvalue decomposition-based methods that permit PSD constraint, it would help a lot if this point is also verified by experimental results/ablation studies.", "limitations": "", "rating": 5, "confidence": 4, "soundness": 3, "presentation": 2, "contribution": 3, "strengths": "1. The motivation of the paper is clear: designing a GNN that fulfills the inductive biases for lattice architected metamaterials.\n\n2. The proposed approach is clean and easy to follow.", "weaknesses": "1. The novelty is somewhat limited in the sense that the proposed model is a direct utilization of MACE, equipped with the proposed layer that satisfies PSD.\n\n2. The dataset is a bit limited and the experiments lack of comparison with other methods. If there are no significant deep learning based (or even GNN-based) methods to compare, it would also help a lot if any FE-based method could be involved. The readers may be curious about how the proposed model can achieve better accuracy/efficiency tradeoff compared with traditional solvers.\n\n3. The implication of the practical usage of the method is limited. Since this paper is developed purely based on practical considerations, it would be better if the paper could show how this model can help in other related/similar tasks or downstream tasks that could potentially benefit from this method., other than just predict the stiffness tensor.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "fcZI13hgnd", "reviewer": "Reviewer_uV3y", "summary": "## Summary\n\nThe paper presents an extension to the [MACE](https://proceedings.neurips.cc/paper_files/paper/2022/file/4a36c3c51af11ed9f34615b81edb5bbc-Paper-Conference.pdf) model by adding a matrix power layer to ensure the positive semidefinite (PSD) nature of the output tensor. Although the idea of maintaining PSD properties might be of value, the paper's contributions are incremental at best. Most notably, the claim about \"high-order rotational equivariance\" is not an original contribution but belongs to the original MACE model. The paper severely suffers from a lack of organization and clarity in writing. Hence, I recommend a rejection of this submission.\n\n## Detailed Comments\n\n### 1. Lack of Definitions and Citations\n\na) The term \"edge shifts\" is used without a definition or reference.\nb) The suffix 'lb' used in the context of training methods is not defined.\nc) Clebsch–Gordan coefficients are mentioned without proper citation.\nd) On page 17, there are numbers in parentheses next to \"Young's modulus,\" \"shear modulus,\" and \"Poisson's ratio,\" without any explanation.\n\n### 2. Unclear Illustrations and Descriptions\n\na) The paper claims that energy conservation equates to a PSD tensor, but this is tucked away in a footnote. This claim needs to be explained in the main text.\nb) The choice of the optimal 've' method, which is squaring matrix A, is not mentioned until late in the results section. This should be mentioned in the methodology section.\nc) Only absolute errors are presented without giving the context of ground truth magnitudes or relative errors. It's unclear for the readers to know if this method is accurate or not at all.\nd) It's unclear if Figure 1b is a plot for ground truth or predictions.\n\n### Flaws in Finite Element Method (FEM) Introduction\n\na) The paper mentions that FEM has \"~10^9 elements\" but does not state the computation time for such a scale. Furthermore, GNNs cannot handle such scales on a single GPU, making the comparison unfair. Instead, try to report a wall time comparison between your method and FEM for the same dataset.\nb) The statement that FEM ensures force equilibrium is incorrect; it is the underlying PDE that ensures this.\nc) Similarly, FEM itself does not ensure PSD properties; this is ensured by the constitutive model.\nd) The paper incorrectly claims that FEM is rotationally equivariant. Special treatments are needed to achieve rotational equivariance in FEM.\n\n## Conclusion\n\nThe paper presents an incremental extension to the MACE model with a focus on preserving the PSD properties of the output tensor. However, the paper lacks clarity in writing and organization, and its contributions are limited. As a result, I recommend a borderline reject for this submission.\n\nThe authors can also consider a workshop or journal submission that focuses on this area.\n\n##\nAfter rebuttal, I think the writing quality and the motivation has been more clear. Hence, decided to improve to 6.", "questions": "See above", "limitations": "", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "See above", "weaknesses": "See above", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Decision: Accept (poster)