#  StrainTensorNet: Predicting crystal structure elastic properties using
SE(3)-equivariant graph neural networks

Teerachote Pakornchote  Annop Ektarawong  Thiparat Chotibut
<EMAIL_ADDRESS>Chula Intelligent and Complex Systems, Department of
Physics, Faculty of Science, Chulalongkorn University, Bangkok, Thailand,
10330

###### Abstract

Accurately predicting the elastic properties of crystalline solids is vital
for computational materials science. However, traditional atomistic scale ab
initio approaches are computationally intensive, especially for studying
complex materials with a large number of atoms in a unit cell. We introduce a
novel data-driven approach to efficiently predict the elastic properties of
crystal structures using SE(3)-equivariant graph neural networks (GNNs). This
approach yields important scalar elastic moduli with the accuracy comparable
to recent data-driven studies. Importantly, our symmetry-aware GNNs model also
enables the prediction of the strain energy density (SED) and the associated
elastic constants, the fundamental tensorial quantities that are significantly
influenced by a material’s crystallographic group. The model consistently
distinguishes independent elements of SED tensors, in accordance with the
symmetry of the crystal structures. Finally, our deep learning model possesses
meaningful latent features, offering an interpretable prediction of the
elastic properties.

Keywords

strain energy tensor, elastic constants, crystallographic group,
equivariant neural networks, density functional theory

Open data

The dataset used in this work has been made available at
https://github.com/trachote/predict_elastic_tensor

##  I Introduction

The elastic properties of crystalline solids such as elastic constants, bulk
modulus, and shear modulus are important macroscopic quantities that determine
materials’ mechanical characteristics. Computational study of which can
provide theoretical guidelines to understand various phenomena in solid
materials, e.g., mechanical stability of crystal structures [1, 2], pressure-
driven deformation and phase transition of materials [3, 4], the response of
materials to sound wave propagation [5, 6], and the hardness of materials [7,
8, 9], to name a few. From an atomistic scale description, an ab initio
approach based on density functional theory has been employed to investigate
the macroscopic elastic properties of materials, yielding, for example,
elastic constants and elastic moduli that are in agreement with experimental
results [10, 11, 12, 13].

Three atomistic scale computational methods are usually adopted to calculate
the elastic constants of crystal structures: an energy-based method [14], a
stress-strain method [10, 15], and a calculation from group velocity [16, 17,
18]. The energy-based and the stress-strain methods are more standard, and
derive the elastic constants using respectively an energy and a stress tensor
acquired from an ab initio calculation. Although both methods yield a
comparable prediction of elastic constants, the energy-based method is
computationally less efficient, involving a larger number of plane-wave basis
set and $k$-point meshes so that the predicted elastic constants converge to
reliable values [19]. Indeed, the more efficient stress-strain method is
commonly employed to determine elastic constants by established softwares such
as VASP, Quantum Espresso and CASTEP [20, 21, 22]. The stress-strain method is
also utilized to create a database for elastic properties of inorganic
compounds [23, 24].

However, atomistic scale simulation can be computationally prohibitive,
especially when the number of atoms in a unit cell grows large. Such
constraint limits an ab initio method’s capability to investigate more complex
crystalline solids. On the other hand, advances in machine learning bring
about alternative approaches to computational materials science. This data-
driven paradigm can reasonably predict the elastic properties of crystal
structures, provided sufficient training data from an ab initio method [25,
26]. Even when the number of atoms in a unit cell is large, such approach can
efficiently predict the bulk and shear moduli of complex materials such as
alloys [27, 28]. Applying machine learning together with ab initio
calculations is also potentially useful in searching for novel superhard
materials [29, 30].

Machine learning models based on graph neural networks (GNNs) have received
increasing attention in studying solid materials and molecules. With GNNs, it
is natural to encode atomistic scaled descriptions of solids into a
computational framework; atomic locations and the associated atomic attributes
can be directly embedded into node attributes of a graph, whereas pairwise
interactions among atoms can be encoded into the edge attributes. Efficient
GNNs training procedures have also been proposed [31, 32, 33], enabling GNNs
to learn the representation of a complex relationship between the input and
its associated prediction [34]. These neural networks can also be endowed with
translation-rotation equivariant property, so that input atomic locations of
molecules (or point clouds) in $\mathbb{R}^{3}$ that differ only in their
orientation or their centroid can be identified. Enforcing an SE(3)
equivariant property helps the networks to extract a more compact (translation
and rotation independent) relationship between the inputs and their
predictions. Due to these appealing features, variations of SE(3) equivariant
GNNs have been developed to study materials science 111Tensor Field Network
(TFN), L1Net, GemNet, and EGNN, have exploited this computational architecture
to make better data-driven predictions on the properties of chemical compounds
[41, 52, 50, 51]. The SE(3)-Transformers model based on TFN with an attention
mechanism also demonstrates a good predictive performance on the standard QM9
dataset [42]. In addition, translation-rotation equivariance concept has also
been incorporated into generative modeling, in an attempt to search for novel
molecules and crystal structures with desirable properties [54, 55]..

In this work, we adopt a data-driven approach using graph neural networks
(GNNs) to predict the elastic properties of crystal structures. Accounting for
the symmetry of crystalline solids, our SE(3) equivariant GNNs take as an
input atomistic descriptions of strained materials, such as atomic locations
and their associated atomic attributes, and predict the strain energy density
(SED) of the materials. The prediction of the SED, which is the energy stored
in a crystal structure when distorted by either tensile or compressive strains
in particular directions, can be obtained efficiently and relatively
accurately, given sufficient training data. The model thus provides an
alternative approach to the standard ab initio energy-based prediction method.

After the SED is computed, we can then calculate the elastic constants and
construct the elastic tensor with a simple analytical expression, see Sec.
II.1. Other macroscopic (scalar) elastic properties including bulk modulus,
shear modulus, Young’s modulus, and Poisson’s ratio immediately follow from
the elastic constants. Sec. IV reports our model prediction results of these
elastic properties. The prediction performances of the scalar elastic
properties are comparable to those of the recent data-driven work [25, 26].
Importantly, however, the data-driven prediction of the SED and the associated
elastic constants, which are fundamental tensorial quantities that depend on
the crystallographic group, are first reported here. The trained model
consistently reveal the independent components of strain energy tensors
dictated by the symmetry of the crystal structures, as evident by the
degeneracy structure of the distance metrics computed in Sec. IV.2. Our
symmetry-aware GNNs architecture as well as the dataset used to train the
model are provided in Sec. III. Also, as opposed to a black-box deep learning
prediction, we show in Sec. IV.3 that the latent features (data
representations) that the model learned are physically meaningful, rendering
our GNNs explainable. We conclude this work in Sec. V, and provide
supplementary results and necessary derivations in the Appendices.

##  II Linear Elasticity Background

###  II.1 The strain energy tensor and the elastic tensor

Assuming linear elastic deformation and no external stress, the strain energy
density (SED) of isotropic materials can be expressed as [36]

$U(\epsilon_{1},\epsilon_{2},...,\epsilon_{6})=\frac{1}{2}\sum_{i,j=1}^{6}C_{ij}\epsilon_{i}\epsilon_{j},$
(1)

where $C_{ij}$ is an elastic constant and $\epsilon_{i}$ is a strain component
$i\in\\{1,2,\dots,6\\}$ in Voigt notation. With the goal of obtaining the
elastic constant, it suffices to focus on the SED tensor of rank 2, which we
now show. Denote $U_{ij}\equiv U(\epsilon_{i},\epsilon_{j})$ for $i\neq j$ and
$U_{ii}\equiv U(\epsilon_{i})$, which are the energy stored in a distorted
crystal structure when the lattice is strained by at most two strain
components (two Voigt indices). In this study, we consider $U_{ij}$ in units
of eV per atom. Since the elastic constant tensor is symmetric, the SED tensor
is also symmetric and can be expressed as

$\displaystyle U_{ij}$
$\displaystyle=C_{ij}\epsilon_{i}\epsilon_{j}+\frac{1}{2}\left(C_{ii}\epsilon_{i}^{2}+C_{jj}\epsilon_{j}^{2}\right)\text{\
for\ }i\neq j\text{ and,}$ (2) $\displaystyle U_{ii}$
$\displaystyle=\frac{1}{2}C_{ii}\epsilon_{i}^{2}.$

Thus, the elastic constants can then be analytically obtained as a function of
the SED:

$C_{ij}=\frac{1}{\epsilon_{i}\epsilon_{j}}\left[(1+\delta_{ij})U_{ij}-\left(1-\delta_{ij}\right)\left(U_{ii}+U_{jj}\right)\right],$
(3)

where $\delta_{ij}$ is the Kronecker delta.

Because the SED tensor of interest and the elastic tensor are symmetric
tensors of rank 2 with $6\times 6$ components, so there are $(6+1)6/2=21$
independent components. In the matrix form, we denote the SED tensor as

$\mathbf{U}=\begin{bmatrix}U_{11}&U_{12}&U_{13}&U_{14}&U_{15}&U_{16}\\\\[3.0pt]
&U_{22}&U_{23}&U_{24}&U_{25}&U_{26}\\\\[3.0pt]
&&U_{33}&U_{34}&U_{35}&U_{36}\\\\[3.0pt] &&&U_{44}&U_{45}&U_{46}\\\\[3.0pt]
&&&&U_{55}&U_{56}\\\\[3.0pt] &&&&&U_{66}\end{bmatrix},$ (4)

where only the upper triangular elements are shown (the lower diagonal
elements are left blank for brevity). For crystalline solids, however, the
upper triangular elements are not completely independent, depending on the
symmetry of the crystal structure. For instance, the SED tensor of a cubic
lattice is

$\mathbf{U}_{cubic}=\begin{bmatrix}U_{11}&U_{12}&U_{12}&U_{14}&U_{14}&U_{14}\\\\[3.0pt]
&U_{11}&U_{12}&U_{14}&U_{14}&U_{14}\\\\[3.0pt]
&&U_{11}&U_{14}&U_{14}&U_{14}\\\\[3.0pt] &&&U_{44}&U_{45}&U_{45}\\\\[3.0pt]
&&&&U_{44}&U_{45}\\\\[3.0pt] &&&&&U_{44}\end{bmatrix}.$ (5)

Note that some components of the elastic constants $C_{ij}$ can be zero in
many crystal structures, e.g. a cubic lattice has 9 zeros out of 21
independent components. However, due to the property of SED from Eq. (2),
$U_{ij}$ is never 0. For the purpose of machine learning regression, working
with SED helps avoid zero inflation problems in the training dataset. Fig. E1
shows the distributions of $C_{ij}$ and $U_{ij}$, which cover different ranges
depending on the indices. Importantly, the distribution of elastic constants
can concentrate around zero, but this is not the case for the SED. As an
illustrative example of how the SED can avoid a zero inflation problem,
consider the elastic constants of a diamond cubic crystal structure. For
diamond, $C_{11}$, $C_{12}$, $C_{14}$, $C_{44}$, and $C_{45}$ are 1054, 126,
0, 562, and 0 GPa, respectively. Then, Eq. (2) gives $U_{11}$, $U_{12}$,
$U_{14}$, $U_{44}$, and $U_{45}$ subjected to 2% strain to be 7.506, 16.807,
11.509, 4.002, and 8.005 meV/atom, respectively. Note that, in this study, we
consider $U_{ij}$ in the unit of eV per atom instead of per volume. The
magnitude of $U_{12}$ is the largest as it is the sum of $C_{12}$ and
$C_{11}$, while $U_{14}$ and $U_{45}$ are smaller (but non-zero) because
$C_{14}$ and $C_{45}$ are zero.

###  II.2 Strain operations

The unit cell of a crystal can be compactly parametrized by a lattice matrix
$\mathbf{L}=[\mathbf{a}\;\mathbf{b}\;\mathbf{c}]$, where $\mathbf{a}$,
$\mathbf{b}$, and $\mathbf{c}$ are vectors of lattice parameters in the
Cartesian coordinate. If the system is applied by strain, the lattice matrix
will be deformed by

$\displaystyle\mathbf{L}^{\prime}(\epsilon_{1},\epsilon_{2},...,\epsilon_{6})$
$\displaystyle=\bm{\upvarepsilon}_{I}\mathbf{L},$ (6)

where $\bm{\upvarepsilon}_{I}$ is the strain matrix

$\bm{\upvarepsilon}_{I}=\begin{bmatrix}\epsilon_{1}&\frac{\epsilon_{6}}{2}&\frac{\epsilon_{5}}{2}\\\\[3.0pt]
\frac{\epsilon_{6}}{2}&\epsilon_{2}&\frac{\epsilon_{4}}{2}\\\\[3.0pt]
\frac{\epsilon_{5}}{2}&\frac{\epsilon_{4}}{2}&\epsilon_{3}\\\
\end{bmatrix}+\mathbf{I}.$ (7)

In this work, we assume the crystal structures will be deformed by at most two
strain components, so that
$\mathbf{L}^{\prime}=\mathbf{L}^{\prime}(\epsilon_{i},\epsilon_{j}).$

Due to an applied strain, the atomic coordinates $\mathbf{r}$ must also be
transformed accordingly as

$\mathbf{r}^{\prime}=\bm{\upvarepsilon}_{I}\mathbf{L}\mathbf{r}_{f},$ (8)

where $\mathbf{r}_{f}$ is a fractional coordinate of an atom in an unstrained
lattice. Noting that $\mathbf{r}_{f}$ typically shifts from its equilibrium
value when the strain is applied (say, in density function theory or in
experiments); however, the atomic relaxation is neglected for simplicity.

Table 1: The minimal set of $ij$ strained component of each crystal system used to generate an initial crystal graph input Crystal systems | # of training indices | Training indices
---|---|---
cubic | 6 | 11, 12, 14, 15, 44, 45
hexagonal | 13 | 11, 12, 13, 14, 15, 16, 33, 34, 36, 44, 45, 46, 66
tetragonal | 14 | 11, 12, 13, 14, 15, 16, 26, 33, 34, 36, 44, 45, 46, 66
trigonal | 16 | 11, 12, 13, 14, 15, 16, 24, 25, 33, 34, 36, 44, 45, 46, 66
orthorhombic | 21 | all
monoclinic | 21 | all
triclinic | 21 | all

##  III Machine Learning with SE(3)-Equivariant Graph Neural Networks

###  III.1 Crystal graphs

A crystal structure can be represented as a multigraph whose nodes and edges,
respectively, encode atoms and their pairwise connections. A pair of nodes
describing an atom of type $m$ and an atom of type $n$ can be connected by
multiple edges, encapsulating the interactions between the atom of type $m$ in
a unit cell and atoms of type $n$ in the unit cell as well as in other
periodic cells of consideration, see Fig 1. An atom of type $m$ can interact
with an atom of the same type in periodic cells if the interaction range of
consideration is large enough, represented by self-loops in the node $m$ of
the crystal graph. Each atom’s location and its atomic features (e.g., atomic
mass, electronegativity, polarizability, atomic radius, and etc.) in a unit
cell is accounted for by the attributes of a single node.

More formally, a crystal graph $\mathcal{G}=(\mathcal{V},\mathcal{E}$)
consists of a set of nodes (vertices) $\mathcal{V}$ and a set of edges
$\mathcal{E}$ defined as

$\displaystyle\mathcal{V}=\\{(\mathbf{f}_{n},\mathbf{r}_{n})\;|\;\mathbf{f}_{n}\in\mathbb{R}^{M},\;\mathbf{r}_{n}=\mathbf{L}\mathbf{r}_{f_{n}}\in\mathbb{R}^{3}\\},$
$\displaystyle\mathcal{E}=\\{\Delta\mathbf{r}_{mn}^{(\mathbf{T})}\;|\;\Delta\mathbf{r}_{mn}^{(\mathbf{T})}=\mathbf{r}_{m}-\mathbf{r}_{n}+\mathbf{T};\
\mathbf{r}_{m},\mathbf{r}_{n}\in\mathbb{R}^{3}\\},$

where $m$ and $n$ are indices of the two atoms in the unit cell,
$\mathbf{f}_{n}$ is a vectorized atomic feature of an atom of type $n$, $M$ is
the number of atomic features, $\mathbf{r}_{f_{n}}$ is the fractional
coordinate of an atom of type $n$, and $\mathbf{T}$ is a translation vector
within the space of interest. Each node embeds two kinds of information:
atomic attributes $\mathbf{f}$ and atomic positions $\mathbf{r}$ in the unit
cell. From this definition, the unit cell and its periodic images is compactly
represented as a multigraph; only a finite number of nodes describing the
atoms in the unit cell is required, and the number of edges joining the two
nodes grows linearly with the number of periodic images of the atoms in the
structure of interest (see Fig. 1).

A crystal graph describing the whole crystal structure will require infinitely
many edges information. In practice, one typically embeds only partial
information of the structure using the description of the unit cell and a
finite number of edges encoding interactions with neighbor atoms [37]. One
method to construct the edges is to connect a target atom to other atoms
within a finite radius; however, still, this method generates an excessive
number of edges, which is computationally infeasible for GNNs machine
learning. Alternatively, some algorithms determine the connection between
atoms using, e.g., the Voronoi tessellation, generating a moderate number of
edges [38, 39]. However, for layered materials, using Voronoi tessellation
causes the edges connecting atoms between layers to be absent. These
interlayer connections are crucial in differentiating the structure strained
in out-of-plane direction from the structures strained in other directions. In
this work, we define an edge between a pair of atoms $\Delta\mathbf{r}$
(dropping the translation vector superscript and the atom indices subscript
for brevity) to be non-zero only if each component $\Delta r_{\beta}$ where
${\beta}\in\\{x,y,z\\}$ satisfies the following spatial cutoff criterion:
$|\Delta r_{\beta}|\leq\min(|a_{\beta}|+|b_{\beta}|+|c_{\beta}|,s)+\delta$
where $a_{\beta}$, $b_{\beta}$, and $c_{\beta}$ are the $\beta$ components of
the corresponding lattice parameters constituting the lattice matrix
$\mathbf{L}$, $s$ is a cutoff distance, and $\delta$ is a small cutoff
extension of our choice. This can reduce the number of edges from the hard
radial cutoff criterion and ensures the connections of atoms between layers
for layered materials by an appropriate choice of $s$ and $\delta$.

###  III.2 Rotational and translational equivariance

If two crystal graph inputs constructed from an identical crystal graph
strained by two different sets of strained components
$(\epsilon_{i},\epsilon_{j})$ and $(\epsilon_{k},\epsilon_{l})$ yielding the
new vertices
$\mathcal{V^{\prime}}=\\{(\mathbf{f}_{n},\mathbf{r^{\prime}}_{n})\\}$ and
$\mathcal{V^{\prime\prime}}=\\{(\mathbf{f}_{n},\mathbf{r^{\prime\prime}}_{n})\\}$,
respectively, are equivalent up to a 3D rotation and a translation, then we
expect our machine learning model to predict the same SED. This symmetry-aware
property can be implemented with geometric deep learning models [40]. We will
equip a latent representation of our model with this symmetry-aware
(equivariant) property, which is defined as follows.

A latent representation $\phi:\mathcal{V}\rightarrow\mathcal{Y}$ is
equivariant if for each transformation
$T_{g}:\mathcal{V}\rightarrow\mathcal{V}$ with $g$ being an element of an
abstract group $G$, there exists a transformation
$S_{g}:\mathcal{Y}\rightarrow\mathcal{Y}$ such that the following condition
holds:

$S_{g}[\phi(v)]=\phi\left(T_{g}[v]\right),$

for all $g\in G,v\in\mathcal{V}.$ In our case, the group $G$ of interest is
SE(3), providing the latent representation with a 3D rotation and translation
equivariant property. This latent feature $\phi$ can be achieved with
SE(3)-equivariant GNNs, known as Tensor Field Network (TFN) [41], and with a
TFN with an appropriate attention mechanism, known as SE(3)-transformers [42],
see more detailed recipes of these two models in Appendix A.

The key idea that enables SE(3)-equivariant property of these two GNNs is that
its message passing kernel is constructed from translational invariant
spherical harmonic bases, with the structure-preserving transformation $S_{g}$
given by the Wigner-D matrices (see Appendix A). Under a 3D rotation, each
spherical harmonic $J$ basis function of the kernel transforms according to

$Y_{J}\left(\mathbf{R}_{g}^{-1}\frac{\Delta\mathbf{r}}{||\Delta\mathbf{r}||}\right)=\mathbf{D}_{J}^{\ast}(g)Y_{J}\left(\frac{\Delta\mathbf{r}}{||\Delta\mathbf{r}||}\right),$

where $\mathbf{D}_{J}(g)$ is the $J^{th}$ Wigner-D matrix, and
$\mathbf{R}_{g}$ is a rotation matrix associated with $g\in\mathrm{SO}(3)$,
making the learned latent representation $\phi$ equivariant [41]. The multi-
head attention mechanism in SE(3)-transformers also induces an equivariant
property on the product rule between each key and value pairs, rendering the
whole message passing algorithm equivariant under SE(3) [42]. In this work, we
use SE(3)-transformers to first build a compressed representation of the
orientations of a strained crystal structure, then such learned latent
representation will be passed to other networks with a high expressive power
to perform the prediction (regression) task of the SED. The following section
contains the entire framework for the SED prediction.

Figure 1:  A 2-dimensional lattice and its strained crystal graph embedding.
Grey dashed lines are periodic boundaries. The whole structure can be spanned
by the 3 atoms in 3 different bases in the unit cell, represented by 3
different nodes. Black lines are (undirected) edges connecting two neighbor
nodes. The strained lattice is represented as a multigraph on the right, which
is the crystal graph input into our SE(3)-equivariant GNN building block. Note
that to describe the crystal graph for the unit cell and the atoms on its
boundary, only three nodes are required. The multi-edges of each pair of nodes
are distinguished by different crystallographic directions.

###  III.3 Model architecture and training procedure

Figure 2:  With a strained crystal graph and the strained component one-hot
vector as the input, StrainTensorNet employs a self-supervised approach
combining both classification and regression networks to predict both the
degeneracy class vector and the $U_{ij}$. The
$\mathbf{f}_{\mathrm{in},n}^{0}\equiv\mathbf{f}_{n}$ denotes the input feature
of the $n^{\textrm{th}}$ node, where a superscript, 0, indicates the input
into the $l=0$ channel of the spherical harmonics in SE(3)-Transformer, which
is a channel whose feature is invariant under rotation (transforms like a
scalar quantity). $\Delta\mathbf{r}^{\prime}_{mn}$ is the edge connecting a
target node $m$ and a neighbor node $n$. The degeneracy class vector from the
classification network will be concatenated with the latent representation of
the strained crystal graph ($\mathbf{f}^{0}_{\mathrm{reg},n}$) for the final
regression network to predict the SED.

We now describe our GNNs-based model that predicts the SED tensor. The
architecture of our model, which we term StrainTensorNet, is illustrated in
Fig. 2. The model is designed to receive two types of input: a crystal graph
representing a strained crystal structure, and a one-hot vector of dimension
$6(6+1)/2=21$ 222Due to the symmetry of SED, we identify only an upper
triangular element indicating the strained components, with a value 1 in the
$ij$ component that the strain operation is applied on and with a value 0
otherwise.

The crystal graph of a strained crystal structure is an input feature of the
SE(3)-equivariant GNNs building block, which generates two latent features.
The first latent feature (${\bf f}^{0}_{\textrm{class},n}$ in Fig. 2) is fed
into a classification neural network that predicts a vector of dimension $21$,
identifying all the upper-triangular elements of the SED tensor whose values
are exactly identical (by symmetry) to the component indicated by the input
one-hot vector. This inherent classification informs the model to discern the
degenerate structure of the SED tensor that depends on the symmetry of an
input crystal graph (see Fig. 3). For example, for a cubic lattice strained in
the direction $11$, the input one-hot vector $(1,0,\dots,0)^{T}$ together with
the latent feature ${\bf f}^{0}_{\textrm{class},n}$ shall give a prediction of
a vector that has the value close to 1 in the indices $ij=11,22,33$, and the
value close to 0 in the other 18 indices. Finally, the predicted degeneracy
class vector (DCV) together with the second latent feature (${\bf
f}^{0}_{\textrm{reg},n}$ in Fig. 2) will be fed into the final neural network
that predicts (regresses on) the SED.

To generate an SE(3)-equivariant latent representation of an input crystal
graph, as alluded to in the previous section, we first used the
SE(3)-Transformer that receives the strained crystal graph input. The material
structure from the database is first transformed according to Eq. (6) by a
strain of a fixed magnitude, which is then converted into the strained crystal
graph input. Noting that in the ab initio calculation, the atomic coordinates
in the strained lattice are optimized to be in the equilibrium positions under
an applied strain, but in this work, the fractional coordinates are kept as in
pre-strain condition. The strained crystal graph input is given by

$\displaystyle\mathcal{V}^{\prime}=\\{(\mathbf{f}_{n},\mathbf{r}^{\prime}_{n})\;|\;\mathbf{f}_{n}\in\mathbb{R}^{M},\;\mathbf{r}^{\prime}_{n}=\bm{\upvarepsilon}_{I}\mathbf{L}\mathbf{r}_{f_{n}}\in\mathbb{R}^{3}\\},$
$\displaystyle\mathcal{E}^{\prime}=\\{\Delta\mathbf{r}_{mn}^{\prime(\mathbf{T}^{\prime})}\;|\;\Delta\mathbf{r}_{mn}^{\prime(\mathbf{T}^{\prime})}=\mathbf{r}^{\prime}_{m}-\mathbf{r}^{\prime}_{n}+\mathbf{T}^{\prime};\
\mathbf{r}^{\prime}_{m},\mathbf{r}^{\prime}_{n}\in\mathbb{R}^{3}\\},$

where the prime notation indicates that the atomic positions and the
translation vector are of the strained lattice.

The training data of the $\\{U_{ij}\\}$ is computed, via Eq. (2), from
$\\{C_{ij}\\}$ extracted from Materials Project database [24] (see Appendix
E). For molecules and materials, if the number of training data is large
enough, only atomic numbers are sufficient for the node (atomic) attributes
[44]. Nevertheless, in this work, the features such as atomic mass,
electronegativity, polarizibility, atomic radius, atomic number, and $ij$
indices, are used as node attributes. Polarizability data is obtained from
[45], whereas other attributes are obtained from Materials Project database.
The atomic number is vectorized through an embedding layer of dimension 512,
which is then concatenated with the other four node attributes (see Fig. 2).
Note that the elastic constants are influenced by the atomic mass since the
phonon frequencies in sound waves, particularly the acoustic modes, are
inversely proportional to the square root of atomic mass. Additionally, the
resistance to changes in bond length correlates with both electronegativity
and polarizability. These three atomic properties are thus incorporated into
node attributes, and they significantly improve the model prediction accuracy.

Since the atomic features are scalar attributes, the input feature vector will
be fed into the $l=0$ channel of the transformer network since such channel
transforms like a scalar (see Appendix. A). We’ll denote each node’s input
feature vector as ${\bf f}^{0}_{\textrm{in},n}={\bf f}_{n}$. For computational
feasibility, we will restrict the output feature vector of the
SE(3)-transformers to consist of 4 channels ${\bf f}^{l}_{\textrm{hid},n}$
with $l\in\\{0,1,2,3\\}$. These latent features are fed into two different
SE(3)-equivariant GNNs without attention (TFNs), outputing a node-wise
classification feature vector ${\bf f}^{0}_{\textrm{class},n}$, and the node-
wise regression feature vector ${\bf f}^{0}_{\textrm{reg},n}$, see Fig. 2.
Note that the attention mechanism is not deployed in the last GNN layer as it
yields a better prediction accuracy.

To train the model to discern different orientations of a strained crystal
graph and classify its appropriate degeneracy class vector, in each epoch, we
draw a new realization of a crystal graph in a different orientation and of a
new one-hot vector in the same class of vector associated with the original
strained component. Specifically, a new crystal graph input is sampled from a
uniformly random orientation (uniformly random Euler angles) around the center
of mass in the unit cell of the original strained crystal graph. A new one-hot
vector of the strained components is chosen such that the element 1 is drawn
uniformly from one of the non-zero components of the DCV. For example, if the
original strained component is $11$ for a cubic material, a new one-hot vector
in each epoch is drawn from the situations where only one of the strained
components 11, 22, or 33 is 1, while the other components are 0.

Our self-supervised approach combines both the regression errors and the
classification errors in the global loss function:

$\mathcal{L}=\frac{1}{N}\sum_{n=1}^{N}|U_{ij}^{(n)}-U_{ij}^{(n),\textrm{pred}}|+\lambda\mathcal{L}_{\textrm{class}},$
(9)

where the SED prediction $U_{ij}^{(n),\textrm{pred}}$ is the output from the
regression network, $\mathcal{L}_{\textrm{class}}$ is the binary cross-entropy
loss function for multi-label classification, $\lambda$ specifies the relative
significance between the classification and the regression errors, and $N$ is
the total number of training samples. We stop the training when the gradient
of the loss is relatively small over multiple epochs. Training was performed
on NVIDIA A100 GPUs.0

###  III.4 Choosing a subset of strained crystal structures for training data

By the virtue of rotational equivariance, we expect StrainTensorNet to
efficiently predict the 21 components of $\mathbf{U}$ using training data
consisting of only non-degenerate components. We have selected the number of
non-degenerate components for training that depends on crystal systems, i.e.
cubic, tetragonal, hexagonal, trigonal, monoclinic, and triclinic. For cubic
materials, only 6 components are used for training (the number of non-
degenerate components for other crystal systems can be found in Table 1).
Recall that, as discussed earlier, the cubic lattice has 5 distinct $U_{ij}$
as shown in Eq. (5); however, the strain tensor transform the cubic structures
into 6 distinct structures (according to the symmetry in Laue class
$m\bar{3}m$) as follows

$\begin{bmatrix}\mathcal{T}_{1}&\mathcal{T}_{2}&\mathcal{T}_{2}&\mathcal{O}_{2}&\mathcal{M}_{2}&\mathcal{M}_{2}\\\
&\mathcal{T}_{1}&\mathcal{T}_{2}&\mathcal{M}_{2}&\mathcal{O}_{2}&\mathcal{M}_{2}\\\
&&\mathcal{T}_{1}&\mathcal{M}_{2}&\mathcal{M}_{2}&\mathcal{O}_{2}\\\
&&&\mathcal{O}_{1}&\mathcal{M}_{1}&\mathcal{M}_{1}\\\
&&&&\mathcal{O}_{1}&\mathcal{M}_{1}\\\ &&&&&\mathcal{O}_{1}\end{bmatrix},$

where each element in the matrix stands for a crystal system of
$\mathbf{L}^{\prime}(\epsilon_{i},\epsilon_{j})$, $\mathcal{T}$,
$\mathcal{O}$, and $\mathcal{M}$ are tetragonal, orthorhombic, and monoclinic
lattices, respectively, and a different number labeling the subscript
indicates a different structure. Despite the fact that $U_{14}$ ($C_{14}$) is
equal to $U_{15}$ ($C_{15}$) by the cubic symmetry,
$\mathbf{L}^{\prime}(\epsilon_{1},\epsilon_{4})$ and
$\mathbf{L}^{\prime}(\epsilon_{1},\epsilon_{5})$ possess different lattice
symmetries which are orthorhombic and monoclinic lattices, respectively. This
is because $\epsilon_{1}$ strains the structure in the $x$ direction, while
$\epsilon_{4}$ and $\epsilon_{5}$ strain the structure in both the $y$ and the
$z$ directions, and both the $x$ and the $y$ directions, respectively. Since
$\mathbf{L}^{\prime}(\epsilon_{1},\epsilon_{4})$ and
$\mathbf{L}^{\prime}(\epsilon_{1},\epsilon_{5})$ are not equivalent up to a
rotation, the SE(3) kernel regards the two inputs as different inputs. By
exploiting SE(3) kernel to identify inputs that are equivalent up to a
rotation, for cubic materials, it suffices to use the training data consisting
of the distinct input structures, strained only in 6 directions, i.e., 11, 12,
14, 15, 44, and 45.

Figure 3:  The distance metrics $d(U_{ij},U_{kl})$ computed from the test set
of the ground-truth SED (top row) and the predicted SED (bottom row). The
columns organize the crystal structure, appearing from left to right is for
cubic, hexagonal, and monoclinic materials respectively. The model prediction
shows an excellent agreement with the ground-truth for crystals with high
symmetry (cubic), and a good agreement with the ground truths for crystals
with lower symmetry (hexagonal and monoclinic).

###  III.5 The distance between SED components

To evaluate our model capability to predict the symmetry-dependent degeneracy
pattern of the SED tensor, we use the Canberra distance between two components
of the SED tensor as a metric:

$d(U_{ij},U_{kl})\equiv\frac{1}{N}\sum_{n=1}^{N}\frac{|U_{ij}^{(n)}-U_{kl}^{(n)}|}{|U_{ij}^{(n)}|+|U_{kl}^{(n)}|},$
(10)

where $N$ represents the number of samples, and $ij$ and $kl$ are Voigt
indices. With this metric, the distance between $U_{ij}$ and $U_{kl}$ is zero
if they belong to the same degeneracy class (their values are identical.) The
top row of Fig.3 shows the ground-truth distance pattern for cubic, hexagonal,
and monoclinic materials. Dark blue and bright green colors indicate the
$d(U_{ij},U_{kl})$ values that are closest to zero and 0.6, respectively. It
is important to note that the degeneracy pattern of $U_{ij}$ depends on the
crystal symmetry, which results in a unique Canberra distance pattern for each
crystal system. For instance, in cubic materials, $U_{11}$, $U_{22}$, and
$U_{33}$ are identical (in the same degeneracy class), so the distance pattern
of cubic materials displays a dark blue box corresponding to indices 11, 22,
and 33. The distance patterns in Figs. 3 and 4 are averaged over the samples
with the same crystal system in the dataset. Hence, the distance pattern of
low-symmetry crystals, such as monoclinic and triclinic, will strongly depend
on the dataset.

##  IV Discussion and Results

###  IV.1 Rationale for the StrainTensorNet model architecture

The underlying principle for using an equivariant network as a core component
is to construct a compact representation of the SED tensor, such that the
crystal structures strained by different sets of $\epsilon_{i}$ and
$\epsilon_{j}$ that are equivalent up to rotation (and translation) possess
the same SED. Such symmetry-dependent SED tensor of specific crystal
structures is demonstrated in the ground-truth distance patterns
$d(U_{ij},U_{kl})$ of Fig. 3 (top), where the StrainTensorNet can well
approximate such symmetry-dependent patterns, see Fig. 3 (bottom).

Fig. G1 shows our earlier model development with the goal to compactly
represent the SED tensor, beginning from the simplest but less expressive to
StrainTensorNet which can compactly represent the SED tensor rather
accurately. All these trial GNN-based models were trained on crystal
structures with fewer than 500 edges. The simplest model (Fig. G1(a)) takes as
an input only the crystal graph with node attributes, i.e., atomic mass,
atomic radius, electronegativity, polarizability, and embedded atomic number.
This minimal equivariant model does not yield a sufficiently accurate
representation of $U_{ij}$; we found that the model is more biased towards
predicting the SED tensor whose distance matrix pattern is more akin to that
of the cubic crystals. This could be because the model is not sufficiently
expressive, and thus attempts to fit the majority of the training dataset
comprising more higher-symmetry structures (see Table 4 for the crystal system
statistics of our curated dataset).

To improve the minimal model’s expressiveness, we introduce the degeneracy
class vector (DCV) as an additional input, to inform the model about the
symmetry-dependent class of the SED tensor. The DCV is first embedded by
fully-connected layers and then concatenated with global max-pooled features
or node attributes, as shown in Fig. G1(b) and (c), respectively. These
symmetry-informed models yield lower MAE and RMSE of the SED tensor compared
to that of the minimal model, while the model in Fig. G1(b) performs slightly
better than the model in Fig. G1(c) (see Table 6). Moreover, the predicted
distance patterns are less biased towards those of the cubic lattice.
Specifically, for crystal structures with lower symmetry such as hexagonal,
tetragonal, trigonal, and orthorhombic, the predicted distance patterns are
relatively similar to their ground-truth distance patterns. As expected,
informing the model about the DCV of the input helps force the distance
between $U_{ij}$ and $U_{kl}$ belonging to the same degeneracy class to be
closer to zero.

While these models with a DCV included as an input give good predictions of
the SED tensor, they require the knowledge of the prestrained crystal symmetry
a priori (to properly assign the value of the DCV). To overcome this
limitation, we use a self-supervised method schematically shown Fig. 2 to
modify the model to predict, rather than to require, the DCV. The
classification neural network that predicts the DCV takes as an input a one-
hot vector of the strained $ij$ components, representing the $\epsilon_{i}$
and $\epsilon_{j}$ components in (7). This self-supervised technique enables
the StrainTensorNet to predict $U_{ij}$ without any prior knowledge of the
prestrained crystal symmetry. Importantly, StrainTensorNet can express the
degeneracy class of the SED tensor relatively well, see Figs. 3 and 4.

###  IV.2 StrainTensorNet’s prediction results

Table 2 summarizes StrainTensorNet’s prediction errors on the test set. To the
best of our knowledge, StrainTensorNet provides the first data-driven
prediction of $U_{ij}$ with a reasonable accuracy. Using these predicted SED
tensors together with Eq. (3), elastic tensors of materials can also be
calculated. Our model’s prediction accuracy significantly depends on the
elemental composition of the compounds. Fig. 5 depicts the MAE of $U_{ij}$ and
$C_{ij}$, averaged over 21 $ij$ components and categorized by elements. It is
interesting to note that the compounds containing period 2 elements (Be, B, C,
N, O, and F), some period 6 transition metals (Ta, W, Re, and Os), and some
actinides (Th, Pa, and Pu), which have high bulk and shear moduli on average,
exhibit a higher MAE for $U_{ij}$ compared to the overall MAE of the entire
test set. As a result, their MAE for $C_{ij}$ are also significantly higher
than the overall MAE of the test set. In contrast, the compounds containing
lanthanides, which exhibit medium bulk and shear moduli on average, except for
Gd, demonstrate a considerably lower MAE than the overall MAE of the test set,
for both $U_{ij}$ and $C_{ij}$.

To further investigate the prediction results of the elastic constants, we
plot the comparisons between the prediction (obtained from the predicted SED
tensors together with Eq. (3)) and the ground truth of the elastic constants
in Fig. 6. Fig. 6 (a) shows the prediction results of every $C_{ij}$
components of all crystal structures, whereas Figs. 6(b)-(f) show the
prediction results of $C_{22}$, $C_{23}$, $C_{26}$, $C_{55}$, and $C_{56}$ of
cubic materials. For cubic materials, the model is trained only on a non-
degenerate subset of 21 $U_{ij}$ components, specifically components 11, 12,
14, 15, 44, and 45. Notably, the model is able to predict other unseen
components reasonably well. However, the prediction of the elastic constant
components that are exactly zero is challenging. For example, the ground-truth
values of $C_{26}$ and $C_{56}$ of cubic materials are exactly 0 GPa, while
the means and standard deviations of these components from our model
predictions are 0.69 and 4.45 GPa, and -0.43 and 4.85 GPa respectively (see
Fig. G3 for the histograms of the prediction results). Although the standard
deviations are relatively large, the means are less than 1 GPa. It is worth
noting that the first-principles calculations also make some errors in these
values, but these are usually not computed since they are known to be exactly
0 GPa.

To see how well StrainTensorNet discerns the geometric relationship between
different strained crystal graph inputs, we plot the Canberra distance metric
of the predicted and ground-truth SED tensors for comparison in Fig. 3. The
predicted distance metric of cubic materials, as shown in the left column of
Fig. 3, closely resembles the ground-truth distance metric, although there are
minor discrepancies in the values that belong to the same degeneracy class.
For instance, the distances between the predicted $U_{11}$, $U_{22}$, and
$U_{33}$ are $4.8-6.2$ $\mu$eV/atom, resulting in the percentage errors of
$1.83-2.43\%$ (computed from
$\frac{1}{N}\sum_{n=1}^{N}|U_{ij}^{(n)}-U_{kl}^{(n)}|/U_{avg}^{(n)}\times
100\%$ where $U_{avg}^{(n)}=(|U_{ij}^{(n)}|+|U_{kl}^{(n)}|)/2$.) For hexagonal
materials, $U_{11}$ are identical to $U_{22}$, but not to $U_{33}$. The
predicted distance metric of hexagonal materials accordingly forms a dark blue
box in the 11 and 22 components. The percentage error between the predicted
$U_{11}$ and $U_{22}$ is 2.22%. For monoclinic materials, the values of each
$U_{ij}$ are not necessarily identical to one another, and the distance metric
is averaged over monoclinic materials data. The resulting distance pattern may
differ significantly depending on the dataset. Notably, these agreements
between the predicted and the ground-truth distance metric in various crystal
structures reveal that StrainTensorNet gives an excellent prediction of the
degenerate structure of SED tensors of strained crystalline materials.

Figure 4:  The distance metric (computed from the test set) of tetragonal,
trigonal, orthorhombic, and triclinic materials (top to bottom respectively)
for (left) the ground-truth and (right) StrainTensorNet’s prediction.

Lastly, the model can also predict other elastic properties, including Voigt’s
bulk modulus ($B_{V}$), Voigt’s shear modulus ($G_{V}$), Young’s modulus
($Y$), and Poisson’s ratio ($\nu$), using the predicted $C_{ij}$ together with
Eqs. (19)$-$(22). Table 2 summarizes the MAE and RMSE of our predicted elastic
properties compared to previous works. The MAE and RMSE of $B_{V}$ are 12.59
and 20.83 GPa, respectively, which are comparable to those of Mazhnik et al.
[25]. The MAE and RMSE of $G_{V}$ are 10.49 and 18.95 GPa, respectively, with
the RMSE of $G_{V}$ being comparable to that reported by Zhao et al. [26]. Our
work also yields smaller errors for $\nu$ than Mazhnik et al. [25] Note that
our dataset (see Sec. E) differs from the dataset in the relevant work of [25,
26]. Fig. G4 presents the MAE of $B_{V}$, $G_{V}$, $Y$ and $\nu$, categorized
by elements. The model produces high MAE for compounds containing period 2
elements, mainly due to their larger MAE in $U_{ij}$ and $C_{ij}$ (see Fig.
5). On the other hand, the model can produce smaller MAE for lanthanide
compounds, as their MAE in $U_{ij}$ and $C_{ij}$ (see Fig. 5) are smaller.

###  IV.3 Interpretability of StrainTensorNet’s latent features

To investigate how latent features facilitate successful prediction of the SED
tensor, we employed the diffusion map to faithfully visualize the three-
dimensional data representation of the high-dimensional latent features of the
entire test set [46, 47]. Figs. 7 and 8 reveal the dimensionality reduction of
the latent features in the diffusion coordinates
$(\varphi_{1},\varphi_{2},\varphi_{3})$. Two different latent features were
considered: the global max-pooling layer from the GNN representation of the
input crystal graph (denoted as A in Fig. 2), and the concatenation between A
and the embedded DCV (denoted as B in Fig. 2). These low-dimensional
representations in Figs. 7 and 8 are colored by the energy scale (left column)
and by the strained component ${ij}$ (middle and right columns).

What does the max-pooled feature of the crystal graph input A represent? Fig.
8(c) shows that the crystal graph latent representations of the same material
that are strained in different $ij$ components are almost identical, in fact
sharing greater resemblance than the representations of different materials
that belong to the same symmetry class. This is consistent with the result of
Bronstein et al. which shows that the latent representation in
SE(3)-equivariant GNNs of the original input is very similar to those of the
mildly spatially distorted inputs, and, interestingly, is less similar to
those of inputs spatially translated further away [40]. Hence, the latent
feature A seems to uniquely encode the information of the input material,
rather than the symmetry of the input. Thus A alone does not suffice to
meaningfully represent the symmetry-dependent SED tensor (Fig. 8(a) and (b)).

On the other hand, with the embedded DCV in the latent feature B, the network
can differentiate different strained components within the same material
class. When strained in 21 different components, the latent features A of
three example materials (trigonal LuTlTe2 phase, tetragonal YMnSi phase, and
cubic Ta3Ru phase) that were not differentiable within the same material class
(Fig. 8(c)) are now clearly segregated in the $\varphi_{2}$ variable through
the help of the embedded DCV, see Fig. 7(c). In fact, the latent feature B
which combines both the input crystal graph representation and the embedded
DCV offers an interpretable latent representation of the SED (Fig. 7), which
we now discuss.

The latent feature B is organized such that the $(\varphi_{1},\varphi_{2})$
coordinate encodes the strain energy density that varies from a smaller value
in the $(+,+)$ quadrant to a larger value in the $(-,-)$ quadrant (see Fig.
7(a)). Since $U_{ij}$ stored by the tensile strain ($i,j\in\\{1,2,3\\}$) is
typically higher than $U_{ij}$ stored by the shear strain
($i,j\in\\{4,5,6\\}$), Figs. 7(a) and (b) consistently reveal that the latent
features of the materials strained by the $11$, $22$, or $33$ component have a
negative $\varphi_{2}$ value, whereas those of the materials strained by the
$44$, $55$, or $66$ components have a positive $\varphi_{2}$ value. In
addition, since $U_{ij}$ for $i\neq j$ is computed from the sum of $C_{ii}$,
$C_{jj}$, and $C_{ij}$, it is typically larger than $U_{ii}$ and $U_{jj}$.
Figs. 7(a) and (b) also consistently reveal that the latent features of the
materials strained by 44, 45, 11, and 12 components are respectively organized
in the $\varphi_{2}$ coordinate from a more positive to a more negative value.
Additionally, Figs 7(a) and (c) show that materials with higher (lower)
average SED will be represented by a larger negative (larger positive)
$\varphi_{1}$ variable. In summary, the max-pooled feature from the graph
neural networks A encodes the material information, and, together with the
information of the DCV, the concatenated latent feature B effectively
encapsulates both the material information and its strained component-
dependent SED.

Figure 5:  The bar charts display the MAE of the predicted $U_{ij}$ and
$C_{ij}$ components, averaged over 21 $ij$ components. The charts are
categorized based on the elements present in the compounds in the test set and
are colored according to the average bulk modulus for (a) and (c), and average
shear modulus for (b) and (d), as indicated in the legend below. The dashed
lines represent the MAE of each property, averaged over the data in the test
set. It can be seen that most elements, whose compounds have high bulk and
shear moduli, exhibit errors of $U_{ij}$ and $C_{ij}$ that exceed the
dataset’s MAE.

Table 2:  Statistics of the dataset and the model prediction Properties | This work | Mazhnik et al. [25] | Zhao et al. [26]
---|---|---|---
Average | MAE | RMSE | Average | MAE | RMSE | RMSE
$U_{ij}$ (meV/atom) | 2.652 | $\mathbf{0.655}$ | $\mathbf{1.288}$ |  |  |  |
$C_{ij}$ (GPa) | 42.92 | $\mathbf{10.37}$ | $\mathbf{17.69}$ |  |  |  |
$B_{V}$ (GPa) | 107.07 | 11.48 | 19.73 | 111.83 | 11.11 | 19.54 | 16.530
$G_{V}$ (GPa) | 50.98 | 9.61 | 17.10 | 54.81 | 8.24 | 11.43 | 15.780
$Y$ | 129.62 | 22.29 | 38.63 | 138.95 | 19.15 | 26.23 |
$\nu$ | 0.401 | $\mathbf{0.037}$ | 0.133 | 0.286 | 0.041 | 0.105 |
Figure 6:  The predicted values (vertical axes) versus the ground-truth values
(horizontal axes) of the elastic constants. (a) All the components of the
elastic constants of the whole test set of all crystal systems. The results
for only cubic materials are shown in (b) $C_{22}$, (c) $C_{23}$, (d)
$C_{26}$, (e) $C_{55}$, and (f) $C_{56}$, where we note that the 22, 23, 26,
55, and 56 indices did not appear in the input of the training set for cubic
materials.  Figure 7:  Three-dimensional diffusion map representation of the
concatenated latent features that are fed into the SED tensor prediction layer
(the latent feature B in Fig. 2). (a) is colored by the SED scale, and (b) is
colored by the strained component $ij$. (c) highlights the latent features of
3 example materials colored by the strained component $ij$ as in (b); the
features of other materials are shown in transparent grey. The 3 materials are
trigonal LuTlTe2 phase ( ), tetragonal YMnSi phase ($\triangle$), and cubic
Ta3Ru phase ($\square$), which have low, medium, and high elastic constants,
respectively. These are all the latent features in the test dataset. More
elaborated discussion can be found in Sec. IV.3.  Figure 8:  Three-dimensional
diffusion map representation of the global max-pooled latent features of the
graph neural networks (i.e., the latent feature A in Fig. 2). (a) is colored
by the SED scale, and (b) is colored by the strained component $ij$. (c)
highlights the latent features of 3 example materials colored by the strained
component $ij$ as in (b); the features of other materials are shown in
transparent grey. The 3 materials are trigonal LuTlTe2 phase ( ), tetragonal
YMnSi phase ($\triangle$), and cubic Ta3Ru phase ($\square$), which have low,
medium, and high elastic constants, respectively. See the discussion and the
interpretation of these plots in Sec. IV.3.

##  V Conclusions

We have demonstrated a novel data-driven framework for predicting the elastic
properties of crystal structures using SE(3)-equivariant graph neural networks
(GNNs). By leveraging SE(3)-equivariant GNNs as building blocks, our self-
supervised deep learning model accurately predicts Voigt’s bulk modulus,
Voigt’s shear modulus, Young’s modulus, and Poisson’s ratio that are
comparable to other non-equivariant deep learning studies [25, 26].

A key contribution is the prediction of the strain energy density (SED) and
its associated elastic constants, which are tensorial quantities that depend
on a material’s crystallographic group. The similarity between the distance
metrics of the SED components between the ground truths and the predictions
demonstrates our model’s capability to identify different symmetry groups of
strained crystal structures. Requiring only a strained crystal graph and the
strained component as the input, our approach offers an efficient alternative
to the standard ab initio method for designing new materials with tailored
elastic properties.

The interpretability of the model is also a notable feature. The learned
latent representations taking into account the degeneracy class vector are
organized in a physically meaningful structure for predicting the SED tensors.
This interpretability aspect enhances the transparency of model prediction,
enabling the justification of whether the prediction is physically relevant.

The combination of interpretability and the consideration of crystallographic
groups sets our model apart from recent data-driven methods for predicting
elastic properties of materials. We hope this work is a stepping stone toward
an efficient data-driven approach for materials discovery and design, and
opens up avenues for approaching more challenging tensor prediction tasks,
e.g., predicting dielectric and piezoelectric tensors, which are second-rank
and third-rank tensors, respectively.

###### Acknowledgements.

This research is supported by the Program Management Unit for Human Resources
and Institutional Development, Research and Innovation (Grant No. B05F650024)
and by Thailand Science research and Innovation Fund Chulalongkorn University
(CU_FRB65_ind (5)_110_23_40). The authors acknowledge high performance
computing resources including NVIDIA A100 GPUs from Chula Intelligent and
Complex Systems Lab, Faculty of Science, and from the Center for AI in
Medicine (CU-AIM), Faculty of Medicine, Chulalongkorn University, Thailand.
The authors also acknowledge the National Science and Technology Development
Agency, National e-Science Infrastructure Consortium, Chulalongkorn University
and the Chulalongkorn Academic Advancement into Its 2nd Century Project
(Thailand) for providing computing infrastructure.

##  Appendix A Ingredients of SE(3)-Transformers

Here we summarize the computational building block of SE(3)-equivariant GNNs
with a built-in attention mechanism (SE(3)-Transformers). More detailed
discussion of this deep learning architecture can be found in the original
proposal [42].

In an SE(3)-equivariant GNN without attention, also known as Tensor Field
Network (TFN) [41], given a multi-channel input feature of a node $j^{th}$,
$\mathbf{f}_{\text{in},j}^{k}$, of a crystal graph, the GNN will update each
crystal graph input feature by the following message-passing convolution

$\mathbf{f}_{\text{out},i}^{l}=\mathbf{W}^{ll}\mathbf{f}_{\text{in},i}^{l}+\sum_{k\geq
0}{\sum_{j\neq
i}^{N}{\mathbf{W}^{lk}\left(\mathbf{r}_{j}-\mathbf{r}_{i}\right)\mathbf{f}_{\text{in},j}^{k}}\
},$ (11)

where $i$ and $j$ are respectively the indices of a target node and its
neighbouring nodes, $l$ and $k$ specify type-$l$ and type-$k$ feature vectors
which respectively transform under SO(3) rotation according to the
$(2l+1)\times(2l+1)$ and $(2k+1)\times(2k+1)$ Wigner-D matrices (type-0 vector
is invariant under SO(3) rotation), $\mathbf{W}^{ll}$ is the diagonal matrix
whose diagonal is identically parametrized by a trainable parameter $w^{ll}$
describing a self-interaction message passing, and $\mathbf{W}^{lk}$ is a
convolution kernel that propagates the channel-$k$ feature from neighboring
nodes to the channel-$l$ feature of the node of interest such that

$\mathbf{W}^{lk}(\mathbf{r}_{j}-\mathbf{r}_{i})=\sum_{J=|k-l|}^{k+l}\varphi_{J}^{lk}\left(||\mathbf{r}_{j}-\mathbf{r}_{i}||\right)\mathbf{W}_{J}^{lk}(\mathbf{r}_{j}-\mathbf{r}_{i}),$
(12)

with

$\mathbf{W}_{J}^{lk}(\mathbf{r}_{j}-\mathbf{r}_{i})=\sum_{m=-J}^{J}Y_{Jm}\left(\frac{\mathbf{r}_{j}-\mathbf{r}_{i}}{||\mathbf{r}_{j}-\mathbf{r}_{i}||}\right)\mathbf{Q}_{Jm}^{lk},$
(13)

where $\varphi_{J}^{lk}\left(||\mathbf{r}||\right)$ is a learnable radial
function, $Y_{Jm}$ is the $m^{\text{th}}$ component ($m$ magnetic quantum
number in quantum mechanics) of the spherical harmonic function
$Y_{J}:\mathbb{R}^{3}\rightarrow\mathbb{R}^{2J+1}$, $\mathbf{Q}_{Jm}^{lk}$ is
a Clebsch-Gordan matrix with dimension $(2l+1)\times(2k+1)$. The
SO(3)-equivariant basis kernel $\mathbf{W}^{lk}_{J}$ imposes a fixed
functional constraint on the convolution kernel along the angular direction,
whereas the radial direction is still trainable via $\varphi_{J}^{lk}$

The SE(3)-Transformer introduces the attention $\alpha_{ij}$ to each edge of a
TFN, so that the message passing rule of Eq. (11) is weighted by the
attentions. Without isolating out the self-interaction term, for brevity, this
weighted message passing rule can be compactly expressed as

$\mathbf{f}_{\text{out},i}^{l}=\sum_{j}\alpha_{ij}\mathbf{v}_{ij}^{l},$ (14)

where the channel-$l$ value message, $\mathbf{v}_{ij}^{l}$, is the same as
that of the above TFN, defined as

$\mathbf{v}_{ij}^{l}=\sum_{k\geq
0}\mathbf{W}_{V}^{lk}(\mathbf{r}_{j}-\mathbf{r}_{i})\mathbf{f}_{\text{in},j}^{k},$
(15)

and the attention on the edge $ij$ (which is invariant under SO(3) [42]) is

$\alpha_{ij}=\frac{\mathrm{exp}(\mathbf{q}_{i}^{\top}\mathbf{k}_{ij})}{\sum_{j^{\prime}}\mathrm{exp}(\mathbf{q}_{i}^{\top}\mathbf{k}_{ij^{\prime}})},$
(16)

with a query vector on a node $i$, $\mathbf{q}_{i}$, and a key vector on the
edge $ij$, $\mathbf{k}_{ij}$, given by, respectively,

$\displaystyle\mathbf{q}_{i}$ $\displaystyle=\bigoplus_{l\geq 0}\sum_{k\geq
0}\mathbf{W}_{Q}^{lk}\mathbf{f}_{\text{in},i}^{k},$ (17)
$\displaystyle\mathbf{k}_{ij}$ $\displaystyle=\bigoplus_{l\geq 0}\sum_{k\geq
0}\mathbf{W}_{K}^{lk}(\mathbf{r}_{j}-\mathbf{r}_{i})\mathbf{f}_{\text{in},j}^{k}.$
(18)

$\bigoplus_{l\geq 0}$ denotes the direct sum over the $l$ channels. Fig. A1
depicts the message passing with the attention mechanism of the
SE(3)-Transformer model. Detailed numerical implementation of these models can
be found in https://github.com/trachote/predict_elastic_tensor.

Figure A1:  A schematic of a GNN message passing rule that aggregates features
from neighboring nodes of a crystal graph with an attention mechanism. A white
square on each edge depicts an attention head. A target node (pink circle)
sends a query vector of the node $i$, ($\mathbf{q}_{i}$, pink arrow), to match
with key vectors ($\mathbf{k}_{ij}$, dotted purple arrow) and the value
message ($\mathbf{v}_{ij}$, purple arrow) from neighboring nodes (blue and
orange circles). The angular momentum channels are not illustrated in the
visualization.

##  Appendix B Macroscopic elastic moduli

Here are the formulae of elastic moduli computable from the elastic constants.

Voigt’s bulk modulus:

$B_{V}=\frac{1}{9}\sum_{i,j=1}^{3}C_{ij}.$ (19)

Voigt’s shear modulus:

$G_{V}=\frac{1}{30}\sum_{i,j=1}^{3}(C_{ii}-C_{ij})+\frac{1}{5}\sum_{i=4}^{6}C_{ii}.$
(20)

Young’s modulus:

$Y=\frac{9BG}{3B+G}.$ (21)

The Possion’s ratio:

$\nu=\frac{3B-2G}{6B+2G}.$ (22)

##  Appendix C Symmetry in the strain energy tensor

The strain energy density (SED) tensor ($\mathbf{U}$) is derived from the
elastic tensor ($\mathbf{C}$). The symmetry in $\mathbf{U}$ follows the
symmetry in $\mathbf{C}$ which depends on the crystal system.

###  C.1 Hexagonal and Tetragonal I
(for Laue class $4/mmm$)

$\mathbf{C}_{\textrm{hex/tet}}=\begin{bmatrix}C_{11}&C_{12}&C_{13}&0&0&0\\\\[3.0pt]
&C_{11}&C_{13}&0&0&0\\\\[3.0pt] &&C_{33}&0&0&0\\\\[3.0pt]
&&&C_{44}&0&0\\\\[3.0pt] &&&&C_{44}&0\\\\[3.0pt] &&&&&C_{66}\end{bmatrix}.$
(23)
$\mathbf{U}_{\textrm{hex/tet}}=\begin{bmatrix}U_{11}&U_{12}&U_{13}&U_{14}&U_{14}&U_{16}\\\\[3.0pt]
&U_{11}&U_{13}&U_{14}&U_{14}&U_{16}\\\\[3.0pt]
&&U_{33}&U_{34}&U_{34}&U_{36}\\\\[3.0pt] &&&U_{44}&U_{45}&U_{46}\\\\[3.0pt]
&&&&U_{44}&U_{46}\\\\[3.0pt] &&&&&U_{66}\end{bmatrix}.$ (24)

###  C.2 Tetragonal II (for Laue class $4/m$)

$\mathbf{C}_{\textrm{tetII}}=\begin{bmatrix}C_{11}&C_{12}&C_{13}&0&0&C_{16}\\\\[3.0pt]
&C_{11}&C_{13}&0&0&-C_{16}\\\\[3.0pt] &&C_{33}&0&0&0\\\\[3.0pt]
&&&C_{44}&0&0\\\\[3.0pt] &&&&C_{44}&0\\\\[3.0pt] &&&&&C_{66}\end{bmatrix}.$
(25)
$\mathbf{U}_{\textrm{tetII}}=\begin{bmatrix}U_{11}&U_{12}&U_{13}&U_{14}&U_{14}&U_{16}\\\\[3.0pt]
&U_{11}&U_{13}&U_{14}&U_{14}&U_{16}^{{}^{\prime}}\\\\[3.0pt]
&&U_{33}&U_{34}&U_{34}&U_{36}\\\\[3.0pt] &&&U_{44}&U_{45}&U_{46}\\\\[3.0pt]
&&&&U_{44}&U_{46}\\\\[3.0pt] &&&&&U_{66}\end{bmatrix},$ (26)

where
$U_{16}^{{}^{\prime}}=-C_{16}\epsilon_{1}\epsilon_{6}+\frac{1}{2}(C_{11}\epsilon_{1}^{2}+C_{66}\epsilon_{6}^{2})$.

###  C.3 Rhombohedral I (for Laue class $\bar{3}m$)

$\mathbf{C}_{\textrm{rhomI}}=\begin{bmatrix}C_{11}&C_{12}&C_{13}&C_{14}&0&0\\\\[3.0pt]
&C_{11}&C_{13}&-C_{14}&0&0\\\\[3.0pt] &&C_{33}&0&0&0\\\\[3.0pt]
&&&C_{44}&0&0\\\\[3.0pt] &&&&C_{44}&C_{14}\\\\[3.0pt]
&&&&&C_{66}\end{bmatrix}.$ (27)
$\mathbf{U}_{\textrm{rhomI}}=\begin{bmatrix}U_{11}&U_{12}&U_{13}&U_{14}&U_{14}^{{}^{\prime\prime}}&U_{16}\\\\[3.0pt]
&U_{11}&U_{13}&U_{14}^{{}^{\prime}}&U_{14}^{{}^{\prime\prime}}&U_{16}\\\\[3.0pt]
&&U_{33}&U_{34}&U_{34}&U_{36}\\\\[3.0pt] &&&U_{44}&U_{45}&U_{46}\\\\[3.0pt]
&&&&U_{44}&U_{46}^{{}^{\prime}}\\\\[3.0pt] &&&&&U_{66}\end{bmatrix},$ (28)

where
$U_{14}^{{}^{\prime}}=-C_{14}\epsilon_{1}\epsilon_{4}+\frac{1}{2}(C_{11}\epsilon_{1}^{2}+C_{44}\epsilon_{4}^{2})$,
$U_{14}^{{}^{\prime\prime}}=\frac{1}{2}(C_{11}\epsilon_{1}^{2}+C_{44}\epsilon_{4}^{2})$,
$U_{46}^{{}^{\prime}}=C_{14}\epsilon_{4}\epsilon_{6}+\frac{1}{2}(C_{44}\epsilon_{4}^{2}+C_{66}\epsilon_{6}^{2})$.

###  C.4 Rhombohedral II (for Laue class $3$)

$\mathbf{C}_{\textrm{rhomII}}=\begin{bmatrix}C_{11}&C_{12}&C_{13}&C_{14}&C_{15}&0\\\\[3.0pt]
&C_{11}&C_{13}&-C_{14}&-C_{15}&0\\\\[3.0pt] &&C_{33}&0&0&0\\\\[3.0pt]
&&&C_{44}&0&-C_{15}\\\\[3.0pt] &&&&C_{44}&C_{14}\\\\[3.0pt]
&&&&&C_{66}\end{bmatrix}.$ (29)
$\mathbf{U}_{\textrm{rhomII}}=\begin{bmatrix}U_{11}&U_{12}&U_{13}&U_{14}&U_{15}&U_{16}\\\\[3.0pt]
&U_{11}&U_{13}&U_{14}^{{}^{\prime}}&U_{15}^{{}^{\prime}}&U_{16}\\\\[3.0pt]
&&U_{33}&U_{34}&U_{34}&U_{36}\\\\[3.0pt]
&&&U_{44}&U_{45}&U_{46}^{{}^{\prime\prime}}\\\\[3.0pt]
&&&&U_{44}&U_{46}^{{}^{\prime}}\\\\[3.0pt] &&&&&U_{66}\end{bmatrix},$ (30)

where
$U_{46}^{{}^{\prime\prime}}=-C_{15}\epsilon_{4}\epsilon_{6}+\frac{1}{2}(C_{44}\epsilon_{4}^{2}+C_{66}\epsilon_{6}^{2})$.

###  C.5 Orthorhombic, Monoclinic, and Triclinic

Strain energy tensors of orthorhombic, monoclinic, and triclinic lattices do
not possess any degenerate structure, and are parametrized by the most general
tensor of the form (4) in the main text.

##  Appendix D The minimal set of strained components for each crystal systems
in a training data set

The equivariant property of the SE(3)-transformer should be able to identify
strained crystal graph inputs that are equivalent up to a rotation, thus to
reduce the number of training data, we require one independent $ij$ component
of the SED tensor for each degeneracy class vector. The other strained
components can be incorporated during the randomized update of the one-hot
strained component vector within the same degeneracy class in each training
epoch, see Sec. III.3. The minimal set of $ij$ indices of each crystal system
can be chosen from the degenerate structure of the SED tensor for each crystal
system in the previous Appendix C. Our chosen set is presented in Table 1.

##  Appendix E Dataset

Material data is pooled from the Materials Project database [24]. There are
12,105 materials in the database that have the elastic tensor. However, only
some materials are used in this work due to memory restriction. The memory
usage increases with the number of edges, so the computationally feasible data
is restricted to the materials with the maximum number of edges of 2,000
constructed by a cell-radii approach. There are 9,296 out of 12,105 materials
from the database with this property (about 76.79%). Fig. E2 shows that the
selected materials have bulk and shear moduli distributing in the range of
$0-600$ GPa, which can be a good representative of the whole database. In
addition, we have found that the validation errors do not decrease
significantly even when the number of edges increases to 3,000. The statistics
of crystal systems curated for training and test sets of these 9,296 materials
is shown is Table 4. The dataset is partitioned into a training set and a
validation set using the stratified method (in scikit-learn) with a random
seed number 0.

Figure E1:  The histograms of (a) $C_{ij}$ and of (b) $U_{ij}$ in the training set and of (c) $C_{ij}$ and of (d) $U_{ij}$ in the test set. For visualization purposes, only some $ij$ components, specifically 11, 12, 14, 33, 44, and 45, are presented. For each $ij$ component, $U_{ij}$ ($C_{ij}$) are distributed over different range of energy (modulus). For $i\neq j$, $U_{ij}$ are distributed over a broader range than $U_{ii}$ are. For $ij=45$, the elastic constant is concentrated around 0 GPa. Importantly, this concentration around zero issue can be ameliorated when the elastic constants are converted to the strain energy density. This conversion to SED reduces the number of zero values in the dataset, which significantly improves model training and model predictive performances. Figure E2: The bulk modulus of materials is plotted against the number of edges created by cell radii method, with the hue indicating the shear modulus. The data below a hard cutoff of 2000 edges can already cover the range of values of bulk and shear moduli of the whole dataset. Table 3: Statistics of number of edges for 12105 materials created by cell radii and CrystalNN methods. | Cell radii | CrystalNN
---|---|---
mean | 2115.16 | 143.47
std | 4665.66 | 129.67
min | 6 | 6
25% | 368 | 58
50% | 1024 | 112
75% | 1776 | 224
max | 75864 | 2272
Table 4: The number of data for training and test sets categorized by crystal systems. Crystal system | Training set | Test set
---|---|---
| Count | % | Count | %
Cubic | 3865 | 46.17 | 441 | 47.83
Hexagonal | 1369 | 16.35 | 146 | 15.84
Tetragonal | 1580 | 18.87 | 163 | 17.68
Trigonal | 530 | 6.33 | 60 | 6.51
Orthorhombic | 731 | 8.73 | 79 | 8.57
Monoclinic | 253 | 3.02 | 27 | 2.93
Triclinic | 43 | 0.51 | 6 | 0.65

##  Appendix F Model parameters

The hyperparameters that significantly affect the model performance are the
number of convolution layers, the number of $l$-channels in the SE(3) kernel,
and the number of hidden units of the learnable radial function in the kernel
$\phi$. The number of convolution layers is 4, where the model performance
only slightly improves with more layers. The learnable radial function $\phi$
is a feedforward network, consisting of an input layer that is fed into a
fully-connected hidden layer with a ReLU activation, which is forwarded to
another fully-connected layer that outputs a hidden feature vector. We use the
number of hidden units in the hidden layer of $\phi$ to be 128 to avoid over-
squashing, and the performance only marginally improves expanding beyond 128
hidden units. The dimension of the hidden feature vectors
$\mathbf{f}^{l}_{\mathrm{hid},n}$ is 32, whereas the dimension of
$\mathbf{f}^{l}_{\mathrm{reg},n}$ and $\mathbf{f}^{l}_{\mathrm{class},n}$ are
128. The dimension of the input atomic number embedding is 512. The maximum
degree of $l$ is 4 where, of course, the model performance also improves if
the degree is higher especially for low-symmetry crystal structures; however,
increasing the maximum degree of $l$ beyond 4 becomes infeasible on our
computational infrastructure, with other hyperparameters fixed. The number of
transformer-head is 2. We use ADAM optimizer and cosine annealing method for a
scheduler with the initial learning rate of $10^{-4}$.

##  Appendix G Trial models and additional results

Table 5: The MAE and RMSE of predicted elastic properties using the model in Fig. G1(b). The dataset for this model consists of materials whose number of edges is less than 2000. Properties | Average | MAE | RMSE
---|---|---|---
$U_{ij}$ (meV/atom) | 2.652 | 0.661 | 1.258
$C_{ij}$ (GPa) | 42.92 | 10.98 | 17.86
$B_{V}$ (GPa) | 107.07 | 11.58 | 19.10
$G_{V}$ (GPa) | 50.98 | 9.91 | 16.31
$Y$ | 129.62 | 22.81 | 37.24
$\nu$ | 0.401 | 0.035 | 0.081
Figure G1:  Earlier trial models with SE(3)-equivariant GNNs as a fundamental building block. (Top) The building block consists of a strained crystal graph input that is fed into SE(3)-equivariant GNNs (SE(3)-transformer followed by a TFN, similar to that of the building block in StrainTensorNet), whose output is the global max-pooled feature from the TFN. (a) A minimal equivariant GNN model for SED regression. (b) An SED regression network that take both the known degeneracy class vector and the minimal equivariant GNN’s representation as an input. (c) The same network as (a) but the input node attributes also incorporates the degeneracy class vector. Figure G2:  From top to bottom rows are, respectively, the distance matric (of the materials whose numbers of edges are less than 500) of the ground truths, of the model in Fig. G1(a), of the model Fig. G1(c), and of the model model in Fig. G1(b). The last row suggests that the concatenation of the strained crystal graph latent representation and the degeneracy class vector significantly increases the expressiveness of the SED regression network, motivating the self-supervised network architecture in the main text.  Table 6: The MAE of predicted elastic properties using the model in Fig. G1(a), (b), and (c). The dataset for these earlier models are the materials whose number of edges is less than 500. Properties | Average | MAE of Fig. G1(a) | MAE of Fig. G1(b) | MAE of Fig. G1(c)
---|---|---|---|---
$U_{ij}$ (meV/atom) | 2.749 | 1.628 | 0.842 | 0.866
$C_{ij}$ (GPa) | 44.10 | 78.90 | 13.29 | 18.66
$B_{V}$ (GPa) | 108.67 | 67.95 | 15.35 | 17.07
$G_{V}$ (GPa) | 53.92 | 79.94 | 12.46 | 14.63
$Y$ | 136.26 | 86.05 | 29.08 | 31.93
$\nu$ | 0.394 | 1.758 | 0.044 | 0.061
Figure G3:  Statistics of the prediction results (vertical axes) vs of the
ground truths (horizontal axes) of the elastic constant components that are
concentrated near 0 GPa. (a) The predicted $C_{26}$ of a cubic system has the
mean and the standard deviation of 0.69 and 4.45 GPa. (b) The predicted
$C_{56}$ of a cubic system has the the mean and the standard deviation of
-0.43 and 4.85 GPa. (c) The predicted $C_{26}$ of all the crystal systems has
the mean and the standard deviation of 0.99 and 6.08 GPa. (d) The predicted
$C_{56}$ of all the crystal systems has the mean and the standard deviation of
-0.20 and 9.22 GPa, respectively. The ground-truth distributions of $C_{26}$
and $C_{56}$ are centered around 0 GPa with vanishing widths. Although some of
the predicted $C_{26}$ and $C_{56}$ are not exactly 0 GPa (for cubic system),
most of the predicted values are near 0 GPa.  Figure G4:  The bar charts
display the MAE of the predicted elastic properties, i.e., (a) $B_{V}$, (b)
$G_{V}$, (c) $Y$, and (d) $\nu$. The charts are categorized based on the
elements present in the compounds in the test set and colored according to the
average bulk modulus, as indicated in the legend below. The dashed lines
represent the MAE of each property, averaged over the data in the test set. It
can be seen that most elements, whose compounds have high bulk modulus,
exhibit errors of $B_{V}$, $G_{V}$, and $Y$ that exceed the dataset’s MAE.

## References

  * Born and Huang [1988] M. Born and K. Huang,  _Dynamical Theory of Crystal Lattices_ , International series of monographs on physics (Clarendon Press, 1988).
  * Mouhat and Coudert [2014] F. Mouhat and F. m. c.-X. Coudert, Necessary and sufficient elastic stability conditions in various crystal systems, Phys. Rev. B 90, 224104 (2014).
  * MACDONALD [1969] J. R. MACDONALD, Review of some experimental and analytical equations of state, Rev. Mod. Phys. 41, 316 (1969).
  * Brazhkin [2007] V. V. Brazhkin, High-pressure synthesized materials: treasures and hints, High Pressure Research 27, 333 (2007).
  * Yamaguchi  _et al._ [1997] M. Yamaguchi, T. Yagi, T. Azuhata, T. Sota, K. Suzuki, S. Chichibu, and S. Nakamura, Brillouin scattering study of gallium nitride: elastic stiffness constants, Journal of Physics: Condensed Matter 9, 241 (1997).
  * Reichmann  _et al._ [2000] H. J. Reichmann, S. D. Jacobsen, S. J. Mackwell, and C. A. McCammon, Sound wave velocities and elastic constants for magnesiowüstite using gigahertz interferometry, Geophysical Research Letters 27, 799 (2000).
  * Teter [1998] D. M. Teter, Computational alchemy: The search for new superhard materials, MRS Bulletin 23, 22 (1998).
  * Chen  _et al._ [2011] X.-Q. Chen, H. Niu, D. Li, and Y. Li, Modeling hardness of polycrystalline materials and bulk metallic glasses, Intermetallics 19, 1275 (2011).
  * Tian  _et al._ [2012] Y. Tian, B. Xu, and Z. Zhao, Microscopic theory of hardness and design of novel superhard crystals, International Journal of Refractory Metals and Hard Materials 33, 93 (2012).
  * Nielsen and Martin [1985] O. H. Nielsen and R. M. Martin, Stresses in semiconductors: Ab initio calculations on si, ge, and gaas, Phys. Rev. B 32, 3792 (1985).
  * Labéguerie  _et al._ [2010] P. Labéguerie, M. Harb, I. Baraille, and M. Rérat, Structural, electronic, elastic, and piezoelectric properties of $\alpha$-quartz and $mx{\text{o}}_{4}$ ($m=\text{Al}$, ga, fe; $x=\text{P}$, as) isomorph compounds: A dft study, Phys. Rev. B 81, 045107 (2010).
  * Wang  _et al._ [2010] Y. Wang, J. J. Wang, H. Zhang, V. R. Manga, S. L. Shang, L.-Q. Chen, and Z.-K. Liu, A first-principles approach to finite temperature elastic constants, Journal of Physics: Condensed Matter 22, 225404 (2010).
  * Ektarawong  _et al._ [2016] A. Ektarawong, S. I. Simak, and B. Alling, Carbon-rich icosahedral boron carbides beyond ${\mathrm{b}}_{4}\mathrm{C}$ and their thermodynamic stabilities at high temperature and pressure from first principles, Phys. Rev. B 94, 054104 (2016).
  * Stadler  _et al._ [1996] R. Stadler, W. Wolf, R. Podloucky, G. Kresse, J. Furthmüller, and J. Hafner, Ab initio calculations of the cohesive, elastic, and dynamical properties of ${\mathrm{cosi}}_{2}$ by pseudopotential and all-electron techniques, Phys. Rev. B 54, 1729 (1996).
  * Le Page and Saxe [2002] Y. Le Page and P. Saxe, Symmetry-general least-squares extraction of elastic data for strained materials from ab initio calculations of stress, Phys. Rev. B 65, 104104 (2002).
  * Sarasamak  _et al._ [2010] K. Sarasamak, S. Limpijumnong, and W. R. L. Lambrecht, Pressure-dependent elastic constants and sound velocities of wurtzite sic, gan, inn, zno, and cdse, and their relation to the high-pressure phase transition: A first-principles study, Phys. Rev. B 82, 035201 (2010).
  * Chernozatonskii  _et al._ [2011] L. A. Chernozatonskii, P. B. Sorokin, A. A. Kuzubov, B. P. Sorokin, A. G. Kvashnin, D. G. Kvashnin, P. V. Avramov, and B. I. Yakobson, Influence of size effect on the electronic and elastic properties of diamond films with nanometer thickness, The Journal of Physical Chemistry C 115, 132 (2011), https://doi.org/10.1021/jp1080687 .
  * Pakornchote  _et al._ [2021] T. Pakornchote, A. Ektarawong, U. Pinsook, and T. Bovornratanaraks, Modifying electronic and elastic properties of 2-dimensional [110] diamond by nitrogen substitution, C 7, 10.3390/c7010008 (2021).
  * Caro  _et al._ [2012] M. A. Caro, S. Schulz, and E. P. O’Reilly, Comparison of stress and total energy methods for calculation of elastic properties of semiconductors, Journal of Physics: Condensed Matter 25, 025803 (2012).
  * Kresse and Furthmüller [1996] G. Kresse and J. Furthmüller, Efficiency of ab-initio total energy calculations for metals and semiconductors using a plane-wave basis set, Computational Materials Science 6, 15 (1996).
  * Giannozzi  _et al._ [2009] P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, G. L. Chiarotti, M. Cococcioni, I. Dabo, A. D. Corso, S. de Gironcoli, S. Fabris, G. Fratesi, R. Gebauer, U. Gerstmann, C. Gougoussis, A. Kokalj, M. Lazzeri, L. Martin-Samos, N. Marzari, F. Mauri, R. Mazzarello, S. Paolini, A. Pasquarello, L. Paulatto, C. Sbraccia, S. Scandolo, G. Sclauzero, A. P. Seitsonen, A. Smogunov, P. Umari, and R. M. Wentzcovitch, Quantum espresso: a modular and open-source software project for quantum simulations of materials, Journal of Physics: Condensed Matter 21, 395502 (2009).
  * Clark  _et al._ [2005] S. J. Clark, M. D. Segall, C. J. Pickard, P. J. Hasnip, M. I. J. Probert, K. Refson, and M. C. Payne, First principles methods using castep, Zeitschrift für Kristallographie - Crystalline Materials 220, 567 (2005).
  * de Jong  _et al._ [2015] M. de Jong, W. Chen, T. Angsten, A. Jain, R. Notestine, A. Gamst, M. Sluiter, C. Krishna Ande, S. van der Zwaag, J. J. Plata, C. Toher, S. Curtarolo, G. Ceder, K. A. Persson, and M. Asta, Charting the complete elastic properties of inorganic crystalline compounds, Scientific Data 2, 150009 (2015).
  * Jain  _et al._ [2013] A. Jain, S. P. Ong, G. Hautier, W. Chen, W. D. Richards, S. Dacek, S. Cholia, D. Gunter, D. Skinner, G. Ceder, and K. a. Persson, The Materials Project: A materials genome approach to accelerating materials innovation, APL Materials 1, 011002 (2013).
  * Mazhnik and Oganov [2020] E. Mazhnik and A. R. Oganov, Application of machine learning methods for predicting new superhard materials, Journal of Applied Physics 128, 075102 (2020).
  * Zhao  _et al._ [2020] Y. Zhao, K. Yuan, Y. Liu, S.-Y. Louis, M. Hu, and J. Hu, Predicting elastic properties of materials from electronic charge density using 3d deep convolutional neural networks, The Journal of Physical Chemistry C 124, 17262 (2020), https://doi.org/10.1021/acs.jpcc.0c02348 .
  * Levämäki  _et al._ [2022] H. Levämäki, F. Tasnádi, D. G. Sangiovanni, L. J. S. Johnson, R. Armiento, and I. A. Abrikosov, Predicting elastic properties of hard-coating alloys using ab-initio and machine learning methods, npj Computational Materials 8, 17 (2022).
  * Linton and Aidhy [2023] N. Linton and D. S. Aidhy, A machine learning framework for elastic constants predictions in multi-principal element alloys, APL Machine Learning 1, 10.1063/5.0129928 (2023), 016109, https://pubs.aip.org/aip/aml/article-pdf/doi/10.1063/5.0129928/16772730/016109_1_online.pdf .
  * Avery  _et al._ [2019] P. Avery, X. Wang, C. Oses, E. Gossett, D. M. Proserpio, C. Toher, S. Curtarolo, and E. Zurek, Predicting superhard materials via a machine learning informed evolutionary structure search, npj Computational Materials 5, 89 (2019).
  * Chen  _et al._ [2021] W.-C. Chen, J. N. Schmidt, D. Yan, Y. K. Vohra, and C.-C. Chen, Machine learning and evolutionary prediction of superhard b-c-n compounds, npj Computational Materials 7, 114 (2021).
  * Kearnes  _et al._ [2016] S. Kearnes, K. McCloskey, M. Berndl, V. Pande, and P. Riley, Molecular graph convolutions: moving beyond fingerprints, Journal of Computer-Aided Molecular Design 30, 595 (2016).
  * Gilmer  _et al._ [2017] J. Gilmer, S. S. Schoenholz, P. F. Riley, O. Vinyals, and G. E. Dahl, Neural message passing for quantum chemistry (2017).
  * Schütt  _et al._ [2018] K. T. Schütt, H. E. Sauceda, P.-J. Kindermans, A. Tkatchenko, and K.-R. Müller, Schnet – a deep learning architecture for molecules and materials, The Journal of Chemical Physics 148, 241722 (2018), https://doi.org/10.1063/1.5019779 .
  * Dym and Maron [2021] N. Dym and H. Maron, On the universality of rotation equivariant point cloud networks, in  _International Conference on Learning Representations_ (2021).
  * Note [1] Tensor Field Network (TFN), L1Net, GemNet, and EGNN, have exploited this computational architecture to make better data-driven predictions on the properties of chemical compounds [41, 52, 50, 51]. The SE(3)-Transformers model based on TFN with an attention mechanism also demonstrates a good predictive performance on the standard QM9 dataset [42]. In addition, translation-rotation equivariance concept has also been incorporated into generative modeling, in an attempt to search for novel molecules and crystal structures with desirable properties [54, 55].
  * Slaughter [2002] W. S. Slaughter,  _The Linearized Theory of Elasticity_ (Birkhäuser Boston, 2002).
  * Xie and Grossman [2018] T. Xie and J. C. Grossman, Crystal graph convolutional neural networks for an accurate and interpretable prediction of material properties, Phys. Rev. Lett. 120, 145301 (2018).
  * Zimmermann  _et al._ [2017] N. E. R. Zimmermann, M. K. Horton, A. Jain, and M. Haranczyk, Assessing local structure motifs using order parameters for motif recognition, interstitial identification, and diffusion path characterization, Frontiers in Materials 4, 10.3389/fmats.2017.00034 (2017).
  * Pan  _et al._ [2021] H. Pan, A. M. Ganose, M. Horton, M. Aykol, K. A. Persson, N. E. R. Zimmermann, and A. Jain, Benchmarking coordination number prediction algorithms on inorganic crystal structures, Inorganic Chemistry 60, 1590 (2021), pMID: 33417450.
  * Bronstein  _et al._ [2021] M. M. Bronstein, J. Bruna, T. Cohen, and P. Velickovic, Geometric deep learning: Grids, groups, graphs, geodesics, and gauges, CoRR abs/2104.13478 (2021), 2104.13478 .
  * Thomas  _et al._ [2018] N. Thomas, T. Smidt, S. Kearnes, L. Yang, L. Li, K. Kohlhoff, and P. Riley, Tensor field networks: Rotation- and translation-equivariant neural networks for 3d point clouds (2018).
  * Fuchs  _et al._ [2020] F. Fuchs, D. Worrall, V. Fischer, and M. Welling, Se(3)-transformers: 3d roto-translation equivariant attention networks, in  _Advances in Neural Information Processing Systems_, Vol. 33, edited by H. Larochelle, M. Ranzato, R. Hadsell, M. Balcan, and H. Lin (Curran Associates, Inc., 2020) pp. 1970–1981.
  * Note [2] Due to the symmetry of SED, we identify only an upper triangular element.
  * Fung  _et al._ [2021] V. Fung, J. Zhang, E. Juarez, and B. G. Sumpter, Benchmarking graph neural networks for materials chemistry, npj Computational Materials 7, 84 (2021).
  * Schwerdtfeger and Nagle [2019] P. Schwerdtfeger and J. K. Nagle, 2018 table of static dipole polarizabilities of the neutral elements in the periodic table, Molecular Physics 117, 1200 (2019).
  * Lafon [2004] S. Lafon,  _Diffusion maps and geometric harmonics_ , Ph.D. thesis, Yale University (2004).
  * Coifman and Lafon [2006] R. R. Coifman and S. Lafon, Diffusion maps, Appl. Comput. Harmon. Anal. 21, 5 (2006).
  * Cousins [1978] C. S. G. Cousins, Inner elasticity, Journal of Physics C: Solid State Physics 11, 4867 (1978).
  * Battaglia  _et al._ [2018] P. W. Battaglia, J. B. Hamrick, V. Bapst, A. Sanchez-Gonzalez, V. Zambaldi, M. Malinowski, A. Tacchetti, D. Raposo, A. Santoro, R. Faulkner, C. Gulcehre, F. Song, A. Ballard, J. Gilmer, G. Dahl, A. Vaswani, K. Allen, C. Nash, V. Langston, C. Dyer, N. Heess, D. Wierstra, P. Kohli, M. Botvinick, O. Vinyals, Y. Li, and R. Pascanu, Relational inductive biases, deep learning, and graph networks (2018).
  * Klicpera  _et al._ [2021] J. Klicpera, F. Becker, and S. Günnemann, Gemnet: Universal directional graph neural networks for molecules, in  _Advances in Neural Information Processing Systems_, edited by A. Beygelzimer, Y. Dauphin, P. Liang, and J. W. Vaughan (2021).
  * Satorras  _et al._ [2021] V. G. Satorras, E. Hoogeboom, and M. Welling, E(n) equivariant graph neural networks, in  _Proceedings of the 38th International Conference on Machine Learning_, Proceedings of Machine Learning Research, Vol. 139, edited by M. Meila and T. Zhang (PMLR, 2021) pp. 9323–9332.
  * Miller  _et al._ [2020] B. K. Miller, M. Geiger, T. E. Smidt, and F. Noé, Relevance of rotationally equivariant convolutions for predicting molecular properties, CoRR abs/2008.08461 (2020), 2008.08461 .
  * Ong  _et al._ [2013] S. P. Ong, W. D. Richards, A. Jain, G. Hautier, M. Kocher, S. Cholia, D. Gunter, V. L. Chevrier, K. A. Persson, and G. Ceder, Python materials genomics (pymatgen): A robust, open-source python library for materials analysis, Computational Materials Science 68, 314 (2013).
  * Xu  _et al._ [2022] M. Xu, L. Yu, Y. Song, C. Shi, S. Ermon, and J. Tang, Geodiff: A geometric diffusion model for molecular conformation generation, in  _International Conference on Learning Representations_ (2022).
  * Xie  _et al._ [2022] T. Xie, X. Fu, O.-E. Ganea, R. Barzilay, and T. S. Jaakkola, Crystal diffusion variational autoencoder for periodic material generation, in  _International Conference on Learning Representations_ (2022).