Accurate Screening of Functional Materials with Machine-Learning Potential and Transfer-Learned Regressions: Heusler Alloy Benchmark

TL;DR

This paper introduces a machine learning-accelerated high-throughput workflow for discovering magnetic materials, specifically Heusler alloys, combining transfer learning and validation with density functional theory to improve prediction accuracy and generalization.

Contribution

It presents a novel ML-HTP workflow for magnetic material discovery using transfer learning and benchmarking of interatomic potentials, advancing high-throughput screening methods.

Findings

High predictive accuracy confirmed by DFT validation

Effective transfer learning improves model generalization

Benchmarking shows performance of different uMLIPs

Abstract

A machine learning-accelerated high-throughput (HTP) workflow for the discovery of magnetic materials is presented. As a test case, we screened quaternary and all-$d$ Heusler compounds for stable compounds with large magnetocrystalline anisotropy energy ($E_{\mathrm{aniso}}$). Structure optimization and evaluation of formation energy and distance to hull convex were performed using the eSEN-30M-OAM interatomic potential, while local magnetic moments, phonon stability, magnetic stability, and $E_{\mathrm{aniso}}$ were predicted by eSEM models trained on our DxMag Heusler database. A frozen transfer learning strategy was employed to improve accuracy. Candidate compounds identified by the ML-HTP workflow were validated with density functional theory, confirming high predictive precision. We also benchmark the performance of different uMLIPs, discuss the fidelity of local magnetic moment prediction, and demonstrate generalization to unseen elements via transfer learning from a universal interatomic potential.

Tables

Table 1. Table 1: Performance comparison of the eSEN uMLIP and MLRM with ALIGNN models and previously reported results. The size of dataset used for MLRM is also listed. 0 0 footnotetext: 1 [ mitraMachineLearningApproach2023 ] Dataset size is about 1000, 2 [ huSearchingHighSpin2020 ] Dataset size is about 65,000, 3 [ luExplainableAttentionCNN2025 ] Dataset size is about 500,000 for a a and m total m_{\mathrm{total}} , and about 450,000 for Δ ​ E \Delta E and Δ ​ H \Delta H , 4 [ miyazakiMachineLearningBased2021 ] Dataset size is 143, 5 [ xieScreeningNewQuaternary2023 ] Dataset size is 16,272, 6 [ kimMachinelearningacceleratedHighthroughputMaterials2018 ] Dataset size is 426,148, 7 [ liuMachineLearningAssisted2025 ] Dataset size is 1153, 8 [ hilgersMachineLearningbasedEstimation2025 ] Dataset size is 408, 9 [ hirohataMachineLearningDevelopment2025 ] Dataset size is 6500.

Property	Metric	eSEN-30M-OAM	ALIGNN-FF	Previous reports
$a$	$R^2$	0.994	0.128	0.80–0.9411footnotemark: 1, 0.9422footnotemark: 2, 0.98733footnotemark: 3, 0.97944footnotemark: 4
	RMSE	0.023	0.330	0.11-0.1255footnotemark: 5
$\Delta E$	$R^2$	0.995	0.453	0.80–0.8811footnotemark: 1, 0.9322footnotemark: 2, 0.98233footnotemark: 3
	RMSE	0.029	0.310	0.11755footnotemark: 5
$\Delta H$	$R^2$	0.980	0.330	0.9166footnotemark: 6, 0.96933footnotemark: 3
Property	Metric	eSEN MLRM	ALIGNN MLRM	Previous reports	Dataset size
$\left\{{𝒎}_i\right\}$	$R^2$	0.989	—	—	27,864
$m_{\mathrm{total}}$	$R^2$	0.986	0.904	0.75–0.8911footnotemark: 1, 0.8277footnotemark: 7, 0.92733footnotemark: 3	27,864
$\sum \left|{𝐦}_i\right|$	$R^2$	0.989	0.891	—	27,864
${\omega }_{\min }$	$R^2$	0.750	0.734	—	8,198
$T_c$	$R^2$	0.910	0.844	0.7688footnotemark: 8, 0.7399footnotemark: 9	2,106
	Accu.	0.910	—	0.7388footnotemark: 8
$E_{\mathrm{aniso}}$	$R^2$	0.680	0.592	—	6,123
\botrule

Table 2. Table 2: Computational cost of each task in the HTP workflow, reported in node-minutes. DFT entries show per-job statistics (mean/median/IQR). ML entries show mean values. For structure relaxation, one job corresponds to a single relaxation from an initial distortion. For structure optimization, the cost per compound should be multiplied by the number of initial structures, which varies across material systems; here, we report the cost for a single structure relaxation. The uMLIP is pre-trained and used without fine-tuning, so its training time is excluded. DFT structure optimizations directly yield local magnetic moments, so no separate timing is reported.

	DFT (mean/median/IQR)	ML training (mean)	ML prediction (mean)
Structure relaxation	5.6/3.8/3.0	—	0.12
$\left\{{𝒎}_i\right\}$	—	876	1e-03
Phonon stability	4542/4008/1442	258	1e-03
$T_c$	70/66/14	78	1e-03
$E_{\mathrm{aniso}}$	142/133/34	192	1e-03

Full text

1 Introduction

The high-throughput (HTP) screening approach has emerged as a powerful strategy for accelerating the discovery of novel materials by systematically exploring large chemical spaces computationally or experimentally [sanvitoAcceleratedDiscoveryNew2017, zhangHighthroughputDesignMagnetic2021, barwalLargeMagnetoresistanceHigh2024]. The density functional theory (DFT)-based HTP workflows have been widely employed to identify materials with target properties. However, as the search space increases, the associated computational cost becomes unsustainable, often restricting screening efforts to a manageable subspace [faleevHeuslerCompoundsPerpendicular2017, huHighthroughputDesignCobased2023]. To address this bottleneck, machine learning (ML) offers a promising route by drastically reducing computational costs. In this work, we demonstrate the robust integration of state-of-the-art ML techniques into the HTP workflow (ML-HTP) through a practical case study focused on screening quaternary and all-$d$ Heusler compounds for stable candidates with strong magnetic anisotropy energy ($E_{\mathrm{aniso}}$).

Initial realizations of the ML-HTP paradigm relied on ML models that utilize compositional descriptors as input features [hilgersMachineLearningbasedEstimation2025, baigutlinMachineLearningAlgorithms2024, mitraMachineLearningApproach2023, liuMachineLearningPredict2021]. These models directly map chemical formulas to target properties, offering efficiency and simplicity. However, composition-based models are inherently unable to distinguish compounds with identical stoichiometry but different atomic arrangements. One workaround involves assigning layer indices to atomic sites, but this approach fixes the number of sites and can yield inconsistent predictions for symmetry-equivalent structures [xieScreeningNewQuaternary2023, luExplainableAttentionCNN2025]. Crystal graph-based models do not have such drawbacks since they explicitly incorporate structure information as input, capturing structure-property relationships more accurately [xieCrystalGraphConvolutional2018a, choudharyAtomisticLineGraph2021]. However, crystal graph-based models introduce an additional computational step, as geometry optimization must precede property prediction.

Although a single DFT optimization typically requires only a few minutes, the cumulative cost of screening a large number of candidate compounds becomes prohibitively expensive, particularly for magnetic systems where multiple magnetic configurations must be considered. A promising solution lies in leveraging universal machine learning interatomic potentials (uMLIPs), which can accelerate structure optimization by several orders of magnitude relative to DFT. The uMLIP field has witnessed rapid advancements in recent years, with many crystal graph-based models proposed. Despite the conceptual appeal, reliable and robust uMLIP-based structure optimization has only become practical recently, especially following the release of the large-scale and diverse Meta Open Materials 2024 Dataset (OMat24) training dataset in 2024 [barroso-luqueOpenMaterials20242024]. This is demonstrated in the current work by benchmarking several uMLIPs, ranging from early-stage implementations to state-of-the-art developments.

With optimized structures, properties can be predicted using machine learning regression model (MLRM). This approach substantially reduces computational cost, particularly for properties that are expensive to compute via DFT, such as phonon spectra, conductivity, magnetic critical temperature ($T_{{\rm c}}$), and $E_{\mathrm{aniso}}$. However, training accurate MLRMs typically requires large, high-quality datasets. To overcome this challenge, transfer learning (TL) techniques can be employed to adapt pretrained ML models to new tasks [yamadaPredictingMaterialsProperties2019, leeTransferLearningMaterials2021a, hoffmannTransferLearningLarge2023]. TL leverages models that have already learned generalizable representations from extensive datasets and fine-tunes them using smaller, task-specific datasets. This strategy enhances predictive accuracy while substantially reducing data requirements.

As a case study, we conducted a ML-HTP screening on Heusler compounds, which have garnered significant attention due to their diverse functional properties, technological potential, and structure complexity [heComputationallyAcceleratedDiscovery2022]. Numerous DFT-HTP screenings have been carried out to identify candidates with different desirable properties [sanvitoAcceleratedDiscoveryNew2017, faleevHeuslerCompoundsPerpendicular2017, nokyGiantAnomalousHall2020, huHighthroughputDesignCobased2023, xingChemicalsubstitutiondrivenGiantAnomalous2024a, xiaoHighthroughputComputationalScreening2025]. In our previous work, we developed DXMag Computational Heusler Database (HeuslerDB), a comprehensive database encompassing nearly all conventional ternary Heusler compounds. The present study significantly extended the search space to include quaternary and all-$d$ Heusler compounds, targeting stability and $E_{\mathrm{aniso}}$ as key screening criteria. In earlier DFT-HTP studies of Heusler compounds, 10 candidates with large $E_{\mathrm{aniso}}$ were identified out of 286 selected compositions, and 15 among 29,784 Co-based structures [faleevHeuslerCompoundsPerpendicular2017, huHighthroughputDesignCobased2023]. This low yield underscores the rarity of such materials and highlights the difficulty of this search problem, making it a stringent test case for ML-HTP approaches.

Previous studies attempted ML approaches to $E_{\mathrm{aniso}}$ in systems outside the Heusler family. They employed early ML models, such as crystal graph convolutional neural networks (CGCNN) and compositional-descriptor models, to predict $E_{\mathrm{aniso}}$ in Fe–Co–N alloys and physics-informed 2D materials, working on the order of hundreds of compounds [xieDataDrivenStudiesMagnetic2021, liaoPredictingMagneticAnisotropy2022, duttaMachineLearningAssisted2022]. Here, we employ state-of-the-art ML methods to extend the scope to hundreds of thousands of compounds with improved accuracy and practicality.

In this work, we demonstrate the use of uMLIP and TL-MLRM s as drop-in replacements for DFT structure optimization and property evaluation within HTP framework, as illustrated in Fig. 1(a). As a practical application, we employed this approach to identify conventional quaternary and all-$d$ Heusler compounds with large $E_{\mathrm{aniso}}$, while simultaneously satisfying thermodynamic, dynamic, and magnetic stability. Structure optimization and thermodynamic stability evaluation were performed using the eSEN-30M-OAM uMLIP [fuLearningSmoothExpressive2025]. The following predictions of local magnetic moment ($\{\boldsymbol{m}_{i}\}$), minimum phonon frequency ($\omega_{\mathrm{min}}$), $T_{{\rm c}}$, and $E_{\mathrm{aniso}}$ were performed using MLRM s. The MLRM s were trained via frozen transfer learning, using eSEN-30M-OAM uMLIP as the base model and fine-tuned using HeuslerDB data and newly computed data. ML-selected candidates were validated through DFT calculations to demonstrate the significant reliability of this ML-HTP approach. We further examine key factors that influence the performance of such ML-HTP workflows, including accuracy of uMLIP-based structure optimization, the magnetic configuration prediction, the performance of frozen transfer learning technique and generalization to unseen elements. The underlying code for the ML-HTP workflow is made available as open-source packages MLIP High-throughput Optimization and Thermodynamics (MLIP-HOT) and MLIP Frozen Transfer Learning (MLIP-FTL), which can be found on our group’s website and Git repository.

2 Results

To accumulate data for training $E_{\mathrm{aniso}}$ MLRM, we first computed the $E_{\mathrm{aniso}}$ of conventional ternary Heusler compounds within HeuslerDB using DFT. The $E_{\mathrm{aniso}}$ of some Heusler compounds were reported in the previous work and the agreement between our DFT results and previous work is demonstrated in Fig. S1 [faleevHeuslerCompoundsPerpendicular2017, huHighthroughputDesignCobased2023]. Among all conventional ternary Heusler compounds, 2190 (7.9 %) exhibit an $E_{\mathrm{aniso}}$ magnitude greater than 1 MJ/m3. When further screened for thermodynamic, dynamical, and magnetic stability, only 135 compounds (0.5 %) meet both the high $E_{\mathrm{aniso}}$ and stability criteria, which are presentated in Table S2. These low percentages underscore the difficulty of identifying stable, high $E_{\mathrm{aniso}}$ compounds and highlight the need for more efficient screening methods demonstrated in current work as a case study. The ML-HTP workflow for this case study is summarized in Fig. 1(b); detailed computational procedures are provided in Sec. 4.

For conventional quaternary compounds, we enumerated all combinations where $X$ and $Y$ are transition metals from the $d$-block (excluding Tc and Hg), and $Z$ is a main-group element from groups 13, 14, or 15 of the $p$-block. In addition, La and Lu were included for $X$ and $Y$ because their 4$f$ orbitals are either empty or fully filled. This exhaustive enumeration, accounting for symmetry constraints, yielded 131,544 unique compositions. For the all-$d$ Heuslers, we extended the screening space to include $d$-block transition metals together with La and Lu across all four sites ($X_{1}$, $X_{2}$, $Y$, and $Z$), resulting in a separate set of 105,763 unique compositions. A schematic of the screened chemical space is presented in Supplementary Information as Fig. S2.

2.1 Validation of ML-HTP selected candidates

Using eSEN-30M-OAM uMLIP in combination with MLRM s, we screened almost all conventional quaternary and all-$d$ Heusler spaces for stable compounds with high $E_{\mathrm{aniso}}$. As a result, 334 and 924 candidates were found, respectively. To evaluate the reliability of this ML workflow, all candidates were validated using DFT calculations. The results are summarized in Fig. 2 and detailed data are provided in Tables S3 and S4. The percentages that DFT results satisfy the screening criteria (i.e., precision) are shown as blue and yellow bars for conventional quaternary and all-$d$ Heusler compounds, respectively. The selection criteria include $c/a$ ratio, $\Delta E$, $\Delta H$, $\{\boldsymbol{m}_{i}\}$, $\omega_{\mathrm{min}}$, $T_{{\rm c}}$, and $E_{\mathrm{aniso}}$. For comparison, the precisions of the ML models, measured on the test set of conventional ternary compounds, are also provided as green bars.

In ML-HTP, the $c/a$ ratio, $\Delta E$, and $\Delta H$ are obtained from the optimized structure and the corresponding energy by eSEN-30M-OAM uMLIP. During the structure optimization process, relaxations were performed starting from multiple initial structures, and the relaxed structure with the lowest energy was selected. Since a non-zero $E_{\mathrm{aniso}}$ requires the Heusler compound to adopt a tetragonal phase, we applied a screening threshold of $|c/a-1|>0.01$ to identify tetragonality. Notably, all ML-selected candidates remained tetragonal in DFT validation, confirming that eSEN-30M-OAM reliably distinguishes between cubic and tetragonal phases. A more detailed discussion of structure optimization and performance for lattice parameters $a$, $c$, and the $c/a$ ratio prediction, along with performance of other uMLIP s, are provided in Sec. 2.3.

Using energies of candidate compounds, elements, and competing phases predicted by eSEN-30M-OAM, the $\Delta E$ and $\Delta H$ were calculated. The criteria of $\Delta E$ $<0$ eV/atom and $\Delta H$ $<0.22$ eV/atom were employed to identify thermodynamically stable candidates, following the thresholds established in our previous DFT-HTP study [xiaoHighthroughputComputationalScreening2025]. Among ML-selected candidates, 99.1% of conventional quaternary and 97.8% of all-$d$ Heusler compounds were validated to have $\Delta E_{\mathrm{DFT}}<0$ eV/atom. Similarly, 96.4% (quaternary) and 98.8% (all-$d$) of the compounds were found to have $\Delta H_{\mathrm{DFT}}<0.22$ eV/atom. These high validation rates demonstrate that state-of-the-art uMLIPs, such as eSEN-30M-OAM, can reliably assess thermodynamic stability.

It is important to note that the eSEN-30M-OAM uMLIP used is a general-purpose, pretrained model without any fine-tuning specific to the Heusler chemical space. Thus, these results highlight its strong generalization, making it an effective drop-in replacement for DFT-based optimization and thermodynamic stability assessment in HTP workflow. Notably, the model achieves strong performance on the studied magnetic systems, despite not explicitly incorporating magnetic moments into its architecture or training. This strong performance and generalization are particularly valuable for the initial screening of novel material systems, where uMLIPs can greatly reduce the search space by rapidly and reliably estimating optimized structures and thermodynamic stability.

The properties $\{\boldsymbol{m}_{i}\}$, $\omega_{\mathrm{min}}$, $T_{{\rm c}}$, and $E_{\mathrm{aniso}}$ were predicted using MLRM s applied to uMLIP-optimized structures. Because $E_{\mathrm{aniso}}$ is a magnetic property, $\{\boldsymbol{m}_{i}\}$ was predicted, and a screening threshold of $\sum|\boldsymbol{m}_{i}|>0.1$ $\mu_{\mathrm{B}}$/f.u was applied to identify magnetic compounds. DFT validations confirmed all ML-selected candidates to be magnetic. Moreover, the $\{\boldsymbol{m}_{i}\}$ MLRM accurately predicts both the magnitude and sign of local moments, as discussed in detail in Sec. 2.5. Magnetic system identification is a critical yet computationally demanding step in DFT-HTP, as multiple initial $\{\boldsymbol{m}_{i}\}$ values must be tested and the low fraction of magnetic systems in some material families can lead to substantial wasted computation. By incorporating the $\{\boldsymbol{m}_{i}\}$ MLRM method as a pre-screening step, the search space can be substantially reduced.

To identify compounds with dynamic stability, magnetic stability, and large $E_{\mathrm{aniso}}$, we applied the criteria $\omega_{\mathrm{min}}$ $>-10$ cm*-1*, $T_{{\rm c}}$ $>300$ K, and |$$E_{\mathrm{aniso}} $|>1$ MJ/m3. Among the ML-selected candidates, 89.2% of conventional quaternary and 93.1% of all-$d$ Heusler compounds were validated to have $\omega_{\mathrm{min}}$ above -10 cm*-1*. For magnetic stability, 81.7% of conventional quaternary and 80.4% of all-$d$ Heusler compounds were validated to have $T_{{\rm c}}$ above 300 K. For the target property $E_{\mathrm{aniso}}$, the validation rates were 82.0% and 68.2%, respectively. To assess the sensitivity to criteria values, we also evaluated the precisions using a range of more stringent thresholds. The results, summarized in Fig. S3, show that the selection precision is not sensitive to the threshold values in the investigated range.

The MLRM s were trained exclusively using train set of conventional ternary Heusler compounds, yet were applied to evaluate quaternary and all-$d$ compositions. By comparing validation rates to the precisions calculated using the test set of conventional ternary compounds, the MLRM s for $\{\boldsymbol{m}_{i}\}$ and $T_{{\rm c}}$ generalize well to these expanded chemical spaces. In contrast, the $E_{\mathrm{aniso}}$ model exhibits lower performance for all-$d$ compounds. This discrepancy can be attributed to the difference in chemical environments: while conventional quaternary compounds retain $Z$-site elements from the $p$-block—consistent with the training set—all-$d$ compounds introduce $Z$ elements from the $d$-block, which were absent during training. $E_{\mathrm{aniso}}$ is a sensitive property, influenced by subtle details of the electronic structure, and thus more susceptible to domain shifts than $\{\boldsymbol{m}_{i}\}$ or $T_{{\rm c}}$.

Relaxing the screening thresholds increases the pool of candidate compounds and might capture promising cases missed initially at the expense of more false positives. Additionally, the curated list of compounds that meet stability criteria, and further magnetic system criterion, serves as an efficient starting point for investigation of other functional properties. For readers interested in exploring an expanded candidate list, the full set of ML-predicted data for 131,544 conventional quaternary and 105,763 all-$d$ Heusler compounds will be accessible through HeuslerDB.

2.2 Distribution of strong $E_{\mathrm{aniso}}$ candidates

In addition to validating predictive precision, it is essential to determine whether the ML models capture known physical trends. A well-established insight is that compounds containing 4d and 5d elements typically exhibit larger $E_{\mathrm{aniso}}$ than those composed of 3d elements, owing to the stronger spin–orbit coupling associated with the heavier atomic nuclei of 4d and 5d elements. This behavior is clearly reproduced in the ML-HTP results, as shown in Fig. 3, which presents the distribution of candidate compounds according to the presence of 4d/5d elements. The figure further highlights the specific 4d/5d elements that appear in the identified compounds. For comparison, Fig. 3 also includes the distribution of compounds with DFT-calculated $E_{\mathrm{aniso}}$ magnitudes exceeding 1 MJ/m3.

2.3 uMLIP optimization performance

In recent years, uMLIPs have advanced rapidly, with numerous new models proposed and trained. To identify the most suitable model for our screening workflow, we benchmarked representative uMLIPs models developed since 2023. These include ALIGNN-FF, CHGNet, SevenNet-l3i5, SevenNet-MF-ompa, HIENet, MatterSim-v1, eqV2-S-OAM, eqV2-M-OAM, eqV2-L-OAM, and the latest eSEN-30M-OAM [choudharyUnifiedGraphNeural2023, dengCHGNetPretrainedUniversal2023, kimDataEfficientMultifidelityTraining2025, yanMaterialsFoundationModel2025, yangMatterSimDeepLearning2024, barroso-luqueOpenMaterials20242024, fuLearningSmoothExpressive2025]. The evaluation focused on structure optimization for 10,000 conventional ternary compounds randomly selected from the ground states in HeuslerDB. To identify the global minimum, 14 initial structures were generated by applying strain to the conventional cell (two formula units) and converting it to the primitive cell (one formula unit). Specifically, the $a$, $b$, and $c$ axes were uniformly scaled by $\pm 10\%$ and $\pm 30\%$, or the $c$-axis alone was varied by $\pm 10\%$, $\pm 20\%$, $\pm 30\%$, $\pm 40\%$, and $\pm 50\%$. The lowest-energy structure from these relaxations was selected as the predicted ground state. Convergence tests for all evaluated models are shown in Fig. S4, demonstrating that the selected ground states are well-converged.

The performance of structure optimization was assessed by comparing the predicted lattice constants $a$ and $c$, and the resulting $c/a$ ratio, with corresponding DFT values. The results are summarized in Fig. 4. The relative error (RE) is defined as the maximum of $\left|\frac{x_{\mathrm{ML}}}{x_{\mathrm{DFT}}}-1\right|$ and $\left|\frac{x_{\mathrm{DFT}}}{x_{\mathrm{ML}}}-1\right|$, where $x_{\mathrm{ML}}$ and $x_{\mathrm{DFT}}$ denote the values predicted by the ML and DFT, respectively. We report the fractions of compounds within 5% and 1% RE tolerances. Among the evaluated models, the eSEN-30M-OAM and eqV2 models achieved the highest accuracy at 5% RE, with eSEN-30M-OAM showing slightly better performance at the stricter 1% RE threshold. A key distinction between the two models is the number of local minima encountered: eqV2-L-OAM identified 91,585 local minima, whereas eSEN-30M-OAM identified 32,606. Both counts are much greater than 10,000, confirming the presence of multiple local minima, but are still significantly less than 10,000×14. This suggests that many different initial distortions ultimately converge to the same local minimum. Importantly, eSEN-30M-OAM found substantially fewer local minima than eqV2-L-OAM. This difference is attributed to the smoother energy landscape of eSEN-30M-OAM. Convergence tests also demonstrate that eSEN-30M-OAM achieves convergence with fewer initial structures, resulting in lower computational cost for HTP purpose.

The predictive performance of total energy ($E$), formation energy ($\Delta E$), and energy above the convex hull ($\Delta H$) using uMLIP was assessed by comparing the uMLIP predictions with DFT values, using the absolute error (AE), $|x_{\mathrm{ML}}-x_{\mathrm{DFT}}|$. The fractions of compounds with AE below 0.01 and 0.05 eV/atom are shown in Fig. 4. Among the benchmarked models, eSEN-30M-OAM and the eqV2 variants showed the highest accuracy for $\Delta E$ and $\Delta H$ at the 0.05 eV/atom threshold, with eSEN-30M-OAM displaying a slight drop in accuracy at the more stringent 0.01 eV/atom level. Predicted total energies from uMLIP were found to be systematically lower than DFT values, reducing direct agreement; however, this offset also applies to elemental references and competing phases, so the relative quantities $\Delta E$ and $\Delta H$ remain in strong agreement with DFT. Given its robust performance in both structure optimization and thermodynamic stability, eSEN-30M-OAM was selected for integration into the ML-HTP workflow.

2.4 Improvements Over Existing Approaches

Previous studies have reported the performance of composition-based models in predicting lattice constants, $\Delta E$, and $\Delta H$ for cubic Heusler compounds. For comparison, we evaluated metrics of eSEN-30M-OAM on the cubic Heusler subset, with results summarized in Table 1. The $R^{2}$ score for the lattice constant $a$ is 0.994, surpassing the previously reported ranges of 0.80–0.94 across different Heusler types and the values of 0.94, 0.979, and 0.987 in other works [mitraMachineLearningApproach2023, huSearchingHighSpin2020, miyazakiMachineLearningBased2021, luExplainableAttentionCNN2025]. Similarly, the $R^{2}$ for $\Delta E$ reaches 0.995, outperforming earlier results of 0.80–0.88 and 0.93, 0.982 [mitraMachineLearningApproach2023, huSearchingHighSpin2020, luExplainableAttentionCNN2025]. The $R^{2}$ for $\Delta H$ is 0.98, exceeding prior values of 0.91 and 0.969 [kimMachinelearningacceleratedHighthroughputMaterials2018, luExplainableAttentionCNN2025]. The root mean squared errors (RMSE) for $a$ and $\Delta E$ are 0.023 Å and 0.029 eV/atom, respectively, which are significantly lower than the 0.11-0.12 Å and 0.117 eV/atom reported in earlier work [xieScreeningNewQuaternary2023].

The MLRM s used in screening were trained on the data from HeuslerDB, supplemented with newly computed $T_{{\rm c}}$ and $E_{\mathrm{aniso}}$ values using optimized structures in HeuslerDB. Test set metrics are summarized in Table 1 and benchmarked against previously reported results. The MLRM for $\{\boldsymbol{m}_{i}\}$ achieved an $R^{2}$ score of 0.989. For comparison with prior studies that used total magnetization ($m_{\mathrm{total}}$) as the target property, our model yielded $R^{2}=0.986$ for $m_{\mathrm{total}}$, exceeding earlier values of 0.75–0.89, 0.82, and 0.927 [mitraMachineLearningApproach2023, liuMachineLearningAssisted2025, luExplainableAttentionCNN2025]. For $T_{{\rm c}}$, the model attained $R^{2}=0.91$ and classification accuracy of 0.91, both substantially higher than previously reported values of $R^{2}=0.76$ and 0.73 and accuracy of 0.73 [hilgersMachineLearningbasedEstimation2025, hirohataMachineLearningDevelopment2025].

To the best of our knowledge, no previous study has predicted phonon stability or $E_{\mathrm{aniso}}$ of Heusler compounds using ML. The effectiveness of the $\omega_{\mathrm{min}}$ and $E_{\mathrm{aniso}}$ models developed here is supported by the validation results in Sec. 2. While phonon stability could also be assessed using uMLIP combined with phonon calculation methods, our regressor-based approach is motivated by both efficiency and accuracy. Only the minimum frequency is needed for phonon stability assessment, thus a regressor is sufficient and much faster than calculating the full spectrum. We also evaluated uMLIP +phonon using CHGNet, MatterSim, and eSEN-30M-OAM on a test set of 1000 conventional ternary compounds, and found stability/instability classification accuracies of 62.5%, 74.9%, and 80.2%, respectively, which are substantially below the regressor’s 93.6%. Notably, CHGNet misclassifies 67.2% of stable compounds as unstable, showing a strong tendency to underestimate stability, while eSEN weakly overestimates it and MatterSim exhibits a more balanced performance. Additional details including example phonon spectra by uMLIP s and DFT are in the Supplementary Information. The $E_{\mathrm{aniso}}$ model achieved an $R^{2}$ of 0.68, lower than those of $\{\boldsymbol{m}_{i}\}$, $\omega_{\mathrm{min}}$, and $T_{{\rm c}}$, highlighting the higher sensitivity and complexity of $E_{\mathrm{aniso}}$ as a target property. Nevertheless, despite the reduced $R^{2}$, its classification accuracy remains satisfactory and sufficient for integration into the ML-HTP workflow.

To benchmark advances in ML techniques since 2023, we applied the ALIGNN-FF uMLIP and ALIGNN MLRM to identify conventional quaternary candidate compounds and evaluated the validation rates of strong $E_{\mathrm{aniso}}$ compounds [choudharyUnifiedGraphNeural2023, choudharyAtomisticLineGraph2021]. In this test, the scalar quantity $\sum|\mathbf{m}_{i}|$ was used directly as the target property rather than being calculated from $\{\boldsymbol{m}_{i}\}$ prediction. The metrics of ALIGNN-FF and ALIGNN MLRM are summarized in Table 1. Using all screening thresholds, only 17 compounds qualified as candidates. To improve statistical robustness, we removed the phonon stability criterion, expanding the candidate list to 107 compounds, of which 26 (24.3%) exhibit |$$E_{\mathrm{aniso}} $|>1\mathrm{MJ/m^{3}}$. While this yield is notably higher than the 7.9% obtained from direct DFT-HTP screening, it remains far below the 82.0% success rate achieved with the eSEN-based ML-HTP workflow. These results highlight the substantial improvements in screening precision enabled by the state-of-the-art eSEN model.

We further tested a hybrid workflow in which structure optimization was performed with eSEN-30M-OAM uMLIP, while property prediction was carried out using ALIGNN MLRMs. This approach yielded 276 candidate compounds, of which 149 (54.0%) were confirmed by DFT to exhibit strong $E_{\mathrm{aniso}}$. The improved yield relative to ALIGNN-FF based optimization underscores the critical importance of accurate structure optimization with eSEN-30M-OAM for enhancing ML-HTP screening. However, the yield still falls short of the 82.0% achieved by the fully eSEN-based workflow, indicating that progress in both the uMLIP and MLRM components is essential for maximal efficiency. We also tested the inverse hybrid configuration, using ALIGNN-FF uMLIP for structure optimization combined with eSEN MLRMs for property prediction. This workflow identified 243 candidates, of which only 76 (31.3%) were validated as strong $E_{\mathrm{aniso}}$ compounds. This marked reduction in performance highlights the pivotal role of selecting an accurate uMLIP for structure optimization.

2.5 Prediction of local magnetic moment

Since the goal of this study is to identify compounds with large $E_{\mathrm{aniso}}$, it is first necessary to determine whether a compound is magnetic. Relying solely on the total magnetization is inadequate, as it cannot capture antiferromagnetic (AFM) or low-moment ferrimagnetic (FiM) compounds. To address this, we employed the total absolute magnetic moment, defined as $\sum_{i}|m_{i}|$, where $m_{i}$ denotes the local magnetic moment at atomic site $i$.

Local magnetic moments were predicted using an MLRM based on the eSEN architecture, trained to output the moment at each atomic site. Restricting to collinear configurations in which all moments are aligned along the $z$-axis, each moment is represented by a scalar whose sign encodes the direction, with $\ell=0$ per site in the output head. To account for the $z$-direction ambiguity—where a configuration and its sign-inverted counterpart (i.e., all local moments flipped) are physically equivalent—we modified the loss function to compute losses for both the predicted $\{\boldsymbol{m}_{i}\}$ and its sign-inverted counterpart, and take the smaller value as the loss. This ensures invariance under global spin inversion.

Fig. 5 (a) shows a scatter plot comparing $\{\boldsymbol{m}_{i}\}$ from the MLRM and DFT for compounds in test set, with histograms along the axes illustrating their distributions. Local moments at the atom site $Z$ are omitted as they are nonmagnetic. Because 73.4% of the test compounds are nonmagnetic, both histograms exhibit a pronounced peak at zero. For magnetic systems, the global sign is adjusted so that the total magnetic moment is positive; since most magnetic compounds are ferromagnetic, positive moments dominate in the distribution. Nearly all points fall along the diagonal and only 1.4% of points lie along either axis, demonstrating that the model accurately predicts both the magnitude and sign of local moments, and reliably distinguishes ferromagnetic (FM) and ferrimagnetic (FiM) systems. For the magnetic/nonmagnetic classification, we evaluated performance using receiver operating characteristic (ROC) and precision-recall (PR) curves, as shown in Fig. S6. The area under the curve (AUC) values are 0.98 and 0.97, respectively, indicating highly accurate classification.

Predicting $\{\boldsymbol{m}_{i}\}$ is a common yet computationally demanding step in HTP studies of magnetic materials. The approach developed in this work achieves high accuracy in $\{\boldsymbol{m}_{i}\}$ prediction and is readily transferable to other systems. A central question, however, is how many training compounds are required to reach satisfactory accuracy. To address this, we performed a learning-curve analysis by training the MLRM on progressively larger subsets of the dataset and evaluating performance on a fixed test set of 5000 compounds, of which 1486 are magnetic. For evaluation, local moments at the $X_{1}$, $X_{2}$, and $Y$ sites were concatenated across samples into a single array, while the $Z$ site was excluded since it is nonmagnetic in conventional Heusler compounds.

The learning curve is shown in Fig. 5 (b), illustrating how model performance improves with increasing training set size. In the top panel, the $R^{2}$ scores for both local moments and their magnitudes are presented. The gap between the two curves indicates that, while the model generally captures the magnitude accurately, it sometimes assigns the incorrect sign. For example, Mn2ScGe with DFT-computed local moments {$2.62$, $3.03$, $-0.29$, $-0.10$} $\mu_{\mathrm{B}}$ is predicted as {$-2.71$, $3.04$, $0.02$, $-0.02$} $\mu_{\mathrm{B}}$ by MLRM. The bottom panel reports two key metrics: (i) classification accuracy for identifying magnetic systems and (ii) the fraction of local moments with absolute error below 0.1 $\mu_{\mathrm{B}}$, evaluated across all compounds and within the magnetic subset. With 5000 training samples, the model achieves a classification accuracy of 0.92, and 90% of all compounds fall within the 0.1 $\mu_{\mathrm{B}}$ error threshold. However, this fraction decreases to 72% when restricted to magnetic compounds subset, indicating that the model identifies whether a site is magnetic with high reliability but remains less accurate in predicting exact $\{\boldsymbol{m}_{i}\}$ values. Increasing the training set to 125,000 samples improves this fraction to 82%, while relaxing the threshold to 0.2 $\mu_{\mathrm{B}}$ further raises it to 92%. Although performance improves with larger datasets, the gains become progressively smaller. These results highlight the critical role of dataset size in improving $\{\boldsymbol{m}_{i}\}$ accuracy and inform the selection of training size in future work focusing on other magnetic systems.

2.6 Frozen transfer learning for MLRM construction

To improve the performance of the MLRM, we employed a frozen transfer learning strategy using the eSEN-30M-OAM uMLIP as the base model. The eSEN-30M-OAM uMLIP was trained on the OMat, MPtrj, and sAlexandria databases, providing a comprehensive training set spanning the periodic table [barroso-luqueOpenMaterials20242024, jainCommentaryMaterialsProject2013, dengCHGNetPretrainedUniversal2023, schmidtMachineLearningAssistedDeterminationGlobal2023, schmidtImprovingMachinelearningModels2024]. Through this training, the embedding and the first several layers learn general chemical and structure representations. To leverage this, we transferred the parameters of the embedding and the first several layers into our MLRM and kept them fixed (frozen layers), while only updating the remaining layers and the output layer (flexible layers). This approach is analogous to a recent work in which the initial layers of ORB, EqV2, or MACE uMLIPs were used to generate feature vectors that were subsequently passed to property prediction models such as MODNet, XGBoost, and MLP [kimLeveragingNeuralNetwork2025a].

The eSEN-30M-OAM model consists of 10 layers. Fig. 6 (a) shows the $R^{2}$ scores of models trained with different numbers of frozen layers, denoted TL-uMLIP-$n$. In the $n=0$ case, the embedding layer is also flexible. Results for $\omega_{\mathrm{min}}$, $T_{{\rm c}}$, and $E_{\mathrm{aniso}}$ are presented. Model performance improves as the number of frozen layers increases up to $n=7$, after which it declines when more layers are frozen. This trend reflects the balance between transferring knowledge from the base model and maintaining sufficient flexibility to adapt to the new task. Based on this analysis, TL-uMLIP-7 was used in the ML-HTP case study for $\omega_{\mathrm{min}}$, $T_{{\rm c}}$, $E_{\mathrm{aniso}}$, and $\{\boldsymbol{m}_{i}\}$ MLRM. For comparison, a model trained from scratch (w/o TL) is included, which yields lower $R^{2}$ scores and highlights the benefit of transfer learning. We also tested a variant initialized from a base model pre-trained on formation energy ($\Delta E$) data from the HeuslerDB database (denoted TL-$\Delta E$-0). TL-uMLIP-0 outperforms TL-$\Delta E$-0, indicating that initialization from the eSEN-30M-OAM model, trained on a comprehensive dataset, is more advantageous.

Besides improving overall performance, transfer learning from uMLIP can significantly enhance model generalization to unseen elements. To demonstrate this, we evaluated model performance using group-wise split analysis. For each group-wise split test, three neighboring elements within the same period were used as the holdout elements. Any compound containing at least one of the holdout elements was reserved as the holdout test data; the remaining compounds were split into train/val with a 9:1 ratio. This setup simulates realistic scenarios where models are applied to material systems containing elements not present in the training data, a typical domain shift in materials science. For direct comparison, we also measured performance using train/val/test sets created using random splits while keeping the counts of train/val/test the same as in the group-wise split. We tested six sets of d-block elements and four sets of p-block elements as holdout sets. We used $T_{{\rm c}}$ as the target property for efficiency. Results obtained using frozen transfer learning and models trained from scratch are shown in Fig. 6 (b).

Using FTL, performance on group-wise splits is generally lower than on random splits, indicating that the model performance on unseen elements is reduced. The drop in performance for p-block holdouts is smaller, which is consistent with the fact that p-block elements are typically nonmagnetic in these compounds and thus less directly related to the target property $T_{{\rm c}}$. For group-wise split tests with models trained from scratch, a much larger performance decrease is observed. This occurs because FTL transfers the general chemical representation learned by uMLIP to property regressors, whereas models trained from scratch do not benefit from this prior knowledge. This highlights the advantage of transfer learning from uMLIP for improving generalization to unseen elements. Variation in random-split performance is due to different train set sizes since the number of holdout compounds varies with the chosen holdout elements.

It should be noted that TL is not always effective. It can be limited when the representations learned by base model do not generalize to the target domain. For example, source and target may differ substantially in material class (e.g. bulk versus surface or molecular systems), or downstream task concerns different physics (e.g. static energetics versus dynamic properties) [yamadaPredictingMaterialsProperties2019a, chenGraphNetworksUniversal2019a, changOvercomingDataScarcity2022]. Furthermore, a recent work, which transfer from model trained using computational data to learn experimental data, revealed that the error of the transfer-learned model decreases according to a power-law as the size of the computational data for the base model increases [minamiScalingLawSim2Real2025]. This suggests that transfer learning may also be less effective when the base model is trained on a small dataset even if the base model performance within its own domain is good.

3 Discussion

We demonstrated the feasibility of combining uMLIP s and MLRM s for HTP screening. As a case study, we identified 334 conventional quaternary Heusler compounds and 924 all-$d$ Heusler compounds that exhibit thermodynamic, dynamical, and magnetic stability, together with large $E_{\mathrm{aniso}}$. The precision of this workflow was confirmed through DFT validation of the candidate list.

For other material systems, if database for a subset is available, the same workflow can be applied to explore the remaining chemical space. For novel materials, uMLIP can be used directly to reduce the search space by filtering for thermodynamic stability, while target properties need to be computed for a representative subset of compounds to train MLRMs for ML-HTP screening of the broader space.

The MLRM s in this study involve a domain shift, as the train dataset does not contain quaternary or all-$d$ Heusler compounds, yet the models are applied to these systems. While performance is satisfactory, further gains are possible using active learning approach. Iteratively refining the models by selecting informative compounds can improve performance and help identify candidates missed by the current one-shot approach. In such iterative framework, current DFT validation results can be added to the training set, forming the first iteration of model refinement. This iterative strategy naturally extends the present work, enabling more thorough exploration of chemical space while progressively improving model accuracy. Such an approach is especially valuable when the target material space differs substantially from the training data.

In the framework of this study, uMLIP perform the critical task of structure optimization, which is traditionally handled by DFT-derived methods, while MLRMs predict the properties of the optimized structures, tasks that are typically carried out using DFT or DFT-derived methods. Together, uMLIP and MLRMs enable a drop-in replacement for DFT in conventional HTP workflows. This replacement is general and can be readily extended to other properties, material classes, and DFT-based HTP pipelines, enabling accelerated HTP screening and discovery of new materials.

4 Methods

4.1 ML-HTP workflow

The ML-HTP workflow is schematically illustrated in Fig. 1(b) and described in detail below.

In step (a) of ML-HTP workflow, the structures were optimized and the $\Delta E$ and $\Delta H$ were calculated using uMLIP. The initial lattice constant was estimated as the average value of known $X_{2}YZ$ Heusler compounds in the HeuslerDB database that share two elemental species with the target composition. A conventional cell in cubic phase with this estimated lattice constant was then constructed. To generate initial structures, the lattice parameters of this cell were systematically varied: the $a$, $b$, and $c$ were uniformly scaled by $\pm 10\%$ or $\pm 30\%$, and, alternatively, the $c$ alone was scaled by $\pm 10\%$ or $\pm 30\%$. All generated structures were subsequently converted to the primitive unit cell and relaxed using the eSEN-30M-OAM uMLIP. The structure with the lowest energy after relaxation was identified as the ground state. The choice of eSEN-30M-OAM was motivated by its superior performance relative to other uMLIPs, as discussed in Sec. 2.3. The selection of initial structures was validated by a convergence test, which is provided in Fig. S4.

For a compound to be thermodynamically stable against decomposition into its constituent elements or competing phases, the formation energy must be negative ($\Delta E<0$) and the distance to the convex hull must be zero ($\Delta H=0$). In practice, however, metastable phases ($\Delta H>0$) at 0 K may become stabilized at finite temperature [sunThermodynamicScaleInorganic2016]. Following our previous work, we adopt a practical stability criterion of $\Delta E<0.0$ eV/atom and $\Delta H<0.22$ eV/atom, which has been shown to effectively capture experimentally accessible compounds [xiaoHighthroughputComputationalScreening2025]. Using the energies of ground state candidates, constituent elements and competing phases predicted by eSEN-30M-OAM, we computed $\Delta E$ and $\Delta H$ to assess thermodynamic stability.

In step (b) of workflow, the $\omega_{\mathrm{min}}$, $\{\boldsymbol{m}_{i}\}$, $T_{{\rm c}}$, and $E_{\mathrm{aniso}}$ were predicted with MLRMs trained on the HeuslerDB and additional computed data. The uMLIP-optimized structures were used as inputs. The construction of these MLRMs is described in Sec. 2.5 and Sec. 2.6. Compounds were classified as dynamically stable for $\omega_{\mathrm{min}}$ $>-10$ cm*-1*, magnetic for $\sum_{i}|\boldsymbol{m}_{i}|>0.1~\mu_{\mathrm{B}}$/f.u., magnetically stable for $T_{{\rm c}}$ $>300$ K, and strong $E_{\mathrm{aniso}}$ candidates for |$$E_{\mathrm{aniso}} $|>1$ MJ/m3. Tetragonal compounds satisfying all of these conditions were designated as promising stable materials with strong $E_{\mathrm{aniso}}$.

In step (c) of workflow, the candidate list was validated with DFT calculations to assess the reliability of the ML-HTP workflow. Structure optimizations were performed starting from various initial spin configurations, consistent with our previous DFT-HTP work. For conventional Heusler compounds, the magnetic moments at the $X_{1}$, $X_{2}$, and $Y$ sites were initialized in configurations where they were either parallel or antiparallel to each other. For all-$d$ Heuslers, spin arrangements on all four sites were considered. To capture possible high-spin and low-spin states, two initial magnitudes of the local moments ($|\boldsymbol{m}_{i}|=1$ and $4~\mu_{\mathrm{B}}$) were tested, along with a nonmagnetic configuration ($|\boldsymbol{m}_{i}|=0$). The uMLIP-optimized structures served as the starting geometries. After structure relaxation, the ground state was identified by comparing total energies. For the resulting ground states, we computed $\Delta E$, $\Delta H$, $E_{\mathrm{aniso}}$, phonon, and $T_{{\rm c}}$ using VASP, OQMD, ALAMODE, and SPRKKR [saalMaterialsDesignDiscovery2013, bahnObjectorientedScriptingInterface2002, tadanoSelfconsistentPhononCalculations2015, AnharmonicForceConstants, ebertCalculatingCondensedMatter2011, liechtensteinLocalSpinDensity1987].

4.2 Computational methods

The uMLIP-based structural optimizations were performed using the Atomic Simulation Environment (ASE) package [hjorthlarsenAtomicSimulationEnvironment2017]. The Fast Inertial Relaxation Engine (FIRE) optimizer was employed, with symmetry constraints enforced throughout the relaxation process [bitzekStructuralRelaxationMade2006]. To ensure consistency in reference energies for computing $\Delta E$ and $\Delta H$, the elemental phases and competing phases were also optimized using the same uMLIP. Their initial geometries were taken from DFT-optimized structures in the Open Quantum Materials Database (OQMD) database [saalMaterialsDesignDiscovery2013, bahnObjectorientedScriptingInterface2002].

The MLRMs were developed to predict $\omega_{\mathrm{min}}$, $\{\boldsymbol{m}_{i}\}$, $T_{{\rm c}}$, and $E_{\mathrm{aniso}}$. To leverage prior knowledge, we employed a frozen transfer learning strategy, as illustrated in Fig. 1(a). Each MLRM was initialized from the pretrained eSEN-30M-OAM uMLIP, with the embedding layer and the first seven message-passing layers kept frozen, while the final three layers and the output layer were fine-tuned. This framework was implemented using a modified version of the FAIRChem package (v1) [shuaibiFAIRChem2025]. For the $\{\boldsymbol{m}_{i}\}$ MLRM training, the loss function was adapted to address the global sign ambiguity of magnetic moments, as described in Sec. 2.5. The training dataset consisted of HeuslerDB together with newly computed $E_{\mathrm{aniso}}$ and $T_{{\rm c}}$ values based on the optimized structures in HeuslerDB. The data were randomly partitioned into training, validation, and test sets in an 8:1:1 ratio. For $\{\boldsymbol{m}_{i}\}$, we used all ground-state entries, resulting in 27,864 data points. The $\omega_{\mathrm{min}}$ data were available for all thermodynamically stable ground states, yielding 8,198 entries. For $T_{{\rm c}}$, 2,106 data points were used, including 750 newly computed values. Since $E_{\mathrm{aniso}}$ data were not included in HeuslerDB, we calculated $E_{\mathrm{aniso}}$ for all magnetic tetragonal ground-state systems, obtaining 6,123 entries.

DFT calculations were performed primarily with Vienna ab initio Simulation Package (VASP) [kresseEfficiencyInitioTotal1996, kresseEfficientIterativeSchemes1996], using the projector augmented wave (PAW) method and the generalized gradient approximation (GGA) with the Perdew–Burke–Ernzerhof (PBE) functional [kresseUltrasoftPseudopotentialsProjector1999, perdewGeneralizedGradientApproximation1996]. $\Delta E$, $\Delta H$, phonon, and $T_{{\rm c}}$ were computed using OQMD, ALAMODE, and SPRKKR following the methodology of our previous DFT-HTP study of ternary Heusler compounds [saalMaterialsDesignDiscovery2013, bahnObjectorientedScriptingInterface2002, tadanoSelfconsistentPhononCalculations2015, AnharmonicForceConstants, ebertCalculatingCondensedMatter2011, liechtensteinLocalSpinDensity1987, xiaoHighthroughputComputationalScreening2025]. The agreement of magnetic moments between VASP and SPRKKR is demonstrated in Fig. S7. Phonon calculations for the all-$d$ compounds MnOsMnRe and MnReMnRu failed due to convergence issues in DFT, and these compounds were treated as unstable in the validation rate analysis. $E_{\mathrm{aniso}}$ was calculated as $E_{\mathrm{aniso}}$ $=E_{\perp}-E_{\|}$ using the force theorem [daalderopFirstprinciplesCalculationMagnetocrystalline1990, xingLatticeDynamicsIts2022, xingFirstprinciplesPredictionPhase2023]. Calculations were performed in the primitive cell with $\boldsymbol{k}$-meshes generated using Python Materials Genomics (pymatgen) at a density of $6000\text{\,}\mathrm{\text{Å}}^{-3}$, and the tetrahedron method with Blöchl corrections was applied [ongPythonMaterialsGenomics2013, blochlImprovedTetrahedronMethod1994]. Input generation, structure manipulation, and symmetry analysis were carried out using pymatgen, ASE, ASE2SPRKKR, and spglib [ongPythonMaterialsGenomics2013, hjorthlarsenAtomicSimulationEnvironment2017, ASE2SPRKKRSoftwarePackage, bahnObjectorientedScriptingInterface2002, togoSpglibSoftwareLibrary2024].

4.3 Computational cost

The computational cost for each task in the HTP workflow is summarized in Table 2. DFT timings are based on validation runs for ML-selected candidates, while ML timings are taken from the ML-HTP screening. For DFT, we report per-job statistics such as mean, median, and interquartile range (IQR). uMLIP structure optimization and MLRM predictions were performed in batches, and mean values are reported as individual per-job timings are not available. MLRM training times correspond to the wall time required to train a single model on one GPU. DFT calculations except phonon were performed on dual-socket Intel Xeon Platinum 8268 (Cascade Lake, 24 cores per CPU, 2.9 GHz, 48 cores per node), phonon calculations were performed on Fujitsu A64FX processors (Armv8.2-A SVE 512 bit, 48 compute cores per CPU, 2.2 GHz), uMLIP structure optimizations were performed on dual-socket Intel Xeon Gold 6230 (Cascade Lake, 20 cores per CPU, 2.1 GHz, 40 cores per node), and MLRMs were trained and applied on NVIDIA RTX 6000 Ada GPUs. All DFT calculations, except for phonon calculations, were performed on 2 nodes; phonon calculations were performed on 6 nodes. The uMLIP structure optimizations were performed on 1 node and MLRM training and prediction were performed on 1 GPU.

In the ML‑HTP case study, uMLIP structure optimization, MLRM training, and MLRM prediction used 1,835 node‑hours, 23 GPU‑hours, and 4 GPU‑hours, respectively. DFT calculations of ML‑selected candidates consumed 1,645 node‑hours for optimizations and 99,680 node‑hours for property evaluations. Ignoring the node difference, the total ML‑HTP cost is 103,160 node‑hours plus 27 GPU‑hours. A DFT-HTP workflow screening the same chemical space would need about 256,160 node‑hours for structure optimizations and 18,802,625 node‑hours for property evaluations, estimated using statistics reported in Table 2. This comparison shows an estimated speed‑up of 185 times for ML‑HTP with DFT validation, or 104 times if used without DFT validation at the end. Please note that these estimates are approximate, and waiting time for job execution, human time spent on debugging and workflow management are not included. The actual speed-ups also vary based on material systems, target properties, and computational resources.

Data availability

The ML-HTP candidate list and DFT validation results are included in the Supplementary Information as Tables S3 and S4. The complete set of all screened compounds, along with ML-predicted properties, will be made available through the HeuslerDB database at https://www.nims.go.jp/group/spintheory/.

Code availability

The developed packages MLIP-HOT and MLIP-FTL will be made available through the Spin Theory Group GitHub repository at https://github.com/nims-spin-theory and our group website at https://www.nims.go.jp/group/spintheory/.

Funding

This study was supported by MEXT Program: Data Creation and Utilization-Type Material Research and Development Project (Digital Transformation Initiative Center for Magnetic Materials) Grant Number JPMXP1122715503 and as “Program for Promoting Researches on the Supercomputer Fugaku” (Data-Driven Research Methods Development and Materials Innovation Led by Computational Materials Science, JPMXP1020230327).

Acknowledgments

This study used computational resources of supercomputer Fugaku provided by the RIKEN Center for Computational Science (Project ID: hp240223), the computer resources provided by ISSP, U-Tokyo under the program of SCCMS, and the computer resources at NIMS Numerical Materials Simulator.

Author contributions

TT conceptualized, designed, and supervised the project; reviewed and edited the manuscript. TT and EX developed the methodology and code implementation; performed the calculations and analysis; EX drafted the manuscript.

Competing interests

All authors declare no financial or non-financial competing interests.