# Prediction and characterization of chemically complex $\sigma$ – phase intermetallics with graph neural network

Wenhao Zhang, Mariano Forti, Runan Xie, Céline Barreteau, Taichi Abe,

Jean-Marc Joubert, Thomas Hammerschmidt, Jean-Claude Crivello

► To cite this version:

Wenhao Zhang, Mariano Forti, Runan Xie, Céline Barreteau, Taichi Abe, et al. Prediction and characterization of chemically complex $\sigma$ – phase intermetallics with graph neural network. Acta Materialia, 2025, 300, pp.121427. 10.1016/j.actamat.2025.121427. hal-05262526

HAL Id: hal-05262526

https://hal.science/hal-05262526v1

Submitted on 16 Sep 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
Graphical Abstract

Prediction and Characterization of Chemically Complex σ-Phase Intermetallics with Graph Neural Network

Wenhao Zhang, Mariano Forti, Runan Xie, Céline Barreteau, Abe Taichi, Jean-Marc Joubert, Thomas Hammerschmidt, Jean-Claude Crivello

DFT database of
intermetallic

<figure><img src="image_2.png" /></figure>

Machine learning
$H_f$ predictions

<figure><img src="image_3.png" /></figure>

<figure><img src="image_4.png" /></figure>

Stability and trends
in σ-phase

<figure><img src="image_5.png" /></figure>

<figure><img src="image_6.png" /></figure>
# Highlights

## Prediction and Characterization of Chemically Complex σ-Phase Intermetallics with Graph Neural Network

Wenhao Zhang, Mariano Forti, Runan Xie, Céline Barreteau, Abe Taichi, Jean-Marc Joubert, Thomas Hammerschmidt, Jean-Claude Crivello

* The formation enthalpy of end-members for the intermetallic σ-phases within a 19-element chemical space, up to quinary systems, is predicted using a graph neural network (GNN) model with a root mean squared error (RMSE) of 10 meV/atom.

* The Gibbs energy and site preferences of σ-phases in all possible ternary, quaternary, and quinary systems are calculated using the Compound Energy Formalism (CEF).

* Bonding analysis are performed to understand the stability and chemical ordering in the σ-structure.
Prediction and Characterization of Chemically Complex
σ-Phase Intermetallics with Graph Neural Network

Wenhao Zhang<sup>a,b</sup>, Mariano Forti<sup>c</sup>, Runan Xie<sup>d</sup>, Céline Barreteau<sup>d</sup>, Abe Taichi<sup>b</sup>, Jean-Marc Joubert<sup>d</sup>, Thomas Hammerschmidt<sup>c</sup>, Jean-Claude Crivello<sup>d,a</sup>

<sup>a</sup> CNRS-Saint-Gobain-NIMS, IRL 3629, Laboratory for Innovative Key Materials and Structures (LINK), 1-1 Namiki, Tsukuba, 305-0044, Ibaraki, Japan

<sup>b</sup> Research Center for Structural Materials, National Institute for Materials Science, 1-2-1 Sengen, Tsukuba, 305-0047, Ibaraki, Japan

<sup>c</sup> Interdisciplinary Centre for Advanced Materials Simulation (ICAMS), Ruhr-University, Universitätsstrasse 150, Bochum, 44801, Bochum, Germany

<sup>d</sup> Univ Paris Est Creteil, CNRS, ICMPE, UMR 7182, 2 rue Henri Dunant, 94320, Thiais, France

Abstract

Determining the thermodynamic stability of intermetallic phases in multi-
component alloys from Compound Energy Formalism (CEF), requires the
formation enthalpy of up to millions of end-member configurations. Here, a
graph neural network (GNN) model, trained from high-throughput DFT cal-
culations, was used to construct such a database of intermetallic σ–phase for
19 metallic elements with a root mean squared error (RMSE) of 10 meV/atom.
We show that a good prediction accuracy of GNN, even without prior knowl-
edge of the relaxed structure, can be achieved by focusing on a single phase
and using a regression model for structure prediction. The resulting 19⁵
formation enthalpy data, corresponding to all end-member configurations
up to quinaries, enables the characterization of the stability and chemical
trends of σ–phase both at 0 K (configurations on the convex hull) and at
high temperature (with disorder using CEF, also up to to quinary systems).
We quantified the important role of entropy in stabilizing the σ–phase in
multicomponent systems and investigated the relationship between chemical

Email addresses: Abe.Taichi@nims.go.jp (Abe Taichi),
thomas.hammerschmidt@rub.de (Thomas Hammerschmidt),
jean-claude.crivello@cnrs.fr (Jean-Claude Crivello)
composition, chemical ordering and enthalpy of formation in the σ−phase. The chemical ordering further is discussed in detail in terms of site and pair preference of different elements. Through electronic structure based bonding analysis, we demonstrate the importance of valence electron concentration (VEC) in explaining and rationalizing the obtained theoretical results.

*Keywords:* Intermetallic, σ−phase, Machine learning, Thermodynamic modeling

**1. Introduction**

The stability of solid solution phases and intermetallic phases across var-
ious compositions and temperatures is one of the keys in alloy design for
diverse applications.[1] As an example, in Ni-based superalloys, avoiding the
formation of brittle intermetallic phases is essential to maintain desirable me-
chanical properties for high-temperature applications.[2] On the other hand,
the intentional inclusion of a secondary phase can also sometimes be benefi-
cial to utilize effects such as precipitation strengthening, which was realized
recently in the case of σ and µ precipitates in high-strength CoCrFeNiMo
high-entropy alloys (HEAs)[3].

The σ−phase is an important intermetallic phase that is observed in vari-
ous binary and ternary systems, as well as in HEAs.[4, 5] It is a topologically
close-packed (TCP) phase with Frank-Kasper polyhedra as basic structural
unit. Its tetragonal crystal structure and the coordination environments of
atoms in the σ−phase are shown in Figure 1. The coordination number
(CN) of the five non-equivalent crystallographic sites in the σ−structures
ranges from 12 to 15, allowing it to accommodate atoms of different sizes.
An important feature of the σ−phases is its wide homogeneity range in the
phase diagram and the associated substitutional disordering. Compared to
other intermetallic phases that tend to be chemically ordered, disordering
can stabilize the σ−phase with increasing temperature due to the large con-
figurational entropy contribution to the Gibbs energy.[6, 7, 8] Therefore, the
σ−phase can be an important competing phase to the high-entropy solid solu-
tion phases.[9, 10, 8] Indeed, a recent survey showed that the σ−phase, along
with the Laves phase, is the most observed intermetallic phase in HEAs.[11]
However, experimental data of σ−phase in multicomponent systems remain
very limited. Therefore, it is important to quantitatively characterize the
phase stability of the σ−phase in multicomponent alloy systems through a
computational approach.

The stability of different phases is given by their relative Gibbs energy at a given temperature and pressure. In the Calphad (Calculation of phase diagram) methodology, the Gibbs energy can be modeled based on thermo-dynamic databases (TDBs) extracted and assessed either from experimental results such as heat capacities, phase equilibria, or from Density Functional Theory (DFT) calculations.[12, 13, 14] Compared to experimental procedures that require many resources, computational methods can be cheaper, and their agreements with experimental findings are often quantitative.[15] Nonetheless, constructing TDB for multicomponent systems solely from first-principles calculations remains challenging due to the combinatorial explosion in the number of energy evaluations required. This challenge arises especially when describing intermetallic systems using sublattices within the Compound Energy Formalism (CEF), where the computational cost scales as $(N_{ele})^M$, where $M$ is the number of sublattices in the phase-model considered and $N_{ele}$ the number of different elements in the database.[16]

Supervised Machine learning (ML) can provide a surrogate model to evaluate the energy of different phases to construct TDBs from a small training database.[17] Since the development of structural-based ML methods, especially the message-passing neural networks on crystal graphs, the error of ML

<figure><img src="image_7.png" /><figcaption>Figure 1: The crystal structure of σ-phase. The space group is *P*4<sub>2</sub>/mnm (136) and the five Wyckoff positions are indicated by different colors. The coordination numbers for each polyhedral are 12 for 2a, 15 for 4f, 14 for 8i1, 12 for 8i2 and 14 for 8j.</figcaption></figure>
prediction on material properties has reduced significantly.[18, 19] According
to recent Matbench results, state-of-the-art graph neural networks (GNNs)
can achieve a mean absolute error (MAE) in energy prediction of around
20 meV/atom, compared to first-principles calculation results.[20] However,
there are still challenges in using these structure-based ML methods for
practical applications. First of all, GNN needs geometric information as
input, and it was shown that well-relaxed crystal structures, typically un-
known at the application stage, are necessary inputs to the GNN models to
achieve the reported high level of accuracy in formation enthalpy.[21, 22, 23]
Secondly, the prediction accuracy of these state-of-the-art models becomes
poor when used on out-of-distribution predictions, suggesting that specif-
ically sampled datasets are necessary to ensure accuracy among different
prediction tasks.[24]

In this work, we first defined a target error of ML prediction and used a
structure-based GNN, combined with regression models for structure param-
eters prediction, to evaluate the formation enthalpy of end-member configura-
tions of σ−phase. A total of 19 elements, including 18 transition metals and
one sp−metal (Al), are considered in order to cover a large chemical space.
Then, prediction results were used to characterize the stability and crystal
chemistry of σ−phase for all systems up to quinary configurations—limited
by the choice of five sublattices that correspond to five different Wyckoff posi-
tions in the σ−phase unit cell. The precipitation of the σ−phase is estimated
in comparison with solid solution phases. In addition, we performed bonding
analysis based on electronic structure calculations to help to rationalize the
chemical trends found from the DFT-ML data.

Related works that used ML methods to study multicomponent alloys in-
clude the prediction of the phases formation (typically using ML classifiers)
from an input chemical composition.[25, 26, 27] However, these predictions
cannot give information about phase fractions or the temperature range in
which secondary phases become stable. On the other hand, there are works
that combine accurate ML energy models with Monte Carlo or molecular
dynamics to perform the atomistic simulations for multicomponent solid so-
lution alloys.[28] However, such applications are usually restricted to a sin-
gle system. This work, focused on a single phase across a large chemical
space, demonstrates that a ML approach can be used to construct accurate,
large scale, thermodynamic databases, which are crucial for thermodynamic
simulation of phase diagrams and phase prediction in multicomponent alloy
design.
## 2. Computational Objectives

### 2.1. Gibbs Energy

Our computational modeling will be based on estimating thermodynamic stability through the Gibbs energy. Since the intermetallic σ-phase features a wide homogeneity range with partial occupancy of atoms on all crystallographic sites, its Gibbs energy as a function of temperature and composition $G(\mathbf{x}, T)$ has to be obtained by optimizing $G$ with respect to site fraction $\mathbf{y}$ given by $y_i^s = N_i^s / \sum_j N_j^s$, where $N_i^s$ is the number of atoms of chemical species $i$ on sublattice $s$. [12, 6] Assuming for simplicity that the vibrational and electronic contributions to the Gibbs energy are similar between the σ-phase and elemental reference structures, the Gibbs energy can be expressed using CEF as:

$$ G(\mathbf{x}, T) = H(\mathbf{x}, T) - TS(\mathbf{x}, T) \\ = \min_{y_i^{2a} \dots y_m^{8j}} \left[ \sum_{i,j,k,l,m} \left[ y_i^{2a} y_j^{4f} y_k^{8i1} y_l^{8i2} y_m^{8j} \right] H_{ijklm}^0 + k_B T \sum_s \frac{m_s}{N} \sum_i y_i^s \ln y_i^s \right] \quad (1) $$

where $S$ is the configurational entropy, $\mathbf{x}$ represents the molar fraction of all the components, $k_B$ is the Boltzmann constant, $m_s$ and $N$ is the multiplicity of site $s$ and the total number of atoms in the unit cell, respectively. $H_{ijklm}^0$ is the formation enthalpy per atom of a structural configuration (end-members) where atoms of species $i, j, k, l, m$ occupy the five sublattices (2a, 4f, 8i1, 8i2, 8j) respectively in the σ-phase unit cell (Figure 1). To evaluate the Gibbs energy in the intermetallic σ-phase, it is therefore required to know the formation enthalpy of all the end-member configurations $H_{ijklm}^0$. 19 elements that are considered in this study give $19^5 \approx 2.5$ million σ-structural configurations (for five sublattices) for which formation enthalpy needs to be known. They were obtained in this work using supervised ML surrogate model.

### 2.2. Target Accuracy

A target accuracy on the predicted formation enthalpy needs to be first identified to ensure that the results are practically useful. This is especially important here since within CEF, crystallographic properties such as site occupancies are also determined from the enthalpy of end-members. These calculated crystallographic properties can later be compared to experimental
observations as a way to assess the accuracy of computational predictions. Previously, Palumbo et al. showed that in the binary σ − system Re–W, an error of 0.5 kJ/mol ( ≈ 5 meV/atom) on end-member energy has almost no effect on calculated thermodynamic and crystallographic properties ( G , y s $_{i}$, and C$_{p}$ ) from CEF, but an error of 2.0 kJ/mol ( ≈ 20 meV/atom) dramatically changes the resulting quantities.[29] Based on this result, it seems suitable to set the target error to be around 1.0 kJ/mol ( ≈ 10 meV/atom) for ML pre- diction, which offers good prediction of thermodynamic properties without being too stringent. Errors are reported mostly as root mean squared error (RMSE) calculated from prediction y and target t as:

$$
\text{RMSE} = \sqrt{\frac{1}{N} \sum_{i} (y_{i} - t_{i})^{2}}
$$

where N is the total number of samples. The mean absolute error (MAE) is typically smaller than RMSE and is calculated as:

$$
\text{MAE} = \frac{1}{N} \sum_{i} |y_{i} - t_{i}|
$$

3. Prediction of Formation Enthalpies

3.1. Database Preparation

Complex intermetallic phases such as the σ − phase are underrepresented in online materials databases, such as the Alexandria Materials Database due to their large unit cell size.[30] Therefore, we constructed a ML training database for the σ − phase. The formation enthalpies of all training structures used in this work were prepared by first-principle calculations using the Pro- jector Augmented Wave pseudopotential method, with PBE functional, as implemented in Vienna Ab initio Simulation Package (VASP).[31, 32] The details of the calculation parameters are provided in Appendix A. Figure 2 shows the workflow to obtain the Gibbs energy for multicomponent σ − phase from the initial data.

The σ − phase configurations up to quaternary systems are systematically sampled to build the training dataset, as shown in Figure 3a. Out of the 19 elements, all the ordered configurations of unary (19) and binary (5130) con- figurations are included in the training datasets, while 4845 ternary and 7752
<figure><img src="image_8.png" /><figcaption>Figure 2: Overall workflow employed in this work, divided into three steps: training of the ML models, prediction of formation enthalpies for end-members, and calculating Gibbs energy of multicomponent σ-phases.</figcaption></figure>
quaternary configurations are equally drawn from each possible chemical sys-
tem. Quinary configurations are not included in the training set since it is dif-
ficult to evenly sample all quinary systems. In total, 17,745 σ−phase config-
urations (0.7% of all possible configurations) with their formation enthalpies
at the relaxed crystal structure are included in the training dataset. Uniform
Manifold Approximation and Projection (UMAP) technique is used to cre-
ate a two dimensional representation, shown in Figure 3b, of all the sampled
σ−configurations (represented by one-hot encoding), and we indeed find an
even distribution of these configurations in the chemical space.[33] For details
of sampling the σ−configurations, see Appendix B. Additional entries were
added successively to augment the relaxed σ−configurations which improve
model accuracy. The final training data includes 108,408 entries enriched
with diverse structural environments that include: 1) 17,745 × 3 = 53,235
DFT calculations of distorted structures (without relaxation) of the same
17,745 σ−configurations, and 2) 37,428 configurations of other common in-
termetallic phases in various structures: Laves, χ and A15 phases (also see
Appendix B).

The errors of formation enthalpy prediction were tested using a test set consisting of 3000 σ−phase configurations from the 19 elements including 1000 ternary, 1000 quaternary and 1000 quinary configurations. These configurations in the test set are randomly sampled from all possible ones that are not included in the training dataset. The even distribution of elements in the test set is illustrated by Figure S4.

3.2. *Formation Enthalpy Model*

First, we prepare the formation enthalpy model which predicts formation
enthalpies from the input crystal structures (how to provide a good input
structure is described in the following section). To obtain high accuracy in
energy predictions, the ALIGNN (Atomistic Line Graph Neural Network)
model is used, which utilizes two graphs to represent a crystal structure
with explicit input of chemical features of the atoms, bond length and bond
angles.[34] ALIGNN is one of the top performing models in the Matbench
formation energy prediction tasks with a small MAE of 21 meV/atom.[35] We
trained ALIGNN model from scratch using the prepared training set, in which
85% were used to train the ALIGNN model while the remaining 15% were
used as validation to select the best performing model after each epoch. The
ALIGNN model was trained with hyperparameters as originally published
<figure><img src="image_9.png" /><figcaption>Figure 3: Statistics of training datasets a) in terms of structure, number of components in the configuration and elements. Inset shows the same figure but in a different scale. b) UMAP representation of the distribution of σ-configurations used in the training set.</figcaption></figure>

but with slightly modified atomic features (see Supplementary Section 2 for detail on model training).

Figure 4a shows the error of enthalpy prediction by the ALIGNN model as a function of increasing number of entries in the training dataset. In these tests, ALIGNN models were provided with the DFT-relaxed structures. Starting with the 17,745 relaxed σ-configurations, a very good test RMSE of about 11 meV/atom (MAE = 7.5 meV/atom) can already be obtained. Since the prediction task is on a single structure type, the level of accuracy achieved here is much better than the reported accuracy in the benchmark, which is trained and tested on various structures from the materials project.[20] With the inclusion of more and more non-relaxed σ-structures and other inter-metallic phases in the training set, the prediction RMSE decreased further to around 7.5 meV/atom (MAE = 5.4 meV/atom) when all training structures are used. According to Figure 4a, the prediction error is likely to reduce further if more training data is available, since there is no observable saturation of prediction accuracy. However, we did not increase the size of the training dataset further since the prediction error is sufficiently small.

Since TCP phases share common Frank-Kasper polyhedra in their crystal structure, it could be wondered whether ML model trained on one structure can be transferrable to another structure. We checked the transferrability
<figure><img src="image_10.png" /><figcaption>Figure 4: a) Prediction error (RMSE) on the σ test set with the ALIGNN model trained on different training sets of various sizes. The three bars at each training set combination indicate the RMSE for ternary, quaternary and quinary systems from left to right. b) Prediction error on other intermetallic phases using the ALIGNN model trained from the σ training sets.</figcaption></figure>

of the ALIGNN model by using the ALIGNN model trained only on the σ-phase to predict the formation energy of other intermetallic phases (A15, Laves and χ) using relaxed structures. Not unexpected, we found that the structure-based ALIGNN model has rather poor extrapolability to unseen structure types. The RMSE plotted in Figure 4b for other intermetallic phases is significantly larger, ranging from 71 meV/atom to 240 meV/atom, compared to about 11 meV/atom in the σ test sets when trained only with relaxed σ-configurations. The differences between the coordination environments present in other intermetallic structures and the training σ-phase may explain these large discrepancies (e.g., there is a large CN16 site in the Laves structure that does not exist in the σ-structure). Interestingly, training with additional distorted σ-phase structures can reduce the prediction error on unseen structure types quite dramatically (as much as 30% for the Laves phase), likely due to the added variety of structural environments with changed bond lengths and angles. Still, these errors are too large to be useful, and it can be concluded that with limited extrapolability, model prediction of GNN could only be reliable if sufficient training data of the corresponding structure type are available.
<figure><img src="image_11.png" /><figcaption>Figure 5: Distribution of error in bond lengths (a), bond angle (b) and energy (c) compared to the result obtained from DFT structure relaxation in 1500 test samples. The root mean squared error (RMSE) in bond length and angle is calculated for all bonds and angles in each of the test samples. The box and whiskers in the box plot indicate data ranges of [25%, 75%] and [5%, 95%], respectively. d) Scatter plot of the prediction, including both the prediction on predicted structures and the prediction on DFT-relaxed structures in the training set of 3000 configurations. Energy differences are shown as a boxplot in the inset. e) Violin plot of prediction error for the seven ternary σ–systems tested. Each ternary system contains 150 ternary configurations. The upper number near the violin plot indicates the root mean squared error of ALIGNN prediction on predicted structures.</figcaption></figure>

3.3. *Estimating Relaxed Structure*

While the trained ALIGNN surrogate model for formation energy has
excellent accuracy, it remains insufficient by itself because relaxed structures
of different configurations are in general unknown before prediction but are
needed as an input. Previous studies showed that prediction with unre-
laxed guessed structures leads to a significant degradation of accuracy.[21, 36]
Several methods were proposed in the literature to overcome this limita-
tion: (i) training the model using only structural prototypes or information
from crystallographic sites without precise atomic positions[37, 38, 39], (ii)
data-augmentation to map unrelaxed structures to the energy of the relaxed
ones.[40], and (iii) optimization of structural parameters based on the trained
surrogate energy model.[41] Nonetheless, the reported prediction errors of
these proposed methods were still substantially larger than the error of a
well-trained model when provided the accurate relaxed structures.

Within a single structure type (σ), we found that direct predictions of
individual structural parameters (volume, c/a ratio and degrees of freedom
in Wyckoff positions) are sufficiently accurate such that the ALIGNN model,
taking these predicted structures, can reliably predict the ground state en-
ergy. Using the structural parameters of well-relaxed σ-configurations as
the training dataset, predictions were done by multi-layer perceptron (MLP)
regressors based on the one-hot encoding (OHE) representation of elements
in the configurations. An independent test set with DFT structure relax-
ation of another 1500 configurations (500 ternary, 500 quaternary and 500
quinary configurations, shown in Figure 3) was used to test the accuracy of
the structural parameter predictions. The mean absolute percentage errors
for these structural parameter predictions ranged from 0.2% to 1.4%, and the
test results are shown, together with the statistics of the relaxed structural
parameters, in Figures S5 and S6.

In closer relationship to the ALIGNN model that take bond length and
bond angle as inputs, the bond length and bond angle differences between
the DFT-relaxed structures and the OHE-predicted structures, within the
unit cell, were directly compared for each of the 1500 test configurations.
The results are shown in Figure 5a and b. We found that the average RMSE
in bond length and bond angles is 0.03 Å and 0.9 degrees, respectively. Two
alternative approaches were also tried out: i) a simple unit cell volume scal-
ing while fixing the c/a ratio and atomic coordinates (atomic volume used for
each elements are shown in Figure S7) and ii) black box Bayesian structural
optimization following Zuo et al.[41] (Details in Supplementary Section 3).
Both alternatives showed larger structural prediction errors. In Figure 5c,
we show that the error in predicting relaxed structures plays a determin-
ing factor in the accuracy of predicted formation enthalpies. The RMSE
of formation enthalpies, predicted by ALIGNN, using structures scaled by
volume is 54 meV/atom, which drops to 23 meV/atom with structures from
Bayesian optimization and around 10 meV/atom using structures from OHE-
MLP predictions. As a result, we used MLP in this work to predict structure
parameters for every end-member configuration, of which the relaxed structure is not known, as structure inputs to the ALIGNN model.

### 3.4. *Model Evaluation*

We used the above workflow to predict the formation enthalpies in the larger test set with 3000 configurations, this time using predicted structures from the aforementioned MLP regressors. Figure 4a shows the learning curve with increasing training data size while the final prediction results are shown in Figure 5d. The RMSE for the final prediction of the formation enthalpy is 10.2 meV/atom (MAE = 6.7 meV/atom) without any prior information on relaxed structure, which is only slightly larger than the prediction error on the DFT-relaxed structure and reaches the desired accuracy. Compared to previous studies that use machine learning to predict formation enthalpy of compounds without prior knowledge of crystal structure (MAE ≈ 30-100 meV/atom), the error obtained in this work is much smaller because we focus on a single phase.[37, 38, 39, 40, 41] The results here also improve significantly upon our previous work that used an MLP regressor to predict the formation enthalpies for the σ-phase (RMSE = 33.6 meV/atom) as well as other related works, due to the use of a large training database and a more complex machine learning model.[17, 42]

In actual applications, all end-member configurations in a given system need to be accessed. To check for any systematic errors in different chemical systems (e.g., systematic underestimation or overestimation in some chemical systems, or error associated with certain elements), we performed additional tests by calculating all end-member σ-configurations from DFT in seven different ternary systems and comparing them with the predictions. Each of these ternary systems contains 150 configurations, and the systems were sampled by simply shuffling all 19 elements under consideration and dividing them into triplets. The result is shown in Figure 5e. In four of the seven ternary systems (Co-Fe-Ru, Os-Re-Ta, Nb-Ti-W, Nb-Pt-W), the predictions are found to be highly accurate with an RMSE error of only around 5 meV/atom. However, the errors are larger for Cr-Cu-Mo, Mn-V-Zr and Al-Ni-Tc due to the overestimation of corresponding formation enthalpy. In these systems, it can be noticed that the error is considerably reduced if relaxed structures are used. This suggests that part of the larger prediction error, especially in Mn-V-Zr, is related to the accuracy in the prediction of relaxed atomic positions. Indeed, the atomic size of Zr is significantly larger than that of Mn and V and may lead to error-prone predicted structures.
Cu and Al may also be considered outliers with the former having filled d-bands and the latter not being a transition metal. In fact, Zr and Cu are traditionally considered not to be σ-phase forming. Thus, from these error characterizations, it can be concluded that while the prediction error is small on average in the whole test set, there can be variations in accuracy depending on the chemical elements. In Table S2, we examined errors related to different elements in the test set by either only including or excluding configurations that contain one type of element. Al, Mn, Cu, Zr and Pt are the elements associated with larger errors on average. In particular, the RMSE is reduced to 7.8 meV/atom when all configurations that contain Zr are excluded from the test set. Additional discussions regarding element-specific ML errors and how they are related to chemical and structural characteristics of certain elements, such as Al and Zr, are given in Section 4 of the Supplementary Material. Based on these observations, it can be expected that active learning—by actively sampling regions in the configuration space where errors are larger—can be useful to achieve consistent error in the entire prediction space, in future works.

In conclusion, an ML process combining both structure prediction and a GNN model is found to give accurate predictions of the formation enthalpy in multicomponent intermetallic σ-phase. With it, formation enthalpy of all σ-phase end-member configurations with 19 elements were evaluated. The prediction results of formation enthalpy are available electronically as a part of the supplementary information.

# 4. Results on the Multicomponent σ-phase

## 4.1. Configurations on the convex hull

While many works used computational methods to characterize σ-phase in different systems, a systematic study of σ-phase in a wide range of systems is not yet available due to computational cost.[6, 29, 8] Each of the ternary (quaternary/quinary) systems requires the knowledge of the formation enthalpy of $3^5 = 243$ ($4^5 = 1024/5^5 = 3125$) end-member configurations—not to mention a chemical space of 19 elements. However, using ML prediction, the formation enthalpies of all 19⁵ end-member configurations can be predicted within a well-defined accuracy (around 10 meV/atom) in this work. From these data, we can develop overall understanding of stability and chemical trends for σ-phase in multicomponent systems.
<figure><img src="image_12.png" /><figcaption>Figure 6: a) The Pettifor chemical scale used to order elements in the following figures, b) the frequency of elements to occupy different sites in σ-configurations on the convex hull and c) a boxplot of the distribution of the formation enthalpies of σ-configurations on the convex hull that are associated with each element.</figcaption></figure>
First, all end-member σ-configurations on the T = 0 K convex hull (up to quinary) were characterized. In Figure 6, the site preferences and the distribution of formation enthalpies associated with each elements are plotted. For each element, its site preference in the σ-structure at 0 K is obtained by counting the number of times it appears at different sites in the configurations on the convex hull, and we have ordered the elements according to the Pettifor chemical scale (illustrated in Figure 6a)[43]. Clear trends emerge in Figure 6b for the site preference. From left to right on the Pettifor scale, the elements progressively prefer to occupy the smaller sites of 2a and 8i2 (both CN12), while they decrease their occupancies on larger sites: 8i1, 8j (both CN14) and 4f (CN15). There is a sharp transition in the region of Cr and Re regarding the preferred occupancy, in line with the fact that these two elements in binary systems can perfer high or low CN depending on the element they are associated with. Especially, from Mn to the right on the Pettifor scale, elements only have limited occurance on the larger sites (CN14 and CN15). As the Pettifor scale orders the considered elements by an increase of the number of valence electrons and an overall decrease of atomic size, the trend shown in Figure 6b is inline with experimental observations.[5, 8] Furthermore, the site preferences between 2a and 8i2 are correlated—the element that prefers to occupy the 2a site equally prefer to occupy the 8i2 site. Similar correlation can be found between the 8i1 and the 8j sites that share the same coordination number (CN14), and to a lesser extent the site 4f (CN15) except for Zr (a non-sigma forming element), suggesting that chemical similarity is shared between sites that has the same coordination number in the σ-structure.

Regarding the values of formation enthalpies of configurations on the convex hull associated with each elements shown in Figure 6c, we found that the absolute value of the formation enthalpies are minimal for Cr, Mo and W while elements with either lower or higher valence electron counts, such as Zr, Ti, Cu and Pt are more exothermic. This trend reveals the importance of the number of valence electrons in determining the formation enthalpies of σ-phase in multicomponent system—the value of formation enthalpy could be associated with the change of the valence electron counts (band filling) upon forming σ-phase with other elements. The chemical bonding and electronic structure of σ-phases will discussed in detail in Section 5.
<figure><img src="image_13.png" /><figcaption>Figure 7: Stability of σ-phase against quinary solid solution phases at equimolar composition. a) Elemental distribution in quinary σ-phase that is more stable than the considered solid solution phase (FCC, BCC and HCP), b) the distribution of energy difference $G_{\sigma} - G_{SS}^{\min}$ (blue) and $H_{\sigma} - G_{SS}^{\min}$ (red) between the σ-phase and the most stable solid solution phase, where $H_{\sigma}$ is the formation enthalpy contribution of Gibbs energy in the σ-phase, and c) the relationship between enthalpy of formation and entropy, for quinary σ-phase that are more stable than solid solution phases.</figcaption></figure>
## 4.2. Multicomponent σ-phase at equimolar compositions

The σ-phase is a very common intermetallic phase found in high-entropy alloys but it is not easy to predict its stability against solid solution phases in multicomponent systems from first-principles calculation due to the large computational cost associated with extensive number of end-member configurations.[11, 44] Using CEF and the machine learning prediction, we are able to characterize the stability of the σ-phase in multicomponent system by comparing its Gibbs energy with the Gibbs energy of BCC, FCC and HCP solid solution phases at equimolar composition. We have considered equimolar composition since i) it is a convenient choice in the discussion with multicomponent solid solution phases and ii) phase diagram calculations required to determine the homogeneity region of σ-phase across all systems, is non-trivial and is beyond the scope of this study. We note that the use of equimolar compositions may overestimate configurational entropy since the composition range of most experimentally observed σ-phase is off-equimolar. Therefore, we provide both a survey of experimental σ-phase composition and an error analysis of the equimolar assumption in Supplementary Section 5. We used $0.6T_m^{\text{ave}}$ for the Gibbs energy calculation, where $T_m^{\text{ave}}$ is the average melting point of constituent elements. This temperature reflects typical σ-phase stabilization at relatively high temperature and aligns with experimental binary σ-phase temperature range (Figure S8).[5]. For solid solution phase, their Gibbs energy is approximated by using binary interaction parameters provided by Chen et al., and is detailed in Appendix C.[45] Quinary system results are shown here; binary, ternary, and quaternary results appear in Supplementary Figures S9-S12.

The Gibbs energy difference $G_\sigma - G_{SS}^{\text{min}}$ between the σ-phase and the most stable solid solution phase at equimolar composition is shown in Figure 7b (blue histogram). Result reveals that the σ-phase is more stable than solid solution in 20.8% of the systems, confirming that the σ-phase is indeed an important competing phase in multicomponent alloys. Considering that many of the experimental σ-phase compounds are found at off-equimolar composition, the σ-phase could appear in even more systems. It is worth to highlight the important role of entropy in stabilizing the σ-phase. In Figure 7b, we compared the value of formation enthalpy of the σ-phase against the Gibbs energy of the solid solution phase ($H_\sigma - G_{SS}^{\text{min}}$, histogram in red) to highlight the fact that almost none of the σ-phase would be stable from only the formation enthalpy. The value of the entropy of these equimolar σ-phase are further shown in Figure 7c. The average value of
configurational entropy is 1.14 $k_B$/atom and can reach as high as 1.4 $k_B$/atom. Experimentally, a configurational entropy value of 1.18 $k_B$ is reported for a quinary $\sigma$-phase Cr$_{46.0}$Mn$_{15.2}$Fe$_{16.3}$Co$_{17.0}$Ni$_{5.5}$, and it agrees well with the range of values that are obtained from the CEF here.[8]

Figure 7a shows the elemental distribution in $\sigma$-phase that is more stable than solid solutions at equimolar composition. While Co, Cr, Ta, W and Mo can indeed be found frequently in the experimental compositions of the $\sigma$-phase, the high frequency of Zr in the elemental distribution is more likely because the omission of other intermetallic phases. Due to the large atomic size of Zr, Laves phase would likely to be preferred.[46] For Cu, its frequent appearance in quinary systems does not correspond to experimental observation in binary or ternary systems where no $\sigma$-phase that contain Cu is observed. However, Cu was found to exist in $\sigma$-phases that are found in high-entropy alloys that include Cu. Our calculation of stability of $\sigma$-phase in binary indeed show that Cu is not preferred to form $\sigma$-phase in binary (see Supplementary Figure S9), similar to experimental results, and we attributes its frequent appearance in multicomponent systems to its role in modulating valence electron concentration, which can lead to the stabilization of the $\sigma$-phase.

Compared to solid solution phases where the value of configurational entropy depends mostly on the number of components and the chemical composition, the configurational entropy in intermetallic phases depends on the chemical composition as well. Figure 7c shows how the calculated configurational entropy in the $\sigma$-phase is correlated with the formation enthalpies—large configurational entropy values are almost exclusively found in systems with positive formation enthalpy, while the $\sigma$-phase is more chemically ordered in systems with large, negative formation enthalpy. The trade-off between enthalpy and entropy is the result of the chemical ordering in the $\sigma$-phase structure and we have characterized the chemical ordering in these equimolar systems by the elemental site preferences, calculated from CEF. The result is shown in Figure 8a. As mentioned already, site mixing is, in general, significant in $\sigma$-phases and it is rare for an element to occupy more than 50% of the same crystallographic site in quinary systems, except in some cases for Ti, Zr and Nb. Larger atoms, which simultaneously have less valence electrons, prefer to occupy the larger sites with higher coordination numbers of 15 and 14 (4f, 8i, 1, 8j) and smaller atoms prefer the lower-coordinated sites of CN12 (2a, 8i, 2)—the same trend that are found for the ground state configurations.
<figure><img src="image_14.png" /><figcaption>Figure 8: a) Statistics of site occupancy $y_i^s$ of each element (i) at the different sites (s) of the σ-phase in all quinary systems at equimolar composition. In the box plot, the bottom and top bar correspond to the data range of 5% and 95%. b) Frequency heat map indicating pair preferences in all equimolar quinary σ-systems. The pair frequencies are given with respect to the baseline of a random solid solution. The heat map is symmetrical by construction.</figcaption></figure>

Based on the optimized site occupancies $y_i^s$, we characterized the tendency for two elements to neighbor each other in multicomponent σ-phase compared to a random solution. Such a “pair preference” is defined here by:

$$N_{ij} = \frac{1}{2} \sum_{\alpha} \sum_{\beta \in NN(\alpha)} (y_i^{\alpha} y_i^{\beta} - c_i c_j)$$

where i, j indicate elements and α, β indicate sites in the unit cell. NN(α) is the set of nearest neighbor atoms of α. Thus the summation is over all bonds between the nearest neighbor atoms in the unit cell. c_i and c_j are the molar fraction of the two elements, corresponding to random solution. The statistical pair preference in equimolar quinary systems were plotted in a heat map between elements in Figure 8b with elements ordered by the Pettifor scale. Since site occupancies are determined by the formation enthalpy of all end-members, the pair frequency should reveal pair interactions between neighboring elements since more frequent pairs should be energetically favorable compared to the less frequent ones. The clear trend in the heat map shows that the most frequent neighbor pairs are those formed between
one earlier atoms and one later atoms on the Pettifor scale (such pairs are located at the upper right and lower left region of the heat map). The least frequent pairs are those between electron-rich elements with small size (from Mn to Cu), which strongly avoid being neighbors with each other. These pair preferences thus determine the value of both enthalpy and entropy through the chemical ordering in the σ-phase. This motivated us to investigate the electronic structure and stability of the σ-phase from a bonding prospective.

**5. Bonding analysis**

In this section, we explain the stability of the σ-phase using bonding analysis in binary systems based on electronic structure calculation, crystal orbital Hamilton population (COHP) and its integrated value (ICOHP). According to its sign, COHP curve between a two-atom pair indicates the bonding/antibonding characteristics of the electron states shared by the two atoms at different energy, while the ICOHP at the Fermi level indicates the total amount of bonding/antibonding between the atom pair.[47] Calculations were performed using the LOBSTER code.[48]

We considered the case of Co forming binary σ-phases with other elements as example, shown in Figure 9a and b. Alloyed with isoelectronic elements, the formation enthalpy of Cr–Co, Mo–Co and W–Co at equimolar composition are very similar, despite the different atomic sizes of Cr, Mo and W (see atomic size used in Figure S7). Ti, Mo and Pt, in contrast, have similar atomic sizes in the σ-structure but different numbers of valence electrons. The formation enthalpy of their binaries with the electron-rich Co differs dramatically. Co–Ti is the most exothermic system with large negative formation enthalpy in almost all binary configurations, followed by the Co–Mo binaries. This observation suggests that electronic factor, rather than atomic size difference, play important role in the formation enthalpy of the σ-phase. At the same time, the energy difference between the different configurations of similar composition is also larger in the Ti–Co binary system, which leads to larger chemical ordering. The reverse is true in the Pt–Co system where the energy of different configurations changes linearly between the two unary end-points with small variance.

In Figure 9c, the COHP curves, summing over all nearest neighbor pairs in the unit cell, are shown for unary Co, unary Ti and the more energetically stable binary configuration CoCoTiCoTi (this elemental sequence specifies the placement of atoms of different kinds in the five sites 2a, 4f, 8i1, 8i2,
<figure><img src="image_15.png" /><figcaption>Figure 9: a) The formation enthalpy of binary configurations in Cr–Co, Mo–Co and W–Co systems, the boxplot on the right indicates the distribution of H<sub>Form</sub> of all configurations. b) The same plot but for Ti–Co, Mo–Co and Pt–Co systems. The arrows from the bottom to the top indicate three configurations: CoCoTiCoTi, CoCoTiTiCo and CoCoPtCoPt. c) Summed crystal orbital hamilton population (COHP) over all nearest neighbor atom pairs for unary Co, unary Ti and binary Co–Ti with configuration CoCoTiCoTi in the σ-structure. d) Integrated COHP (ICOHP) per bond (averaged) for Co–Co, Ti–Ti and Co–Ti nearest neighbor pairs in the configuration CoCoTiCoTi. The gray dashed curve shows the average ICOHP of Co–Co and Ti–Ti e) and f) display the same results as c) and d) but for the case of Pt.</figcaption></figure>

8j in the σ-phase structure). The COHP curves reveal that the position of the Fermi level is crucially related to the formation enthalpies in these three example configurations. In unary Co with a large valence electron concentration (VEC), antibonding states (with positive COHP value) are occupied, while the bonding states (with negative COHP value) are not fully filled in unary Ti in the σ-structure. However, in binary configuration CoCoTiCoTi, the Fermi level is just below the antibonding states, leading to a much stable compound. Pt–Co binary shown in Figure 9e is a counterexample to the
Ti–Co binary. Due to high VEC in both Co and Pt, electrons occupy the same amount of antibonding states without noticeable difference at different compositions in the Co–Pt binary system. Thus, the formation enthalpy of configurations in Pt–Co binary simply interpolates between unary Co and Pt as shown in Figure 9b. Therefore, COHP bonding analysis supports with previous empirical observations, structure maps and coarse-grained electronic structure calculations[5, 49, 50, 46, 51] and the optimal VEC could be approximately 6.5.

Just as the summed COHP can explain the formation enthalpies difference between different configurations, the atom-pair-wise COHP are found to explain the chemical ordering in the σ-structure. This is shown in Figure 9d and f with the average integrated COHP (ICOHP) curve for different types of nearest neighbor pairs in the two separate binary configurations: CoCoTiCoTi and CoCoPtCoPt. These ICOHP value at the Fermi level should be a good indicator of average bond energy between nearest neighbor pair. Consider the configuration CoCoPtCoPt first, the ICOHP values at the fermi level for three different types of nearest neighbor bonds are almost the same and, as a result, there is only very weak chemical ordering in the binary Co–Pt system. The situation is opposite in the CoCoTiCoTi configuration where the Co–Co nearest neighbor bonds are not favored energetically due to overfilling of electron states in the antibonding energy region. Therefore, strong chemical ordering can be expected in the Co–Ti binaries since bonds between Co and Ti are energetically favored over a mixture of Co–Co and Ti–Ti bonds (indicated by the gray dashed line). The relationship between bond energy and chemical ordering clearly explains the trend in pair preferences in Figure 8b where high VEC atoms (associated with band overfilling) tend to avoid neighboring each other.

## 6. Conclusion

In conclusion, we demonstrated the capability to construct a thermodynamic database for intermetallics using the machine learning method. All formation enthalpy of end-member configurations in the σ-structure, within a 19-element chemical space, can be predicted by the ALIGNN model with an accuracy (RMSE) of 10 meV/atom. Without prior DFT calculations, the relaxed crystal structures of the end-members are not known in advance. This is a bottleneck for achieving good prediction accuracy in structure-based ML models using GNN. Here, we shown that within a single complex structural
type, a regression model using OHE of configurations can give sufficiently accurate structure estimates leading to small prediction errors in formation enthalpy. The high accuracy is partly due to the use of the state-of-the-art ALIGNNN model, but also from training the model on various structural environments, including distorted σ-structure and additional TCP structures. Using the fully predicted database (about 2.5 million configurations), we systematically characterized both the structural configurations on the 0 K convex hull as well as the σ-phase at equimilar compositions considering substitutional disorder. The σ-phase is a common secondary phase in multicomponent high-entropy alloy system and it is important to determine the stabilization mechanism of the σ-phase compared to solid solutions. To this end, we have quantified the role of disorder in the stabilization of multicomponent σ-systems, and demonstrated a trade-off relationship between the enthalpy and entropy of the σ-phase. Electronic structure analysis using COHP shows that the values of both formation enthalpy and entropy are related to the chemical composition by valence electron count band filling, which helps us to rationalize the chemical trends observed from the theoretical characterization and predictions.

While our machine learning prediction of formation enthalpy is rather accurate, there are still several limitations that affect the accuracy of thermodynamic calculations for the σ-phase using the results presented here. First of all, we have limited ourselves to the compound energy formalism within the Bragg-Williams approximation, ignoring any excess Gibbs energy due to mixing on the same sublattice and the possibility of short-range ordering. Moreover, the entropy of formation from lattice vibrational and electronic excitation are also ignored, due to the limitation of this machine learning workflow and the expensive computational cost associated to evaluate vibrational spectrum from first-principles calculation. Dynamic stability, an important criterion for phase stability, is also not considered.[52] These are all directions to improve the thermodynamic description of not only the σ-phase considered here but other intermetallic phases in future works.

## Data Availability

The training and test dataset and scripts, along with the predicted formation enthalpy for all the σ-phase end-member configurations will be made publicly available on zenodo repository upon the acceptance of the manuscript.
## Acknowledgements

This work is a part of project "Artificial Intelligence for Intermetallic Materials"-AIIM (ref. ANR-22-CE92-0036) funded and supported by the French National Research Agency (ANR) and the German Deutsche Forschungsgemeinschaft (DFG No. 505643559). First-principles calculations were performed using HPC resources from GENCI-CINES (No. A0060906175) and NIMS Numerical Materials Simulator (Tsukuba, Japan). In addition, WZ acknowledges the financial support from Saint-Gobain (SGKK).

## Supporting Informations

The following supplementary files are available:

Section S1 Minimization of the Gibbs energy in the compound energy formalism.

Section S2 Details of the machine learning model training.

Section S3 Details of using bayesian optimization for structure relaxation.

Section S4 Prediction errors related to different elements.

Section S5 Survey of experimental compositions of the σ-phase.

Table S1 Training parameters of the multi-layer perception used in structure parameters prediction

Table S2 Prediction error of associated with different elements.

Table S3 Survey of representative compositions of the σ-phase in the literature.

Figure S1 ALIGNN predicted energy volume curves for selected σ-configurations

Figure S2 Distribution of elemental features used for machine learning and the structural deviation parameters in the σ-phase for structures in the training data.

Figure S3 Comparison of entropy term $-TS$ in the Gibbs energy between off-equimolar composition and equimolar composition.

Figure S4 The distribution of elements in the two test sets containing 3000 and 1500 σ phase configurations, respectively.

Figure S5 Statistics of relaxed structural parameters for all training configurations of σ phase.

Figure S6 The prediction results of cell volume and structure parameters using one-hot encoding and MLP regressor.

Figure S7 Atomic volumes in σ-phase obtained from linear regression of unit cell volume.
Figure S8 Comparison and distribution of temperature $0.6T_m^{\text{ave}}$ used in the compound energy formalism calculation.

Figure S9 Stability comparison between the σ-phase and the most stable solid solution phase in binary, ternary and quaternary systems.

Figure S10 Statistics of site occupancy and neighboring preferences in equimolar binary σ-phase.

Figure S11 Statistics of site occupancy and neighboring preferences in equimolar ternary σ-phase.

Figure S12 Statistics of site occupancy and neighboring preferences in equimolar quaternary σ-phase.
## Appendix A. DFT calculations

Density functional theory (DFT) calculations are performed using the Vienna Ab initio Simulation Package (VASP) version 6.4.1 using PBE functionals with an energy cutoff of 400 eV.[31, 32] The Brillouin zone (BZ) is sampled with a Γ centered regular *k*-grid with a total number of *k*-points 3000/*N*<sub>atom</sub> where *N*<sub>atom</sub> is the number of atoms in the calculation unit cell, giving well converged total energy. Specifically, for *σ*-phase, this condition gives a *k* grid of 8 × 8 × 4 in the reciprocal cell. Gaussian smearing with a smearing width of 0.1 eV is used for BZ summation. Tetrahedron summation should yield a more accurate total energy, but the difference compared to Gaussian smearing is negligible. Structures are set to relax until all forces on atoms are smaller than 0.01 eV/Å.

The DFT calculations are performed without spin-polarization, which saves considerable computational time, not only because the number of electronic states is halved, but also because determining global energy minimum of the magnetic ground state with certainty is very challenging by high-throughput calculations. Moreover, while the spin-polarized calculation could give a more accurate description of formation enthalpy of *σ*-phase at 0 K if the specific system is magnetic, it could be less relevant at higher temperature (> 1000 K) that are of more practical interest. Previous and recent literature reported that the magnetic transition temperatures in Fe, Ni and Co containing binary and ternary systems are in general below 100 K due to itinerant electron magnetism, with the only exception of Fe-V system which show relatively high $T_c$ = 315 K.[53, 4, 54, 55] On the other hand, the effect of spin-polarization on the formation enthalpy could later be added by extrapolating from spin-polarized calculations on a smaller scale, or using Calphad method to treat magnetic contribution to the Gibbs energy separately, as previously shown.[56, 57]

Formation enthalpies of the intermetallic phases are calculated by subtracting the energy per atom at their standard reference structures. The energy of the elemental reference structures are calculated with spin-polarization considered.
# Appendix B. Dataset Preparation

A total of 19 elements are included in our datasets: Al, Ti, V, Cr, Mn, Fe, Co, Ni, Cu, Zr, Nb, Mo, Tc, Ru, Ta, W, Re, Os and Pt.

Within the compound energy formalism, the end-member configurations are obtained by assigning atoms of different species to a set of distinct sublattices in a unit cell. In the case of the σ-structure, the 30 atoms in the unit cell are divided into the five sublattices given by non-equivalent Wyckoff positions (2a, 4f, 8i1, 8i2, 8j). Such a partition enables the consideration of up to quinary systems within the σ-phase structure and results in 2.5 million possible end-member configurations. In order to obtain a set of training configurations without potential bias, a systematic sampling strategy is used as follows:

1. All unary (19 configurations) and binary configurations (171 binary combinations of elements) are included (5149 in total).

2. For a total of 969 ternary combinations of 19 elements, we randomly select 5 from 150 possible σ-phase configurations for each ternary system (4845 in total).

3. For a total of 3876 quaternary combinations of 19 elements, two are randomly selected from 240 different σ-phase configurations for each quaternary system (7752 in total).

17746 σ-end-member configurations (0.7% of all possible configurations) make up the training set. The distribution of elements and chemical systems in the training configurations is balanced due to its construction.

Secondly, distorted structures of those training σ-configurations are included to provide more training examples with even more different bond lengths and bond angles. These distorted structures are generated by adding small Gaussian-like deviations to both the relaxed cell volume and atomic positions. Symmetry is not considered in this process. Each relaxed structure is distorted randomly three times, and DFT calculation without structure relaxation is used to evaluate their energy, resulting in an additional 17746 × 3 = 53238 structure-energy pairs.

Apart from the σ-structure, other intermetallic structures were used in the training set as well to reduce prediction error. A15, C15 (MgCu₂), C14 (MgZn₂), C36 (MgNi₂) and χ-structures are included. The sampling strategies are the same as one used for the σ-structures. In detail:
1. A15 structure has eight atoms among two non-equivalent positions (2a and 6c) in the unit cell. All 19 unary and 342 binary configurations are considered, resulting in a total of 361 configurations.

2. C15 (MgCu₂) structure has 24 atoms among two non-equivalent positions (8a and 16d) in the unit cell. All 19 unary and 342 binary configurations are considered, resulting in a total of 361 configurations.

3. C14 (MgZn₂) structure has 12 atoms in a hexagonal unit cell distributed in three non-equivalent sites (2a, 4f and 6h). All 19 unary and 1026 binary configurations are considered, together with three configurations sampled from each of the 969 ternary systems, resulting in a total of 3952 configurations.

4. C36 (MgNi₂) structure has five non-equivalent positions in a hexagonal cell (4e, 4f, 6g, 6h and 4f). The amount of configurations from unary to quaternary systems is exactly the same as for σ-structure.

5. χ-structure has four non-equivalent positions in the 58-atom body-centered cubic cell (2a, 8c, 24g₁ and 24g₂). All 19 unary configurations and all 2394 binary configurations are included. Five configurations are sampled from each of the 969 ternary systems (4845), and two configurations from each of the 3876 quaternary systems are sampled (7752).

There are in total 38114 configurations sampled from other intermetallic phases. Their relaxed structures and corresponding formation enthalpy are used for training the ALIGNN model.

Two test sets of σ-configurations without overlap with each other and the training set are constructed to test and select the structure prediction method (1500 configurations) and the overall prediction accuracy (3000 configurations). The smaller/larger training set contains 500/1000 ternary, 500/1000 quaternary, and 500/1000 quinary systems. All of them are sampled randomly. For these training and test configurations, DFT calculation is performed to relax the structure and obtain the corresponding formation energy.

Finally, to test the ML prediction of formation energy for all configurations in specific chemical systems, we divided all the 19 elements (by random shuffle and a split into triplets) to form seven ternary systems: Nb–W–Ti, Cr–Cu–Mn, Ru–Fe–Co, Ta–Re–Os, Tc–Al–Ni, Zr–Mn–V and Pt–Nb–W. The formation energy of all possible ternary configurations from these seven systems (150 × 7 = 1050) are obtained from DFT structure relaxation and calculation.
## Appendix C. Gibbs energy of the solid solution phases

The Gibbs energy of the solid solution phases are estimated using binary interaction model with interaction parameters from Chen et al.[45] For the BCC, FCC and HCP solid solutions, their Gibbs energy is expressed as:

$$G_{\alpha}(x, T) = \sum_{i} x_{i} H_{i}^{\alpha} + \sum_{i} \sum_{j>i} \Omega_{ij}^{\alpha} x_{i} x_{j} + k_{B} T \sum_{i} x_{i} \ln x_{i} + G_{\text{mag}}$$

where $x_i$ is the composition of element i in the solid solution phase and $H_i^\alpha$ is the formation enthalpy of solid solution phase ($\alpha$ = BCC/FCC/HCP) of element i in the paramagnetic state. This is obtained by subtracting magnetic ordering energy from the formation enthalpy calculated using DFT with spin-polarization, following Wang et al.[57] $\Omega_{ij}$ is the binary interaction parameter between element i and element j in the solid solution phase taken from Chen et al.[45] The use of only binary interaction to model multicomponent solid solution phase is an approximation. However, according to Bokas et al., the error of this approximation is relatively small with a root mean squared error of 16 meV/atom.[58] The third term on the right hand side is the configurational entropy of the solid solution and $G_{\text{mag}}$ is the temperature dependent contribution to the Gibbs energy due to magnetic ordering. The magnetic contribution to the Gibbs energy is necessary to account for the energy difference between magnetically ordered ground state in DFT calculation and paramagnetic state at high temperature. The magnetic moments and transition temperature of the solid solution phase are estimated by linear interpolation of the corresponding parameters in elemental phases (taken from SGTE unary database) to calculate $G_{\text{mag}}$.[59]

## References

[1] O. Senkov, J. Miller, D. Miracle, C. Woodward, Accelerated exploration of multi-principal element alloys with solid solution phases, Nat Commun 6 (1) (2015) 6529. doi:10.1038/ncomms7529.

[2] C. Rae, R. Reed, The precipitation of topologically close-packed phases in rhenium-containing superalloys, Acta Materialia 49 (19) (2001) 4113-4125. doi:10.1016/S1359-6454(01)00265-8.

[3] W. Liu, Z. Lu, J. He, J. Luan, Z. Wang, B. Liu, Y. Liu, M. Chen, C. Liu, Ductile CoCrFeNiMo$_x$ high entropy alloys strengthened by hard
intermetallic phases, Acta Materialia 116 (2016) 332-342. doi:10.1016/
j.actamat.2016.06.063.

[4] A. K. Sinha, Topologically close-packed structures of transition metal
alloys, Progress in Materials Science 15 (2) (1972) 81-185. doi:10.
1016/0079-6425(72)90002-3.

[5] J.-M. Joubert, Crystal chemistry and Calphad modeling of the σ phase,
Progress in Materials Science 53 (3) (2008) 528-583. doi:10.1016/j.
pmatsci.2007.04.001.

[6] J.-C. Crivello, M. Palumbo, T. Abe, J.-M. Joubert, Ab initio ternary -
phase diagram: The Cr-Mo-Re system, Calphad 34 (4) (2010) 487-494.
doi:10.1016/j.calphad.2010.09.002.

[7] K. Han, R. Sahara, T. Abe, K. Oikawa, N. Ueshima, I. Ohnuma, Phase equilibria of the Co-Cr-Mn ternary system at 700 °C, Journal of Alloys and Compounds 965 (2023) 171315. doi:10.1016/j.jallcom.2023.171315.

[8] J.-M. Joubert, Y. Kalchev, A. Fantin, J.-C. Crivello, R. Zehl, E. Elkaim,
G. Laplanche, Site occupancies in a chemically complex σ-phase from
the high-entropy Cr-Mn-Fe-Co-Ni system, Acta Materialia 259 (2023)
119277. doi:10.1016/j.actamat.2023.119277.

[9] D. Miracle, O. Senkov, A critical review of high entropy alloys and
related concepts, Acta Materialia 122 (2017) 448-511. doi:10.1016/j.
actamat.2016.08.081.

[10] T. Abe, K. Han, Y. Goto, I. Ohnuma, T. Koyama, Stabilization of Equiatomic Solutions Due to High-Entropy Effect, Mater. Trans. 64 (4) (2023) 877-884. doi:10.2320/matertrans.MT-M2022167.

[11] M.-H. Tsai, R.-C. Tsai, T. Chang, W.-F. Huang, Intermetallic Phases in High-Entropy Alloys: Statistical Analysis of their Prevalence and Structural Inheritance, Metals 9 (2) (2019) 247. doi:10.3390/met9020247.

[12] N. Saunders, A. P. Miodownik, CALPHAD (Calculation of Phase Diagrams): A Comprehensive Guide, Elsevier, Burlington, 1998.
[13] A. Walle, G. Ceder, Automating first-principles phase diagram calculations, JPE 23 (4) (2002) 348–359. doi:10.1361/105497102770331596.

[14] S. Shang, Y. Wang, Z.-K. Liu, First-principles calculations of phonon and thermodynamic properties in the boron-alkaline earth metal binary systems: B-Ca, B-Sr, and B-Ba, Phys. Rev. B 75 (2) (2007) 024302. doi:10.1103/PhysRevB.75.024302.

[15] P. Y. Chew, A. Reinhardt, Phase diagrams—Why they matter and how to predict them, The Journal of Chemical Physics 158 (3) (2023) 030902. doi:10.1063/5.0131028.

[16] B. Sundman, J. Ågren, A regular solution model for phases with several components and sublattices, suitable for computer applications, Journal of Physics and Chemistry of Solids 42 (4) (1981) 297–301. doi:10.1016/0022-3697(81)90144-X.

[17] J.-C. Crivello, J.-M. Joubert, N. Sokolovska, Supervised deep learning prediction of the formation enthalpy of complex phases using a DFT database: The σ - phase as an example, Computational Materials Science 201 (2022) 110864. doi:10.1016/j.commatsci.2021.110864.

[18] J. Gilmer, S. S. Schoenholz, P. F. Riley, O. Vinyals, G. E. Dahl, Neural Message Passing for Quantum Chemistry, arXiv:1704.01212 [cs] (Jun. 2017). arXiv:1704.01212.

[19] T. Xie, J. C. Grossman, Crystal Graph Convolutional Neural Networks for an Accurate and Interpretable Prediction of Material Properties, Phys. Rev. Lett. 120 (14) (2018) 145301. doi:10.1103/PhysRevLett. 120.145301.

[20] A. Dunn, Q. Wang, A. Ganose, D. Dopp, A. Jain, Benchmarking materials property prediction methods: The Matbench test set and Automatminer reference algorithm, npj Comput Mater 6 (1) (2020) 138. doi:10.1038/s41524-020-00406-3.

[21] C. W. Park, C. Wolverton, Developing an improved crystal graph convolutional neural network framework for accelerated materials discovery, Phys. Rev. Materials 4 (6) (2020) 063801. doi:10.1103/PhysRevMaterials.4.063801.
[22] A. Kolluru, M. Shuaibi, A. Palizhati, N. Shoghi, A. Das, B. Wood, C. L. Zitnick, J. R. Kitchin, Z. W. Ulissi, Open Challenges in Developing Generalizable Large-Scale Machine-Learning Models for Catalyst Discovery, ACS Catal. 12 (14) (2022) 8572–8581. doi:10.1021/acscatal.2c02291.

[23] J. Damewood, J. Karaguesian, J. R. Lunger, A. R. Tan, M. Xie, J. Peng, R. Gómez-Bombarelli, Representations of Materials for Machine Learning, Annu. Rev. Mater. Res. 53 (1) (2023) 399–426. doi:10.1146/annurev-matsci-080921-085947.

[24] S. S. Omee, N. Fu, R. Dong, M. Hu, J. Hu, Structure-based out-of-distribution (OOD) materials property prediction: A benchmark study, npj Comput Mater 10 (1) (2024) 144. doi:10.1038/s41524-024-01316-4.

[25] N. Islam, W. Huang, H. L. Zhuang, Machine learning for phase selection in multi-principal element alloys, Computational Materials Science 150 (2018) 230–235. doi:10.1016/j.commatsci.2018.04.003.

[26] Y. Zhang, C. Wen, C. Wang, S. Antonov, D. Xue, Y. Bai, Y. Su, Phase prediction in high entropy alloys with a rational selection of materials descriptors and machine learning models, Acta Materialia 185 (2020) 528–539. doi:10.1016/j.actamat.2019.11.067.

[27] R. Sikdar, N. Roy, B. Selvaratnam, V. Mishra, D. Mumbaraddi, A. Mondal, K. Buxi, P. P. Jana, A. Mar, Ahead by a Century: Discovery of Laves Phases Assisted by Machine Learning, Inorg. Chem. 63 (13) (2024) 5972–5981. doi:10.1021/acs.inorgchem.3c04647.

[28] X. Liu, J. Zhang, Z. Pei, Machine learning for high-entropy alloys: Progress, challenges and opportunities, Progress in Materials Science 131 (2023) 101018. doi:10.1016/j.pmatsci.2022.101018.

[29] M. Palumbo, S. G. Fries, T. Hammerschmidt, T. Abe, J.-C. Crivello, A. A. H. Breidi, J.-M. Joubert, R. Drautz, First-principles-based phase diagrams and thermodynamic properties of TCP phases in Re-X systems (X=Ta, V, W), Computational Materials Science 81 (2014) 433–445. doi:10.1016/j.commatsci.2013.08.051.
[30] J. Schmidt, N. Hoffmann, H.-C. Wang, P. Borlido, P. J. M. A. Carriço, T. F. T. Cerqueira, S. Botti, M. A. L. Marques, Machine-Learning-Assisted Determination of the Global Zero-Temperature Phase Diagram of Materials, Advanced Materials 35 (22) (2023) 2210788. doi:10.1002/adma.202210788.

[31] G. Kresse, J. Furthmüller, Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set, Phys. Rev. B 54 (16) (1996) 11169-11186. doi:10.1103/PhysRevB.54.11169.

[32] J. P. Perdew, K. Burke, M. Ernzerhof, Generalized Gradient Approximation Made Simple, Phys. Rev. Lett. 77 (18) (1996) 3865-3868. doi:10.1103/PhysRevLett.77.3865.

[33] E. Becht, L. McInnes, J. Healy, C.-A. Dutertre, I. W. H. Kwok, L. G. Ng, F. Ginhoux, E. W. Newell, Dimensionality reduction for visualizing single-cell data using UMAP, Nat Biotechnol 37 (1) (2019) 38-44. doi:10.1038/nbt.4314.

[34] K. Choudhary, B. DeCost, Atomistic Line Graph Neural Network for improved materials property predictions, npj Comput Mater 7 (1) (2021) 185. doi:10.1038/s41524-021-00650-1.

[35] R. Ruff, P. Reiser, J. Stühmer, P. Friederich, Connectivity optimized nested line graph networks for crystal structures, Digital Discovery 3 (3) (2024) 594-601. doi:10.1039/D4DD00018H.

[36] J. Noh, G. H. Gu, S. Kim, Y. Jung, Uncertainty-Quantified Hybrid Machine Learning/Density Functional Theory High Throughput Screening Method for Crystals, J. Chem. Inf. Model. 60 (4) (2020) 1996-2003. doi:10.1021/acs.jcim.0c00003.

[37] A. Jain, T. Bligaard, Atomic-position independent descriptor for machine learning of material properties, Phys. Rev. B 98 (21) (2018) 214112. doi:10.1103/PhysRevB.98.214112.

[38] J. Schmidt, L. Pettersson, C. Verdozzi, S. Botti, M. A. L. Marques, Crystal graph attention networks for the prediction of stable materials, Sci. Adv. 7 (49) (2021) eabi7948. doi:10.1126/sciadv.abi7948.
[39] R. E. A. Goodall, A. S. Parackal, F. A. Faber, R. Armiento, A. A. Lee, Rapid discovery of stable materials by coordinate-free coarse graining, Sci. Adv. 8 (30) (2022) eabn4117. doi:10.1126/sciadv.abn4117.

[40] J. Gibson, A. Hire, R. G. Hennig, Data-augmentation for graph neural network learning of the relaxed energies of unrelaxed structures, npj Comput Mater 8 (1) (2022) 211. doi:10.1038/s41524-022-00891-8.

[41] Y. Zuo, M. Qin, C. Chen, W. Ye, X. Li, J. Luo, S. P. Ong, Accelerating materials discovery with Bayesian optimization and graph deep learning, Materials Today 51 (2021) 126-135. doi:10.1016/j.mattod.2021.08.012.

[42] Y. Su, J. Wang, Machine learning aided high-throughput first-principles calculations to predict the formation enthalpy of σ phase, Calphad 82 (2023) 102599. doi:10.1016/j.calphad.2023.102599.

[43] D. G. Pettifor, Bonding and Structure of Molecules and Solids, Clarendon Press; Oxford University Press, Oxford: New York, 1995.

[44] R.-C. Tsai, A.-C. Fan, W.-F. Huang, Y.-W. Shi, D. B. Miracle, W.-Y. Lai, M.-H. Tsai, Composition of Intermetallic Phases in High-Entropy Alloys: Compositional Inheritance of High-Entropy Sigma Phases from Binary Counterparts, High Entropy Alloys & Materials 2 (1) (2024) 33-40. doi:10.1007/s44210-023-00027-7.

[45] W. Chen, A. Hilhorst, G. Bokas, S. Gorsse, P. J. Jacques, G. Hautier, A map of single-phase high-entropy alloys, Nat Commun 14 (1) (2023) 2856. doi:10.1038/s41467-023-38423-7.

[46] B. Seiser, R. Drautz, D. Pettifor, TCP phase predictions in Ni-based superalloys: Structure maps revisited, Acta Materialia 59 (2) (2011) 749-763. doi:10.1016/j.actamat.2010.10.013.

[47] V. L. Deringer, A. L. Tchougréeff, R. Dronskowski, Crystal Orbital Hamilton Population (COHP) Analysis As Projected from Plane-Wave Basis Sets, J. Phys. Chem. A 115 (21) (2011) 5461-5466. doi:10.1021/jp202489s.

[48] R. Nelson, C. Ertural, J. George, V. L. Deringer, G. Hautier, R. Dronskowski, LOBSTER: Local orbital projections, atomic charges,
and chemical-bonding analysis from projector-augmented-wave-based
density-functional theory, J Comput Chem 41 (21) (2020) 1931-1940.
doi:10.1002/jcc.26353.

[49] M.-H. Tsai, K.-Y. Tsai, C.-W. Tsai, C. Lee, C.-C. Juan, J.-W. Yeh, Criterion for Sigma Phase Formation in Cr- and V-Containing High-Entropy Alloys, Materials Research Letters 1 (4) (2013) 207-212. doi: 10.1080/21663831.2013.831382.

[50] T. Hammerschmidt, A. F. Bialon, D. G. Pettifor, R. Drautz, Topologically close-packed phases in binary transition-metal compounds: Matching high-throughput ab initio calculations to an empirical structure map, New J. Phys. 15 (11) (2013) 115016. doi:10.1088/1367-2630/15/11/115016.

[51] B. Seiser, T. Hammerschmidt, A. N. Kolmogorov, R. Drautz, D. G. Pettifor, Theory of structural trends within 4d and 5d transition metal topologically close-packed phases, Phys. Rev. B 83 (22) (2011) 224116. doi:10.1103/PhysRevB.83.224116.

[52] A. Vesti, D. Music, T. Hammerschmidt, M. Palumbo, P. A. T. Olsson, First-principles study of dynamic instability and phase transformation in the WRe and WOs sigma phase, Phys. Rev. Materials 9 (6) (2025) 063602. doi:10.1103/PhysRevMaterials.9.063602.

[53] D. Read, E. Thomas, J. Forsythe, Evidence of itinerant electron ferro-magnetism in sigma phase alloys, Journal of Physics and Chemistry of Solids 29 (9) (1968) 1569-1572. doi:10.1016/0022-3697(68)90099-1.

[54] J. Ciesłak, B. F. Costa, S. M. Dubiel, M. Reissner, W. Steiner, Magnetic ordering above room temperature in the sigma-phase of Fe₆₆V₃₄, Journal of Magnetism and Magnetic Materials 321 (14) (2009) 2160-2165. doi: 10.1016/j.jmmm.2009.01.005.

[55] J. Cieslak, J. Tobola, M. Reissner, Magnetic properties of sigma-phase FeCrX (X=Co, Ni) alloys: Experimental and theoretical study, Acta Materialia 123 (2017) 35-43. doi:10.1016/j.actamat.2016.10.010.

[56] A. Breidi, M. Andasmas, J.-C. Crivello, N. Dupin, J.-M. Joubert, Experimental and computed phase diagrams of the Fe-Re system, J. Phys.:
Condens. Matter 26 (48) (2014) 485402. doi:10.1088/0953-8984/26/48/485402.

[57] P. Wang, T. Hammerschmidt, U. R. Kattner, G. B. Olson, Structural stability of Co–V intermetallic phases and thermodynamic description of the Co–V system, Calphad 68 (2020) 101729. doi:10.1016/j.calphad.2019.101729.

[58] G. Bokas, W. Chen, A. Hilhorst, P. Jacques, S. Gorsse, G. Hautier, Unveiling the thermodynamic driving forces for high entropy alloys formation through big data ab initio analysis, Scripta Materialia 202 (2021) 114000. doi:10.1016/j.scriptamat.2021.114000.

[59] A. Dinsdale, SGTE data for pure elements, Calphad 15 (4) (1991) 317–425. doi:10.1016/0364-5916(91)90030-N.
Supplementary Information for:
Prediction and Characterization of Chemically Complex
σ-Phase Intermetallics with Graph Neural Network

Wenhao Zhang$^{a,b}$, Mariano Forti$^{c}$, Runan Xie$^{d}$, Céline Barreteau$^{d}$,
Abe Taichi$^{b,*}$, Jean-Marc Joubert$^{d}$, Thomas Hammerschmidt$^{c,*}$, Jean-Claude Crivello$^{a,*}$

a. CNRS-Saint-Gobain-NIMS, IRL 3629, Laboratory for Innovative Key Materials and Structures (LINK), 1-1 Namiki, Tsukuba, Japan

b. Research Center for Structural Materials, National Institute for Materials Science, 1-2-1 Sen-
gen, Tsukuba, Japan

c. Interdisciplinary Centre for Advanced Materials Simulation (ICAMS), Ruhr-University, Uni-
versitätsstrasse 150, Bochum, Germany

d. Univ Paris Est Creteil, CNRS, ICMPE, UMR 7182, 2 rue Henri Dunant, Thiais, France

Abe.Taichi@nims.go.jp; thomas.hammerschmidt@rub.de; jean-claude.crivello@cnrs.fr

<footnote>\* Corresponding authors:</footnote>
## S1 Minimization of the Gibbs Energy

To determine the Gibbs energy curve at each composition from the compound energy formalism, it is necessary to minimize $G(x, T)$ with respect to internal degrees of freedom $0 \le y_i^s \le 1$ subject to the following constraints:

$$ \sum_{i} y_{i}^{s} = 1 \\ \sum_{s} m_{s} y_{i}^{s} = N_{\text{atom}} $$

In this work, numerical minimizations of $G(x, T)$ are done with both equality and inequality constraints utilizing the gradient:

$$ \frac{\partial G}{\partial y_i^s} = \sum_{j,k,l,m,n} \left[ y_j^{2a} y_k^{4f} y_l^{8i1} y_m^{8i2} y_n^{8j} \right] \frac{\delta_i^s}{y_i^s} H_{kjlmn}^0 + k_B T \sum_s m_s (1 + \ln y_i^s) \\ \delta_i^s = \begin{cases} 1 & \text{if } (i,s) \in \{(j,2a), (k,4f), (l,8i1), (m,8i2), (n,8j)\} \\ 0 & \text{otherwise} \end{cases} $$
## S2 Details of ML Model Training

Multi-layer perceptron (MLP) combined with one-hot encoding (OHE) representation of σ-phase configurations is used to predict nine structure parameters: volume, c/a ratio, $x_{4f}$, $x_{8i(1)}$, $y_{8i(1)}$, $x_{8i(2)}$, $y_{8i(2)}$, $x_{8j}$, $z_{8j}$. The learning rate and $\alpha$ are tuned by 8-fold CV for each of these structure parameters separately while the network hidden layer is fixed to be one with 150 units. These hyperparameters are listed in Table S1. The scikit-learn python library is used for this learning task.

<table><caption>Table S1: Learning Rates and Weight Decay Parameters of the Multi-Layer Perception used in Structure Parameters Prediction</caption><thead><tr><th>Internal parameter</th><th>Learning rate</th><th>α</th><th>Parameter</th><th>Learning rate</th><th>α</th></tr></thead><tbody><tr><td>Volume</td><td>1 × 10<sup>-3</sup></td><td>0.05</td><td>c/a</td><td>5 × 10<sup>-3</sup></td><td>0.05</td></tr><tr><td>x<sub>4f</sub></td><td>1 × 10<sup>-3</sup></td><td>0.1</td><td>/</td><td>/</td><td>/</td></tr><tr><td>x<sub>8i(1)</sub></td><td>5 × 10<sup>-3</sup></td><td>0.1</td><td>y<sub>8i(1)</sub></td><td>1 × 10<sup>-3</sup></td><td>0.1</td></tr><tr><td>x<sub>8i(2)</sub></td><td>1 × 10<sup>-3</sup></td><td>0.1</td><td>y<sub>8i(2)</sub></td><td>5 × 10<sup>-3</sup></td><td>0.05</td></tr><tr><td>x<sub>8j</sub></td><td>5 × 10<sup>-3</sup></td><td>0.05</td><td>z<sub>8j</sub></td><td>5 × 10<sup>-3</sup></td><td>0.1</td></tr></tbody></table>

We used the ALIGNN model (network structure and hyperparameters) as is published, except that the number of nearest neighbors is increased to 16 for dealing with the closed-packed intermetallic structures. Furthermore, initial node features that should encode atomic properties are slightly tuned by 1) removing the atomic volume, which should correlate with covalent radius, and 2) changing the classes in the “block” feature from (s, p, d, f) to a more detailed classification of (alkali metal, alkaline earth, transition metal, ..., lanthanide, actinide). All of the atomic properties are one-hot encoded and concatenated together. The size of the initial feature vector is 93, almost the same as the size of the original feature vector (92). The model is trained using a GPU on a workstation configured with NVIDIA RTX A5500 graphics card. Training of the final model with around 108,000 entries in the training set took roughly one day to complete with 150 epochs. For the final model, the training and validation errors (RMSE) are 2.70 and 6.62 meV/atom, respectively.
## S3 Structure Relaxation using Bayesian Optimization

The problem of structural relaxation is the minimization of the internal energy of the crystal with respect to the atomic positions. In DFT calculations, the minimization of the atomic positions relies on the calculated atomic forces. In this work, we have also tried to use the well-trained ALIGNN model to relax the structure by minimizing the predicted energies with respect to degrees of freedom in the σ-structure.

Different from DFT calculation, only energies are used to guide the structure relaxation process. The starting structures are the ones obtained by fitting the volume. Initially, iterative gradient (from finite difference method) based minimization (BFGS) is tried. However, during test it is found that energies would not reduce and convergence was not reached. This is likely to be caused by noisy gradient obtained by surrogating the ALIGNN model, which is shown below by an energy-volume curve of a few selected configurations. Training the ALIGNN model to predict force on atoms could enable gradient based structure relaxation, but it was not tried since most of the training entries are well relaxed.

<figure><img src="image_16.png" /><figcaption>Figure S1: Energy predicted by the ALIGNN model at volumes differing from the relaxed volumes in five σ-phase configurations taken from the test set.</figcaption></figure>

While the gradient curve is found to be noisy, clear minimal of energies was predicted for volumes near the relaxed ones. Thus, a global, black-box minimization method, such as Bayesian optimization, could be effective to locate the energy minimal that correspond to a relaxed structure. For each structure, Bayesian optimization was done to subsequently optimize cell parameters (volume and c/a ratio) and internal coordinates. The two steps were then repeated one more time to better locate the energy minimal. In each of these optimization processes, 200 points were initially explored before the optimization process start and 150 points were subsequently
explored using upper confidence bound (UCB) as acquisition function. Regarding Bayesian op-
timization, it is worth mentioning the underlying Gaussian process regression, scaling with the
number of data points as $O(n^3)$, is more time consuming than energy prediction from ALIGNN
model. The computational cost associated with Bayesian optimization needs to be solved if such
minimization to be applied in a larger scale. As discussed in the article, the accuracy of formation
enthalpies prediction is much improved (23.3 meV/atom) using Bayesian optimization, compared
to structures from statistic average, with an overall underestimation of around -12 meV/atom.
The bond length and bond angles are also more accurately predicted. However, for only a single
structural type, Bayesian optimization is not as accurate compared to one-hot encoding structure
prediction in both final structures and energies.
## S4 Prediction errors related to different elements

<table><caption>Table S2: Prediction error of formation enthalpies from ALIGNN model and predicted structure for subsets within the 3000 σ phase test set. “w/” indicate the subset of configurations that contain the element while “w/o” indicate the subset of configurations without the corresponding element. All elements in the test sets are evenly distributed as plotted in Figure S7.</caption><thead><tr><th>Elements</th><th>Al</th><th>Ti</th><th>V</th><th>Cr</th><th>Mn</th><th>Fe</th><th>Co</th><th>Ni</th><th>Cu</th></tr></thead><tbody><tr><td>RMSE w/ (meV/atom)</td><td>12.4</td><td>9.9</td><td>8.9</td><td>9.8</td><td>11.4</td><td>10.5</td><td>11.0</td><td>9.8</td><td>11.0</td></tr><tr><td>RMSE w/o (meV/atom)</td><td>10.3</td><td>10.3</td><td>10.6</td><td>10.3</td><td>9.9</td><td>10.2</td><td>10.0</td><td>10.3</td><td>10.0</td></tr></tbody><thead><tr><th>Elements</th><th>Zr</th><th>Nb</th><th>Mo</th><th>Tc</th><th>Ru</th><th>Ta</th><th>W</th><th>Re</th><th>Os</th><th>Pt</th></tr></thead><tbody><tr><td>RMSE w (meV/atom)</td><td>16.3</td><td>8.2</td><td>8.8</td><td>7.3</td><td>9.1</td><td>10.1</td><td>9.4</td><td>8.3</td><td>8.9</td><td>11.4</td></tr><tr><td>RMSE w/o (meV/atom)</td><td>7.8</td><td>10.7</td><td>10.6</td><td>10.9</td><td>10.5</td><td>10.3</td><td>10.4</td><td>10.7</td><td>10.6</td><td>9.9</td></tr></tbody></table>

Table S2 reports the prediction error associated with each elements by comparing prediction errors in configurations that contain/not contain a given elements. The notable trends that can be observed from Table S2 regarding the errors are:

1. Al and Zr are associated with notable increase of error.

2. Electron-rich/small 3d elements (Mn, Fe, Co, Cu) are associated slightly larger error compared to the more electron-poor elements such as V, Nb, Mo and Re.

In general, machine learning prediction can be accurate for a structural configuration only if the machine learning model is trained with a sufficient number of structural configurations that are similar to it. This can help to understand the larger error in the case of configurations with Al and Zr. Different from all other elements in the training data, Al is the only element that belongs to group IIIA in the third row of the periodic table and the only element which does not have *d* electrons. Zr, on the other hand, has a significantly larger atomic radius. The elemental features of Zr and Al are therefore relatively different from other elements that we considered in machine learning. The difference in elemental features between Zr, Al and the rest of the elements is shown in the Figure S2 where the input atomic features of different elements are represented in a two dimensional representation given by UMAP technique.

Similarly, machine learning prediction error can also be large if bond length or bond angle in a configuration deviate quite largely from the statistical average, i.e., the crystal structure is quite different from a typically σ-phase structure. The right hand side of Figure S2 shows the statistics of deviation parameter Δ(x) for each structure in the training set from the average crystal structure of the σ-phase, measured by atomic positions in the unit cell. While most of the configurations are not far from the statistical average structure, some configurations deviate significantly from the average (indicated by orange color). The inset of the figure shows the
<figure><img src="image_17.png" /><figcaption>Figure S2: Left: distribution of elemental features for each elements used for the machine learning using UMAP. Right: deviation Δ(x) of fractional atomic position from the average value and the elemental distribution in the structures with large Δ(x).</figcaption></figure>

distribution of elements in the configurations that deviate the most from the average, and it is clear that configurations that contain Zr lead to very large structural deviations because of its large atomic size. This structural deviation from average could also contribute to the very large error associated with the element Zr.

Finally, the slightly larger error associated with the electron-rich elements (Mn, Fe, Co, Cu), compared to the electron-poor ones, could be due to the fact that there are more electron-poor elements in the training dataset. Secondly, these errors could be related to differences in electronic structure. Through bonding analysis, we have shown through COHP that an average electron concentration of around 6-7 is preferred for the stability of the σ-phase. It could be possible that the electronic structure (therefore bonding) is slightly different when the configuration is more electron-rich compared to electron-poor, leading to errors in enthalpy prediction for configurations with elements Mn, Fe, Co, Cu.
S5 Survey of experimental compositions of the σ–phase

Although equimolar composition is used to characterize σ–phase in multicomponent systems in the manuscript, this choice could inherently lead to an overestimation of entropy of the σ–phase since many of the experimentally observed σ–phase exist at an off-equimolar composition. Here, we provide a comparison between entropy values evaluated at equimolar and off-equimolar compositions in seom selected systems. First, we surveyed literature for the composition of σ–phase that are observed experimentally. The homogeneity region or composition of σ–phase are taken from three literature sources: i) (Joubert 2008) for binary systems, ii) (Seiser et al. 2011) for ternary σ–phase and iii) (Tsai et al. 2024) for σ–phase in high-entropy alloy systems. We only included systems that contain only the elements considered in this study. For binaries, the homogeneity regions are extracted from Figure 2 in (Joubert 2008), and we selected composition at the center of the homogeneity region. For ternary systems listed in (Seiser et al. 2011), a representative composition is chosen at the approximate center of the entire homogeneity in the isothermal section of phase diagrams, at temperature close to 0.6*T<sub>m</sub><sup>ave</sup>, where T<sub>m</sub><sup>ave</sup> is the mean value of melting temperature of all elements. Lastly, the composition of σ–phase in high-entropy alloy systems are directly taken from the supplementary materials of (Tsai et al. 2024), but we have limited the phase to contain at most 5 elements. The representative compositions are listed in the following Table S3.

<figure><img src="image_18.png" /><figcaption>Figure S3: Comparison of entropy term −*TS* in the Gibbs energy between off-equimolar compositions listed in Table S3 and equimolar composition in the same system.</figcaption></figure>

Figure S3 shows the comparison between compound energy formalism calculation done at equimolar composition and off-equimolar composition at 0.6*T<sub>m</sub><sup>ave</sup>* of the system surveyed. In majority of the binary and ternary systems, the values of the entropy −*TS* term in off-equimolar compositions are not too far away from the corresponding values in equimolar compositions (within and around 20 meV/atom). The difference is larger for *σ*–phase composition found in high-entropy-alloys (up to 40 meV/atom). Some of the cases in which the values of −*TS* differs
<table><caption>Table S3: Survey of representative composition in the literature.</caption><thead><tr><th>Type</th><th>Systems</th></tr></thead><tbody><tr><td>Binary</td><td>V<sub>19</sub>Mn<sub>81</sub>, Cr<sub>22</sub>Mn<sub>78</sub>, V<sub>24</sub>Re<sub>76</sub>, W<sub>29</sub>Tc<sub>71</sub>, Mo<sub>30</sub>Tc<sub>70</sub>, Cr<sub>32</sub>Tc<sub>68</sub>, Mo<sub>35</sub>Mn<sub>65</sub>, Mo<sub>39</sub>Re<sub>61</sub>, Re<sub>40</sub>Fe<sub>60</sub>, Ta<sub>41</sub>Re<sub>59</sub>, Tc<sub>41</sub>Fe<sub>59</sub>, W<sub>42</sub>Re<sub>58</sub>, Nb<sub>46</sub>Re<sub>54</sub>, Cr<sub>47</sub>Fe<sub>53</sub>, V<sub>49</sub>Fe<sub>51</sub>, Mo<sub>50</sub>Fe<sub>50</sub>, Re<sub>51</sub>Mn<sub>49</sub>, Cr<sub>57</sub>Co<sub>43</sub>, V<sub>57</sub>Co<sub>43</sub>, Tc<sub>60</sub>Mn<sub>40</sub>, Re<sub>61</sub>Cr<sub>39</sub>, Mo<sub>62</sub>Co<sub>38</sub>, Mo<sub>62</sub>Ru<sub>38</sub>, Nb<sub>63</sub>Al<sub>37</sub>, Nb<sub>64</sub>Os<sub>36</sub>, W<sub>63</sub>Ru<sub>37</sub>, V<sub>65</sub>Ni<sub>35</sub>, Mo<sub>65</sub>Os<sub>35</sub>, Cr<sub>66</sub>Ru<sub>34</sub>, Nb<sub>66</sub>Pt<sub>34</sub>, Ta<sub>68</sub>Al<sub>32</sub>, Cr<sub>67</sub>Os<sub>33</sub>, Ta<sub>67</sub>Os<sub>33</sub>, W<sub>72</sub>Os<sub>28</sub>, Ta<sub>73</sub>Pt<sub>27</sub></td></tr><tr><td>Ternary</td><td>Co<sub>30</sub>Cr<sub>55</sub>Fe<sub>15</sub>, Co<sub>35</sub>Cr<sub>50</sub>Mo<sub>15</sub>, Co<sub>35</sub>Cr<sub>60</sub>Mo<sub>5</sub>, Co<sub>25</sub>Cr<sub>50</sub>Re<sub>25</sub>, Co<sub>40</sub>Cr<sub>55</sub>Ti<sub>5</sub>, Co<sub>40</sub>Cr<sub>55</sub>W<sub>5</sub>, Co<sub>5</sub>Fe<sub>55</sub>Re<sub>40</sub>, Co<sub>30</sub>Fe<sub>35</sub>V<sub>35</sub>, Co<sub>33</sub>Mn<sub>33</sub>V<sub>34</sub>, Co<sub>20</sub>Ni<sub>20</sub>V<sub>60</sub>, Cr<sub>33</sub>Fe<sub>33</sub>Mn<sub>34</sub>, Cr<sub>35</sub>Fe<sub>50</sub>Mo<sub>15</sub>, Cr<sub>20</sub>Fe<sub>30</sub>Re<sub>50</sub>, Cr<sub>20</sub>Fe<sub>55</sub>V<sub>25</sub>, Cr<sub>45</sub>Fe<sub>45</sub>W<sub>10</sub>, Cr<sub>35</sub>Mo<sub>35</sub>Ni<sub>30</sub>, Cr<sub>40</sub>Ni<sub>10</sub>Re<sub>50</sub>, Cr<sub>60</sub>Ni<sub>30</sub>W<sub>10</sub>, Cr<sub>10</sub>Re<sub>50</sub>W<sub>40</sub>, Fe<sub>25</sub>Mn<sub>35</sub>Mo<sub>40</sub>, Fe<sub>33</sub>Mn<sub>33</sub>V<sub>34</sub>, Fe<sub>45</sub>Mo<sub>30</sub>V<sub>25</sub>, Fe<sub>50</sub>Ni<sub>5</sub>Re<sub>45</sub>, Fe<sub>30</sub>Ni<sub>20</sub>V<sub>50</sub>, Mn<sub>35</sub>Ni<sub>25</sub>V<sub>50</sub>, Mo<sub>40</sub>Nb<sub>5</sub>Re<sub>55</sub>, Mo<sub>40</sub>Ni<sub>35</sub>Re<sub>25</sub>, Mo<sub>45</sub>Re<sub>45</sub>Ru<sub>10</sub>, Mo<sub>25</sub>Re<sub>55</sub>Ta<sub>20</sub>, Mo<sub>20</sub>Re<sub>65</sub>V<sub>15</sub>, Mo<sub>22</sub>Re<sub>56</sub>W<sub>22</sub>, Mo<sub>20</sub>Re<sub>70</sub>Zr<sub>10</sub>, Mo<sub>55</sub>Ru<sub>40</sub>Ta<sub>5</sub>, Mo<sub>55</sub>Ru<sub>35</sub>W<sub>10</sub>, Nb<sub>15</sub>Re<sub>50</sub>W<sub>35</sub>, Ni<sub>15</sub>Re<sub>63</sub>V<sub>22</sub>, Os<sub>33</sub>Ta<sub>33</sub>W<sub>34</sub>, Re<sub>30</sub>Ru<sub>20</sub>W<sub>50</sub>, Re<sub>45</sub>Ta<sub>20</sub>W<sub>35</sub>, Re<sub>55</sub>V<sub>5</sub>W<sub>40</sub></td></tr><tr><td>HEA</td><td>Cr<sub>48.4</sub>Fe<sub>18.6</sub>Co<sub>16.6</sub>Ni<sub>11.5</sub>Al<sub>4.6</sub>, Cr<sub>32.6</sub>Mn<sub>28.1</sub>Fe<sub>32.6</sub>Ni<sub>6.0</sub>Al<sub>0.7</sub>, Cr<sub>40.7</sub>Mo<sub>13.4</sub>Fe<sub>38.3</sub>Ni<sub>4.1</sub>Al<sub>3.5</sub>, V<sub>37.7</sub>Fe<sub>31.9</sub>Ni<sub>22.4</sub>Al<sub>2.2</sub>Cu<sub>5.8</sub>, Cr<sub>40.5</sub>Mn<sub>12.2</sub>Fe<sub>24.6</sub>Co<sub>21.9</sub>Cu<sub>0.7</sub>, Cr<sub>35.0</sub>Mn<sub>22.7</sub>Fe<sub>23.2</sub>Co<sub>19.2</sub>, Cr<sub>46.8</sub>Mn<sub>16.9</sub>Fe<sub>21.3</sub>Co<sub>4.9</sub>Ni<sub>9.9</sub>, Cr<sub>24.1</sub>Mo<sub>25.6</sub>Fe<sub>16.9</sub>Co<sub>17.8</sub>Ni<sub>15.6</sub>, Cr<sub>50.3</sub>Fe<sub>25.9</sub>Co<sub>17.5</sub>Ni<sub>6.3</sub>, Cr<sub>45.8</sub>Fe<sub>27.5</sub>Co<sub>14.9</sub>Ni<sub>10.1</sub>Ta<sub>1.7</sub>, Cr<sub>24.1</sub>V<sub>21.9</sub>Fe<sub>19.4</sub>Co<sub>19.3</sub>Ni<sub>15.3</sub>, Cr<sub>55.5</sub>Fe<sub>29.1</sub>Co<sub>12.4</sub>Ni<sub>2.6</sub>Zr<sub>0.4</sub>, Mo<sub>22.1</sub>V<sub>23.2</sub>Fe<sub>19.9</sub>Co<sub>19.9</sub>Ni<sub>14.9</sub>, Cr<sub>44.3</sub>Mn<sub>18.9</sub>Fe<sub>26.4</sub>Ni<sub>10.5</sub>, Cr<sub>20.7</sub>Mo<sub>39.9</sub>Fe<sub>16.2</sub>Ni<sub>23.1</sub></td></tr></tbody></table>

most significantly between off-equimolar composition and equimolar composition are labeled in the figure, and they correspond to cases where the σ-phase from a lower order system does not extend very much to a higher order system. For example, the composition Re<sub>55</sub>W<sub>40</sub>V<sub>5</sub> in the ternary case.
Supplementary Figures and Tables

<figure><img src="image_19.png" /><figcaption>Figure S4: The distribution of elements in the two test sets containing 3000 and 1500 σ phase configurations, respectively.</figcaption></figure>
<figure><img src="image_20.png" /></figure>

Figure S5: Relaxed structure parameter statistics for all the training configurations of σ phase. For unit cell volume, the relaxed volume is compared to a sum of elemental volume obtained by linear regression, which is listed in Figure S7.
<figure><img src="image_21.png" /><figcaption>Figure S6: The prediction of cell volume and structure parameters using one-hot encoding and MLP regressor. The R² scores in the test set, containing 1500 relaxed structure, for each of the structure parameters (from top to bottom, left to right) are 0.998, 0.869, 0.841, 0.589, 0.903, 0.919, 0.655, 0.635 and 0.441.</figcaption></figure>
<figure><img src="image_22.png" /><figcaption>Figure S7: Atomic volume in σ-phase obtained from linear regression of unit cell volume with respect to the number of atoms of different species in the unit cell.</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>Figure S8: Left: comparison of temperatures used for the compound energy formalism calculation ($0.6T_m^{ave}$) with the experimental temperature range for the σ-phase in binary systems. Right: statistics of $0.6T_m^{ave}$ in all ternary, quaternary and quinary systems.</figcaption></figure>
<figure><img src="image_24.png" /><figcaption>Figure S9: Relative stability of the σ-phase in binary, ternary and quaternary systems compared to solid solution phase at equimolar compositions. First row: elemental distribution in σ-phase that are more stable than solid solution. Second row: the distribution of Gibbs energy difference between σ-phase and the solid solution phase, and third row: the relationship between entropy S and formation enthalpy H for the σ-phase that are more stable than solid solution. The σ-phase is more stable than solid solution phase in 18.6%, 27.2% and 22.4% of the systems from binary to quaternary.</figcaption></figure>
<figure><img src="image_25.png" /><figcaption>Figure S10: Statistics of site occupancy and neighboring preference in equimolar binary σ phase.</figcaption></figure>

<figure><img src="image_26.png" /><figcaption>Figure S11: Statistics of site occupancy and neighboring preference in equimolar ternary σ phase.</figcaption></figure>
<figure><img src="image_27.png" /><figcaption>Figure S12: Statistics of site occupancy and neighboring preference in equimolar quaternary σ phase.</figcaption></figure>