Arsen Sultanov
email: arsen.sultanov@cnrs.fr

Jean-Claude Crivello
email: jean-claude.crivello@cnrs.fr

Tabea Rebafka
email: tabea.rebafka@sorbonne-universite.fr

Nataliya Sokolovska
email: nataliya.sokolovska@sorbonne-universite.fr

Data-Driven Score-Based Models for Generating Stable Structures with Adaptive Crystal Cells

The discovery of new functional and stable materials is a big challenge due to its complexity. This work aims at the generation of new crystal structures with desired properties, such as chemical stability and specied chemical composition, by using machine learning generative models. Compared to the generation of molecules, crystal structures pose new diculties arising from the periodic nature of the crystal and from the specic symmetry constraints related to the space group. In this work, score-based 1 probabilistic models based on annealed Langevin dynamics, which have shown excellent performance in various applications, are adapted to the task of crystal generation. The novelty of the presented approach resides in the fact that the lattice of the crystal cell is not xed. During the training of the model, the lattice is learned from the available data, whereas during the sampling of a new chemical structure, two denoising processes are used in parallel to generate the lattice along the generation of the atomic positions. A multigraph crystal representation is introduced that respects symmetry constraints, yielding computational advantages and a better quality of the sampled structures. We show that our model is capable of generating new candidate structures in any chosen chemical system and crystal group without any additional training. To illustrate the functionality of the proposed method, a comparison of our model to other recent generative models, based on descriptor-based metrics, is provided.

2

Introduction

The discovery of new materials with desired properties is challenging. We are particularly interested in the search for novel stable crystal structures. Discovering new stable crystal structures is of paramount signicance, since there is an acute need to nd new materials with certain characteristics. From the practical point of view, it is dicult to discover a new material with given properties. Traditionally, new materials are found through laborious experimental research, where human experts apply their knowledge, intuition, and understanding of mechanisms together with some simple heuristics to discover new chemical compounds. In recent decades with the advent of powerful computers and the density functional theory (DFT), it has become possible to estimate the stability of hypothetical crystal structures using high-throughput numerical simulations. This, in turn, allows new materials to be discovered through extensive brute-force screening. However, this approach is still far from being ecient, as the computation is extremely expensive in terms of CPU resources.

Machine learning oers a completely dierent approach. In essence, machine learning can be thought of as automated learning from data. Methods can be designed to make predictions or decisions by learning structures and discovering patterns in existing data. More recent models are trained to generate new instances, with similar properties to those in the training data. Applied in the context of materials discovery, generative machine learning methods can provide fast and viable solutions to problems that are otherwise expensive and dicult to solve.

Background

Given the huge search space in which new materials could be found, the application of machine learning methods can be of great help. However, this type of approach requires a large amount of data to explore the search space. In the case of the crystal structures, the available data is very limited. Another major problem is the representation of chemical crystal structures. Careful engineering of feature sets and ecient data representation are crucial for building scalable and eective machine-learning models.

Genetic algorithms are among the oldest statistical methods used in materials science.

They have been considered for searching the complex chemical space, but their application is computationally too expensive for systems with more than two dierent elements. [START_REF] Le | Discovery and Optimization of Materials Using Evolutionary Approaches[END_REF][START_REF] Henson | Designing Algorithms To Aid Discovery by Chemical Robots[END_REF] Recent advances in machine learning methods propose surrogate algorithms for materials science applications that are cost eective compared to genetic algorithms. This research direction seems to be fruitful. [START_REF] Xie | Machine learning assisted synthesis of metal-organic nanocapsules[END_REF] Some promising results have been reported on the prediction of chemical and physical properties. [START_REF] Ward | A general-purpose machine learning framework for predicting properties of inorganic materials[END_REF] The relationships between structure and properties have also been investigated. [START_REF] Dong | Bandgap prediction by deep learning in congurationally hybridized graphene and boron nitride[END_REF][START_REF] Oliynyk | Disentangling Structural Confusion through Machine Learning: Structure Prediction and Polymorphism of Equiatomic Ternary Phases ABC[END_REF] Note that there is no single optimal machine learning method for all materials science applications, so Support Vector Machines (SVM) were used to predict reaction outcomes of syntheses of templated vanadium selenites. [START_REF] Raccuglia | Machine-learning-assisted materials discovery using failed experiments[END_REF] Random Forests (RF) were successfully applied to predict catalyst performance in C-N coupling reactions. [START_REF] Ahneman | Predicting reaction performance in CN cross-coupling using machine learning[END_REF] Deep learning has been used to predict the reactivity of reaction combinations with high accuracy. [START_REF] Granda | Controlling an organic synthesis robot with machine learning to search for new reactivity[END_REF] A number of papers report promising results on machine learning for chemical synthesis. 1012 The XGBoost model was tested on a trial-error synthesis process of single crystals to predict their crystallisation propensity. [START_REF] Xie | Machine learning assisted synthesis of metal-organic nanocapsules[END_REF] Note, however, that it has also been reported that data mining methods for chemical applications, e.g. for the synthesis of inorganic-organic hybrid materials, are still rather limited. [START_REF] Raccuglia | Machine-learning-assisted materials discovery using failed experiments[END_REF][START_REF] Moosavi | Capturing chemical intuition in synthesis of metal-organic frameworks[END_REF] An important part of machine learning algorithms are generative models, which aim to produce new data points that resemble the data shown during training. Best known for their ability to generate realistic images, generative models are capable of mimicking many other types of data, including text, video and, most notably, chemical structures.

Finding and extrapolating relevant patterns in data is at the heart of new data generation.

Generative models have shown their eciency and some promise for the generation of new materials. Most generative models are based on the construction of a material latent space to encode the information of the data set. This continuous material vector space is then used to generate new data points, such as crystal structures. It is also possible to introduce a mapping between a generated material and the desired property, leading to data generation with given properties. [START_REF] Kim | Generative adversarial networks for crystal structure prediction[END_REF] Recently, various generative deep learning models, have been adapted to the general task of suggesting new stable material candidates, most of which are based on either variational autoencoders (VAE) 1417 or generative adversarial networks (GAN). 1820 Notably, a recent method explores the latent space of a trained VAE generative model aiming at targeted exploration of chosen chemical systems. 21 Generative machine learning models are even more developed in the domain of molecular generation, 22,[START_REF] Kang | Conditional Molecular Design with Deep Generative Models[END_REF] with many examples relying on string representations of molecules. [START_REF] Bagal | Molecular Generation Using a Transformer-Decoder Model[END_REF][START_REF] Sattarov | De Novo Molecular Design by Combining Deep Autoencoder Recurrent Neural Networks with Generative Topographic Mapping[END_REF] A number of deep learning architectures for materials science applications have been developed. 2630 Some of the rst attempts to generate new candidate structures using machine learning methods used 3D voxelised image representations, 1416,[START_REF] Long | Constrained crystals deep convolutional generative adversarial network for the inverse design of crystal structures[END_REF] since the generative models were originally applied in the computer vision and image processing communities. Other approaches have used atomic coordinates and lattice parameters directly. [START_REF] Nouira | CrystalGAN: Learning to Discover Crystallographic Structures with Generative Adversarial Networks[END_REF][START_REF] Kim | Generative Adversarial Networks for Crystal Structure Prediction[END_REF] However, an important problem is that these representations are not invariant to the isometries of the space. Graph neural networks are widely and successfully exploited in molecular chemistry [START_REF] Shi | Learning Gradient Fields for Molecular Conformation Generation[END_REF][START_REF] Hoogeboom | Equivariant Diusion for Molecule Generation in 3D[END_REF] and also in crystal property prediction 3336 and do not suer from this drawback. In addition, they can capture the local structure quite accurately. On the downside, they often struggle to capture the long-range dependencies that are crucial in periodic materials. Score-based models 3739 have been reported to achieve the state-of-the-art generative performance. These models are motivated by non-equilibrium thermodynamics and the main idea is to learn a process that gradually corrupts an input by adding Gaussian noise.

The transformations used to reconstruct the input are learned. After a number of iterations, a score model is able to reconstruct the input from nothing but Gaussian noise.

The score-based generative models have recently been developed for the generation of 3D molecules. [START_REF] Shi | Learning Gradient Fields for Molecular Conformation Generation[END_REF][START_REF] Hoogeboom | Equivariant Diusion for Molecule Generation in 3D[END_REF][START_REF] Xu | GeoDi: a Geometric Diusion Model for Molecular Conformation Generation[END_REF]

Problem denition

In this work the goal is to develop a model for the generation of new stable crystal structures that incorporates known symmetry constraints of chemical structures, while being more exible and expressive than existing methods. Furthermore, we seek a method that generates structures with user-specied chemical composition (up to equivalent symmetric positions), so that it can be applied to any chemical system of interest. There are several major hurdles to the application of generative models in materials science. Firstly, the data available is very limited: the open-access Materials Project 41 database contains just over 33, 000 stable structures. This is much less than used for the training of generative models in other elds. In image processing, for example, there are databases containing several million observations or more. Secondly, the nature of the data is quite specic, as the structure should ideally be represented in a way that is invariant to the rigid transformations of space such as rotations and translations. Furthermore, an important characteristic of a crystal structure which makes it particularly dicult to generate but has to be taken into account, is its periodic nature. This challenge has recently been addressed by using graph representations of chemical compounds. 4245 However, the proposed methods still lead to suboptimal performance.

Our contributions

We challenge to overcome the drawbacks of the state-of-the-art methods, and our contribution is manifold: First, we propose a score-based model for chemical structures with periodic nature, namely crystals. To the best of our knowledge, we are the rst to propose a method that applies the score-based model not only to the atomic site coordinates, but also to the lattice itself. This double denoising process results in a better exploration of the search space and leads to a much more meaningful generation of candidate structures. Second, our model takes the number of symmetrically non-equivalent atoms per element as an input. Therefore it can be directly used for the task of targeted exploration of systems of any desired composition without additional training. Third, our model respects the given constraints of a chosen space group, by using a novel symmetry-consistent multigraph description. This helps to capture long-range patterns and improves the quality of the generated structure.

Finally, the proposed approach is validated on the stable structures of the Materials Project database. Our numerical experiments consist of three parts. Firstly, the performance of our model is compared to that of alternative state-of-the art models using statistical metrics. Secondly, it is shown that our score-based model is able to generate existing crystal structures that were not seen by the algorithm during training. Lastly, examining structures that are generated by our model by computing their electronic and mechanical stability by DFT/phonon calculations, we identify original unknown structures of great relevance.

Organization of the paper

Section 2 introduces the crystal structure notations used in our paper. Section 3 provides a general formulation of the Score Matching with Langevin Dynamics (SMLD) model for new data generation. In Section 4 we adapt this approach to the specic problem of generating new crystal structures. Finally, Section 5 presents the results of numerical experiments carried out using our model.

Notations

A crystal structure is described by a unit cell (arrangement of atoms) and its Bravais lattice, which denes the periodic conditions. The unit cell is a region of space bounded by a parallelepiped dened by three lattice vectors and the number, nature and position of the atoms contained in that region of space. A crystal structure is then generated by replicating and translating the unit cell according to its lattice vectors. Formally, we dene a crystal structure as M = (A, L, X, H) 2. L = (l 1 , l 2 , l 3 ) ∈ R 3×3 denotes the lattice matrix dening the form of the unit cell. The whole structure is invariant to translations of the form

∈ A N × R N ×3 × R 3×3 × H,
k 1 l 1 +k 2 l 2 +k 3 l 3 , (k 1 , k 2 , k 3 ) ∈ Z 3 .
Alternatively, the lattice can be dened by the lengths a, b, c of its vectors and the angles α, β, γ between them, denoted by L 6 = (a, b, c, α, β, γ). This representation is invariant to the choice of the coordinate system. When angles are xed, we also use the representation given by L 3 = (φ(a), φ(b), φ(c)), where φ(x) = k log(exp(x/k) + 1)

for some positive real number k. Note that L 3 takes its values in R 3 , while the lengths (a, b, c) are necessarily positive real numbers.

3. X ∈ R N ×3 is the matrix of atomic Cartesian coordinates of the atomic positions that are not equivalent with respect to the symmetry operations. Therefore, XL -1 corresponds to the matrix of the atomic coordinates relative to the unit cell. Note that X is usually much smaller than the full matrix of all atomic coordinates, since there are a lot of equivalent atomic positions in most structures.

4.

H ∈ H = {H 1 , ..., H 230 } denotes the crystallographic space group dening symmetry constraints for the structure. [START_REF] Aroyo | Bilbao Crystallographic Server: I. Databases and crystallographic computing programs[END_REF] A space group H contains the symmetry operations

H = {h i , i ∈ 1, ..., |H|}.
For a given space group, it holds that X and h i (X) dene the same atomic positions, up to permutation, for all i ∈ 1, ..., |H|. Moreover, by the symmetry constraints of H, the position of one atom may determine the position of some of the other atoms in the unit cell. Hence, storing all atomic positions is redundant, and in our model for every set of equivalent positions, we only keep track of the position of one of the atoms. Some of the symmetry groups also put constraints on the lattice.

Depending on the space group H, the lattice is bound to be either triclinic (H

dX(t) = s (X(t)) dt + σdB(t), (1)
where X is the process of interest, t is the time, B is a Wiener process and s is the Stein score of the distribution p dened by s(x) = ∇ x log p(x). At x, the score points in the direction of higher density. Intuitively, the process X(t) tends to push towards the high density regions of p. More formally, one can show that, with σ = √ 2, the invariant distribution of this process is p, i.e., as time t becomes large, the distribution of X(t) tends to the distribution p. Thus, the process X can be used as a sampling method, since simulating the process X(t) yields a distribution close to p, which corresponds to our goal.

From a physical point of view, -log p(x) can be seen as an the energy function U (x)

of a system in state x and p(x) is the corresponding Gibbs distribution for a given choice of the temperature. In this case, the score s in the rst term of the stochastic dierential equation ( 1) corresponds to -∇ x U (x), which pushes the system towards a minimum of its energy function, while the second term represents the external forces.

Our generative model has two main parts: training and sampling. First, during the training, the score s has to be learned from the available data set. Typically, a deep neural network parameterized by θ is trained to give an estimate s θ . Second, new data points are sampled from p using score matching and annealed Langevin dynamics. An Euler scheme is used to obtain a numerical approximation of the Langevin process.

Score matching

Let D = {x 1 , ..., x n } be a dataset of sample size n with unknown distribution p and p 0 (x) =

1 n x∈D δ x (x) the associated empirical distribution, where δ x is the Kronecker delta. As p 0 is not dierentiable, it cannot be used to construct an estimate of the score s. However, a smooth density estimate is given by the Gaussian kernel density estimator pσ dened by

pσ (x) = 1 n x∈D f N dim(x) (x,σ 2 I) (x),
where f N dim(x) (x,σ 2 I) denotes the multivariate Gaussian density of dimension dim(x) with mean x and covariance matrix σ 2 I, where I is the identity matrix and dim(x) the dimension of x. Denote p σ the limit of pσ when n → ∞, which can be expressed as

p σ (x) = f N dim(x) (x,σ 2 I) (x)p(x)dx. (2)
Now, for training the score, say s θ (•, σ), which may be parameterized by a neural network, the straightforward objective function is the mean squared loss between the score estimate and the score of pσ given by

E x∼pσ ∥s θ (x, σ) -∇ x log pσ (x)∥ 2 2 . (3)
simpler objective given by

ℓ σ (θ) = E x∼p 0 ,x|x∼N dim(x) (x,σ 2 I) s θ (x, σ) -∇ x log f N dim(x) (x,σ 2 I) (x) 2 2 . ( 4
)
It has been shown that these two objectives, equations 3 and 4, are equivalent, [START_REF] Vincent | A Connection Between Score Matching and Denoising Autoencoders[END_REF] where two objectives ℓ 1 and ℓ 2 are said to be equivalent if ℓ 1 = αℓ 2 + β for some α > 0, β ∈ R. In particular that means that gradient descent optimization will yield the same result when performed on either objective.

The gradient on the right-hand side of (4) has a simple analytical expression, which is convenient for its numerical evaluation. However, the expectation is not explicit, but can be approximated by Monte-Carlo simulations. In our case, we use minibatch sampling to estimate ℓ σ (θ). Let {x 1 , ..., x m } be a minibatch of size m, then for every

x i , i ∈ 1, ..., m we sample one data point xi ∼ N dim(x i ) (x i , σ 2 I) that is used to compute an empirical counterpart of ℓ σ (θ).
That is, the objective becomes

lσ (θ) = 1 m m i=1 s θ (x i , σ) -∇ xi log f N dim(x i ) (x i ,σ 2 I) (x i ) 2 2 . ( 5
)
When the sample size n tends to innity, the minimizer s θ * (x, σ) of ( 5) equals the score of the distribution p σ almost surely.

Clearly, the choice of the value of σ has an impact on the result. If σ is small, pσ (and p σ ) is a multimodal density, where all mass is concentrated on small neighbourhoods of the data points. As a consequence, an iterative sampling procedure is likely to get stuck in a mode with no means to escape from it. On the contrary, a large value of σ implies a at density pσ , which favours the exploration of the whole data space during sampling, but which is a rather biased estimate of the target distribution p.

To address this issue, one may use Langevin dynamics, [START_REF] Song | Generative Modeling by Estimating Gradients of the Data Distribution[END_REF] that consists in combining the solutions for a predened set of decreasing noise levels σ 1 > ... > σ F > 0 with F being the number of noise levels. This approach provides stable score estimates, and during sampling all modes of the distribution may be reached. The score s θ is trained by minimizing the following objective, which is a weighted average of the loss functions dened above, given by

L(θ) = F k=1 w k lσ k (θ). (6)
By choosing the weights as w k = σ 2 k , the terms in the sum in ( 6) are of nearly the same order of magnitude. [START_REF] Song | Generative Modeling by Estimating Gradients of the Data Distribution[END_REF] As such, no particular noise level is prioritized over others. For computational reasons, we use a unique neural network s θ to estimate the score for all noise levels. To distinguish the scores corresponding to dierent noise levels, we pass the noise level σ k as one of the inputs to the neural network that provides the estimate of the score of p σ k .

Sampling

The trained estimator s θ is used for sampling. The Langevin dynamics process in (1) can be approximated by the Euler scheme

x t = x t-1 + ϵ 2 s(x t-1 ) + √ ϵz t , (7)
where z t has standard normal distribution N dim(xt) (0, I) and ϵ is the step size. Under convenient assumptions, the scheme converges to the stationary distribution of the Langevin dynamics Stochastic Dierential Equation (SDE). Using this scheme and starting from some randomly initialized data point x 1,0 , T iterations according to equation ( 7) are performed.

The nal value x 1,T can be considered as a realization of p σ 1 and is used as the initial value for the sampling at the next noise level, that is, we set x 2,0 = x 1,T . The noise parameter σ 1 should be large enough, so that the process is unlikely to get stuck in a local maximum of the p σ 1 , and the data space can be explored rapidly. The next Langevin simulation samples from p σ 2 and the result is used as the initialization for the sampling from p σ 3 and so on, until we obtain a sample from p σ F . Since we need the nal distribution p σ F to be close to p, one should choose σ F to be small enough. We restrict our attention to space groups whose lattices are neither triclinic or monoclinic.

For all these crystal systems, i.e., for H ∈ {H 16 , ..., H 230 }, the angles between the lattice vectors are xed in the conventional lattice representation. Thus, it is sucient to consider the generation of the lengths of the lattice vectors. As mentioned above, the transformed vector lengths L 3 are real-valued, which is convenient for adding Gaussian noise to the denoising process.

Now, in our model, the Langevin dynamics process is applied to both the atomic coordinates X and the lattice parameters L 3 , as shown in Figure 1. Note that X and L 3 are objects of dierent nature which is new compared to existing score-based models, where the denoising process typically acts only on a single type of object. To handle this issue, we consider the processes on X and L 3 separately. That is, the score s is separated into two parts, one related to the atomic positions, and another corresponding to the lattice. We note To obtain a exible model, we choose dierent step sizes and noise levels for each of the Langevin dynamics processes in our model. Let (ϵ X,k ) k and (ϵ L,k ) k be the step size sequences for the atomic coordinates and the lattice parameters, respectively. They are dened as

ϵ X,k = γ X σ 2 X,k /σ 2 X,F , ϵ L,k = γ L σ 2 L,k /σ 2 L,F ,
where {σ X,i } F i=1 , γ X and {σ L,i } F i=1 , γ L are two dierent sets of parameters. For the sake of faster convergence we choose low temperatures for the lattice Langevin process by additionally multiplying the noise term by a factor of 0.15. The resulting lattices are likely to be rather close to a local maximum of the associated probability density function. For the atomic positions, we keep the normal temperature for a better exploration of the sample space. As illustrated in Figure 2, the sampling process starts o with a randomly initialized structure, with lattice parameters drawn from a Gaussian distribution and atomic coordinates sampled uniformly in the unit cell. Then, one iteration of the equation ( 7) is done in two steps. First, the atomic positions are updated, and then a step over the lattice parameters L 3 is done. Our learning procedure is drafted as Algorithm 1.

Note that the actual number of atoms in the structure is not xed due to the symmetry constraints of the chosen space group H. For example, if an atomic position x coincides with all of its symmetric positions, i.e., h(x) = x for h ∈ H, then there is a single atomic site corresponding to x in A. However, if all symmetric positions are dierent, that is, h(x) ̸ = h(x) for all h, h ∈ H and h ̸ = h, then there are |H| atomic sites in the structure corresponding to x. Therefore, as the atomic positions are modied during sampling, it may happen that the number of atoms changes. Algorithm

for k ← 1 to F do ϵ X,k ← γ X σ 2 X,k /σ 2 X,F ▷ Atomic positions step ϵ L,k ← γ L σ 2 L,k /σ 2 L,F ▷ Lattice parameters step for t ← 1 to T do Draw z X,t ∼ N dim(X) (0, I) Draw z L 3 ,t ∼ N 3 (0, I) X k,t ← X k,t-1 + 1 2 ϵ X,k s θ|X (X k,t-1 , L t-1 , σ X,k ) + √ ϵ X,k z t L k,t ← L k,t-1 + 1 2 ϵ L,k s θ|L (x k,t-1 , L t-1 , σ L,k ) + 0.15 √ ϵ L,k z l,t end for X k+1,0 ← X k,T L k+1,0 ← L k,T
end for

Multi-graph representation

The data points (A, L, X, H) are complex objects that are not conventional for machine learning. In addition, using the atomic coordinates X is prohibitive because this representation is not invariant to rotations and translations. Thus, the choice of an appropriate data representation is an important issue.

First, we denote by X ′ the coordinates of all positions obtained by applying symmetry operations h ∈ H to the atomic positions in X, which contain only one position for each set of equivalent positions in the real structure. Note that the number of coordinates in X ′ can be larger than the number of actual atoms in a structure, since some of the positions in X ′ may correspond to coinciding symmetric positions.

As in some other material generation models, 17,31 we adopt a multi-graph crystal representation (V, E, W ), which has the desired invariance properties, while preserving all structural information. In our graph representation, however, the graph vertices V correspond to atomic positions in X, i.e., non-equivalent atomic positions within the unit cell. The edges E correspond to the interatomic bonds, such that one end of the bond is in X and the other is in X ′ . However, since there is a single node for each set of equivalent positions in X ′ , these bonds create multiple edges and self-loops in the graph representation. Another source of multiple edges and self-loops is the periodicity of the structure since each atom can be connected to multiple copies of each atom obtained by the periodicity. W are the angles between these edges, i.e., angles W (i, j, k) = α ijk formed by triplets of nodes (i, j, k) such that (i, j) ∈ E and (j, k) ∈ E.

Introducing symmetry relations into the multigraph helps to account for the long-range atomic relationships within the structure.

Additional details can be found in the Supplementary Information. Here, we describe the estimation of the score s |X , while the case of s |L is presented in Section 4.3. Recall that for a given pair of vertices, say i and j, there may be multiple edges e ijk due to the periodicity of the structure. For each edge e ijk , denote by dx(e ijk ) ∈ R 3 the vector of the dierences between the coordinates of the associated atomic sites. Let w ijk be the real output provided by the GNN for edge e ijk . Then, the score of the atomic coordinates of a given node i with coordinates x i is obtained as a weighted sum of normalized vectors dx(e ijk ) over all edges originating at this node, that is,

s θ|X (x i ) = j∈N (i) k:e ijk ∈E(i,j)
w ijk dx(e ijk ) ∥dx(e ijk )∥ ,

where N (i) is the set of neighbours of node i and E(i, j) is the set of multiple edges between i and j. Since we want a unique score vector for each set of symmetry equivalent atomic positions, the scores for nodes x i that were added by symmetries are not evaluated. In other words, the score is calculated only for the positions in X ′ that are also present in X.

A nal post-sampling step is added to the procedure, which removes some of the atoms in X ′ that we consider to come from symmetry operations applied to special positions in the unit cell. More precisely, any two atoms obtained by applying symmetry operations to an atomic site are considered to be the same atom if the distance between them is less than a given threshold, which we x to 0.3 Å. Thus, it is the number of non-equivalent atomic positions that is xed, not the exact composition.

Our implementation uses a Dimenet++ graph neural network [START_REF] Gasteiger | Fast and Uncertainty-Aware Directional Message Passing for Non-Equilibrium Molecules[END_REF] as the estimator for our model, which uses Bessel basis representation for the interatomic distances and Fourrier-Bessel basis representation for the angles between bonds in its message passing function.

The original Dimenet++ aggregates the messages in the atoms and performs pooling over the node features at the end, producing an output in the form of a vector for the whole structure. However, our method requires an output for each edge in the graph, so we modied the output blocks by replacing the dense layers on the atomic features with dense layers on the edge features and removing the pooling.

Crystal cell generation

The lattice, which denes the periodicity of the crystal structure, is not unique. Usually, either the conventional unit cell or the primitive unit cell is used. The former can have up to 4 times more atoms in face-centered cubic unit cell, but the constraints it has to respect are much simpler than those of the primitive cell. There is an invertible mapping between the two types of the unit cell, so we assume that for each primitive unit cell there is a corresponding conventional unit cell and vice versa.

As the symmetry constraints are easier to handle in the conventional unit cell, we consider a Langevin process acting on the lattice parameters of the conventional unit cell and aim to learn to denoise them. However, a large number of atoms in the unit cell implies more stringent requirements in terms of memory and the execution time for the GNN. Therefore, , a colinear vector with the gradient that changes the edge length by the predicted amount. We then take the mean of these estimates for all edges e, yielding the nal estimate of the noise of the lattice parameter given by

to
∆L 3,θ = 1 |E| e∈E d e,θ G e ∥G e ∥ 2 2
.

The lattice score estimate s θ|L is then given by s θ|L = ∆L 3,θ σ k .

20

Numerical results

In this section the performance of our model is assessed. First, it is compared to other recently developed generative models via structure descriptor-based metrics on a set of structures sampled by the dierent models. Then, the model's capacity of generating known and unknown structures is investigated for a variety of chemical systems.

Materials project data

Our study uses most of the nearly stable structures (energy above hull is below 0.1eV /atom) available in the Materials Project [START_REF] Jain | Commentary: The Materials Project: A materials genome approach to accelerating materials innovation[END_REF] (we accessed the data in February 2022). There are 86, 508 structures in total in our experiments.

Each structure is stored in a CIF le containing a matrix L dening the lattice, and a matrix X dening the atomic positions and the composition A. The space groups H were identied using the Phonopy package. [START_REF] Togo | First principles phonon calculations in materials science[END_REF] It is important to emphasise that we did not use any additional information about the chemical properties of the structures other than the atomic number to identify the element species. As shown in Figure 4, 30, 358 structures belong to the monoclinic and triclinic crystal families. While the atomic position score s θ|X is learned on the entire dataset, the lattice score s θ|L 3 does not use monoclinic and triclinic crystals in training, as we do not consider crystal families with variable angles in conventional lattice representation. Furthermore, the distribution of the space groups of the structures, displayed in Figure 5, is far from uniform, with a considerable number of space groups not being represented in the data at all. We also see that the dataset is heavily biased towards the oxygen element, which occupies 35%

of the atomic sites and with 40, 282 structures containing this atom (Figure 6). As our model allows the user to specify the desired number of non-equivalent atomic sites, we examine the distribution of the number of non-equivalent atomic sites per element in the database. As shown in Figure 7, in most cases there are very few non-equivalent sites per element in any structure, mostly only one or two per element. Based on this observation, it is reasonable to expect that new stable compounds will also be of rather simple nature.

Comparison with other models

For a comparative study, we consider two recent state-of-the-art crystal generative models as a baseline. First, the Crystal Diusion Variational Autoencoder (CDVAE) 17 is a variational autoencoder that uses a diusion process to denoise atomic coordinates from a latent representation. In that respect, this model is the most similar to ours. The second one is the Physics Guided Crystal Generative Model (PGCGM), [START_REF] Zhao | Physics guided deep learning for generative design of crystal materials with symmetry constraints[END_REF] which uses a Generative Adversarial Network to sample crystal structures within a given set of 20 space groups.

For a fair comparison, all models should be trained on the same dataset. We trained CDVAE and our model on the so-called MP-20 dataset proposed in the CDVAE paper, 17 which consists of the Materials Project structures that have at most 20 atomic sites in its unit cell. We also removed the structures belonging to monoclinic and triclinic space groups, thus reducing the number of materials in the database from [START_REF] Belsky | New Developments in the Inorganic Crystal Structure Database (ICSD): Accessibility in Support of Materials Research and Design[END_REF] and Open Quantum Materials Database (OQMD). [START_REF] Kirklin | The Open Quantum Materials Database (OQMD): assessing the accuracy of DFT formation energies[END_REF] Concerning sampling, a major dierence between our model and the other two is that our model takes the chemical composition (up to the number of symmetrically equivalent atomic sites) and the space group operations as input. Therefore, when sampling with our model, we feed the model with the composition and space group of one of the materials in the test set.

For each model, 1000 structures are generated and we compute the corresponding coveragerecall (COV-R) and coverage-precision (COV-P) metrics proposed in the CDVAE paper, 17 that are based on the comparison of the CrystalNN and Magpie ngerprints of the generated structures and the structures in the test set. The larger the metrics, the better generated structures imitate the test data. More precisely, COV-P is dened as the proportion of generated structures that have close neighbours in the feature space among the test structures.

Inversely, COV-R is the proportion of test data structures that are close to the generated structures. Notably, COV-P heavily depends on the size of the test dataset, while COV-R heavily depends on the number of sampled structures.

Inspired by the standard for the image generation task, Fréchet inception distance metric, [START_REF] Heusel | Trained by a Two Time-Scale Update Rule Converge to a Local Nash Equilibrium[END_REF] we propose to consider two similar metrics for crystal structures. The metric assumes that the structure descriptors are normally distributed for both the generated and the test data and computes the Fréchet distance between those distributions. The lower the distance, the more similar the test and generated data are. We consider CrystalNN as the rst descriptor, and the activations of the second to last layer of the Crystal Graph Convolutional Neural Network (CGCNN) [START_REF] Xie | Crystal Graph Convolutional Neural Networks for an Accurate and Interpretable Prediction of Material Properties[END_REF] as the second. The distributions were calculated based on 1000 random samples for each model.

The results are displayed in the Table 1, where we note these metrics FD for Fréchet Distance. In case of the PGCGM model, since the model generates structures from only 20

space groups, we also compute the comparison metrics with the test data containing only structures from the same 20 space groups, the result shown in numbers between parentheses.

The sharp drop in COV-P can be explained by the dependence of the metric on the size of the test data set. We see that our model performs well across the board outperforming or nearly matching the performance of the other models. Finally, to examine the validity of the proposed metrics, we also evaluate these metrics for 1000 stable structures randomly chosen from the training dataset. Since we can assume that these structures come from the same distribution as the test structures, the corresponding structure descriptors should also have close distributions, resulting in close to perfect metric values. We observe from Table 1 that the results are matching these expectations, supporting the validity of the chosen metrics.

Sampling known stable structures

The goal here is to assess the ability of our model to generate stable structures, which are not used during training. Here and further below we used our model trained on the complete Materials Project dataset as described in the Section 5.1. First we did a 9:1 train-test split, using the test data to track the loss and choose the stopping point. In total, we were able to reconstruct 11 of the 36 structures, while others are close but have dierent atomic arrangements. Concerning the lattice parameters in the reconstructed structures they may be slightly dierent from those in the original structures, but these deviations can simply be corrected by applying the DFT relaxation scheme.

Thus, we conclude that our model is able to reproduce known structures with a large number of dierent space groups across all targeted crystal families. In particular, the lattice denoising is indeed able to capture the desired lattice shape, while adapting to the constraints imposed by the chosen space group.

Generating new stable structures

To test the ability to generate new stable crystal structures, we attempt to sample new stable ternary compounds. We chose several chemical systems with no known stable ternary compounds in most of them. As for the space group, we prioritize 32 of the most common (Figure 5) space groups, sampling 60 structured per space group per composition. Taking into account the distribution of the number of non-equivalent atomic sites per element (Fig-

ure 7), which has to be xed before sampling, we chose one or two non-equivalent positions for each of the three elements in the composition, with a total of at most four non-equivalent sites per structure.

We have tested our model for several compositions in several systems. Table 3 displays those that are considered to be the most stable by DFT and we also validated their mechanical stability by additional phonon calculations. It turns out that with the estimated formation energy, these structures lie close to the convex hull formed by other known binary congurations composed of the same elements and calculated under the same DFT conditions. Although no structure seems to support the known ground state, it appears that the generated phases are close to stable and remind us of the diculty of discovering higher-order compounds (more than binaries, such as ternary compound). The method demonstrates its One of the challenges to be addressed in the near future is to speed up the generative process to explore the chemical composition space more extensively. However, one has to keep in mind that the verication procedure by DFT calculation remains costly but necessary, as an accurate regression is still not precise enough.

29

Conclusion and outlook

Our model successfully applies a score-based generative method to the task of generating new materials. It compares favourably with existing generative models, showing matching or superior results in terms of the studied metrics.

The model is able to reconstruct some of the existing stable structures, in particular it succeeds in denoising the lattice parameters together with the atomic positions. The method performs well across a wide range of space groups and lattice systems, and is not restricted to compositions of a certain number of chemical elements.

We have also shown that the model is capable of generating entirely new stable structures.

The model generates structures with a chosen space group and composition, which can be potentially useful in determining the structure of experimentally known materials with unknown unit cells. More generally, the model may be combined with any random sampler of compositions and space groups.

It is noteworthy that only the atomic numbers have been entered to indicate the nature of the selected elements, without any additional information such as chemical descriptors that would provide explicit correlations between the properties of dierent elements. An interesting line of research involves the integration of such prior chemical knowledge to further enhance the model's performance.

At present, sampling a single structure takes several minutes, which in turn makes validation of the model cumbersome. The main bottleneck of the current implementation is the construction of the multigraph, which is performed on CPU and is not yet parallelized.

Another possible valuable extension of the model consists in conditioning the model on the exact composition of the desired structure. At present, the number of non-equivalent sites is xed in advance, but the nal ratio of total sites per element still varies widely, which is not always desired. 30

where the following notations are used: 1 .

 1 A ∈ A N denotes the composition of a structure with N atoms per unit cell. It indicates the chemical element of every atom in the unit cell and A is the set of all chemical elements.

4

  New crystal score-based modelHere we develop a score-based model specically tailored for the generation of crystal formations. Our approach is such that the user can choose the elements and their composition in a new structure that is generated by the model. This oers the possibility of a targeted search taking into account prior expert knowledge or the practitioner's intuition. More precisely, to sample a new structure M = (A, L, X, H) the user rst chooses the composition A and also the symmetry space group H, which helps the algorithm to sample realistic structures.

Figure 1 :

 1 Figure 1: Model training. While learning the positions score s θ|X , the lattice remains unchanged. While learning the lattice score s θ|L , both the lattice and the atomic positions are perturbed. Orange arrows correspond to the score estimates given by the model as shown on Figure 3.

Figure 2 :

 2 Figure 2: Sampling with annealed Langevin dynamics. Orange arrows correspond to the steps of the Langevin dynamics simulations, as shown in Figure 3. The sample from the most perturbed distribution, corresponding to σ 1 is used as initialization for the subsequent dynamics, denoted as σ 2 .

4. 2

 2 Graph neural networkWith the multigraph representation at hand, it is natural to use Graph Neural Networks (GNNs), a class of deep neural networks that operate on graph data. The forward pass of a GNN generally consists of iteratively updating the multigraph representation by passing messages from each node to its neighbours. In the context of crystal generative models, a GNN learns an appropriate edge output to estimate the score function, as elaborated further below. All the steps of the score evaluation are illustrated on Figure3and are represented by the orange arrow in Figures1 and 2.

Figure 3 :

 3 Figure 3: Score estimation, representing the orange arrow in Figures1 and 2. The crystal structure is dened by its lattice parameters L 3 and the atomic positions X. All symmetric positions X ′ are added. A multigraph is then constructed and passed to a graph neural network estimator that outputs the score.

  estimate the score of the perturbed conventional unit cell parameters, we still use the multigraph representation associated with the corresponding primitive cell which we feed to the network. During training, the noise ∆L 3 , which is added to the lattice representation L 3 , has to be computed. The aim of lattice score training is such that the GNN computes an estimate of s |L 3 = ∆L 3 σ k . Note that changing the lattice also changes the Cartesian coordinates of the atomic sites, and the estimate of ∆L 3 is based on the changes in relative distances between atoms induced by the deformation of the lattice. More specically, given a multigraph representation (V, E, A) of the structure, for each edge e ∈ E the network provides the amount of contraction/dilatation ∆d e of the distance d e between the atoms associated with edge e due to the lattice deformation. Independently, the gradient G e = ∇ L 3 d e of d e with respect to the conventional lattice parameters L 3 is approximated numerically. Indeed, for a suciently small lattice noise ∆L 3 , it holds that ⟨G e , ∆L 3 ⟩ ≈ ∆d e , where ⟨a, b⟩ denotes the scalar product between a and b. So, ∆L 3 is naturally approximated by ∆d e,θ G e ∥G e ∥ 2 2

Figure 4 :

 4 Figure 4: Distribution of crystal families in the dataset used in training. Monoclinic and triclinic systems are excluded from the dataset when training s θ|L 3 .

Figure 5 :

 5 Figure 5: Distribution of the space groups of the structures in the dataset.

Figure 6 :

 6 Figure 6: Number of non-equivalent atomic sites by element in the data set.

Figure 7 :

 7 Figure 7: Distribution of the number of non-equivalent atomic sites per element in the data set.

Figure 8 :

 8 Figure 8: Example of new generated stable compounds veried by rst-principles calculations. The phonon dispersion curves indicate no imaginary frequency. The electron localization function (ELF) is mapped onto the range 0 (blue) to 1 (red) and represents the probability of nding an electron in the neighbourhood space of a reference electron located at a given point.

  1 Sampling with Langevin dynamics Require: {σ X,k } F

	k=1 , γ X	▷ Position denoising process coecients
	{σ L,k } F k=1 , γ L	▷ Lattice denoising process coecients
	T	▷ Step count per noise level
	Output: X, L	▷ Atomic positions, lattice parameters
	Initialize X 1,0 , L 1,0	▷ Initial atomic coordinates and lattice parameters

  45, 229 to 36, 334. We used the same 60:20:20 split of the data into training, test and validations sets to train the CDVAE and our model. The PGCGM material generator is publicly available albeit lacking the code to train on new data. It is known that PGCGM is trained on data collected from Material Project, Inorganic Crystal Structure Database (ICSD)

Table 1 :

 1 Comparative study of our model to PGCGM and CDVAE using various metrics. Values of the metrics for 1000 randomly picked stable structures.

	Model	COV-P	COV-R	FD CrystalNN	FD CGCNN
	This work	0.972	0.965	0.102	0.682
	PGCGM	0.944 (0.864)	0.76 (0.765)	0.225 (0.303)	4.359 (4.209)
	CDVAE	0.982	0.862	0.199	0.703
	1000 stable	0.998	0.968	0.012	0.044
	structures

  Some of the experimentally observed structures are not included in the training data set, in particular those structures that are unstable at 0K and no pressure. Taking advantage of this fact, we attempt to sample these structures with our model by providing the correct space group and elemental composition, that is, the number of non-equivalent positions for each element. This way, we hand-picked 36 structures that have been observed experimentally in the Materials Project, but are not included in the training set. These 36 structures have been chosen to provide a set of samples as representative as possible of the whole chemical space (dierent crystal families and dierent types of chemical bonds). For each structure, we then generate 50 samples with the same number of non-equivalent positions, with the space group symmetries xed in advance, and perform DFT relaxation on the output structures. This results in several original structures being found in the generated samples, some of them are shown in the Table2. The complete list of tested structures with their material project IDs can be found in the Supplementary Material.

Table 2 :

 2 Some of the experimentally known structures that were successfully generated by the model when the correct composition and the space group were given as input.

	Elements	Crystal family	Generated	Lattice	Lattice
		and space	structure	parameters of	parameters
		group number	before DFT	the sampled	according to
			relaxation	structure, Å	Materials
					project, Å
	Li-Mg-P	Cubic, 216		4.43	4.08
	Ta-Mn-O	Hexagonal, 191		5.207, 3.12	5.1, 3.07
	Cr-Co-P	Orthorhombic,		3.386, 5.889,	3.52, 5.73, 6.68
		62		6.949
	Tl-Cd-S	Hexagonal, 187		4.81, 6.33	3.89, 7.08
	Ce-In-Cu	Hexagonal, 194		5.398, 7.592	4.74, 6.73
	Mn-Ca-Sn	Tetragonal, 129		4.471, 7.387	4.53, 7.35
	Ba-Ti-O	Tetragonal, 123		4.326, 5.587	4.11, 5.04
	Pu-I-O	Tetragonal, 129		4.0, 6.7	4.01, 9.47

Table 3 :

 3 Samples for compositions with no known stable ternary compounds.

	Chemical	Crystal family	Generated	Heat of	Energy from
	system	and space	composition	Formation,	the Hull,
		group		kJ/mol	kJ/mol
	Cu-Si-K	Orthorhombic,	Cu 2 KSi	-5.29	+0.61
		P nma
	Fe-Ta-Sn	Orthorhombic,	Fe 2 TaSn	-8.96	+4.14
		F mmm
	Sb-Ta-V	Hexagonal,	SbTa 2 V 3	-12.55	+1.57
		P 6 3/mmc
	Al-K-Pd	Hexagonal,	Al 3 KPd 6	-55.77	+15.38
		P 6/mmm

Acknowledgements

Calculations were performed using HPC resources from GENCI-CINES (No. A0060906175)

In addition, we acknowledge the nancial support from the CNRS-MITI (CRISTAT project).

Data and Software Availability

The implemented code along with the used data is available online at https://github.com/ findooshka/diffusion-atoms.