Substitutional Alloying Using Crystal Graph Neural Networks
Materials discovery, especially for applications that require extreme
operating conditions, requires extensive testing that naturally limits the
ability to inquire the wealth of possible compositions. Machine Learning (ML)
has nowadays a well established role in facilitating this effort in systematic
ways. The increasing amount of available accurate DFT data represents a solid
basis upon which new ML models can be trained and tested. While conventional
models rely on static descriptors, generally suitable for a limited class of
systems, the flexibility of Graph Neural Networks (GNNs) allows for direct
learning representations on graphs, such as the ones formed by crystals. We
utilize crystal graph neural networks (CGNN) to predict crystal properties with
DFT level accuracy, through graphs with encoding of the atomic (node/vertex),
bond (edge), and global state attributes. In this work, we aim at testing the
ability of the CGNN MegNet framework in predicting a number of properties of
systems previously unseen from the model, obtained by adding a substitutional
defect in bulk crystals that are included in the training set. We perform DFT
validation to assess the accuracy in the prediction of formation energies and
structural features (such as elastic moduli). Using CGNNs, one may identify
promising paths in alloy discovery.
Dario Massa, Daniel Cieśliński, Amirhossein Naghdi, Stefanos Papanikolaou
# Substitutional Alloying Using Crystal Graph Neural Networks

###### Abstract

Materials discovery, especially for applications that require extreme operating conditions, requires extensive testing that naturally limits the ability to inquire the wealth of possible compositions. Machine Learning (ML) has nowadays a well established role in facilitating this effort in systematic ways. The increasing amount of available accurate DFT data represents a solid basis upon which new ML models can be trained and tested. While conventional models rely on static descriptors, generally suitable for a limited class of systems, the flexibility of Graph Neural Networks (GNNs) allows for direct learning representations on graphs, such as the ones formed by crystals. We utilize crystal graph neural networks (CGNN) to predict crystal properties with DFT level accuracy, through graphs with encoding of the atomic (node/vertex), bond (edge), and global state attributes. In this work, we aim at testing the ability of the CGNN MegNet framework in predicting a number of properties of systems previously unseen from the model, obtained by adding a substitutional defect in bulk crystals that are included in the training set. We perform DFT validation to assess the accuracy in the prediction of formation energies and structural features (such as elastic moduli). Using CGNNs, one may identify promising paths in alloy discovery.

**Keywords:** GNN, CGNN, Graphs, Substitutional Alloys, Materials

Discovery, Neural Networks, Machine Learning
## 1 Introduction

The use of machine learning (ML) Michalski et al (2013); LeCun et al (2015) methods in material science to accelerate materials discoveryCurtarolo et al (2013) is at the base of the so-called material informatics (MI) Ramakrishna et al (2019); Ramprasad et al (2017); Takahashi and Tanaka (2016); L. Ward (2017); Rajan (2005). By training ML models on large databases, such as OQMD or the Materials Project high-throughput electronic structure calculation databases Saal et al (2013); Jain et al (2013); Curtarolo (2012); Hachmann et al (2011); NOMAD ([https://nomad-coe.eu](https://nomad-coe.eu)), the goal is to achieve predictions of material properties with quantum accuracy.

As in statistical mechanics with the need for identifying appropriate order parameters of novel phases and structures, the key challenge in ML algorithms is to identify effective system descriptors that can function as structure identifiers. A large variety of descriptors have been proposed, including fixed-length feature vectors of material elemental or electronic properties Seko et al (2015); Xue et al (2016); Isayev et al (2017), as well as structural descriptors, based on rotational and traslational invariant transformations of atomic coordinates, like the Coulomb matrix Rupp et al (2012), atom-centered symmetry functions (ACSFs) Behler (2011), social permutation invariant coordintes (SPRINT) Pietrucci and Andreoni (2011), smooth overlap of atomic positions (SOAP) De et al (2016) and global minimum of root mean-square distance Sadeghi et al (2013). However, these solutions are often system-specific, and are not suitable for vast compositional and structural space exploration.

For this reason, a topic of fervent interest in the materials science community is the use of graph neural networks (GNNs) Zhou et al (2020); Wu et al (2021), which allow to learn representations directly and in a flexible way,
focused on molecular systems Jorgensen et al (2018); Schutt et al (2017); Duvenaud et al (2015); Wu et al (2018); Kearnes et al (2016); Coley et al (2017), surfaces Back et al (2019); Palizhati et al (2019); Gu et al (2020) and periodic crystalsSchutt et al (2017); Xie and Grossman (2018, 2018); Chen et al (2019); Dunn et al (2020); Louis et al (2020); Park and Wolverton (2020); Karamad et al (2020); Chen et al (2019). GNNs can be regarded as the generalization of convolutional neural networks (CNN) to graph-structured data, from which the internal materials representations can be learned and used for prediction of target properties Reiser et al (2022); even though larger amounts of data is required with respect to conventional ML models, GNNs take advantage of the unambiguous physics-guided real-space local associations between the system's degrees of freedom, hence for any type of atomic crystalline structure Gong and Yan (2021). The common idea of GNN-based models is to represent atoms as nodes (V) and their chemical bonds as edges (E) in a graph G(V,E), which can be fed to a trained neural network to create node-level embeddings (learned representations of each atom in its individual chemical environment) through convolutions with neighbouring nodes and edges Fung et al (2021). Therefore, given a set of learnable weights (W), and (y) a target material property, the GNN model reformulates the prediction task as the mapping \(f(G:W)\rightarrow\) y.

A direct benefit of the crystal material GNN-converted graph encoding is the naturally derived vector characterization of the atoms and edges Louis et al (2020). The work of Xie _et al_. Xie and Grossman (2018) represents the pioneering example of a crystal graph convolutional neural network (CGCNN) architecture, which has been later extended in the iCGCNN by Park _et al_. Park and Wolverton (2020) to include 3-body correlations on neighbouring atoms, information on the Voronoi tasselated structure and an optimized chemical representation of interatomic bonds in the crystal graphs.
For the discovery of new materials, one may take various exploring paths, involving high-throughput computational Curtarolo et al (2013) and experimental Liu et al (2019) methods. However, the combined approach of machine-learning methods and compositional manipulation, has very quickly acquired a well established role in materials science, and it is applied in a wide range of property optimization searches like for zinc blende semiconductors Mannodi-Kanakkithodi et al (2022), perovskites Zhai et al (2022); Balachandran et al (2018); Ye et al (2018); Sharma et al (2020); Klug et al (2017); Sampson et al (2017) and others Guan (2019); Ning et al (2017); Oba et al (2018); Deml et al (2015, 2014); Wan et al (2021); Varley et al (2017); Mannodi-Kanakkithodi et al (2020).

In this work, we utilize a particular improvement of the originally proposed Xie and Grossman (2018) CGCNN model, the MatErials Graph Network (MEGNet) model from Chen, Ye and coworkers Chen et al (2019), introduced in Sec.(2.1.1), that has the merit of being developed and tested both on molecules and crystals, with the possibility of defining global state attributes including temperature, pressure and entropy. We test the capabilities of graph networks to predict the properties of single-atom substitutionally defected crystals with the MEGNet model. After considering a pre-trained model on the Materials Project (MP) database (Sec.(2.1.3)), we focus on the formation energies, bulk and shear moduli predictions, both comparing the results obtained in datasets of similarly defected structures (Sec.(3.1)) and the effects of almost all the possible single-atom defects in the same matrices (Sec.(3.2)). To validate the predictions, as described in Sec.(2.2) and Sec.(3.3), we perform Density Functional Theory (DFT) calculations, and we find CGNNs have both a great potential, but also limitations in predicting properties of defected bulk crystals, and promote materials discovery.
## 2 Methods

### Machine Learning framework

#### MEGNet description

In the present work, we utilize the MEGNet modelChen et al (2019). The reasons for this choice lie in the structure and performance of the model:

1. It is characterized by a low number of attributes, one for the atom (atomic number) and one for the bond (spatial distance), but MEGNet outperforms previous graph-based modelsChen et al (2019), as CGCNNXie and Grossman (2018) and MPNNJorgensen et al (2018), with higher number of attributes, as well as SchNetSchutt et al (2017), with a similar low number.
2. the MEGNet framework includes a global state attribute, essential for state-property relationship predictions in materials,
3. The graph network construction of MEGNet has been developed and tested for both molecules and crystals

Here, we limit ourselves to present the main features of the model, but for a more exhaustive explanation we recommend the reader to the original work of Chen _et al._Chen et al (2019) and references there-in. In particular, given a graph \(G(E,V,\mathbf{u})\), where

\begin{table}
\begin{tabular}{|l|r|r|} \hline Parameter & Value & Short Description \\ \hline \hline nfeat\_node & 94 & number of atom features \\ \hline nfeat\_global & 2 & number of state features \\ \hline ngauss\_centers & 110 & number of gaussians \\ \hline converter\_cutoff & 4 & cutoff radius \\ \hline megnet\_blocks & 3 & number of MEGNetLayer blocks \\ \hline optimizer & Adam & optimizer of the model weights \\ \hline lr & 1e-3 & learning rate \\ \hline n1 & 64 & number of hidden units in layer 1 \\ \hline n2 & 32 & number of hidden units in layer 2 \\ \hline n3 & 16 & number of hidden units in layer 3 \\ \hline \end{tabular}
\end{table}
Table 1: Parameters from the pre-trained MEGNet model.
* \(V\) is the set of \(N^{v}\) atomic attribute vectors \(\mathbf{v}_{i}\);
* \(E=\left\{\left(\mathbf{e}_{k},\mathbf{r}_{k},\mathbf{s}_{k}\right)\right\}_{k=1 \dots N}^{e}\) is the set of \(N^{e}\) bond attribute vectors, with \(\mathbf{r}_{k}\) and \(\mathbf{s}\) being the indexes of the connected atoms;
* \(\mathbf{u}\) is the global state attribute vector;

the role of graph network is to recursively update a input graph \(G(E,V,\mathbf{u})\) to an output graph \(G(E^{\prime},V^{\prime},\mathbf{u}^{\prime})\), with a progressive and inclusive information flow going from bonds to atoms, and finally to the global state. In particular, first the attributes of each bond are updated through a function \(\phi_{e}\), applied on the concatenation of the self-attributes, the ones of the connecting \(\mathbf{v}_{s_{k}}\) and \(\mathbf{v}_{r_{k}}\) atoms, and of the global state \(\mathbf{u}\), as in

\[\mathbf{e}_{k}^{\prime}=\phi_{e}\big{(}\mathbf{v}_{s_{k}}\bigoplus\mathbf{v}_ {r_{k}}\bigoplus\mathbf{e}_{k}\bigoplus\mathbf{u}\big{)} \tag{1}\]

The update of atomic attributes involves the average over i-th atom connecting bonds \(\bar{\mathbf{v}}_{i}^{e}=\dfrac{1}{N^{e}}\sum_{k=1}^{N_{i}^{e}}\left\{ \mathbf{e}_{k}^{\prime}\right\}_{r_{k}=i}\), the i-th atom self-attributes \(\mathbf{v}_{i}\) and the global state ones \(\mathbf{u}\), as in

\[\mathbf{v}_{i}^{\prime}=\phi_{v}\big{(}\bar{\mathbf{v}}_{i}^{e}\bigoplus \mathbf{v}_{i}\bigoplus\mathbf{u}\big{)} \tag{2}\]

Finally, an information flow from all three attribute groups is involved in the update of the global state attibutes, as in

\[\mathbf{u}^{\prime}=\phi_{u}\big{(}\bar{\mathbf{u}}^{e}\bigoplus\bar{ \mathbf{u}}^{v}\bigoplus\mathbf{u}\big{)} \tag{3}\]

where \(\bar{\mathbf{u}}^{e}=\dfrac{1}{N^{e}}\sum_{k=1}^{N^{e}}\left\{\mathbf{e}_{k}^ {\prime}\right\}\) and \(\bar{\mathbf{u}}^{v}=\dfrac{1}{N^{v}}\sum_{i=1}^{N^{v}}\left\{\mathbf{v}_{i}^ {\prime}\right\}\).

As mentioned before, for our systems of interest, namely periodic crystals, the atomic number is the only node attribute. For bonds, the spatial distance
is expanded in a Gaussian basis set, centered at a linearly spaced \(r_{0}\) locations between \(r_{0}=0\) and \(r_{0}=r_{\rm cut}\), and characterized by a given width \(\sigma^{1}\). Finally, the global state is simply a two zeros placeholder for global information exchange.

#### Data collection

We consider crystal structures collected through the Python Materials Genomics interface (pymatgen) Ong et al (2013) to the Materials Application Programming Interface from Materials Project Jain et al (2013). When creating the dataset, there were 126301 structures in the database, that had formation energy (E\({}_{\rm form}\)) property, and 13102 structures that had bulk modulus (K\({}_{\rm VRH}\)) and shear modulus (G\({}_{\rm VRH}\)) properties.

#### Pre-trained model

Our focus in this work is on the prediction capabilities of MEGNet for _minimally_ defected systems, that are clearly not in the training database, given the training of a dataset of undefected structures. In order to do so, we consider the substitution of a single atom in a supercell, hoping that CGNN

Figure 1: Parity plots for the pre-trained model on the MP dataset. The plots involve the predictions on the bulk modulus (K\({}_{\rm VRH}\)), the shear modulus (G\({}_{\rm VRH}\)) and the formation energy (E\({}_{\rm form}\)).
training captures atomic similarities, based on combinations of atomic radii, valence electrons, and other atomic properties. Table(1) shows some of the paramaters of the model, and a more complete list can be found at the default implementation of the class Chen et al (2022)

We report parity plots of Fig.(1) for all three properties of interest in this study: bulk modulus (K\({}_{\rm VRH}\)), shear modulus (G\({}_{\rm VRH}\)) and formation energy (E\({}_{\rm form}\)). To evaluate the model accuracy in predicting the properties of interest for the present study, the mean-absolute error (MAE) is used as the evaluation metric. Table(2) presents the MAE values for each predicted property over the dataset, which provides insights into the pre-trained model performance.

### Validation with DFT

We verify the accuracy of the model's predictions for system properties such as bulk modulus K\({}_{\rm VRH}\), shear modulus G\({}_{\rm VRH}\) and formation energy E\({}_{\rm form}\), after single-atom substitution is implemented in \(2\times 2\times 2\) supercells. We perform DFT calculations with Quantum Espresso (QE) Giannozzi et al (2009, 2017, 2020) and its THERMO_PW Dal Corso (2023) driver for the calculation of structural properties. Pseudo-Potentials for all involved atomic species are Ultra-Soft and with Perdew-Burke-Ernzerhof (PBE) Perdew et al (1996) functional. The Methfessel-Paxton smearing Methfessel and Paxton (1989) has been introduced to correctly investigate metallic systems, and the calculations have been set as spin-polarized, for possible non-zero magnetization effects.

\begin{table}
\begin{tabular}{|l|r|} \hline Property & MAE \\ \hline \hline K\({}_{\rm VRH}\) & 6.143 GPa \\ \hline G\({}_{\rm VRH}\) & 10.489 GPa \\ \hline E\({}_{\rm f}\) & 0.029 eV/atom \\ \hline \end{tabular}
\end{table}
Table 2: MAEs of the model for the prediction of the bulk modulus (K\({}_{\rm VRH}\)), shear modulus (G\({}_{\rm VRH}\)) and formation energy (E\({}_{\rm form}\)).
Convergence is checked on the number of k-points, plane-wave cutoff energy, and also, the energy smearing spreading (degauss parameter in QE), for each case, after a preliminary variable-cell relaxation of pure crystals, and then, further fixed-cell relaxation with optimal parameters, for final equilibrium bulk structures. The common acceptance threshold in the variation of the total energy upon parameter change is set at \(10^{-5}\) Ry. Forces and total energy convergence thresholds for ionic minimization are set to a common value of \(10^{-5}\)a.u. and \(10^{-6}\)a.u. respectively. After optimization of pure crystals, fixed-cell relaxation is performed on supercells with single-atom substitutions, and its structural properties are then extracted through the THERMO_PW driver.

The computation of formation energies for validation purposes is performed on the case of single-atom substitutional defects (D) applied on pure bulk crystals (M), as follows,

\[\mathrm{E_{form}}\big{(}\mathrm{M_{1-x}D_{x}}\big{)}=\mathrm{E}\big{(} \mathrm{M_{1-x}D_{x}}\big{)}-\big{(}1-\mathrm{x}\big{)}\mathrm{E}\big{(} \mathrm{M}\big{)}-\mathrm{x}\mathrm{E}\big{(}\mathrm{D}\big{)} \tag{4}\]

where x is the atomic fraction of substitutional defects, \(\mathrm{E}\big{(}\mathrm{M_{1-x}D_{x}}\big{)}\) is the total energy per atom of the compound, while \(\mathrm{E}\big{(}\mathrm{M}\big{)}\) and \(\mathrm{E}\big{(}\mathrm{D}\big{)}\) are the total energies per atom of the precursor species that compose it; the latter energies are obtained by relaxing the ground-state lattices of these species using the same aforementioned QE parameters as for the compound, as necessary to ensure computational consistency in the evaluation of accurate values for the defects formation energies.
## 3 Results: Property predictions for

substitutional alloying with CGNN

There is a large variety of ways to systematically evaluate the effect of substitutional alloying on the properties of crystals Kaxiras (2003). Here, we focus on two key questions:

1. Data Science of Defects: What is the effect of substituting the same defect in a large variety of systems?
2. Is there qualitative and quantitative effects from atom-substituting various elemental defects in the same host crystalline matrix?

While it has not yet been possible to perform an exhaustive search of the kind, in this work, for the first part, the basin of host systems is represented by the Materials Project crystals dataset introduced in the previous sections Jain et al (2013), while for the latter part, the host systems are pure metallic bulk crystalline supercells of Al, Ni, Mo and Au.

### An elemental defect seeing a wealth of different crystalline environments

First, we focus on the prediction of system properties, where the systematic substitutional defective process is applied using the same replacement atom, on randomly selected sites of crystals that exist in the the Materials Project crystals dataset Jain et al (2013). The aim is to explore how material properties change after single-atom substitution, evaluating deviations from original pure crystalline predictions, and how they depend on the specific substitution. For we consider three elemental cases: Rb, Mn and H, as the key replacement atoms that we will mutually compare. The reason for the choice lies in the drastic elemental differences among their unique characteristics and the assumption
that, on the base of these, we might be able to gain deeper understanding on how the model behaves in the predictions of previously unseen defect-induced changes in the system properties, basing on its learned notion of local environment. Table(3) reports, as an example, the values of the atomic weight, radius and electronic configuration of each considered replacement atom.

In Fig.(2), we display the predictions of MEGNet for the bulk moduli of Rb-defected systems with respect to the original, pure ones. As contained in Table(4), this kind of substitutional defect causes the largest root-mean-squared (RMSD) deviations among the three considered atomic species, with a value of approximately 6.2 GPa. In this section, we only report plots referring to the largest RMSD cases, but analogous plots can be found in the Appendix Sec.(A).

\begin{table}
\begin{tabular}{|c|c|c|c|} \hline  & Weight(u) & Radius (pm) & El. Configuration \\ \hline \hline H & 1.008 & 53 & 1s\({}^{1}\) \\ \hline Mn & 54.938 & 161 & [Ar] 4s\({}^{2}\) 3d\({}^{5}\) \\ \hline Rb & 85.468 & 265 & [Kr] 5s\({}^{1}\) \\ \hline \end{tabular}
\end{table}
Table 3: Differences in the atomic properties of the three elements considered for the single-atom substitutional process. As an example, we report here the atomic weight, radius and electronic configuration.
Also the shear moduli predictions show the largest RMSD for the case of a Rb-defect, with a value of 0.0628 log(GPa)2. The results are shown in Fig.(3) and on Table(5). According to the model, it seems that both K\({}_{\rm VRH}\) and G\({}_{\rm VRH}\) upon substitution of a Rb atom, are showing a tendency to decrease their value, respectively, implying an increase in their compressibility and decrease in hardness. It is worth noticing that even though the only physical atomic

\begin{table}
\begin{tabular}{|l|r|} \hline Defect & RMSD (GPa) \\ \hline \hline Rb & 6.1645 \\ \hline Mn & 3.1274 \\ \hline H & 1.2975 \\ \hline \end{tabular}
\end{table}
Table 4: RMSD (GPa) for the prediction of the bulk modulus (K\({}_{\rm VRH}\)) in Rb, Mn, and H single-atom substitutionally defected systems with respect to the non-defected ones.

Figure 2: Predicted bulk moduli in the Rb-defected MP crystals with respect to their prediction in non-defected ones. An example structure from the MP dataset is shown, CaN2, in which one atom of Ca has been replaced with a Rb atom. Here only the case of Rb-defected systems is shown, due to its largest RMSD among the set of considered defects. Similar plots for the Mn- and H-defected systems are contained in the Supplementary Materials.
feature that the model exploits is the atomic number, the prediction of larger changes in the structural properties for involved defects with larger radii can be regarded as a reasonable one.

Furthermore, in Fig.(4), we show the predicted deviations in the formation energy of the MP-defected samples for the case of elemental H, and also for all three elements in Tab.(6). As it can be seen, the deviations are all in the order of 0.02 eV/atom, and are progressively less prominent going from H, to Rb and Mn.

### A crystalline matrix seeing a wealth of different elemental defects

The case of systematically changing the single-atom substitutional defect species in the same crystalline matrix is complementary to the prior addressed
one. We focus on pure metallic bulk host matrices, namely Al, Ni, Mo and Au, even though the test could be done on any crystalline matrix. In Fig.(5), we show the key aspects of the performed calculation, with the supercell shown on the left (a), and the elements considered on the right (b). The compositional space considered for the single-atom substitution covers almost the entire periodic table, and also comprises the species highlighted as host matrices when they are not selected as such. Even though our focus is on simple systems, the process can be insightful on the capabilities of the model to learn and predict, with minimal input, physio-chemical trends throughout the periodic table.
We show periodic table plots for the properties of interest, K\({}_{\rm VRH}\), G\({}_{\rm VRH}\) and E\({}_{\rm form}\), for the host matrix which displays the largest deviations in the properties, this is the case of Mo, but analogous plots for the other matrices are included in the Supplementary Material (SM). With the help of the visualization style, we characterize the predictions of the effect on the chemical distance between the substitutional defects and the host matrix on the properties of interest, and how it correlates with the well known trends along the periodic table. In Fig.(6), we show the prediction of the host supercell bulk modulus variation when defected with one of the elements from the periodic table. In particular, in the plot we refer to

\[\rm K^{\prime}{}_{\rm VRH}=K_{VRH}^{Mo(X)}-K_{VRH}^{Mo} \tag{5}\]

where K\({}_{\rm VRH}^{Mo(X)}\) is the bulk modulus of the Mo matrix defected by atom X, and K\({}_{\rm VRH}^{Mo}\) its value for pure Mo (labelled with a red flag in the figure).

Even though Sr represents an outstanding outlier in decreasing the bulk modulus of the Mo crystal, at this scale it is still possible to appreciate how the

Figure 5: Periodic table plot (b) explaining the process of selecting a set of host matrices (in light blue) and substitutional defects (in green). To a given a selected host matrix, the non-selected ones represent defects too. The host matrices are \(3\times 3\times 3\) supercells of the highlighted species (a).
variation happens along the periods: for the 3d, 4d and 5d transition metals from 3rd to 12th group, the defect-induced variations are mainly small, while a tendency to increase is observed, in modulus, in the post-transition metals and, remarkably, in the alkali and alkaline-earth. This behaviour can be interpreted in terms of the well known variations of the bulk modulus along the periodic table of elements, which seem correlated with respect to the defect-induced ones contained in this plot. The effect of substitutional alloying species which, in their bulk and standard temperature and pressure (STP) conditions, show a lower bulk modulus, is, as a tendency, to lower the bulk modulus of their host system.

Due to the strong effect induced by the Sr defect, we further focus on transition metals to evaluate how defect-induced variations fluctuate in the

Figure 6: Predicted bulk modulus variation (K’\({}_{\mathrm{VRH}}\)) for a single-atom substitutionally defected Mo supercell with respect to the undefected one, for each possible defect atomic specie from the provided periodic table. Similar plots for Al, Au and Ni supercells are provided in the Supplementary Materials. The red flag highlights the zero relative difference, meaning the pure Mo matrix selection.
compositional vicinity of the host matrix. Fig.(7) shows the results along the list of 3d, 4d and 5d transition elements, where the variations are harder to distinguish from the prior periodic table plot. We can recognize a decreasing trend that is supportive of our interpretation, but in the compositional vicinity of the host matrix the fluctuations in the defect-induced variations are comparable with the MAE on bulk moduli prediction, just as the differences between the curves.

However, in view of a trend analysis of the variations along the periodic table, Alkali (like K), Alkaline-earth (like Sr and Ba) and post-transition metals (like Se and Te), can support the given interpretation in terms of correlation with the bulk modulus of the impurity, given their associated fluctuations are much larger than the model predictions errors.

Figure 7: Predicted bulk modulus variation (K’\({}_{\rm VRH}\)) for a single-atom substitutionally defected Mo supercell with respect to the undefected one, along the 3d, 4d and 5d series of the periodic table. The black dashed line highlights the pure Mo case.
We may perform analogous investigations for the shear modulus of a pure host matrix, like Mo, and how it gets influenced by single-atom substitutional alloying, spanning over almost the entire periodic table, as shown in Fig.(8). We follow the same protocol and consider the variation \(\rm G^{`VRH}=G_{VRH}^{Mo(X)}-G_{VRH}^{Mo}\). In this case, the colormap reveals Si as an outlier towards hardening of the Mo matrix. This is a feature that cannot be explained by a correlation in defect- and defect-induced properties trends, since the predicted value is larger that the model prediction error on the shear modulus. We believe that the power of this investigation method is in paving the way to an efficient exploration of substitutional alloying, with the twofold possibility of looking for comparable performances (discovery of alternatives) or outstanding ones (discovery of exceptionals). Similarly to the investigation of the bulk modulus, the Alkali and Alkaline-earth metals like K, Rb, Sr and Ba are among the substitutional species providing the largest decreases in the shear modulus.

Overall, the effects and fluctuations caused by the substitutional defects on the defected host can always be highlighted, but it is not the aim of this work to find an exhaustive explanation for the existing predictions: The reasons for such trends may be due to any of the input parameters, such as the atomic number and bond lengths, or an abstract notion of local environment which is good enough to show reasonable correlations with existing alternative descriptors (i.e. atomic properties). The variations along the 3d, 4d and 5d transition metals are, for most of the cases, below the model MAE for \(\rm G_{VRH}\) predictions, therefore no meaningful extrapolation is possible, but we report the plot in Fig.(B9) of the Appendix.
For what concerns the formation energy, by definition of the latter, the results in Fig.(9) can be interpreted as the gain or loss in stability after the single-atom substitutional alloying has taken place. Fluorine represents a strong outlier, raising the formation energy by 0.14 eV/atom with respect to the pure Mo matrix. As shown in Fig.(10), it is evident that there is a trend that spans over the periods, and an overall interesting correlation could be found with the known trends for the atomic electronegativity along the periodic table, suggesting that the higher the latter for the substitutional defect in the Mo matrix, the higher is the resulting formation energy.

Figure 8: Predicted shear modulus variation (G’\({}_{\rm VRH}\)) for a single-atom substitutionally defected Mo supercell with respect to the undefected one, for each possible defect atomic specie from the provided periodic table. Similar plots for Al, Au and Ni supercells are provided in the Supplementary Materials. The red flag highlights the zero relative difference, meaning the pure Mo matrix selection.
Figure 9: Predicted formation energy variation (E’\({}_{\rm form}\)) for a single-atom substitutionally defected Mo supercell with respect to the undefected one, for each possible defect atomic specie from the provided periodic table. Similar plots for Al, Au and Ni supercells are provided in the Supplementary Materials. The red flag highlights the zero relative difference, meaning the pure Mo matrix selection.
### Property prediction - DFT validation

The validation of artificial intelligence methods predictions is a crucial step in order to quantify their quality. Even though the proposed graph network based method has been already validated for its predictions in bulk crystals, this work aims at testing its capabilities in the presence of single-atom substitutional defects. As explained in Sec.2.2, we compare the model predictions of K\({}_{\mathrm{VRH}}\), G\({}_{\mathrm{VRH}}\) and E\({}_{\mathrm{form}}\) with the DFT based ones, obtained with the THERMO_PW driver of QE. Table(7) shows the results for the validation on the properties in the case of an Al matrix and single-atom substitutional defects, including B, C, I, Ni and Zr.

Comparing the MAE values of our DFT calculations, with the ones of the model for non-defected systems in Table(2), we find good performances of the

Figure 10: Predicted formation energy variation (E’\({}_{\mathrm{form}}\)) for a single-atom substitutionally defected Mo supercell with respect to the undefected one, along the 3d, 4d and 5d series of the periodic table. The black dashed line highlights the pure Mo case.
model with respect to K\({}_{\rm VRH}\) and G\({}_{\rm VRH}\) predictions, but large errors when it comes to the formation energies: the first may be regarded as a success, given that the model is predicting properties of a new class of systems, and given the computational limitations; the second one is a negative signature, even though the (i) defect dependent nature of the error order of magnitude opens up a deeper window of investigation on its reasons, and (ii) the validation set for defected systems is extremely small compared to the undefected MP dataset, on which initial MAEs have been evaluated.

### Size effects

In substitutional alloying, the defect atomic specie is usually present in a dilute concentration, in the range of 0.1-10%. For this reason, it is of interest the study of how the property predictions vary with defect concentration, that we propose in the saturation plots of Fig.(11), Fig.(12) and Fig.(13). Our example system follows the choice of the Mo host matrix, with the H, Mn and Rb substitutional defects we focused on, respectively in the second and first part of the previously reported results.
asymptotic over-dilute level (\(<0.1\%\)). Moreover, while the defect formation energy, as expected, shows a common descent to zero-level for all the defect species, the predicted structural properties seem to be sensitive on them, with the Rb-defected Mo crystal conserving a visible difference in the property value even at the limit of 0.1% concentration, both for K\({}_{\rm VRH}\) and G\({}_{\rm VRH}\). Even though it is out of the present work aims to validate the asymptotic-size behaviour of the predictions with accurate but expensive DFT calculations, we believe these results to stand in favour of the positive model understanding of a crystalline defected system environment.

Figure 11: Saturation plot of a Mo host matrix K\({}_{\rm VRH}\) when substitutionally defected with H, Mn or Rb, for different supercell sizes.
Figure 12: Saturation plot of a Mo host matrix \(\rm G_{VRH}\) when substitutionally defected with H, Mn or Rb, for different supercell sizes.

Figure 13: Saturation plot of a Mo host matrix \(\rm E_{form}\) when substitutionally defected with H, Mn or Rb, for different supercell sizes.
Let us continue along the previously paved path of analysis which also looks into the effects of systematically changing the substitutional atomic specie among the large variety contained in the periodic table, as previously done in Fig.(6), with the only exception of selecting the smallest and largest supercell sizes from the previous saturation plots of Fig.(11), respectively \(2\times 2\times 2\) (1% concentration, 16 atoms cell) and \(8\times 8\times 8\) (10% concentration, 1024 atoms cell). In Fig.(14), focusing on the first row of plots which deal with the K'\({}_{\rm vRH}\) variation, and comparing with the previously investigated supercell case of size \(3\times 3\times 3\) (2% concentration) of Fig.(6), one can notice that the scale of property variations changes accordingly: smaller defect concentrations lead to smaller defect-induced effects, and viceversa, as expected. In particular, in the largest supercell case the induced variations range is reduced to a 3% of the \(3\times 3\times 3\) system's one. However, the composition map outliers are left unchanged. The shear modulus variations, in the second row of tables, sees the emergence of new outliers in the chemical neighbourhoods of the previously obtained ones in \(3\times 3\times 3\) supercells, while the formation energy variations undergo both a change in scale and a complete change in the map outliers. Following the brief assessment of predictions quality performed for each of the interesting properties in the previous section, we expect the formation energy variations to suffer by non-negligible fluctuations, leading to the possible need of further validation. However, the main aim of the present discussion is to underline the power of the overall approach in highlighting the path towards composition search in substitutional alloying, which can effectively drive to specific desired emerging properties.
## 4 Conclusions

In this work, we utilized a convolutional graph-neural network, based on the MEGNet architecture Chen et al (2022), in order to attempt the design of novel alloys. Alloying involves substitutional and interstitial alloying at relatively low concentrations, thus, single-defect properties shall be informative on overall designing capabilities and guidance. For this purpose, we utilized the MP database, and we focused on the prediction of the properties of single-atom

Figure 14: Predicted variations K’\({}_{\mathrm{VRH}}\) (first row plots), G’\({}_{\mathrm{VRH}}\) (second row plots), E’\({}_{\mathrm{form}}\) (third row plots) for a single-atom substitutionally defected Mo supercell with respect to the undefected one, for each possible defect atomic specie from the provided periodic table. In plots of column a) a \(2\times 2\times 2\) supercell is considered, while in column b) a \(8\times 8\times 8\) supercell.
substitutionally defected bulk crystals,in the context of both (i) systematic substitution with a specific set of species in a wide variety of crystals from the entire Materials Project dataset and (ii) systematic substitution with a variety of atomic species in a specific set of pure bulk crystals. We also validated some of the results with our own DFT calculations. We believe that such approaches might provide novel insights into alloy design, especially if predictions include extended lattice defects such as dislocations and/or grain boundaries.

Acknowledgments.We would like to thank J. Llorca for insights and fruitful discussions and suggestions. This research was funded in part by the European Union Horizon 2020 research and innovation program under Grant Agreement No. 857470 and from the European Regional Development Fund via Foundation for Polish Science International Research Agenda PLUS program Grant No. MAB PLUS/2018/8.

### Data availability

Requests for additional data and availability can be directed to the provided lead contact. The employed DFT data, ML model and database can be found on the GitHub: [https://github.com/danielcieslinski/suballoy](https://github.com/danielcieslinski/suballoy)

## Declarations

Declarations
