An equivariant graph neural network for the elasticity tensors of all
  seven crystal systems
The elasticity tensor that describes the elastic response of a material to
external forces is among the most fundamental properties of materials. The
availability of full elasticity tensors for inorganic crystalline compounds,
however, is limited due to experimental and computational challenges. Here, we
report the materials tensor (MatTen) model for rapid and accurate estimation of
the full fourth-rank elasticity tensors of crystals. Based on equivariant graph
neural networks, MatTen satisfies the two essential requirements for elasticity
tensors: independence of the frame of reference and preservation of material
symmetry. Consequently, it provides a unified treatment of elasticity tensors
for all seven crystal systems across diverse chemical spaces, without the need
to deal with each separately.. MatTen was trained on a dataset of
first-principles elasticity tensors garnered by the Materials Project over the
past several years (we are releasing the data herein) and has broad
applications in predicting the isotropic elastic properties of polycrystalline
materials, examining the anisotropic behavior of single crystals, and
discovering new materials with exceptional mechanical properties. Using MatTen,
we have discovered a hundred new crystals with extremely large maximum
directional Young's modulus and eleven polymorphs of elemental cubic metals
with unconventional spatial orientation of Young's modulus.
Mingjian Wen, Matthew K. Horton, Jason M. Munro, Patrick Huck, Kristin A. Persson
# A Universal Equivariant Graph Neural Network for the Elasticity Tensors of Any Crystal System

###### Abstract

The elasticity tensor that describes the elastic response of a material to external forces is among the most fundamental properties of materials. The availability of full elasticity tensors for inorganic crystalline compounds, however, is limited due to experimental and computational challenges. Here, we report the materials tensor (MatTen) model for rapid and accurate estimation of the full fourth-rank elasticity tensors of crystals. Based on equivariant graph neural networks, MatTen satisfies the two essential requirements for elasticity tensors: independence of the frame of reference and preservation of material symmetry. Consequently, it provides a universal treatment of elasticity tensors for all crystal systems across diverse chemical spaces. MatTen was trained on a dataset of first-principles elasticity tensors garnered by the Materials Project over the past several years (we are releasing the data herein) and has broad applications in predicting the isotropic elastic properties of polycrystalline materials, examining the anisotropic behavior of single crystals, and discovering new materials with exceptional mechanical properties. Using MatTen, we have discovered a hundred new crystals with extremely large maximum directional Young's modulus and eleven polymorphs of elemental cubic metals with unconventional spatial orientation of Young's modulus.

## 1 Introduction

All materials exhibit elastic behavior under small external loads and return to their original shape upon the release of these loads [1]. The elasticity tensor provides a fundamental and complete description of a material's response to such loads. It offers a lens to examine the inherent strength of the bonding in a material and enables the understanding, analyzing, and designing of many macroscopic physical properties of materials. Commonly employed mechanical properties (for instance, Young's modulus and Poisson's ratio), thermal properties (for instance, thermal conductivity) and acoustic properties (for instance, sound velocity) can be mathematically derived from the elasticity tensor. These properties have long been leveraged, for example, by materials scientists to search for ultrahard materials [2, 3] and by geophysicists to interpret seismic data [4, 5]. More recently, the anisotropic elastic behavior of inorganic solid electrolytes has been found to play a decisive role in determining the stability of electrodeposition at the interfaces in solid-state batteries [6, 7]. Moreover, in solid-state synthesis, one would want to utilize the elasticity tensor to determine the local stability of a material, so as to avoid unsuccessful synthesis of materials that spontaneously transform into different structures [8, 9].
In spite of the importance, elasticity tensor data from experimental measurement is limited to a small number of materials. For example, for inorganic crystalline compounds, experimental data is only on the order of a few hundred, considering entries both tabulated in handbooks and scattered in individual papers [10]. The major difficulty lies in preparing large enough single crystals for accurate experimental measurement using current techniques such as resonant acoustic spectroscopy [11]. In the past decade, efficient and reliable electronic structure calculation methods such as density functional theory (DFT) [12] with automation frameworks [13, 14, 15] have enabled high-throughput computational investigation of materials. Using this approach in the Materials Project [16], we produced an initial dataset of elasticity tensors for 1181 crystals in 2015 [10], which has been expanded over time to 10276, which we now release as a new dataset with this work. Nevertheless, this only accounts for 6.6% of the more than 154000 crystals in the Materials Project, let along the even greater number of crystals recorded in the Inorganic Crystal Structure Database (ICSD) [17] and predicted by universal interatomic potentials [18].

It is therefore no surprise that machine learning (ML) has gathered substantial interest as a means to develop efficient surrogate models for the prediction of elastic properties. In a nutshell, state-of-the-art ML models for elastic properties encode compositional information [19, 20, 21] and/or structural information [20, 21, 22, 23] in a material as feature vectors and then map them to a target using some regression algorithms. This approach is limited to derived scalar elastic properties such as bulk modulus and shear modulus, and separate models have to be built for each derived property. Ideally, one would hope to predict the full elasticity tensor, from which all other derived elastic properties can be obtained. Along this direction, there have been attempts to predict individual tensor components [24, 25]. These models are great first steps, but essentially they still predict separate scalars in a specific coordinate system and thus unavoidably violate the transformation rules of tensors. Geometric deep learning [26] such as equivariant graph neural networks (GNNs) [27, 28, 29, 30] can directly operate in the space of high-rank tensors and adhere to their transformation rules. Although the major use case is still for scalar molecular and materials properties, a handful of works have explored the application in predicting vectors such as the molecular dipole moment [31, 32], demonstrating the viability of direct deep learning of tensorial properties.

In this work, we develop the Materials Tensor (MatTen) model for rapid and accurate estimation of the fourth-rank elasticity tensors of inorganic compounds. Our model, based on equivariant GNNs, takes a crystal structure as input and outputs its full elasticity tensor with all transformation rules satisfied. Other elastic properties such as bulk modulus and shear modulus can then be computed from the full elasticity tensor. The capabilities of MatTen are demonstrated via the study of both isotropic and anisotropic elastic properties. Using MatTen, we screened the Materials Project database for the discovery of new materials with large maximum directional Young's modulus. On average, the values of the newly found materials are more than three times larger than existing data, as verified by first-principles calculations. In addition, we have identified 11 unconventional polymorphs of elemental cubic metals whose maximum directional Young's moduli are in the \(\langle 100\rangle\) directions.

## 2 Results

### Symmetry and irreducible representation of the elasticity tensor

The elasticity tensor \(\mathbf{C}\) is a fourth-rank tensor that fully characterizes the elastic behavior of a material. Given that it is the second derivative of the total elastic energy with respect to the strain tensor and that the strain tensor is symmetric [33, 34], the elasticity tensor possesses major symmetry \(C_{ijkl}=C_{klij}\) and minor symmetry \(C_{ijkl}=C_{jikl}=C_{ijlk}\) (in indicial notation, where \(i,j,k,l\in\{1,2,3\}\)). Consequently, only 21 of the 81 components of \(\mathbf{C}\) are independent. It is convenient to write the elasticity tensor in a contracted matrix \(\mathbf{c}\) (\(c_{\alpha\beta}\), where \(\alpha,\beta\in\{1,2,\ldots,6\}\)) with a pairs of indices \(ij\) in the tensor notation replaced with a single index \(\alpha\) in the matrix notation: \(11\to 1\); \(22\to 2\); \(33\to 3\); \(23,32\to 4\); \(13,31\to 5\); and \(12,21\to 6\). This Voigt matrix
[35] is a \(6\times 6\) matrix symmetric about the main diagonal, reflecting the fact that the elasticity tensor has 21 independent components.

The intrinsic material symmetry of a crystal can further reduce the number of independent components [8, 34]. For example, copper is a cubic crystal with mirror planes and three-fold rotation axes (point group m\(\bar{3}\)m), and such symmetry results in a number of only three independent components. Formally, the material symmetry imposes the following constraints on the elasticity tensor [36, 37]:

\[C_{ijkl}=Q_{ip}Q_{jq}Q_{kr}Q_{ls}C_{pqrs}, \tag{1}\]

where \(\mathbf{Q}\in G\subset SO(3)\), and \(G\) is the material symmetry group of the crystal, which is a subgroup of the rotation group \(SO(3)\). An interesting question is: how many unique symmetry classes exist and what is the number of independent components in each class?

It turns out that there exists only eight distinct classes (Fig. 1), proved via a purely algebraic approach by directly identifying the equivalence classes corresponding to Eq. (1) [36]. Of the eight classes, one is for isotropic materials, and each of the other seven corresponds to a crystal system [38, 37]. In our opinion, there is still significant confusion on this topic. For example, the categorization by Wallace [39] and populated by Nye [34], which incorrectly gives two unique classes for each of the tetragonal and trigonal cases (Fig. S1 in the Supplementary Information (SI)), is still widely cited in recent works [40, 41, 42]. We refer to Ref. [37] for a historical note on the development of the categorization.

The Voigt matrix provides a visual way to observe the material symmetry of a crystal reflected in the elas

Figure 1: Symmetry classes and independent components of the elasticity tensor. The value in the parentheses after the name indicates the number of independent components. All matrices are symmetric about the main diagonal, with the lower triangular part omitted in the depiction. For single crystals considered in this work, the isotropic case does not apply.
ticity tensor (Fig. 1). The values of the matrix components, however, depend on the choice of the coordinate system and do not show any systematic pattern upon coordinate transformation [43, 44], making it difficult to build predictive models for elasticity tensors. This can be overcome by the _harmonic decomposition_[45], where the space of all elasticity tensors is factored into the direct sum of irreducible representations of \(SO(3)\). Consequently, any elasticity tensor can be written in the form,

\[\mathbf{C}=h_{1}(\lambda)+h_{2}(\eta)+h_{3}(\mathbf{A})+h_{4}(\mathbf{B})+h_{5}(\mathbf{H}), \tag{2}\]

where \(\lambda\) and \(\eta\) (scalars) are the _isotropic_ part, \(\mathbf{A}\) and \(\mathbf{B}\) (second-rank traceless symmetric tensors) are the _deviatoric_ part, and \(\mathbf{H}\) (fourth-rank fully symmetric tensor with a zero trace) is the _harmonic_ part [45]. The harmonic decomposition has two advantageous characteristics. First, for a given \(\mathbf{C}\), the values of \(\lambda,\eta,\mathbf{A},\mathbf{B},\mathbf{H}\) and the functions \(h_{1},\dots,h_{5}\) are uniquely determined [36, 45] (see SI for their expressions). Second, each part in Eq. (2) transforms in a known manner with respect to \(SO(3)\) operations, enabling the construction of predictive models that leverage recent advancements in geometric deep learning.

### Equivariant graph neural networks for high-rank tensors

The MatTen model captures the structure-property relationship of crystalline materials. It takes a crystal structure as input, represents it as a three-dimensional crystal graph, performs equivariant feature update on the crystal graph, and finally outputs a tensor property of the material built according to the reverse process of the harmonic decomposition in Eq. (2).

In the GNN model (Fig. 2), the input crystal is represented as a graph \(\mathcal{G}(V,E)\), with atoms as the nodes and bonds as the edges. The feature \(F_{i}\in V\) characterizes atom \(i\), and \(\vec{\mathbf{r}}_{ij}\in E\) denotes the distance vector between a pair of atoms \(i\) and \(j\). The embedding modules operate on the graph representation of the crystal. The atom embedding module encodes the atomic number \(Z_{i}\) using a one-hot scheme, which becomes the initial atom feature \(F_{i}\). The distance vector \(\vec{\mathbf{r}}_{ij}\) is separated into two parts: the unit vector \(\hat{\mathbf{r}}_{ij}\) from atom \(i\) to atom \(j\) and the scalar distance \(r_{ij}\) between them. The former is expanded on real spherical harmonics \(Y^{l}_{m}\), and the latter is expanded on the Bessel radial basis functions [46]. These embedding modules extract structural information from the crystal and provide them to the interaction blocks.

The interaction blocks iteratively refine the atom features via convolution operations. The architecture of the interaction block follows the design of Tensor Field Network [27] and NequIP [29]. Unlike many existing GNNs for molecules and crystals [47, 48, 49, 22] that utilize scalar features, here, the atom feature \(F_{i}\) is a set of scalars, vectors, and high-rank tensors. Formally, it is a geometric object consisting of a direct sum of irreducible representations of the \(SO(3)\) rotation group [27, 29]. There are two major benefits of using such geometric features. First, they are incorporated as inductive bias which can improve model accuracy and reduce the amount of training data. Second, from them one can easily construct other physical tensors such as the elasticity tensor in this work.

The convolution on these geometric objects is an equivariant function, meaning that if the input feature \(F_{i}\) to the convolution is transformed under a rotation in \(SO(3)\), the output is transformed accordingly. This is achieved via the tensor product convolution,

\[F_{i}=\sum_{j\in\mathcal{N}_{i}}R(r_{ij})Y(\hat{\mathbf{r}}_{ij})\otimes F_{j}, \tag{3}\]

where \(\mathcal{N}_{i}\) denotes the set of neighboring atoms for atom \(i\), \(R\) indicates a multilayer perceptron (MLP) on the radial basis expansion of \(r_{ij}\), and \(Y\) indicates the spherical harmonics expansion of \(\hat{\mathbf{r}}_{ij}\). The tensor product \(\otimes\) between \(Y\) and \(F_{j}\) ensures that the operation is equivariant [27], and the use an MLP makes the convolution learnable. After a skip connection [50], the feature \(F_{i}\) is passed through a nonlinear activation function and finally normalized using an equivariant normalization function [51].
The output head maps the refined features from the interaction blocks to the target materials tensor of interest. First, the features of all atoms are aggregated to obtain a representation of the crystal. For intensive properties such as the elasticity tensor, we adopt the mean pooling by averaging the features such that the representation of the crystal is independent of the number of atoms. Next, an appropriate subset of the pooled irreducible representations that correspond to the target tensor of interest is selected and then assembled as the model output. For the elasticity tensor, the selected ones are two scalars, two second-rank traceless symmetric tensors, and a fourth-rank harmonic tensor. They are assembled to an elasticity tensor according to Eq. (2).

Overall, MatTen is a function \(\mathbf{C}=f(x)\) that maps a crystal structure \(x\) to its elasticity tensor \(\mathbf{C}\). The function \(f\) is equivariant to the \(SO(3)\) group transformation, that is, for any \(x\) and \(g\in SO(3)\), we have \(D_{y}(g)f(x)=f(D_{x}(g)x)\), where \(D_{x}(g)\) and \(D_{y}(g)\) are rotation matrices parameterized by \(g\) for the crystal structure and the elasticity tensor, respectively. The equivariance ensures that the model can produce an elasticity tensor \(\mathbf{C}\) that respects the orientation of the input crystal structure. In other words, the choice of a

Figure 2: Schematic overview of the MatTen graph neural network model. (a) The model takes a crystal structure as input, performs message passing with equivariant graph neural network (GNN) layers, and outputs the full elasticity tensor satisfying all symmetry requirements. (b) Architecture of the model. The model consists of four modules: the position embedding module converts an input displacement vector \(\vec{\mathbf{r}}_{ij}\) between atoms \(i\) and \(j\) to radial (\(R\)) and spherical (\(Y\)) expansions; the atom embedding module encodes the atomic number \(Z_{i}\) as irreducible atom features \(F_{i}\) using a one-hot encoding; interaction blocks iteratively refine the features using convolutions based on tensor product; and the output head selects the appropriate irreducible features corresponding to the elasticity tensor and assembles them as the output tensor \(\mathbf{C}\). Following Ref. [44], the fourth-rank harmonic part of the elasticity tensor is depicted as a generic matrix. The \(\oplus\) denotes addition.
specific coordinate system does not affect the model. This _independence of the frame of reference_ characteristic is an indispensable property for models that predict tensors. In addition, any such model should also _preserve the material symmetry_ of the crystal. By construction, MatTen guarantees the material symmetry reflected in the elasticity tensor. Concretely, if the predicted elasticity tensor is represented as a Voigt matrix, the symmetry and number of independent components in Fig. 1 are automatically maintained for any crystal system (proof in the SI).

### Elastic properties of polycrystals

The MatTen model directly outputs the full elasticity tensor. To assess its performance, we computed several commonly used elastic properties for polycrystals from the elasticity tensor. Fig. 3 illustrates the outcomes on bulk modulus, shear modulus, and Young's modulus obtained using the Hill average scheme [52] (Section 4). The predictions of MatTen closely align with the DFT reference values, achieving mean absolute errors (MAEs) of 7.37 GPa, 8.38 GPa, and 20.59 GPa for bulk modulus, shear modulus, and Young's modulus, respectively. For comparison, we trained two additional models. The first is a variant of MatTen, where the tensor output head of MatTen is replaced with a scalar output head, referred to as MatSca hereafter. The second is the AutoMatminer algorithm, an automated machine learning pipeline designed for predicting scalar materials properties [20]. We evaluated their performance in predicting the elastic moduli, and the results are listed in Table 1. Both MatTen and MatSca have smaller MAEs than AutoMatminer across all three moduli, owning to the effectiveness of the underlying equivariant neural networks in learning materials properties from structures. The performance between MatTen and MatSca are comparable. However, it is worth highlighting that while an individual MatSca model was trained for each modulus, a single

Figure 3: Performance of MatTen on various elastic properties. (a) Bulk modulus, (b) shear modulus, and (c) Young’s modulus predicted by MatTen versus DFT reference values. (d) Scaled error by crystal systems. The error bar indicates the standard deviation obtained from an ensemble of five models trained with different initialization. MAE: mean absolute error; MAD: mean absolute deviation.
MatTen model successfully produced all the elastic moduli, demonstrating the versatility and universal applicability of the MatTen model.

Upon closer examination, we have identified some unphysical predictions. Although the theory of elasticity allows negative bulk and Young's moduli in general, under the thermodynamic conditions considered here, all moduli should be positive [53]. The predicted moduli, however, occasionally yield negative values, despite all crystals in the training set having positive moduli. This is an inherent challenge faced by machine learning regression models in general, albeit with physical inductive biases embedded in the model such as the symmetry requirements in MatTen. Nevertheless, this is not a concern in practical use. One can filter out the unphysical ones predicted by MatTen. Moreover, the number of crystals with negative predicted moduli remains minimal, accounting for only 3, 2, and 2 out of the 1021 test data for bulk, shear, and Young's moduli, respectively.

To assess how MatTen performs for different elastic properties as well as for different crystal systems, we computed the scaled error, \(\mathrm{SE}=\mathrm{MAE}/\mathrm{MAD}\), in which MAE and MAD are the mean absolute error and mean absolute deviation, respectively (see Section 4). MAD quantifies the distance of reference values to their mean, and larger MAD means the reference values are more scattered. A model that makes accurate predictions for each data point will have an SE of 0, and a model that always predicts the mean of the dataset will have an SE of exactly 1. Comparing between properties, we see from Fig. 3d that the SE of bulk modulus is smaller than these of shear modulus and Young's modulus across all crystal systems. This suggests that bulk modulus is easier to predict, in agreement with existing observations [19, 20, 21, 22].

Next, we compare between crystal systems. The dataset used for model development has an uneven distribution in terms of the number of materials in each crystal system (Fig. S2 in the SI). For example, it contains 4217 cubic crystals, fewer than 800 trigonal and monoclinic crystals, and only 60 triclinic crystals. As a result, the SE and the error bar of the predicted moduli are larger for trigonal, monoclinic, and triclinic crystals in general (Fig. 3d). Despite the slightly higher errors, it is notable that the model can still perform well for the crystal systems with a low presence in the training data, particularly for triclinic crystals. This is primarily because MatTen internally treats all crystals the same, enabling crystal systems with fewer data to leverage the abundant data from other crystal systems and acquire enhanced representations.

The elastic moduli have values across different orders from near zero to hundreds (Figs. S3-S6 in the SI). To mitigate the challenge of learning values across a broad spectrum, some existing models [19, 20, 21, 22] adopt the approach of predicting the logarithm of the moduli. Unlike these models, MatTen directly predicts the full tensor without any data transformation, and all other elastic properties (including their logarithms) can be computed from it. The logarithms of the bulk, shear, and Young's moduli obtained from MatTen are comparable to those that learn in the logarithmic space (Table S1 in the SI). Further, the performance of MatTen can be further improved with additional training data. The MAE almost decreases linearly with the logarithm of the number of data used to train the model (Fig. S8 in the SI).

\begin{table}
\begin{tabular}{c c c c c c c} \hline \hline  & \multicolumn{2}{c}{\(K\) (GPa)} & \multicolumn{2}{c}{\(G\) (GPa)} & \multicolumn{2}{c}{\(E\) (GPa)} \\ \cline{2-7}  & MAE & \(\mathrm{MAE}/\mathrm{MAD}\) & MAE & \(\mathrm{MAE}/\mathrm{MAD}\) & MAE & \(\mathrm{MAE}/\mathrm{MAD}\) \\ \hline MatTen & 7.37 (0.10) & 0.130 (0.002) & 8.38 (0.16) & 0.280 (0.005) & 20.59 (0.35) & 0.275 (0.005) \\ MatSca & 7.32 (0.09) & 0.129 (0.002) & 8.63 (0.07) & 0.288 (0.002) & 19.87 (0.43) & 0.265 (0.006) \\ AutoMatminer & 9.84 (0.34) & 0.174 (0.006) & 9.27 (0.32) & 0.309 (0.011) & 22.10 (0.77) & 0.295 (0.024) \\ \hline \hline \end{tabular}
\end{table}
Table 1: Performance of the models in predicting the bulk, shear, and Young’s moduli. The value in a pair of parentheses is the standard deviation obtained from an ensemble of five models trained with different initialization. MAE: mean absolute error; MAD: mean absolute deviation.
### Anisotropic elastic properties

Crystals are inherently anisotropic, and thus their elastic properties can vary depending on the direction of measurement. This anisotropy arises from a crystal's structure, including the symmetry of the lattice and the arrangement of the atoms. MatTen predicts the full elasticity tensor, and, for the first time with a machine learning model, we are able to investigate the anisotropic elastic behavior of crystals. We focus our discussion on the directional dependence of Young's modulus (further results on shear modulus given in the SI).

Young's modulus \(E\) discussed in Section 2.3 is an averaged property for isotropic polycrystals. But for single crystals, Young's modulus depends on the direction along which the strain/stress is applied and measured. Given the elasticity tensor \(C_{ijkl}\) (equivalently, the compliance tensor \(S_{ijkl}\), see Section 4), the directional Young's modulus can be computed as [34, 42]

\[E_{\text{d}}(\mathbf{n})=(n_{i}n_{j}n_{k}n_{l}S_{ijlk})^{-1}, \tag{4}\]

where \(\mathbf{n}\) is a unit vector that specifies the direction of measurement. The direction dependence of Young's modulus can be visualized with a three-dimensional plot (Fig. 4a). Interactive visualization can be obtained via, for example, the elate package [54]. Alternatively, via a spherical coordinates transformation: \(\mathbf{n}=[\sin\theta\cdot\cos\varphi,\sin\theta\cdot\sin\varphi,\cos\theta]\) (Fig. 4c), it can be represented in two dimensions (Fig. 4d, with a Robinson map projection [55]). Such plots make it easier to visually investigate the anisotropic characteristics

Figure 4: Directional Young’s modulus \(E_{\text{d}}\). (a) Three-dimensional representation of Young’s modulus predicted by MatTen, for (b) the CaS rocksalt crystal. With (c) a spherical coordinate system and the Robinson map projection, (d) \(E_{\text{d}}\) is represented in two dimensions. (e) Distribution of the prediction error in \(E_{\text{d}}\); also included is the prediction error in the isotropic Young’s modulus from Hill average.
of Young's modulus. For example, for the cubic rocksalt CaS crystal (Fig. 4b), the maximum directional Young's modulus \(E_{\text{d}}^{\text{max}}\) is along the \(\langle 100\rangle\) directions (for instance, \(\theta=90^{\circ}\) and \(\varphi=0^{\circ}\)), while the minimum \(E_{\text{d}}^{\text{min}}\) is along the the \(\langle 111\rangle\) directions (for instance, \(\theta=54.7^{\circ}\) and \(\varphi=45^{\circ}\)). In fact, for cubic crystals such as CaS, the extreme values of \(E_{\text{d}}\) are guaranteed to occur in these two directions [34]. Eq. (4) can be simplified as \(E_{\text{d}}(\mathbf{n})=[S_{1111}-2(S_{1111}-S_{1122}-2S_{2323})(n_{1}^{2}n_{2}^{ 2}+n_{2}^{2}n_{3}^{2}+n_{3}^{2}n_{1}^{2})]^{-1}\) for cubic crystals, expressed in terms of their three independent elasticity tensor components. It can be mathematically shown that, if

\[S_{1111}-S_{1122}-2S_{2323}<0, \tag{5}\]

Young's modulus achieves its maximum \(E_{\text{d}}^{\text{max}}\) in the \(\langle 100\rangle\) directions and minimum in the \(\langle 111\rangle\) directions; otherwise, if \(S_{1111}-S_{1122}-2S_{2323}>0\), the two extremes switch directions (derivation given in the SI). Eq. (5) is satisfied by CaS, and thus we observe the maximum and minimum in the \(\langle 100\rangle\) and \(\langle 100\rangle\) directions, respectively. In addition, \(E_{\text{d}}\) of CaS possesses symmetry (for example, 3-fold rotational axis along the cube diagonals) consistent with a cubic crystal, further confirming that the predicted elasticity tensor _preserves material symmetry_.

To quantitatively assess the ability of MatIen in predicting anisotropic elastic properties, we measured the error between the model predicted directional Young's modulus \(E_{\text{d}}^{\text{pred}}\) and the DFT reference \(E_{\text{d}}^{\text{ref}}\). For CaS, MatIen prediction closely follows DFT reference, with a maximum under-prediction of \(8.7\) GPa along the \(\langle 111\rangle\) directions and a maximum over-prediction of 9.4 GPa along the \(\langle 100\rangle\) directions (Fig. S9 in the SI). In addition to this example crystal, we calculated the error for the entire test set. This is achieved by computing \(D=\frac{1}{M}\int_{\theta}\int_{\varphi}|E_{\text{d}}^{\text{pred}}(\theta, \varphi)-E_{\text{d}}^{\text{ref}}(\theta,\varphi)|\,\text{d}\theta\text{d}\varphi\) for each crystal, where the integrations were performed by uniformly sampling 100 values for each of \(\theta\) and \(\varphi\), and \(M=10000\) is a normalization factor resulted from the integration. The distribution of the prediction error \(D\) of \(E_{\text{d}}\) for the test set is plotted in 4e. It has a Gaussian-like shape, and it almost overlaps with that of isotropic Young's modulus obtained using the Hill average. Similar behaviors are observed for shear modulus and MAEs by crystal systems (Figs. S10 and S11 in the SI). These observations suggest that, on average, MatIen performs equally well for the anisotropic and averaged isotropic elastic properties. Accurate prediction of anisotropic properties offers a comprehensive understanding of the elastic characteristics exhibited by crystals, enabling the discovery of new materials through the utilization of these predictive capabilities.

### Discovering new crystals with extreme properties

The maximum of the directional Young's modulus, \(E_{\text{d}}^{\text{max}}\), characterizes the smallest possible deformation of a crystal under applied external loading. It helps in the selection and orientation of materials to minimize shape change to guarantee the reliability of precision devices such as micro-electromechanical systems [56].

We have also applied the MatIen model to discover new crystals with large \(E_{\text{d}}^{\text{max}}\). The candidate crystals are identified via a screening process. We first filtered crystals from the Materials Project database based on their energy above the convex hull values, selecting those with a value of \(\leq 50\) meV/atom. This energy determines the thermodynamic stability of a crystal and has been shown to correlate with the synthesizability of crystals [57, 58]. We further narrowed down the selection to crystals with fewer than 50 atomic sites to reduce computational cost and remove crystals already present in our dataset. This resulted in 53480 crystals for further analysis. Next, we employed MatTen to predict the full elasticity tensors for these crystals and compute their \(E_{\text{d}}^{\text{max}}\). The top 100 crystals with the highest \(E_{\text{d}}^{\text{max}}\) were then supplied for DFT computation to obtain their elasticity tensors. Fig. 5 presents a histogram of \(E_{\text{d}}^{\text{max}}\) for the newly discovered 100 crystals. Their \(E_{\text{d}}^{\text{max}}\) values all fall at the tail of the distribution of \(E_{\text{d}}^{\text{max}}\) from the training data. Quantitatively, the mean of \(E_{\text{d}}^{\text{max}}\) for the newly discovered crystals is 606 GPa, while that for the training data is 174 GPa, corresponding to the expected value for a randomly selected material. This demonstrates the effectiveness of leveraging MatTen to discover materials with extreme elastic properties.

In addition, we have identified 11 unconventional polymorphs of elemental cubic metals regarding the direction of the extreme values of Young's modulus. Five of them are already experimentally synthesized
(four stable ground-state polymorphs and one metastable polymorph [58]), and the other six are metastable polymorphs that have not been experimentally observed yet. It is believed that \(E_{\text{d}}^{\text{max}}\) is along the \(\langle 111\rangle\) directions (meaning Eq. (5) is not satisfied) for all cubic metals except Mo [34]. We suspect there exists other polymorphs of elemental cubic metals satisfying the criterion in Eq. (5). To test this, we performed a screening of the Materials Project database. From the dataset used for developing MatTen, we have identified six such crystals. Mo is among them, and the other five are V (one polymorph), Cr (two polymorphs), and W (two polymorphs). For crystals not in the dataset, we first used MatTen to predict their elasticity tensors and then selected the 18 crystals satisfying the criterion for verification using DFT. DFT predicted 12 crystals satisfying the criterion; six are elastically unstable structures whose Voigt matrix has negative eigenvalues [8, 9], and the other six successful ones are Mn, Na, K, Cs, Rh, and Tl (one polymorph for each). Among the 11 newly identified crystals, polymorphs of alkali metals Na, K, and Cs have DFT calculated \(S_{1111}-S_{1122}-2S_{2323}\) values much more negative than that of Mo, but they are all hypothetical crystals that have not been experimentally synthesized yet. Five polymorphs are indeed experimental observed, and they are all neighbors of Mo in the periodic table, namely V, Cr, Mn, and W. Among them, four are thermodynamically stable ground-state polymorphs, and one polymorph of Cr is metastable. Crystal structures, ground-state information, and elasticity tensors of the 11 unconventional polymorphs are provided in the SI.

## 3 Discussion

A universal model such as MatTen that can predict the full elasticity tensors of inorganic crystalline compounds across crystal systems and chemical species brings new possibilities to probe and design materials with targeted mechanical properties. Notably, the full elasticity tensor enables the exploration of any anisotropic elastic property. Besides the directional Young's modulus studied in this work, many other elastic properties possess anisotropic behavior, such as linear compressibility, shear modulus, and Poisson's ratio. It should be noted, however, that a crucial aspect regarding the practical use of the model relies on the robustness of the input structure. Given two structures of a crystal where the atomic coordinates are slightly different, we would expect the elastic properties to be similar. If, otherwise, the model is extremely sensitive to the input, then it is not ideal for practical application. Extra work such as highly accurate DFT structure relaxation is needed before applying the model for predictions. We tested MatTen by using

Figure 5: Discovery of new crystals with large maximum directional Young’s modulus. The training data and new discoveries are separately normalized such that each sums to 100 percent.
structures directly queried from the Materials Project database and structures with tighter geometry optimization. The MAE of \(E_{\text{d}}^{\text{max}}\) between using the two types of structures is 6.55 GPa (Fig. S12 in the SI), more than three times smaller than the MAE (22.36 GPa) of \(E_{\text{d}}^{\text{max}}\) between model predictions and DFT references. This suggests that MatTen is robust enough to its input, and reasonably optimized structures (for example, from online databases) would not introduce extra error larger than the intrinsic error in the model.

The MatTen model is not limited to inorganic crystalline compounds and even elasticity tensors in general. The elastic behavior of other classes of materials such as two-dimensional layered materials and molecular crystals play a significant role in determining their functionality, and MatTen can be directly applied to model their elasticity tensors. Of course, a curated dataset of reference elasticity tensors is needed. Such data already exists, for example, in the Computational 2D Materials Database (C2DB) [59]. Moreover, besides elasticity tensor, MatTen can be applied to other tensorial properties of materials. These can be broadly categorized into two classes: material-level property and atom-level property. While the former means a single tensor for each crystal, the latter means a separate tensor for each atom in the crystal. Other material-level properties such as piezoelectric and dielectric tensors can be modeled by updating the output head as in Fig. 2 to use the corresponding irreducible representations of the tensor of interest (for example, a single second-rank symmetric matrix for the dielectric tensor). For atom-level properties such as the neutron magnet resonance (NMR) tensor, instead of using a mean pooling to aggregate atom features, one can directly map the features of an atom to a tensor for that atom. Using MatTen, we have conducted such an analysis for NMR tensors of silicon oxides and found that MatTen significantly outperforms both historic analytical models and other machine learning models by more than 50% for isotropic and anisotropic NMR chemical shift [60].

One potential limitation of the proposed approach is the reliance on a relatively large dataset to develop the model. We have curated a dataset of 10276 elasticity tensors which took millions of CPU hours to obtain. Such large datasets for other tensorial properties may not be readily available, but they begin to emerge. For example, the Materials Project has about 3000 piezoelectric [61] tensors and 7000 dielectric [62] tensors. This amount of data might still be a good start to training faithful models given that piezoelectric and dielectric are third- and second-rank tensors, respectively, which are much simpler than the fourth-rank elasticity tensor. In fact, for the second-rank NMR tensor, we only used a dataset of 421 crystals to obtain the best-performing model [60]. Another possibility is to apply a transfer learning approach, where the model is first trained on a different property with large data (for instance, elasticity tensor) and then finetuned on the target property of interest (for instance, piezoelectric tensor).

## 4 Methods

### Data generation

The elasticity tensors were computed by a liner fitting of the stresses and strains obtained from DFT calculations using the Vienna Ab Initio Simulation Package (VASP) [63]. The calculations follow the same procedures discussed in Ref. [10] with two general improvements. First, to get more precise stresses for calculating the elasticity tensor, the projection operators in VASP are evaluated in the reciprocal space, that is, the setting LREAL=False was adopted. Second, to reduce numerical error in the calculations, the stresses are symmetrized according to the crystal symmetry. The entire workflow was implemented in the open-source atomate package [15].
## Model architecture

A crystal structure is converted to a graph using a distance-based approach, where an edge is created between a pair of atoms if their distance is smaller than a cutoff radius \(r_{\text{cut}}\). Periodic boundary conditions are considered in the graph construction.

For atom \(i\), its atomic number \(Z_{i}\) is embedded as a vector with \(c\) components using a one-hot encoding to obtain the initial atom feature \(h_{i}\). The unit vector \(\hat{r}_{ij}\) from atom \(i\) to atom \(j\) is expanded using a spherical harmonics basis consisting up to a degree of \(l=4\). (Explicitly, this corresponds to the "\(0\text{e}+1\text{o}+2\text{e}+3\text{o}+4\text{e}\)" irreducible representations in e3nn notation [51]). The distance \(r_{ij}\) between atoms \(i\) and \(j\) is expanded into a vector \(R\) using the radial basis functions [46],

\[\text{RBF}_{n}(d)=\sqrt{\frac{2}{r_{\text{cut}}}}\frac{\sin(\frac{n\pi}{r_{ \text{cut}}}d)}{d} \tag{6}\]

where \(n=1,2,\ldots,\) is an index of the radial basis.

With the atom features \(F\), the spherical harmonics expansion of \(\hat{r}_{ij}\), and the radial basis expansion \(R\) of \(r_{ij}\) obtained from the embedding layers, the interaction block performs tensor product-based convolution to refine the atom features. This is achieved via Eq. (3), more specifically [27],

\[\mathcal{L}_{acm_{o}}^{l_{o}}(\{\vec{\mathbf{r}}_{ab}\},\{F_{bcm_{i}}^{l_{i}}\})= \sum_{m_{i},m_{f}}C_{(l_{f},m_{f})(l_{i},m_{i})}^{(l_{o},m_{o})}\sum_{b\in \mathcal{N}_{a}}R_{c}^{(l_{f},l_{i})}(r_{ab})Y_{l_{f}}^{m_{f}}(\hat{\mathbf{r}}_{ ab})F_{bcm_{i}}^{l_{i}}, \tag{7}\]

where \(a\) denotes the center atom, \(b\) denotes all its neighbors \(\mathcal{N}_{b}\) within the cutoff \(r_{\text{cut}}\); \(l\) is an integer indicating the degree of the spherical harmonic function, and \(m=-l,\ldots,l\); the subscripts \(i\), \(o\), and \(f\) indicate input, output, and filter, respectively; \(c\) is the channel index (for example, for the embedding layer, it indicates the components of the one-hot encoding); \(C\) denotes the Clebsch-Gordan coefficients; and finally, \(R_{c}^{(l_{f},l_{i})}\) are learnable multilayer perceptrons (MLPs), which take the RBF expansion as the input and contain most of the parameters of the model. Essentially, this combines the atom features of neighbors \(b\) to be the new atom features of the center atom \(a\), in the same spirit of a message-passing graph neural network. A major characteristic of Eq. (7) is that the use of the spherical harmonics and the Clebsch-Gordan coefficients together imply that convolutions are equivariant [27].

The atom features \(V\) is also passed through a self-interaction,

\[\mathcal{S}_{acm}^{l}(V_{ac^{\prime}m}^{l})=\sum_{c^{\prime}}W_{cc^{\prime}}^ {l}V_{ac^{\prime}m}^{l}, \tag{8}\]

where \(W_{cc^{\prime}}^{l}\) are learnable model parameters. The updated atom features is then obtained as

\[V^{k+1}=\mathcal{L}(V^{k})+\mathcal{S}(V^{k}), \tag{9}\]

where \(V^{k}\) denotes the features in interaction block \(k\), and \(V^{k+1}\) the features in the next interaction block. Indeed, \(V^{k+1}\) are further passed through a nonlinearity and a normalization functions. For each scalar part \(s\) in \(V\), the nonlinearity is chosen to be the SiLU function [64],

\[\text{SiLU}(s)=\sigma(s)s, \tag{10}\]

and for each non-scalar part \(\mathbf{t}\) in \(V\), the gated nonlinearity [65] is adopted,

\[G(\mathbf{t})=\sigma(x)\mathbf{t}, \tag{11}\]

where \(\sigma\) is the sigmoid function, and \(x\) is a scalar obtained from Eq. (7) by setting \(l_{o}=0\) and \(m_{o}=0\). Finally, the equivariant batch normalization [51] is applied to the features to avoid gradient vanishing or exploding.
The readout head aggregates the features of individual atoms to obtain a representation of the material via a mean pooling,

\[V_{\text{mat}}=\frac{1}{N}\sum_{a}^{N}V_{a}, \tag{12}\]

where \(V_{a}\) denotes the features of atom \(a\) for a crystal of a total number of \(N\) atoms. From \(V_{\text{mat}}\), the appropriate irreducible representations ("2x0e + 2x2e + 4e" in e3nn notation) that correspond to the elasticity tensor (two scalars, two second-rank traceless symmetric tensors, and one fourth-rank traceless symmetric tensor) are selected to construct the elasticity tensor according to Eq. (2).

### Model training

The dataset of 10276 elasticity tensors is split into three subsets for training, validation, and testing with a split ratio of 8:1:1. The model parameters are optimized using the training set, model hyperparameters are determined based on model performance on the validation set, and error analysis is performed using the test set unless otherwise stated. We train the model with the Adam optimizer to minimize a mean-squared-error loss function \(L=\sum_{i}^{B}\|\mathbf{C}_{i}-\mathbf{C}_{i}^{\text{ref}}\|^{2}\) with a mini-batch size \(B\) of 32. Note, \(\mathbf{C}_{i}\) denotes the irreducible representation of the model predicted elasticity tensor with 21 components (see Eq. (2)), but not the Cartesian tensor with 81 components. Similarly, \(\mathbf{C}_{i}^{\text{ref}}\) denotes the corresponding reference DFT values. The learning rate is set to 0.01, and a reduce-on-plateau learning rate scheduler is used, which decreases the learning rate by a factor of 0.5 if the validation error does not decrease for 50 epochs. The training stops when the validation error does not decrease for 150 consecutive epochs, and a maximum of 1000 epochs are allowed for the optimization. We performed a grid search to obtain model hyperparameters such as the \(r_{\text{cut}}\) and \(c\), and their optimal values are listed in Table S3 in the SI.

### Compliance tensor

The compliance tensor \(\mathbf{S}\) is a fourth-rank tensor defined from the inverse stress-strain relation \(\mathbf{\epsilon}=\mathbf{S}\mathbf{\sigma}\), where \(\mathbf{\epsilon}\) and \(\mathbf{\sigma}\) are the second-rank strain tensor and stress tensor, respectively. The compliance tensor in Voigt notation \(\mathbf{s}\) can be obtained as the inverse of the elasticity tensor Voigt matrix,

\[\mathbf{s}=\mathbf{c}^{-1}, \tag{13}\]

which is a 6 by 6 symmetric matrix. The relationships between the components of the full compliance tensor \(S_{ijkl}\) and the Vogit matrix \(s_{\alpha\beta}\) are [34]:

\[S_{ijkl} =s_{\alpha\beta},\text{ when }\alpha\text{ and }\beta\text{ are 1, 2, or 3,}\] \[2S_{ijkl} =s_{\alpha\beta},\text{ when either }\alpha\text{ or }\beta\text{ are 4, 5, or 6,} \tag{14}\] \[4S_{ijkl} =s_{\alpha\beta},\text{ when both }\alpha\text{ and }\beta\text{ are 4, 5, or 6.}\]

### Averaged elastic moduli of polycrystals

Given the elasticity tensor of a single crystal, the averaged bulk, shear, and Young's moduli of polycrystalline materials can be computed using different average schemes. The Voigt average assumes that the strain is the same in each grain, equal to the macroscopically-applied strain [35]. The Voigt bulk modulus is

\[K_{V}=[(c_{11}+c_{22}+c_{33})+2(c_{12}+c_{23}+c_{31})]/9, \tag{15}\]
and the Voigt shear modulus is

\[G_{V}=[(c_{11}+c_{22}+c_{33})-(c_{12}+c_{23}+c_{31})+3(c_{44}+c_{55}+c_{66})]/15. \tag{16}\]

The Reuss average assumes that the stress is the same in each grain, equal to the macroscopically applied stress [66]. The Reuss bulk modulus is

\[K_{R}=1/[(s_{11}+s_{22}+s_{33})+2(s_{12}+s_{23}+s_{31})], \tag{17}\]

and Reuss the shear modulus is

\[G_{R}=15/[4(s_{11}+s_{22}+s_{33})-4(s_{12}+s_{23}+s_{31})+3(s_{44}+s_{55}+s_{66 })]. \tag{18}\]

The Voigt and Reuss averages are the two extreme cases. The Hill average takes their arithmetic mean and is considered the most accurate in a wide range of experimental conditions [52]. The Hill bulk modulus is

\[K_{H}=(K_{V}+K_{R})/2, \tag{19}\]

and the Hill shear modulus is

\[G_{H}=(G_{V}+G_{R})/2. \tag{20}\]

Given the bulk modulus and the shear modulus (from any of the Voigt, Reuss, and Hill schemes), Young's modulus [67] can be computed as

\[E=9KG/(3K+G). \tag{21}\]

In this work, we report the bulk modulus, shear modulus, and Young's modulus from the Hill scheme.

### Scaled error

The mean absolute error (MAE) and mean absolute deviation (MAD) are defined as \(\text{MAE}=\frac{1}{N_{1}}\sum_{i}^{N_{1}}|y_{i}-y_{i}^{\text{pred}}|\) and \(\text{MAD}=\frac{1}{N_{2}}\sum_{i}^{N_{2}}|(y_{i}-\bar{y})|\), in which \(y_{i}\) is the reference value of data point \(i\), \(y_{i}^{\text{pred}}\) is the model prediction for the data point, and \(\bar{y}\) is the average of all reference values. The numbers \(N_{1}\) and \(N_{2}\) do not necessarily need to be the same. This is the case in Fig. 3 and Table 1, where \(N_{1}\) is the number of crystals in a specific crystal system and \(N_{2}\) is the total size of the test set. The scaled error (SE) is then computed as \(\text{SE}=\text{MAE}/\text{MAD}\).

### Software implementation

The MatTen model was implemented using the e3nn package [51] built on top of PyTorch [68], and the training was performed using Pytorch-Lightning [69]. The DFT calculations were performed using the atomate workflow [15] and all crystal structure processing was performed using the Python Materials Genomics (pymatgen) [13]. Directional Young's modulus was obtained using the elate package [54].

## Data Availability

The elasticity tensors used for model development, the 100 new crystals with large maximum directional Young's modulus, and the elemental cubic metals are available at [https://doi.org/10.5281/zenodo.8190849](https://doi.org/10.5281/zenodo.8190849) The elasticity tensors are also available from the Materials Project database via the web interface at [https://materialsproject.org](https://materialsproject.org) or the API at [https://api.materialsproject.org](https://api.materialsproject.org).
## Code Availability

The MatTen model and training scripts are released as an open-source repository at [https://github.com/wengroup/matten](https://github.com/wengroup/matten).

## Acknowledgements

The method development was supported by the startup funds from the Presidential Frontier Faculty Program at the University of Houston. Support for software and data infrastructure development was provided by the U.S. Department of Energy, Office of Science, Office of Basic Energy Sciences, Materials Sciences and Engineering Division under contract No. DE-AC02-05-CH11231 (Materials Project program KC23MP). This work used computational resources provided by the Research Computing Data Core at the University of Houston, the National Energy Research Scientific Computing Center (NERSC), a U.S. Department of Energy Office of Science User Facility operated under contract No. DE-AC02-05CH11231, and the Lawrencium computational cluster resource provided by the IT Division at the Lawrence Berkeley National Laboratory (Supported by the Director, Office of Science, Office of Basic Energy Sciences, of the U.S. Department of Energy under contract No. DE-AC02-05CH11231).

## Author contributions

Conceptualization, investigation, and writing - original draft: M.W.; data curation and formal analysis: M.W., M.K.H., J.M.M., and P.H.; writing - review & editing: M.W., M.K.H., J.M.M., P.H., and K.A.P.; project administration: M.W.; funding acquisition: M.W. and K.A.P.

## Competing interests

The authors declare no competing interests.

## References

* Hetnarski and Ignaczak [2016] Richard B Hetnarski and Jozef Ignaczak. _The mathematical theory of elasticity_. CRC Press, 2016. doi: 10.1115/1.1849176.
* Kaner et al. [2005] Richard B Kaner, John J Gilman, and Sarah H Tolbert. Designing superhard materials. _Science_, 308(5726):1268-1269, 2005. doi: 10.1126/science.1109830.
* Tehrani et al. [2018] Aria Mansouri Tehrani, Anton O Oliynyk, Marcus Parry, Zeshan Rizvi, Samantha Couper, Feng Lin, Lowell Miyagi, Taylor D Sparks, and Jakoah Brgoch. Machine learning directed search for ultraincompressible, superhard materials. _Journal of the American Chemical Society_, 140(31):9844-9853, 2018. doi: 10.1021/jacs.8b02717.
* Anderson et al. [1968] Orson L Anderson, Edward Schreiber, Robert C Liebermann, and Naohiro Soga. Some elastic constant data on minerals relevant to geophysics. _Reviews of Geophysics_, 6(4):491-524, 1968. doi: 10.1029/SP026p0237.
* Kark et al. [2001] Bijaya B Kark, Lars Stixrude, and Renata M Wentzcovitch. High-pressure elastic properties of major materials of earth's mantle from first principles. _Reviews of Geophysics_, 39(4):507-534, 2001. doi: 10.1029/2000RG000088.
* Monroe and Newman [2004] Charles Monroe and John Newman. The effect of interfacial deformation on electrodeposition kinetics. _Journal of The Electrochemical Society_, 151(6):A880, 2004. doi: 10.1149/1.1710893.
* Ahmad and Viswanathan [2017] Zeeshan Ahmad and Venkatasubramanian Viswanathan. Stability of electrodeposition at solid-solid interfaces and implications for metal anodes. _Physical review letters_, 119(5):056003, 2017. doi: 10.1103/PhysRevLett.119.056003.
* Mouhat and Coudert [2014] Felix Mouhat and Francois-Xavier Coudert. Necessary and sufficient elastic stability conditions in various crystal systems. _Physical review B_, 90(22):224104, 2014. doi: 10.1103/PhysRevB.90.224104.
* Tolborg et al. [2022] Kasper Tolborg, Johan Klarbring, Alex M Ganose, and Aron Walsh. Free energy predictions for crystal stability and synthesisability. _Digital Discovery_, 1(5):586-595, 2022. doi: 10.1039/D2DD00050D.
* De Jong et al. [2015] Maarten De Jong, Wei Chen, Thomas Angsten, Anubhav Jain, Randy Notestine, Anthony Gamst, Marcel Sluiter, Chaitanya Krishna Ande, Sybrand Van Der Zwaag, Jose J Plata, et al. Charting the complete elastic properties of inorganic crystalline compounds. _Scientific data_, 2(1):1-13, 2015. doi: 10.1038/sdata.2015.9.
* Du and Zhao [2017] Xinpeng Du and Ji-Cheng Zhao. Facile measurement of single-crystal elastic constants from polycrystalline samples. _npj Computational Materials_, 3(1):17, 2017. doi: 10.1038/s41524-017-0019-x.
* Lejaeghere et al. [2016] Kurt Lejaeghere, Gustav Bihlmayer, Torbjorn Bjorkman, Peter Blaha, Stefan Blugel, Volker Blum, Damien Caliste, Ivano E Castelli, Stewart J Clark, Andrea Dal Corso, et al. Reproducibility in density functional theory calculations of solids. _Science_, 351(6280):aad3000, 2016. doi: 10.1126/science.aad3000.
* Ong et al. [2013] Shyue Ping Ong, William Davidson Richards, Anubhav Jain, Geoffroy Hautier, Michael Kocher, Shreyas Cholia, Dan Gunter, Vincent L Chevrier, Kristin A Persson, and Gerbrand Ceder. Python materials genomics (pymatgen): A robust, open-source python library for materials analysis. _Computational Materials Science_, 68:314-319, 2013. doi: 10.1016/j.commatsci.2012.10.028.
* Jain et al. [2015] Anubhav Jain, Shyue Ping Ong, Wei Chen, Bharat Medasani, Xiaohui Qu, Michael Kocher, Miriam Brafman, Guido Petretto, Gian-Marco Rignanese, Geoffroy Hautier, et al. Fireworks: a dynamic workflow system designed for high-throughput applications. _Concurrency and Computation: Practice and Experience_, 27(17):5037-5059, 2015. doi: 10.1002/cpe.3505.
* Mathew et al. [2017] Kiran Mathew, Joseph H Montoya, Alireza Faghaninia, Shyam Dwarakanath, Muratahan Aykol, Hanmei Tang, Iek-heng Chu, Tess Smidt, Brandon Bocklund, Matthew Horton, et al. Atomate: A high-level interface to generate, execute, and analyze computational materials science workflows. _Computational Materials Science_, 139:140-152, 2017. doi: 10.1016/j.commatsci.2017.07.030.
* Jain et al. [2013] Anubhav Jain, Shyue Ping Ong, Geoffroy Hautier, Wei Chen, William Davidson Richards, Stephen Dacek, Shreyas Cholia, Dan Gunter, David Skinner, Gerbrand Ceder, et al. Commentary: The materials project: A materials genome approach to accelerating materials innovation. _APL materials_, 1(1):011002, 2013. doi: 10.1063/1.4812323.
* Hellenbrandt [2004] Mariette Hellenbrandt. The inorganic crystal structure database (icsd)--present and future. _Crystallography Reviews_, 10(1):17-22, 2004.
* Chen and Ong [2022] Chi Chen and Shyue Ping Ong. A universal graph deep learning interatomic potential for the periodic table. _Nature Computational Science_, 2(11):718-728, 2022. doi: 10.1038/s43588-022-00349-3.
* Wang et al. [2021] Anthony Yu-Tung Wang, Steven K Kauwe, Ryan J Murdock, and Taylor D Sparks. Compositionally restricted attention-based network for materials property predictions. _Npj Computational Materials_, 7(1):77, 2021. doi: 10.1038/s41524-021-00545-1.
* Dunn et al. [2020] Alexander Dunn, Qi Wang, Alex Ganose, Daniel Dopp, and Anubhav Jain. Benchmarking materials property prediction methods: the matbench test set and automatminer reference algorithm. _npj Computational Materials_, 6(138), 2020. doi: 10.1038/s41524-020-00433-0.
* De Breuck et al. [2021] Pierre-Paul De Breuck, Geoffroy Hautier, and Gian-Marco Rignanese. Materials property prediction for limited datasets enabled by feature selection and joint learning with modnet. _npj Computational Materials_, 7(1):83, 2021. doi: 10.1038/s41524-021-00552-2.
* Chen et al. [2019] Chi Chen, Weike Ye, Yunxing Zuo, Chen Zheng, and Shyue Ping Ong. Graph networks as a universal machine learning framework for molecules and crystals. _Chemistry of Materials_, 31(9):3564-3572, 2019. doi: 10.1021/acs.chemmater.9b01294.
* Choudhary and DeCost [2021] Kamal Choudhary and Brian DeCost. Atomistic line graph neural network for improved materials property predictions. _npj Computational Materials_, 7(1):185, 2021. doi: 10.1038/s41524-021-00650-1.
* Ahmad et al. [2018] Zeeshan Ahmad, Tian Xie, Chinmay Maheshwari, Jeffrey C Grossman, and Venkatasubramanian Viswanathan. Machine learning enabled computational screening of inorganic solid electrolytes for suppression of dendrite formation in lithium metal anodes. _ACS central science_, 4(8):996-1006, 2018. doi: 10.1021/acscentsci.8b00229.
* Revi et al. [2021] Vivek Revi, Saurabh Kasodariya, Anjana Talapatra, Ghanshyam Pilania, and Alankar Alankar. Machine learning elastic constants of multi-component alloys. _Computational Materials Science_, 198:110671, 2021. doi: 10.1016/j.commatsci.2021.110671.
* Bronstein et al. [2021] Michael M Bronstein, Joan Bruna, Taco Cohen, and Petar Velickovic. Geometric deep learning: Grids, groups, graphs, geodesics, and gauges. _arXiv preprint arXiv:2104.13478_, 2021.
* Thomas et al. [2018] Nathaniel Thomas, Tess Smidt, Steven Kearnes, Lusann Yang, Li Li, Kai Kohlhoff, and Patrick Riley. Tensor field networks: Rotation-and translation-equivariant neural networks for 3d point clouds. _arXiv preprint arXiv:1802.08219_, 2018.
* Satorras et al. [2021] Victor Garcia Satorras, Emiel Hoogeboom, and Max Welling. E (n) equivariant graph neural networks. In _International conference on machine learning_, pages 9323-9332. PMLR, 2021.
* Batzner et al. [2022] Simon Batzner, Albert Musaelian, Lixin Sun, Mario Geiger, Jonathan P Mailoa, Mordechai Kornbluth, 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, 2022. doi: 10.1038/s41467-022-29939-5.
* Takamoto et al. [2022] So Takamoto, Satoshi Izumi, and Ju Li. Teanet: Universal neural network interatomic potential inspired by iterative electronic relaxations. _Computational Materials Science_, 207:111280, 2022. doi: 10.1016/j.commatsci.2022.111280.
* Schutt et al. [2021] Kristof Schutt, Oliver Unke, and Michael Gastegger. Equivariant message passing for the prediction of tensorial properties and molecular spectra. In Marina Meila and Tong Zhang, editors, _Proceedings of the 38th International Conference on Machine Learning_, volume 139 of _Proceedings of Machine Learning Research_, pages 9377-9388. PMLR, 18-24 Jul 2021.
* Nguyen and Lunghi [2022] Vu Ha Anh Nguyen and Alessandro Lunghi. Predicting tensorial molecular properties with equivariant machine learning models. _Physical Review B_, 105(16):165131, 2022. doi: 10.1103/PhysRevB.105.165131.
* Fedorov [1968] Fedor I Fedorov. _Theory of elastic waves in crystals_. Springer, 1968. doi: 10.1007/978-1-4757-1275-9.
* Nye [1985] John Frederick Nye. _Physical properties of crystals: their representation by tensors and matrices_. Oxford university press, 1985. doi: 10.1063/1.3060200.
* Voigt [1966] Woldemar Voigt. _Leltrbuch der kristallphysik (mit Ausschluss der Kristalloptik)_. Vieweg+Teubner Verlag Wiesbade, 1966. doi: 10.1007/978-3-663-15884-4.
* Forte and Vianello [1996] Sandra Forte and Maurizio Vianello. Symmetry classes for elasticity tensors. _Journal of Elasticity_, 43:81-108, 1996. doi: 10.1007/BF00042505.
* Tadmor et al. [2012] Ellad B Tadmor, Ronald E Miller, and Ryan S Elliott. _Continuum mechanics and thermodynamics: from fundamental concepts to governing equations_. Cambridge University Press, 2012. doi: 10.1017/CBO9781139017657.
* Chadwick et al. [2001] Peter Chadwick, Maurizio Vianello, and Stephen C Cowin. A new proof that the number of linear elastic symmetries is eight. _Journal of the Mechanics and Physics of Solids_, 49(11):2471-2492, 2001. doi: 10.1016/S0022-5096(01)00064-3.
* Wallace [1972] Duane C. Wallace. _Thermodynamics of Crystals_. Wiley & Sons, 1972. doi: 10.1119/1.1987046.
* Singh et al. [2021] Sobhit Singh, Logan Lang, Viviana Dovale-Farelo, Uthpala Herath, Pedram Tavadze, Francois-Xavier Coudert, and Aldo H Romero. Mechanistic: A python library for analysis of mechanical and elastic properties of bulk and 2d materials. _Computer Physics Communications_, 267:108068, 2021. doi: 10.1016/j.cpc.2021.108068.
* Li et al. [2022] Yunguo Li, Lidunka Vocadlo, and John P Brodholt. Elast: A toolkit for thermoelastic calculations. _Computer Physics Communications_, 273:108280, 2022. doi: 10.1016/j.cpc.2021.108280.
* Ran et al. [2023] Zheng Ran, Chunming Zou, Zunjie Wei, and Hongwei Wang. Velas: An open-source toolbox for visualization and analysis of elastic anisotropy. _Computer Physics Communications_, 283:108540, 2023. doi: 10.1016/j.cpc.2022.108540.
* Itin and Hehl [2013] Yakov Itin and Friedrich W Hehl. The constitutive tensor of linear elasticity: its decompositions, cauchy relations, null lagrangians, and wave propagation. _Journal of Mathematical Physics_, 54(4):042903, 2013. doi: 10.1063/1.4801859.
* Itin [2020] Yakov Itin. Irreducible matrix resolution for symmetry classes of elasticity tensors. _Mathematics and Mechanics of Solids_, 25(10):1873-1895, 2020. doi: 10.1177/1081286520913596.
* Backus [1970] George Backus. A geometrical picture of anisotropic elastic tensors. _Reviews of geophysics_, 8(3):633-671, 1970. doi: 10.1029/RG008i003p00633.
* Gasteiger et al. [2020] Johannes Gasteiger, Janek Gross, and Stephan Gunnemann. Directional message passing for molecular graphs. In _International Conference on Learning Representations (ICLR)_, 2020.
* Xie and Grossman [2018] Tian Xie and Jeffrey C Grossman. Crystal graph convolutional neural networks for an accurate and interpretable prediction of material properties. _Physical review letters_, 120(14):145301, 2018. doi: 10.1103/PhysRevLett.120.145301.
* Wen et al. [2020] Mingjian Wen, Samuel M Blau, Evan Walter Clark Spotte-Smith, Shyam Dwaraknath, and Kristin A Persson. Bondnet: a graph neural network for the prediction of bond dissociation energies for charged molecules. _Chemical Science_, 12(5):1858-1868, 2020. doi: 10.1039/D0SC05251E.
* Wen et al. [2022] Mingjian Wen, Samuel M. Blau, Xiaowei Xie, Shyam Dwaraknath, and Kristin A. Persson. Improving machine learning performance on small chemical reaction data with unsupervised contrastive pre-training. _Chemical Science_, 13:1446-1458, January 2022. doi: 10.1039/D1SC06515G.
* He et al. [2016] Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun. Deep residual learning for image recognition. In _2016 IEEE Conference on Computer Vision and Pattern Recognition (CVPR)_, pages 770-778, 2016. doi: 10.1109/CVPR.2016.90.
* Geiger and Smidt [2022] Mario Geiger and Tess Smidt. e3nn: Euclidean neural networks. _arXiv preprint arXiv:2207.09453_, 2022.
* Hill [1952] Richard Hill. The elastic behaviour of a crystalline aggregate. _Proceedings of the Physical Society. Section A_, 65(5):349, 1952. doi: 10.1088/0370-1298/65/5/307.
* Lakes and Wojciechowski [2008] Rod Lakes and KW Wojciechowski. Negative compressibility, negative poisson's ratio, and stability. _physica status solidi (b)_, 245(3):545-551, 2008. doi: 10.1002/pssb.200777708.
* Gaillac et al. [2016] Romain Gaillac, Pluton Pullumbi, and Francois-Xavier Coudert. Elate: an open-source online application for analysis and visualization of elastic tensors. _Journal of Physics: Condensed Matter_, 28(27):275201, 2016. doi: 10.1088/0953-8984/28/27/275201.
* Robinson [1974] Author H Robinson. A new map projection: Its development and characteristics. _International yearbook of cartography_, 14(1974):145-155, 1974.
* Huang et al. [2012] Yunhan Huang, Arvind Sai Sarathi Vasan, Ravi Doraiswami, Michael Osterman, and Michael Pecht. Mems reliability review. _IEEE Transactions on Device and Materials Reliability_, 12(2):482-493, 2012. doi: 10.1109/TDMR.2012.2191291.
* Sun et al. [2016] Wenhao Sun, Stephen T Dacek, Shyue Ping Ong, Geoffroy Hautier, Anubhav Jain, William D Richards, Anthony C Gamst, Kristin A Persson, and Gerbrand Ceder. The thermodynamic scale of inorganic crystalline metastability. _Science advances_, 2(11):e1600225, 2016. doi: 10.1126/sciadv.1600225.
* Bartel [2022] Christopher J Bartel. Review of computational approaches to predict the thermodynamic stability of inorganic solids. _Journal of Materials Science_, 57(23):10475-10498, 2022. doi: 10.1007/s10853-022-06915-4.
* Gjerding et al. [2021] Morten Niklas Gjerding, Alireza Taghizadeh, Asbjorn Rasmussen, Sajid Ali, Fabian Bertoldo, Thorsten Deilmann, Nikolaj Rorbaek Knesgaard, Mads Kruse, Ask Hjorth Larsen, Simone Manti, et al. Recent progress of the computational 2d materials database (c2db). _2D Materials_, 8(4):044002, 2021. doi: 10.1088/2053-1583/ac1059.
* Venetos et al. [2023] Maxwell C Venetos, Mingjian Wen, and Kristin A Persson. Machine learning full nmr chemical shift tensors of silicon oxides with equivariant graph neural networks. _The Journal of Physical Chemistry A_, 127(10):2388-2398, 2023. doi: 10.1021/acs.jpca.2c07530.
* De Jong et al. [2015] Maarten De Jong, Wei Chen, Henry Geerlings, Mark Asta, and Kristin Aslaug Persson. A database to enable discovery and design of piezoelectric materials. _Scientific data_, 2(1):1-13, 2015. doi: 10.1038/sdata.2015.53.
* Petousis et al. [2017] Ioannis Petousis, David Mrdjenovich, Eric Ballouz, Miao Liu, Donald Winston, Wei Chen, Tanja Graf, Thomas D Schladt, Kristin A Persson, and Fritz B Prinz. High-throughput screening of inorganic compounds for the discovery of novel dielectric and optical materials. _Scientific data_, 4(1):1-12, 2017. doi: 10.1038/sdata.2016.134.
* Kresse and Hafner [1993] Georg Kresse and Jurgen Hafner. Ab initio molecular dynamics for liquid metals. _Physical review B_, 47(1):558, 1993. doi: 10.1016/0022-3093(95)00355-X.
* Hendrycks and Gimpel [2016] Dan Hendrycks and Kevin Gimpel. Gaussian error linear units (gelus). _arXiv preprint arXiv:1606.08415_, 2016.
* Weiler et al. [2018] Maurice Weiler, Mario Geiger, Max Welling, Wouter Boomsma, and Taco S Cohen. 3d steerable cnns: Learning rotationally equivariant features in volumetric data. _Advances in Neural Information Processing Systems_, 31, 2018.
* Reuss [1929] Andras Reuss. Berechnung der fliesgrenze von mischkristallen auf grund der plastizitatsbedingung fur einkristalle. _ZAMM-Journal of Applied Mathematics and Mechanics/Zeitschrift fur Angewandte Mathematik und Mechanik_, 9(1):49-58, 1929.
* Anand et al. [2022] Lallit Anand, Ken Kamrin, and Sanjay Govindjee. _Introduction to mechanics of solid materials_. Oxford University Press, 2022. doi: 10.1093/oso/9780192866073.002.0004.
* Paszke et al. [2019] Adam Paszke, Sam Gross, Francisco Massa, Adam Lerer, James Bradbury, Gregory Chanan, Trevor Killeen, Zeming Lin, Natalia Gimelshein, Luca Antiga, et al. Pytorch: An imperative style, high-performance deep learning library. _Advances in neural information processing systems_, 32, 2019.
* [69] William Falcon. Pytorch lightning. [https://github.com/Lightning-AI/lightning](https://github.com/Lightning-AI/lightning), 2023. accessed 2023-06-11.
Supplementary Information:

A Universal Equivariant Graph Neural Network for the Elasticity Tensors of Any Crystal System

Mingjian Wen\({}^{1}\)1, Matthew K. Horton\({}^{2,3}\), Jason M. Munro\({}^{2}\), Patrick Huck\({}^{4}\), Kristin A. Persson\({}^{5,6}\)

Footnote 1: Email: mjwen@uh.edu

\({}^{1}\) Chemical and Biomolecular Engineering, University of Houston, Houston, 77204, TX, USA

\({}^{2}\) Materials Sciences Division, Lawrence Berkeley National Laboratory, Berkeley, 94720, CA, USA

\({}^{3}\) Microsoft Research, Redmond, 98052, WA, USA

\({}^{4}\) Energy Technologies Area, Lawrence Berkeley National Laboratory, Berkeley, 94720, CA, USA

\({}^{5}\) Molecular Foundry, Lawrence Berkeley National Laboratory, Berkeley, 94720, CA, USA

\({}^{6}\) Department of Materials Science and Engineering, University of California, Berkeley, Berkeley, 94720, CA, USA

## The ten symmetry classes of elasticity tensors

Early approaches take inspiration from crystallography. Out of the 32 distinct crystallographic point groups, only 11 are centrosymmetric (meaning the point group contains an inversion center as one of its symmetry elements), each forming a unique diffraction pattern. The diffraction patterns of other noncentrosymmetric crystals each is the same as one of the 11 centrosymmetric crystals. Based on the diffraction patterns, the 32 distinct point groups can be categorized into 11 classes, called the Laue groups [3]. According to the Laue groups, Wallace [1] classifies the elastic tensors into 12 classes (the additional 1 being the isotropic class that does not apply for single crystals), and they reduce to 10 classes considering the number of independent components (Fig. S1). The results are widely cited, including the classical book on the subject by Nye [2] and many recent papers [4, 5, 6]. This crystallographic approach seems reasonable; however, the conclusions are incorrect. The tetragonal and trigonal systems are each divided into two symmetry classes, but the distinctions can be eliminated by a different choice of the coordinate system [7, 8]. Then each of the tetragonal and trigonal systems will have 6 independent components.

## Harmonic decomposition of the elasticity tensor

In the harmonic decomposition, the elasticity tensor can be written as

\[\mathbf{C}=h_{1}(\lambda)+h_{2}(\eta)+h_{3}(\mathbf{A})+h_{4}(\mathbf{B})+h_{5}(\mathbf{H}). \tag{1}\]

The appropriate values for each of the term is as follows:

\[\lambda=[2C_{ppmm}-C_{pmpm}]/15,\] \[\eta=[-3C_{ppmm}-C_{pmpm}]/90,\]
\[A_{ij} =[15C_{ijmm}-12C_{imjm}-5\delta_{ij}C_{ppmm}+4\delta_{ij}C_{pmpm}]/21,\] \[B_{ij} =[-6C_{ijmn}+9C_{imjm}+2\delta_{ij}C_{ppmm}-3\delta_{ij}C_{pmpm}]/21,\] \[H_{ijkl}= (C_{ijkl}+C_{iklj}+C_{iljk})/3\] \[-[(C_{ijmm}+2C_{imjm})\delta_{kl}+(C_{klm}+2C_{kmlm})\delta_{ij}\] \[+(C_{ikmm}+2C_{imkm})\delta_{jl}+(C_{jlmn}+2C_{jmlm})\delta_{ik}\] \[+(C_{ilmm}+2C_{imlm})\delta_{jk}+(C_{jkmn}+2C_{jmkm})\delta_{il}]/21\] \[+(C_{ppmm}+2C_{pmpm})(\delta_{ij}\delta_{kl}+\delta_{ik}\delta_{ jl}+\delta_{il}\delta_{jk})/105,\]

and

\[h_{1}(\lambda)=\delta_{ij}\delta_{kl}\lambda,\]

\[h_{2}(\eta)=(\delta_{ik}\delta_{jl}+\delta_{il}\delta_{jk})\eta,\]

\[h_{3}(\mathbf{A})=\delta_{ij}A_{kl}+\delta_{kl}A_{ij},\]

\[h_{4}(\mathbf{B})=\delta_{ik}B_{jl}+\delta_{jl}B_{ik}+\delta_{il}B_{jk}+\delta_{ jk}B_{il},\]

\[h_{5}(\mathbf{H})=H_{ijkl},\]

where \(\delta_{ij}\) is the Kronecker delta.
This decomposition follows Ref. [9], and as mentioned there that "... other forms of harmonic decomposition are possible: It suffices to use invertible linear combinations of \(\mathbf{A}\) and \(\mathbf{B}\) and, analogously, invertible linear combinations of \(\lambda\) and \(\eta\)." See Ref. [10, 11] for such examples. Nevertheless, harmonic decomposition is unique to linear combinations.

## Proof of MatTen satisfying material symmetry

The MatTen model \(\mathbf{C}=f(x)\) is equivariant to \(SO(3)\) transformations, satisfying

\[D_{y}(g)f(x)=f(D_{x}(g)x). \tag{2}\]

This comes from the fact that each layer of MatTen is equivariant, and the composition of such layers is also equivariant. We refer to Ref. [12] for proof of the equivariance of the layers. The representation \(D_{x}(g)\) in the space of crystal structures can be written as \(D_{x}(g)=R_{ip}\), and the representation \(D_{y}(g)\) in the space of stiffness tensors can be written as \(D_{g}(g)=R_{ip}R_{jq}R_{kr}R_{ls}\), where \(R\in SO(3)\) is a rotation matrix.

Let \(Q\in P\), where \(P\) denotes the set of rotations in the point group of a crystal, we will have \(P\subset SO(3)\). Therefore, for \(R=Q\), Eq. S(2) is satisfied, i.e.,

\[Q_{ip}Q_{jq}Q_{kr}Q_{ls}C_{prqs}=f(Qx) \tag{3}\]

Owning to material symmetry, we have \(Qx=x\), that is, the crystal structure \(x\) is indistinguishable before and after the transformation. Thus, \(f(Qx)=f(x)\). Plugging it into Eq. S(2), we have

\[Q_{ip}Q_{jq}Q_{kr}Q_{ls}C_{prqs}=f(x)=C_{ijkl}, \tag{4}\]

which is Eq.(1) in the main text. Once this is satisfied, the material symmetry will be reflected in the stiffness tensor as discussed in the main text and proved in Ref. [9].
Figure S3: Distribution of the bulk, shear, and Young's moduli in the dataset.
Figure S5: Distribution of shear modulus in the dataset by crystal system.
## Additional results on isotropic properties
## Directional Young's modulus

Here we prove that, for cubic crystals,

\[\text{if}\ S_{1111}-S_{1122}-2S_{2323}<0,E_{\text{d}}^{\text{max}}\text{is along}\ \langle 100\rangle\ \text{and}\ E_{\text{d}}^{\text{max}}\text{is along}\ \langle 111\rangle, \tag{5}\]

otherwise,

\[\text{if}\ S_{1111}-S_{1122}-2S_{2323}>0,E_{\text{d}}^{\text{max}}\text{is along}\ \langle 111\rangle\ \text{and}\ E_{\text{d}}^{\text{max}}\text{is along}\ \langle 100\rangle. \tag{6}\]

and

\[\text{if}\ S_{1111}-S_{1122}-2S_{2323}=0,\text{ the materials is isotropic regarding Young's modulus.} \tag{7}\]

The inverse of the directional Young's modulus is

\[E_{\text{d}}(\mathbf{n})^{-1}=n_{i}n_{j}n_{k}n_{l}S_{ijlk}, \tag{8}\]

where \(S_{ijkl}\) is the compliance tensor and \(\mathbf{n}\) is an unit direction vector. For a cubic crystal, the 21 non-zero components can be classified into three groups [2]:

* \(S_{1111}=S_{2222}=S_{3333}\)
* \(S_{1122}=S_{2211}=S_{2233}=S_{3322}=S_{3311}=S_{1133}\)
* \(S_{2323}=S_{2332}=S_{3223}=S_{3232}=S_{1212}=S_{1221}=S_{2112}=S_{1313}=S_{13 31}=S_{3113}\).

Substituting these into Eq. S(8), we have,

\[E_{\text{d}}(\mathbf{n})^{-1} =S_{1111}(n_{1}^{4}+n_{2}^{4}+n_{3}^{4})+2S_{1122}(n_{1}^{2}n_{2}^ {2}+n_{2}^{2}n_{3}^{2}+n_{3}^{2}n_{1}^{2})+4S_{2323}(n_{1}^{2}n_{2}^{2}+n_{2}^ {2}n_{3}^{2}+n_{3}^{2}n_{1}^{2})\] \[=S_{1111}[1-2(n_{1}^{2}n_{2}^{2}+n_{2}^{2}n_{3}^{2}+n_{3}^{2}n_{1} ^{2})]+2S_{1122}(n_{1}^{2}n_{2}^{2}+n_{2}^{2}n_{3}^{2}+n_{3}^{2}n_{1}^{2})+4S_ {2323}(n_{1}^{2}n_{2}^{2}+n_{2}^{2}n_{3}^{2}+n_{3}^{2}n_{1}^{2})\] \[=S_{1111}-2(S_{1111}-S_{1122}-2S_{2323})(n_{1}^{2}n_{2}^{2}+n_{2} ^{2}n_{3}^{2}+n_{3}^{2}n_{1}^{2})\] \[=S_{1111}-2(S_{1111}-S_{1122}-2S_{2323})f. \tag{9}\]

In the second equality, we have used \(n_{1}^{4}+n_{2}^{4}+n_{3}^{4}=(n_{1}^{2}+n_{2}^{2}+n_{3}^{2})^{2}-2(n_{1}^{2}n_ {2}^{2}+n_{2}^{2}n_{3}^{2}+n_{3}^{2}n_{1}^{2})=1-2(n_{1}^{2}n_{2}^{2}+n_{2}^{2 }n_{3}^{2}+n_{3}^{2}n_{1}^{2})\), in which \((n_{1}^{2}+n_{2}^{2}+n_{3}^{2})^{2}=1\), because \(\mathbf{n}\) is a unit vector. In the last equality, we have defined \(f:=n_{1}^{2}n_{2}^{2}+n_{2}^{2}n_{3}^{2}+n_{3}^{2}n_{1}^{2}\).

From Eq. S(9), it seen that Eq. S(7) is valid.

In fact, \(f\) has its maximum value of \(1/3\) along the \(\langle 111\rangle\) directions, and the minimum of \(f\) is 0 along the \(\langle 100\rangle\) directions (derived below). As a result, Eq. S(5) and Eq. S(6) are valid. (Note that Eq. S(9) gives the inverse of the directional Young's modulus.)

\begin{table}
\begin{tabular}{c c c c c c c} \hline \hline  & \multicolumn{2}{c}{\(\log_{10}(K)\)} & \multicolumn{2}{c}{\(\log_{10}(G)\)} & \multicolumn{2}{c}{\(\log_{10}(E)\)} \\ \cline{2-7}  & MAE & MAE/MAD & MAE & MAE/MAD & MAE & MAE/MAD \\ \hline Matflen & 0.046 (0.002) & 0.166 (0.006) & 0.094 (0.002) & 0.331 (0.010) & 0.087 (0.002) & 0.309 (0.018) \\ AutoMatminer & 0.050 (0.002) & 0.187 (0.009) & 0.090 (0.002) & 0.307 (0.006) & 0.086 (0.002) & 0.301 (0.009) \\ \hline \hline \end{tabular}
\end{table}
Table S1: Prediction of the bulk modulus \(K\), shear modulus \(G\), and Young’s modulus \(E\) in logarithmic space. \(K\), \(G\), and \(E\) are in the units of GPa. The results for MatTen are calculated from a single model, while a separate AutoMatminer model is trained for each property. The value in a pair of parentheses is the standard deviation from an ensemble of five models trained with different initialization. MAE: mean absolute error; MAD: mean absolute deviation.
Below, we show that the maximum of \(f\) is \(1/3\) along the \(\langle 111\rangle\) directions, and the minimum of \(f\) is 0 along the \(\langle 100\rangle\) directions.

Let \(n_{1}^{2}=a,n_{2}^{2}=b\) and \(n_{3}^{2}=c\), we have \(a+b+c=1\) because \(\mathbf{n}\) is a unit vector. Thus,

\[f=ab+bc+ca=ab+c(b+a)=ab+(1-a-b)(b+a)=a+b-ab-a^{2}-b^{2}. \tag{10}\]

Let

\[\frac{\partial f}{\partial a} =1-b-2a=0 \tag{11}\] \[\frac{\partial f}{\partial b} =1-a-2b=0,\]

and solve the equations, we have \(a=b=c=1/3\), i.e. \(n_{1}^{2}=n_{2}^{2}=n_{3}^{3}=1/3\), At these values, \(f=1/3\) and we can verify that it is a maximum. This also suggests \(\mathbf{n}\) is along the \(\langle 111\rangle\) family of directions.

The other extreme values of \(f\) are located at the boundaries of \(a\) (or \(b\) or \(c\)). Since \(n_{1}\) is a component of the unit vector, then \(n_{1}\in[-1,1]\), i.e. \(a\in[0,1]\). So, the extreme value is obtained when

* \(n_{1}=0\), \(n_{2}=\pm 1\), \(n_{3}=0\)
* \(n_{1}=0\), \(n_{2}=0\), \(n_{3}=\pm 1\)
* \(n_{1}=\pm 1\), \(n_{2}=0\), \(n_{3}=0\).

These are the \(\langle 100\rangle\) directions, at which the minimum is \(f=0\).
Figure S10: Distribution of the prediction error for shear modulus. The directional shear modulus can be computed as \(G_{\mathrm{d}}(\mathbf{n},\mathbf{m})=(n_{i}m_{j}n_{k}m_{l}S_{ijkl})^{-1}\), where \(\mathbf{n}\) and \(\mathbf{m}\) are two direction unit vectors, and \(\mathbf{S}\) is the compliance tensor [6]. The data of \(G_{\mathrm{d}}\) is obtained by sampling in a way similar to \(E_{\mathrm{d}}\) discussed in the main text. The prediction error is the difference between MatIen prediction and DFT reference.
Fig. S12 shows the \(E_{\text{d}}^{\text{max}}\) for the 100 new crystals. The MAE between predicted \(E_{\text{d}}^{\text{max}}\) with initial structure and predicted \(E_{\text{d}}^{\text{max}}\) with further relaxed structure is 6.55 GPa. It is much smaller than the MAE (22.36 GPa) between MatTen prediction and DFT reference for the test set. This demonstrates the robustness of MatTen with respect to the structure of the input crystal as discussed in the main text. As shown in Fig. S12 b, if we consider the 100 new crystals instead of the test set, the MAE between MatTen prediction and DFT is much higher, with a value of 48.69 GPa. This is expected, since, for the 100 new crystals, we are probing extreme values at the edge of the training data distribution, while the test set follows the same distribution of the training data. This signifies the importance of further confirmation with more accurate computation (DFT in this case) and even experiments once the search space has been narrowed down via the screening using the model.

\begin{table}
\begin{tabular}{c c c c c c} \hline Materials Project ID & Formula & \(\Delta S_{\text{DFT}}\) & \(\Delta S_{\text{MatTen}}\) & Experimentally observed & Ground-state polymorph \\ \hline mp-129 & Mo & -0.00150 & -0.00191 & Yes & Yes \\ mp-146 & V & -0.00994 & -0.01006 & Yes & Yes \\ mp-17 & Cr & -0.00267 & -0.00253 & Yes & No \\ mp-90 & Cr & -0.00369 & 0.00055 & Yes & Yes \\ mp-91 & W & -0.00056 & -0.00042 & Yes & Yes \\ mp-11334 & W & -0.00285 & -0.00331 & No & No \\ mp-35 & Mn & -0.02022 & -0.00277 & Yes & Yes \\ mp-1186040 & Na & -0.15435 & -0.32852 & No & No \\ mp-1184808 & K & -0.34397 & -0.12334 & No & No \\ mp-949029 & Cs & -0.53668 & -4.79427 & No & No \\ mp-1239193 & Rh & -0.03371 & -0.06590 & No & No \\ mp-1187790 & TI & -0.09708 & -0.05340 & No & No \\ \hline \end{tabular}
\end{table}
Table S2: Polymorphs of elemental cubic metal with \(E_{\text{d}}^{\text{max}}\) along \(\langle 100\rangle\) directions and \(E_{\text{d}}^{\text{min}}\) along \(\langle 111\rangle\) directions. \(\Delta S=S_{1111}-S_{1122}-2S_{2323}\). Among the crystal structures with the same composition, the one having the lowest energy is called the ground-state polymorph and is stable with respect to phase transition into other structures [13]. The crystal structures and the elasticity tensors of these metals are provided as well. See Data Availability in the main text.
## Model hyperparameters