Unified Model for Crystalline Material Generation

1 Introduction

One of the fundamental challenges of materials science is to obtain new thermodynamically stable materials that satisfy certain desirable properties. Recently, considerable attention has been devoted to the generation of crystalline (crystal) materials due to their wide range of usage in our modern society, e.g. metal alloys or semiconductors. As an example of application, we are interested in generating new crystal materials for developing new solar panels with a band gap enabling hydrolyse allowing hydrogen production from H2O. Another application consists in discovering metal hydride for hydrogen storage applications which helps solve problems related to clean energy production and storage, which is one of the major challenges facing our society.

To discover new materials with desirable properties, highthroughput screening based on Machine Learning (ML) models is the most widely employed technique. Despite several success stories, in particular, for organic molecule generations, more progress is required for crystal generation. Crystals are three-dimensional periodic structures composed of a wide variety of chemical bonds and atoms which are often represented as a parallelepiped tiling, a.k.a crystal lattice or unit cell. The periodic structure of crystals makes it difficult to process when training generative models. Recently, graphbased representation models based on geometrically equivariant ML techniques have led to substantial performance increases across a wide range of supervised tasks such as property predictions and classification [Sch¨utt et al., 2017; Jørgensen et al., 2018; Gasteiger et al., 2020b; Gasteiger et al., 2020a; Chen et al., 2019; Choudhary and DeCost, 2021; Klicpera et al., 2021] and material generation [Satorras et al., 2021; Xie et al., 2022; Long et al., 2021; Ekstr¨om Kelvinius et al., 2022; Gibson et al., 2022] for crystals. In [Klipfel et al., 2023], a general framework that formulates an Equivariant Message Passing Neural Network (EMPNN) on periodic structures acting on crystal lattice without any label from the interaction forces and stress tensors has been proposed. This model enforces a structuring bias for crystals using group actions with respect to the Euclidean Euc(3) and SL3(Z) groups by means of an equivariance property of Message Passing Neural Network (MPNN) layers. Results on the denoising tasks show the capability of the EMPNN model to perform arbitrary crystal lattice deformation by improving the total energy of the structure (lowering the total energy of a structure leads to obtaining thermodynamic stability). In a similar vein, [Xie et al., 2022] proposed a model adapted from GemNet [Klicpera et al., 2021], which is an equivariant GNN model to rotation and translation groups for organic chemistry, to periodic structures. Contrarily to the model from [Klipfel et al., 2023], the model from [Xie et al., 2022] is capable to act on the atomic positions, but not capable of modifying the geometry (the shape of crystal lattice).

In this paper, we propose and analyse two unified models that process both the geometry and atomic positions of crystal materials at the same time. This is important for generation since the model will be able to modify the whole geometry during the inference and hence reach an optimum one. We first propose an Equivariant Graph Neural Network (EGNN) model as an extension of the EMPNN model that adds action on the atomic positions in addition to the shape of the lattice. While outperforming the EMPNN model, the architecture of EGNN remains suboptimal as it does not allow dealing with chemical structures. Namely, EGNN relies on a very basic embedding to represent angles whereas the message-passing schema only relies on the edges, but not the whole geometric information of the triplets. To this end, we propose GemsNet, a geometry architecture inspired by GemNet [Klicpera et al., 2021] that is capable to modify both atomic positions and lattice geometry at the same time within the same Equivariant GNN by exploiting all the geometric information of the triplets. To the best of our knowledge, GemsNet and EGNN are the first unified models that act simultaneously on crystal lattice geometry and properties. The quantitative and qualitative analyses of our models show the effectiveness of such a unified approach.

2 RelatedWorks

Within the area of discovering new organic or crystal materials, we can distinguish different classes of related works according to representing technique used to encode molecules, e.g. Fingerprint representation [Kim et al., 2020b; Nouira et al., 2018] or voxel representation [Court et al., 2020; Kim et al., 2020a; Noh et al., 2019; Long et al., 2021]. From organic chemistry to material science, models based on graph representation are by far the most used in a variety of tasks such as classification [Sch¨utt et al., 2017; Jørgensen et al., 2018; Gasteiger et al., 2020b; Gasteiger et al., 2020a; Chen et al., 2019; Choudhary and DeCost, 2021; Klicpera et al., 2021] and generation [Satorras et al., 2021; Xie et al., 2022; Long et al., 2021; Ekstr¨om Kelvinius et al., 2022; Gibson et al., 2022].

The graph-based representation of materials has the advantage of representing the periodicity of structures as well as the local environment of each atom. The graph-based representation allows using Graphical Neural Networks (GNN) architectures to manipulate materials. GNNs are capable to process sparse data and take into account invariance or equivariance to many action groups that act on molecules to deform them. Existing works (e.g. [Klicpera et al., 2021]) are equivariant to SO(3) due to a spherical basis which allows predicting lattice properties and performing simulations. However, this is not enough to deform crystal lattices when the shape of the lattice is unknown in advance. In [Klipfel et al., 2023], group actions with respect to the Euclidean group Euc(3) and SL3(Z) group are incorporated by the equivariance property of MPNN layers to act on crystal lattices. This has the advantage of modifying the material regardless of its representation, i.e. the modification applied by the MPNN layer is independent of the orientation or the way the material is paved. Inspired by GemNet [Klicpera et al., 2021] an equivariant GNN model to rotation SO(3) and translation groups for organic chemistry, [Xie et al., 2022] proposed a model for periodic structures that contrarily to [Klipfel et al., 2023] allows to act on atomic positions, but not capable of modifying the shape of crystal lattice. While existing models have shown interesting results on denoising tasks, i.e. performing arbitrary crystal lattice deformation by improving the total energy of the structure, acting only on atomic positions or the crystal geometry is not enough to perform generation tasks. Our models perform simultaneously on both atomic positions and the crystal geometry within the same equivariant architecture.

Recently, several molecular dynamics models for generations have been proposed by approximating DFT simulation with GNN [Pickard and Needs, 2011; Ekstr¨om Kelvinius et al., 2022; Gibson et al., 2022; Cheon et al., 2020]. These methods mainly learn interaction forces and stress tensors to lower the total energy of a structure with methods analogue to DFT calculation. While relying on self-simulations to gather data is interesting, it only concerns a small set of specific structures. To discover new materials, however, a lot of additional information about interaction forces is required, which is not always available. In addition, we cannot rely on randomly generated structures, as they lead in general to unstable structures. Our models do not require any additional labels.

3 Background

A crystalline (crystal) structure can be defined as a cloud of atoms and a repetition pattern. The repeated pattern represents periodicity and is often described as a parallelepiped called a cell. The repetition of the cell is called a lattice. The periodic structure is obtained with the tiling of the space by the crystal cell. In particular, a given atom inside of the cell is repeated in multiple positions because of the tiling in space. As a consequence, the local environment of an atom can overlap with adjacent repetitions. We follow [Klipfel et al., 2023] to represent a crystal material as a graph. It is an oriented graph where each edge is represented by triplets containing the index of the source node, the index of the destination node and the relative cell coordinate of the destination node. Notice that the definition of the graph provided in [Klipfel et al., 2023] generalizes most of the existing graph definitions [Jørgensen et al., 2018; Chen et al., 2019; Satorras et al., 2021].

Definition 1. The representation space of featured materials MF is the disjoint union ∐n∈N MFn where:

Chemical materials are represented in M = MN, with atomic numbers as feature sequence z.

GLd(R) defines the shape of the lattice, i.e. the periodicity, where d stands for the dimension of the material. F is the feature space that encodes chemical information such as atomic number or charge. Our aim is to consider networks capable of deforming the geometry of a structure in order to minimize the total energy and hence obtain a stable structure. To deform the geometry of a structure, we consider a neural network that predicts a collection {y, y1, · · · , yn}, where y ∈ GL3(R) acts on the material lattice ρ resulting in the updated lattice ρ′ where yi ∈ R3 acts on the atomic positions as follows:

The atomic positions are brought back into the crystal lattice by truncation. The equivariance property ensures that the actions performed on a given material are independent of any choice of representation. Namely, a given material has multiple equivalent associated representations, each of which

with depends on the ordering of the atoms, the orientation of the crystal and the tailing of the space. Our goal is then to process materials independently from their representations. To this end, we use group actions to determine the relationship between different representations. Group actions are the mathematical tools that enable one to pass from one choice of representation to another. The group Euc(3) of translations and rotations keeps track of the orientation and the position, while SL3(Z) keeps track of the choice of a cell. This last group is needed because crystals are infinite collections of atoms while computers are only capable to fit a finite description. A natural approach to obtain such a finite description is to only work with a finite piece of a given crystal that is big enough to contain all the information of the structure. However, we need to ensure that the obtained results are independent of the choice of the piece of material we are working with. Since one can go from a description based on one choice of a piece of crystal to another using the SL3(Z) action, results produced by models equivariant to this action will not depend on a particular choice of a representative piece of material. More formally, equivariance is defined as follows:

For more details about how an element g of the group acts on the structure M, we refer to Proposition 1 and 2 in [Klipfel et al., 2023].

We now recall the EMPNN model proposed in [Klipfel et al., 2023] which performs arbitrary deformation by acting on relative atomic distances and angles. In this model, the spatial equivariance is enforced by an MPNN layer. Intuitively the EMPNNmodel takes advantage of the local invariance (input quantities are themselves invariant: distance, angle, etc.) and equivariance of the graph-based representation of materials to define equivariant actions on crystal lattices.

Definition 2. Let M = (ρ, x, z) in MFn be a material and ci > 0 for 1 ≤ i ≤ n denotes cutoff distances. We define a directed 2-graph Γ = ΓM,c by the graded components:

• Γ0 = {1, . . . , n}

• Γ2 = {(γ, γ ′ ) ∈ Γ1 × Γ1 ∣∣ src(γ) = src(γ ′ )}

The following definition introduces the notations needed for our models.

Definition 3. Let consider M = (ρ, x, z) ∈ MF and Γ = ΓM,c, we introduce the following notations:

• eijτ = (xj - xi + τ) for edge vector in lattice coordinates,
• vijτ = ρeijτ for the edge vector in physical space,
• rijτ = ||vijτ || for the physical edge length,
• θijτkτ ′ as the unoriented angle between vijτ and vikτ ′

Let us also write eγ, vγ, rγ, θγγ ′ for the same quantities when we do not need to make vertices explicit. Note that rγ and θγγ ′ are natural Euclid invariants.

4 UnifiedModels for Crystal Material

In this section, we propose two unified models that are capable to act on both atomic positions and crystal lattice. The first model called EGNN, which starts from the EMPNN model proposed in [Klipfel et al., 2023] which already acts on crystal lattices. The improvement consists in carefully adding a new component that modifies atomic positions without losing the equivariance property. The other way around, we introduce a second model called GemsNet which is based on GemNet [Klicpera et al., 2021; Xie et al., 2022] which acts on atomic positions. In GemsNet, we add components allowing to deform the crystal lattices within the same equivariant GNN architecture. In the following, we describe these two models.

4.1 Equivariant Graph Neural Network

The GNN model uses a message-passing schema with the formation of messages from the chemical and geometric information of the edges (Equation 2) and updates the chemical features of the nodes with a GRU (Equation 3). To generate a deformation of a crystal lattice, we first predict the weighting of the vector fields generated from the edges and the triplets (Equation 4). Then, we aggregate the weighting of all the vector fields to deform the lattice (Equation 5). Figure 1 provides an overview of our EGNN model.

The EGNN model modifies atomic positions by updating them according to the following equations:

With wlijτ are the weights at the layer l for the set of edges N(i) connected to the atom i and eijτ is the vector corresponding to the edge (i, j, τ). Intuitively, the GNN architecture predicts the atom displacement along the edges. We compute an average of these predictions. The atoms displacement is computed as the sum of the displacement predicted by all the layers of the network. We then obtain the weightings wlijτ from a prediction on the chemical embedding of the atoms and the embedding of the interatomic distances according to the setup given in [Sch¨utt et al., 2017], which is expressed by Equation 9.

We use c = 10 ˚A and µk = k0.1 ˚A with 0 ≤ µk < c. Notice that the size of the hidden layer (from Equation 8) is 256.

4.2 GemsNet

For GemsNet model, we start from the work of [Klicpera et al., 2021] on organic chemistry and its adaptation to crystal materials proposed by [Xie et al., 2022]. The overview of our proposed architecture is shown in Figure 2. More precisely, we add a module to predict the actions on the lattice at the same time with atomic positions based only on triplets from Γ2. To perform crystal lattice deformation, we simply add components to the "output block" of GemNet to be able to predict deformations from the triplets based on the chemical embeddings. Recall that the "output block" is responsible for making predictions about the structures. Namely, it predicts the energy (when it is needed) and the interaction forces, which are used to generate a trajectory and to update the atomic positions). We add a feed-forward network to weight vector fields with wlγγ ′ ∈ RNfields and obtain actions on the lattice.

with φρθ is an MLP and || is concatenation. θγγ ′ is the concatenation of geometric features of GemNet triplet, so one RBF per edge and one CBF per angle such that: embγγ ′ = eRBF(rγ )||eRBF(rγ ′ )||eCBF(rγ , rγγ ′ ).

With the final representations of RBF and CBF are given as follows:

Where u(r) is used to remove the interactions having distances greater than the cutoff distance. Following the recommendation given by [Gasteiger et al., 2020a], we chose p = 6.

The hidden layers (Equation 10) are of size 64. Notice that mγγ ′ is invariant. As such, we can express ρ′ as follows:

with λγγ ′ ∈ RNfields×3×3 are the vector fields that define how the GNN acts on the crystal lattice. We use gradients of the geometric invariant such as:

5 Experiments

In this section, we experimentally compare the two proposed models with a number of baseline methods with respect to reconstruction and denoising tasks 1 . The denoising task (or geometry optimization task) aims to evaluate the capability of our model to perform arbitrary crystal lattice deformation by improving the total energy of the structure. To perform denoising, we apply a small random deformation to a stable structure leading to a less stable one with a high energy level (as energy increases in all directions locally). We can therefore generate pairs of stable and less stable structures that can be used to teach our models how to deform the less stable structure to reach stability. Making a structure more stable comes down to lowering the energy level (thermodynamic stability). The reconstruction task (or lattice reconstruction task) aims to restore a given crystal lattice from a non-noisy atomic coordinate following the same settings reported in [Klipfel et al., 2023]. More specifically, the reconstruction aims to build a crystal lattice from scratch. We start from the point cloud as if it was in a cubic lattice of one angstrom on a side. The main working hypothesis is that there is a single stable lattice which corresponds to the starting atomic positions. Notice that it is impossible to perform such a task with chemical simulation techniques such as DFT. Interestingly enough, with reconstruction capabilities, a model can be easily integrated into VAE or GANs architectures to perform generation.

1 Code and data are available at https://github.com/aklipf/ GemsNet

Datasets. We considered datasets of stable crystals where each structure is in local minima of formation energy. First, we use Perov-5 [Castelli et al., 2012a; Castelli et al., 2012b] which contains perovskite (cubic) structures with highly uniform shapes, but different chemical compositions between structures. Second, we consider Carbon-24 [Pickard, 2020] which is composed of carbon atoms having a large variety of shapes. Finally, we conduct experiments on Mp-20 [Jain et al., 2013] which is a subset of the materials project proposed in [Xie et al., 2022] that has a large sample of shapes and chemical compositions. It is the most representative of ordinary structures. We used the same training, validation and test splits as reported in [Xie et al., 2022].

Random deformation. As mentioned above, the random deformation aims to introduce noisy structures needed to assess the capability of our models to reach stability. To this end, we sample noises and add them to the atomic positions. The lattice parameters are removed and replaced by a cubic lattice of dimension 1. The following equation gives precision about how to generate the noise.

Loss function. As loss function, we aim to determine the optimal trajectory to denoise a structure. Indeed, as the crystals are periodic, several trajectories are possible to obtain the same result. The optimal trajectory is then defined as the trajectory that minimizes the distance we are moving through.

More formally, we have:

We denote by p : GL3(R) → R6 an application that associates the lattice parameters (a, b, c, α, β, γ) to the lattice ρ. We recall that yi is defined in 1. We can then define a reconstruction loss with L1 norm as follows:

Evaluation metrics. As evaluation metrics, we rely on angle and length metrics which are respectively Mean Absolute Error (MAE) metrics on the lattice parameters (α, β, γ) and (a, b, c). The positioning metric is the distance between the optimal position and the actual position. We take into account the periodicity of the material for the atomic positions. If the atom is close to an edge and the optimal position is close to the opposite edge, then the distance to be covered is shorter than the distance that crosses the entire lattice. The distance metric is therefore the shortest distance between the position of an atom and one of its possible optimal positions. N is the number of structures and M the number of atoms.

Training. We train our two neural network models for 128 epochs. We use hard training from 45 min for the smallest dataset such as Perov-5 to 5 hours for the biggest datasets such as Mp-20. Our models are trained on an NVIDIA RTX A6000 GPU. The different hyperparameters are given in Table 3.

Lattice reconstruction task. Following the same setting proposed in [Klipfel et al., 2023], we use Equation 18 without applying noise to xi, i.e. ϵ := 0. We only use the loss term on the Equation 20 of lattice parameters. Table 2 reports the obtained results. We can see that GemsNet has better performances than existing models on Carbon-24 and Mp-20 datasets. In particular, there is an improvement in the length parameters. However, we obtain identical or slightly worse results for the angular parameters. The results are more variable on Mp-20 than on Carbon-24. This is expected since Carbon-24 is a database composed only of pure carbon structures under pressure, while Mp-20 contains many different chemical species and forms. It is therefore possible that the training is less stable on Mp-20 because there is more variability from one structure to another. Overall, GemsNet achieves good performance. This is interesting because GemsNet has more than 2 times fewer parameters and it is much faster in terms of inference compared to existing models.

Geometry Optimization Task

In this experiment, we use our EGNN as a baseline since to the best of our knowledge, there is no equivariant GNN that is capable of acting on both atomic positions and lattice shape at the same time. EGNN can be seen as a hybridisation between a model which only allows acting on the geometry (lattice shape) in an equivariant way (e.g. [Klipfel et al., 2023]) and a model which only allows acting on the atomic positions (e.g. [Satorras et al., 2021]). As can be seen in Table 1, our two models achieve good performances and are able to find the original geometry of the structures in the database. Notice that GemsNet outperforms EGNN for both atomic positions and lattice parameters.

For GemsNet, the reconstruction error on the Perov-5 database is 61% lower for angles and 20% lower for the lengths of the sides of the lattice. The positioning error is greater for the atomic positions, but the positioning metrics remain low in both cases. Recall that Perov-5 is a perovskite base where structures have very similar geometry. As such, it is the least difficult and least interesting dataset. For Carbon24 the structures are all in pure carbon, but there is more variety in the geometry. GemsNet allows a reconstruction error reduction of 28% for the lattice parameters and a 20% reduction in the reconstruction error of the atomic positions. Finally, Mp-20 is the most representative and varied database. On this database, GemsNet allows an average reduction of 10% of the error of reconstruction of the angular parameters of the lattice; a reduction of the error of almost 40% of the parameters of lengths and a reduction of 13% of atomic position error. Finally, we can see that the accuracy of the reconstruction is still very good, especially for the databases which are the most natural, namely with Mp-20 which is a subset of Materials Project.

Qualitative Analysis

Figure 3 depicts some examples of structures fromMaterials Project Mp-20 that have been denoised by GemsNet. We can see that the initial structures are too noisy and it is not always easy to reach the stability. One can notice that the majority of the structures are very well reconstructed by GemsNet. Indeed, one can sometimes notice minimal differences in the alignments of the atoms. We can also notice that the shape of the lattices is very close. It is important to note that even crystal structures which are composed of 4 elements are generally well reconstructed whereas many generative models in the literature are limited to 2 or 3 different elements.

However, there are some structures which are visibly poorly reconstructed as we can see in the top right corner of Figure 3. Several elements are remarkable in these structures. The atomic positions look realistic compared to the stable structures, but it is the lattice that makes the shape of the crystals different. Moreover, the reconstruction error is generally big. Errors produced by GemsNet are therefore rare but they are much more important when they occur. The reconstruction error could therefore come in large part from crystals which are outliers. Indeed, if we take a closer look at the badly reconstructed structure at the top right of Figure 3, we can notice that this structure is composed of 4 elements. It is therefore a difficult structure to study. Generative models generally struggle to deal with ternary structures (composed of 3 elements) even when they contain a light element such as oxygen. Our structure contains 4 metals including a heavy atom of potentially radioactive uranium, which is uncommon compared to the other elements in the Mp-20 dataset. We can see that GemsNet manages to solve the denoising task when there is uranium in some cases. Moreover, GemsNet handles structures with a few different chemical species, e.g. the binary structure at the bottom left of Figure 3. In general, crystals with heavy atoms are notoriously more difficult to process but GemsNet manages to process some of them. Overall, GemsNet handles a large part of the dataset structures very well but fails on some outliers. These failures can be nuanced because some structures are really more difficult to analyse, even with ab-initio simulation.

6 Conclusion

We have proposed two unified equivariant models which act simultaneously on the lattice of a crystal and atomic positions. This is essential for generating new crystal materials since our models are capable to modify the whole geometry during inference and reach the optimum state. While our models provide good performances, GemsNet allows a high level of accuracy for the predicted lattice parameters and for the generated atomic positions. In addition, GemsNet is faster in inference time and contains more than 2 times fewer parameters than existing models. This equivariant architecture opens the way for other applications such as the generation of stable crystalline structures which is a challenging problem.

Acknowledgments

This work has been supported by ANR-22-CE23-0002 ERIANA, ANR-20-THIA-0004 and by HPC resources from GENCI-IDRIS (Grant 2022-[AD011013338]).