We first described the workflow of our ML code, as shown in Fig. 1. Stable 1617 samples were extracted from the C2DB dataset, of which 1292 samples were labeled as nonmagnetic 2D materials and 325 as magnetic. This imbalanced dataset could not be used to train the classification models. Therefore, 800 2D magnetic materials were randomly selected from the 2DMP database. The new dataset included 2417 samples that were used to pre-train the classification model for magnetic/nonmagnetic 2D materials. After obtaining model N-M, the 615 samples with magnetic moments of > 1μ_B were fed to high-throughput calculations to train the FM-AFM classification model. In this step, an ferromagnetic/antiferromagnetic (FM/AFM) classification model and a dataset of 2D materials with precise magnetic moments were obtained. In the next step, the AFM 2D materials were discarded and a regression model of magnetization (model-R) was trained. 545 stable magnetic samples were selected from the recently updated C2DB and duplicate samples removed from the previous database. Non-prior samples were adopted as a validation set to predict new materials. Finally, 2D materials with magnetic moments of > 7μ_B were selected and defined as strongly magnetized materials. Our hierarchical filtering process does not depend on empirical parameters or physical intuition, and the database is updated in real time. Finally, the interpretability of the model was analyzed. Our ML predictions and DFT calculations did not target the Curie point since Kabiraj et al. have already calculated this through graphics processing unit (GPU) acceleration²⁸. Additionally, the C2DB database does not contain the Curie point, and no experimental values have been reported for the magnetic transition critical temperatures of most 2D materials.

The ultimate goal of the current work was to develop novel 2DFM materials with strong magnetization. Therefore, a two-step ML model was constructed. The first step of ML was to classify 2D materials into magnetic and nonmagnetic groups according to their class labels of magnetic states in the C2DB and 2DMP databases, and data-M-N was adopted for this task. Model M-N performs binary classification to distinguish magnetic and nonmagnetic 2D materials from the dataset data M-N. The following four algorithms were selected to train this binary model: Bernoulli Naive Bayes (BNB), K-Nearest Neighbor (KNN), Support Vector Machine Classifier (SVC), and Random Forest Classifier (RFC). All the ML algorithms used in the present work were implemented through the Scikit-learn library (version 1.1.3)⁴⁴. Since the dataset used in the first model was of medium size, the training and test sets were randomly split at a ratio of 7 : 3. The ML models were trained with the adoption of 2D material data in the training set, and the generalization ability and accuracy of the ML model were tested using 2D material data in the test set. A 10-fold cross-validation approach was used to optimize the hyperparameters in all the models. The performance metrics of the test set for all the four models are presented in Table 1. As shown, KNN and BNB exhibit similar metrics. SVC with the radial basis function (rbf) kernel performed better than KNN and BNB. The RFC classifier achieved the best performance for all the metrics. The receiver operating characteristic (ROC) curves and typical confusion matrices for the four classification tasks are shown in Fig. 2. The optimized hyperparameters for all the models are presented in Table S2, and the Supplemental Information (SI) details the classification techniques utilized. The classification area under the curve (AUC) value of the RFC model for the non-magnetic/magnetic classifications exceeded 90% after re-examination of distinguishable samples. The confusion matrix for each classifier presents the counts of the classes predicted using the four algorithms versus the true classes of the test set. Excellent performance was achieved, particularly for the RFC classification, with an AUC value of 0.94.

To further verify the generalizability of the model, a dataset containing 10 2D materials was constructed: CrI₃, CrGeTe₃, PrTe₃, TiPb₉O₁₁, TaSe, InCl₃, TiGeTe₆, KB(CO₂)₄, Ag₂WS₄, and NaFeAs. In this small validation set, CrI₃ and CrGeTe₃ are known 2D magnetic materials, and the other eight materials are labeled as nonmagnetic materials and are taken from 2DEP. The test results are presented in Table S3, where the boolean values “True/False” represent magnetic and nonmagnetic materials, respectively. The four algorithms accurately classified CrI₃ and CrGeTe₃ as magnetic materials⁴⁵. However, the most accurate classification model remained the RFC algorithm. Once again, our goal was to predict strongly magnetized ferromagnetic materials. Therefore, at this stage, DFT calculations were not conducted to determine whether the materials that have not been reported have magnetism or not.

In the second step, the FM/AFM classification of the selected magnetic candidates was performed with the adoption of model-F-A with four ML algorithms. The second step included two sub-steps: high-throughput DFT calculation and FM/AFM classification. After simple screening of the old C2DB database, 615 samples exhibit thermodynamic and dynamic stability, along with magnetic moments > 1μ_B. However, as mentioned previously, the 615 samples were all labeled as FM. The calculation results in the C2DB database were based on a unit cell, and most of them were obtained using Generalized Gradient Approximation (GGA) and Heyd Scuseria Ernzerhof (HSE) approximations. This resulted in the fact that all materials being incorrectly calculated as ferromagnetic, whereas the antiferromagnetic materials were completely ignored. In addition, for magnetic systems, the GGA may underestimate the band gap and magnetic moment of 2D strong-correlation materials⁴⁶. To obtain an accurately labeled dataset of FM/AFM, the 615 2D candidates were re-calculated using the GGA+U approach. The detailed Hubbard U values of the TM are presented in Table S1, where only the one on element Tc is null. A 2 × 1 × 1 super-cell was built for the unit cell, which included only one magnetic atom. In our calculations, AFM represented the spin-antiparallel arrangement of two magnetic atoms, and more complex magnetic structures were not taken into consideration. Structural relaxation was performed on all materials, and the free energies of FM and AFM materials were calculated. The energy differences between the FM and AFM are presented in Table S4. The calculation results indicated that 357 materials are ferromagnetic, 223 antiferromagnetic, 26 exhibit ferromagnetic energy equal to the antiferromagnetic energy, and nine materials are quinary compounds that were discarded in our ML model. 26 materials with equal ferromagnetic and antiferromagnetic energies were divided into antiferromagnetic sets. Therefore, the ratio of the number of samples in the FM and AFM sets was 1.43 : 1. The 2DFM materials were further classified into diverse prototypes according to their structural symmetries and space groups. Among the 357 2DFM materials, the 10 most common ones, representing all 300 structures, are shown in Fig. 3. The most common space group was that of P-3m1, which corresponded to 77 similar structures including many MXenes and TM dichalcogenides. It is found that most 2DFM materials have not been reported, and their structural prototypes are novel and compelling⁴⁷. Moreover, it is emphasized that the Hf₂Br₆ with P-62m space group and the Os₂Cl₆ with P6/mmm space group both have a TM in the sixth period; they have not been reported before and are rare in 2D ferromagnetism.

After high-throughput DFT calculations, the 606 samples obtained were used to train the FM and AFM classification model F-A. Since these dataset were small and the performance of the generalization ability of training sets on test sets may be affected by a particular random training-test set split, the train and test set were randomly split at an ratio of 8:2. Hyperparameter optimization for the models involved a 10-fold cross-validation technique. The main metrics for the FM/AFM classification of the test set are presented in Table 2. The KNN and RFC classifiers performed significantly better than the SVC and Gaussian Naive Bayes (GNB) classifiers. To increase the classification accuracy of the Naive Bayesian algorithm, the Bernoulli Naive Bayesian was replaced with the Gaussian Naive Bayesian; however, the accuracy did not increase. Again, the RFC classifier achieved an excellent test-set accuracy of 0.918. For the magnetic/nonmagnetic classification models, the four metrics were relatively average. However, for the FM/AFM classifier, the precision was higher than the accuracy, and the recall was lower than the precision for the KNN and RFC classifiers. The ROC curves and typical confusion matrices for the four classification tasks are shown in Fig. 4. The optimized hyperparameters for all the models are presented in Table S5. The classification AUC value of the RFC model for the FM/AFM classifications exceeded 95% after re-examination of distinguishable samples. The confusion matrix for each classifier presents the counts of the classes predicted using the four algorithms versus the true classes of the test set. Excellent performance was achieved, especially for RFC classification, with an AUC value of 0.97.

In orderto verify the generalization ability of model F-A and quickly predict new ferromagnetic materials, 542 samples were selected and labeled as magnetic in the updated C2DB database and they were screened to obtain 463 2D materials with a magnetic moment of > 1μ_B. Compared with the previous dataset of 615 samples, duplicate compounds and quinary compounds were removed, resulting in 167 candidates. This small dataset were used as a validation set and imported into the RFC model F-A for classification. It is predicted that 60 of these compounds are ferromagnetic, as shown in Table S6.

Feature selection is critical for ML models. Fig. 5 shows the sorting results for the material features of the RFC model. As indicated by Fig. 5a, the following were among the top 10 features in the dataset examined: (i) the electronegativity, covalent radius, Mendeleev number, and melting (“mode electronegativity”, “mean electronegativity”, “Covalent radius”, “Minimum MendeleevN”, “Range MeltingT”, and “Maximum MeltingT”), which are correlated to the elemental properties; (ii) the number of the space group (SGnumber), which is linked to the structural symmetry; and (iii) GSvolume_pa, NpUnfilled, and NUnfilled, which affect the contributions of different elements in the materials to the top of the valence band, thus affecting the electronic band structure and magnetic properties of the compounds. As shown in Fig. 5a, a magnetic/nonmagnetic classification model (model-M-N) with high generalization ability and scores can be obtained with the adoption of only element and component features, even without considering structural features as descriptors. However, after testing, the same set of features was used to train the AFM/FM classification models, which resulted in extremely low scores for each metric. This is originated from the fact that the feature vector of a single TM element cannot determine the magnetic order of the system, and both ferromagnetic and antiferromagnetic interactions are associated with the exchange interactions of the adjacent magnetic atoms. The types of TM ions, the corresponding ligand atoms, and their ionic radii lead to the formation of different types of crystal fields in the compounds. As reported by Wang et al. the crystal fields are an important feature used to describe magnetic ordering²⁶articularly the ability to classify FM and AFM materials. Therefore, during the training of model F-A, the properties of the ligand atoms were introduced as feature vectors. We also introduced structural features as well as 20 atom centered symmetry function (ACSF) features associated with structures. Fig. 5b shows the sorting results for the material features of the RFC model F-A. All the elemental features and their descriptions (spin magnetic moments, ionic radius, covalent radius, dipole polarization, formation heat, electronegativity, and vdw radius) were obtained from the Python Mendeleev package 0.11.0³⁴. A simple summation process was performed on the elemental features of TM atoms and their coordinate atoms, and the d orbital electron number of the TM was added to the feature sets. As shown in Fig. 5b, the sum of the ionic radii of the TM and coordination atoms and the ionic radii of the coordinate atoms were the most important features. In contrast, the ACSF did not have a significant impact on the nature of the target materials.

According to our high-throughput calculations, most 2D magnetic materials have not been reported, and a special one is Cr₂NF₂, which has a negative energy difference between FM and AFM states, indicating that the ferromagnetism is stable. However, according to our recent work, Mo₂NCl₂ is a zigzag-type AFM^48,49. This is ascribed to the fact that, in our high-throughput calculations, the FM and Néel AFM orders were simply considered, whereas the ground-state magnetic structure of this functionalized MXene material is a zigzag-type AFM state. Therefore, three magnetic structures were selected, as shown in Fig. S1, and it is found that Cr₂NF₂ tends to be a potential antiferromagnetic/ferroelectric multiferroic material. To estimate the Nèel temperature (T_N) of the ferrimagnetic Cr₂NF₂, Monte Carlo (MC) simulations based on the 2D Heisenberg Hamiltonian model were performed. Details regarding the MC simulations are provided in the SI. From the calculated curve of C_v, as shown in Fig. 6a, it is found that Cr₂NF₂ had a fairly high T_N of 169 K. The non-zero magnetic moment at 0 K reflects the antiferromagnetic behavior of the ground state, as shown in Fig. 6b. Cr₂NF₂ has a higher magnetic transition temperature than CrI₃ (T_C = 45 K) and Cr₂Ge₂Te₆ (T_C = 30 K)^50,51, suggesting its potential applications in nano-spintronic devices. In addition, the ferroelectric polarization switching process was also simulated with the adoption of the nudged elastic band (NEB) approach, as shown in Fig. 6c. As the paraelectric state is dynamically unstable, the energy profile is a double-well curve. Two degenerate polarized states pass through a saddle point (0.276 eV).

Most reports on magnetic moment prediction lay great emphasis on the atomic magnetic moment of each TM⁵², or the magnetic moment per atom⁵³. So it is hard to believe that it fully reflects the total magnetization of the system. According to model F-A and the high-throughput DFT results, there were 383 candidates, including 357 ferromagnetic and 26 energy-degenerate FM/AFM states, which were used to train the regression model of the net magnetic moments of the unit cell. A material information platform (Matminer) was used to describe elemental features of the materials in the dataset⁵⁴. Elemental properties mainly include elemental information, the electronic configuration, and the material composition. In addition, we added 20 sine Coulomb matrix eigenvalues and 116 SOAP features associated with the structure. Thus, 273 features were included in the dataset. The input of high-dimensional features not only leads to the risk of overfitting but also models inefficiency. The most important features were evaluated and retained using 10-fold cross-validation recursive feature elimination (RFECV). Through this procedure, 24 features were obtained, which were used for the subsequent model training and testing. In Table S8, 24 selected significant features are shown, and the physical explanation of each feature is presented in Table S8. In addition, the number of unpaired d orbital electrons of the TM elements (N-e-unpaired) and the magnetic moment of a single TM atom calculated using DFT (M-single atom) were important features for training the regression model. Unlike classification models, during the training of regression models, three ensemble algorithms were selected, Random Forest Regression (RFR), Gradient Boosting Decision Tree (GBDT), eXtreme Gradient Boosting (XGboost), and an Artificial Neural Network (ANN) algorithm⁵⁵. The selected hyperparameters for the four models are presented in Table S7.

Fig. 7 shows the results of the training and testing of the net magnetic moment of the unit cell for the four regression models. As indicated by Fig. 7a to d, the training scores of all four models were high, indicating that the 25 features (24 optimized features and N-e-unpaired) are closely related to the net magnetic moment. The R² scores of the three ensemble and ANN models were all > 0.9, and the GBDT model exhibited the best performance, with R² = 0.945. It also exhibited the smallest Root Mean Squared Error (RMSE) and Mean Absolute Error (MAE) of 0.649μ_B and 0.430μ_B, respectively. The test results indicated that all the four models exhibit excellent generalization abilities, as shown in Fig. 7e to h. Among the four ML regression models, the highest performance was achieved by the GBDT model with the smallest MAE of 0.556μ_B, exhibiting superior prediction accuracy. The regression models trained using 25 features (24 optimized features and M-single atoms) are shown in Fig. S2. The training accuracy of the model containing N-e-unpaired features was similar to that of the model containing M-single-atom features. However, expensive and accurate DFT calculations are required by the M-single-atom features. Therefore, in our subsequent prediction of 2D materials with strong magnetization, the ML model used was a regression model with N-e-unpaired features.

Shapley Additive exPlanations (SHAP) value analysis is an approach for evaluating ML model interpretability. It can provide a clear graph of how diverse features compete with each other and determine the target property⁵⁶. Specifically, it can reveal the quantitative local contribution of each feature to the prediction target property of a single sample, which is difficult to explain by feature importance. Therefore, SHAP value analysis was implemented on the optimized dataset of the magnetic moment regression model. The SHAP values of the 25 most important elemental and magnetic features of the magnetic moment regression model are shown in Fig. 8. The exchange interaction between the nearest-neighbor (nn) and next-nearest-neighbor (nnn) atoms in 2D monolayer materials has a more significant effect on the magnetic moment than the elemental or atomic properties⁵⁷. Taking Cr₂C as an example, the crystal field affects the electronic structure, which in turn affects the magnetic interaction of magnetic systems, as shown in Fig. 8a. The sine Coulomb matrix, volume per atom (VPA), is related to the structural properties, which are also vital for the prediction of net magnetization, as shown in Fig. 8b. The local influence of the optimal representation set was also analyzed. Cr₂C, Cr₂CO₂, and Cr₃C₂O₂ were selected as representative materials, and their SHAP values for the preferred features are shown in Fig. 8c. The net magnetic moment is 6.92μ_B, 6.42μ_B, and 5.49μ_B per supercell for Cr₂C, Cr₂CO₂, and Cr₃C₂O₂, respectively. Among the preferred features for enhancing the net magnetic moment, the N-e-unpaired is the most positive contribution. Conversely, the “sine Coulomb matrix eig 3”, and mean number of unfilled electrons are the features with the most negative contributions. The N-e-unpaired of pristine and O₂ functionalized Cr₂C are similar in contribution. The contribution of N-e-unpaired increases from two to three layers of TM (Cr₂C to Cr₃C₂). Additionally, the proportion of VPA in Cr₂CO₂ is relatively high, and the maximum covalent radii of Cr₂CO₂ and Cr₃C₂O₂ gradually increase. This leads to a shorter distance between the two nn Cr³⁺ ions. Thus, the Cr-C-Cr FM super-exchange interaction in monolayer MXenes tends to be stronger. Owing to the complexity of 2D structures and the various factors affecting their magnetism, the SHAP value provides only limited information on the physical mechanism. Nevertheless, the SHAP value analysis makes our ML model for predicting 2D materials with strong magnetization more interpretable. Significantly, the SHAP method was employed for a quantitative analysis of the impact of ML-selected features on magnetic properties, thereby revealing the underlying physical insights of our models. This interpretable framework has the potential to unlock the “black box” of ML, which could lead to groundbreaking ML-aided material design advancements.

The aim of this section is to design 2D materials with large net magnetic moments for the aforementioned reasons. 60 novel 2DFM materials were predicted with the adoption of the ML classification model F-A, as shown in Table S6. We input these 60 2DFM materials into model-R for regression prediction of the magnetic moment. After the regression calculations, screening was performed to obtain 13 2DFM materials with magnetic moments > 7μ_B, as shown in Table 3. Thus, starting from an initial dataset of 60 ferromagnetic candidates, a small dataset of 13 samples was obtained, which was calculated to be stable and exhibit net magnetic moments of > 7μ_B. It is emphasized that the goal was achieved without high-throughput DFT calculations. Here, the anomalous point appearing in the regression model was re-examined, i.e., the magnetic moment prediction of Fe₄S₈, as shown in Fig. 9a.

An anomalous point also occurs in the regression model with a feature containing the magnetic moment of a single atom (M-single atom) such as the yellow dot in the rectangular box in Fig. S3. In the C2DB database, Fe₄S₈ is labeled as an FM state with a net magnetic moment of 3.759μ_B, which is clearly the magnetic moment of a single Fe atom. This is consistent with the GGA+U calculations. The unit cell used in the DFT calculations contained four Fe atoms, the calculated net magnetic moment was 15.253μ_B, with an average magnetic moment of 3.813μ_B for each Fe atom. However, the crystal symmetry of Fe₄S₈ is relatively low, and the crystal field is not a conventional octahedral, tetrahedral, or triangular prism. From a purely structural perspective, the nearest-neighbor Fe-Fe interaction is a direct exchange interaction. According to the interpretability analysis of the regression model, a shorter nearest neighbor direct exchange tends to exert a more significant impact on the net magnetic moment, resulting in a smaller net magnetic moment per unit cell. This is why anomalous prediction points appear in the ML regression models. In addition, the dynamic stability of Fe₄S₈ is relatively low, and further experimental confirmation of its magnetic properties is required.

In the predicted set, most of the 2D materials with large net magnetic moments contained Mn atoms. Os₂Cl₈ and Fe₃GeSe₂ are exceptions; their structures are shown in Fig. 9b and c, respectively. Os₂Cl₈ was classified as a magnetic state with a magnetic moment of 7.668μ_B in the C2DB dataset, which was smaller than the predicted value of 8.435μ_B. Os₂Cl₈ exhibits high thermodynamic and dynamic stability, which was also reported for the C2DB dataset. In addition, there are no relevant reports on Os₂Cl₈. Fe₃GeSe₂ has the same structural prototype as Fe₃GeTe₂, which was reported to be a ferromagnetic monolayer by Li et al⁵⁸. The net magnetic moment of Fe₃GeSe₂ is 6.4297μ_B, which is smaller than the predicted value of 7.3087μ_B. It is also found that the S-functionalized Mn-based MXene Mn₄N₃S₂ has the largest net magnetic moment 12.431μ_B among the sub-dataset, which is larger than that of C2DB (10.776μ_B). Both Fe₃GeTe₂ and Mn₄N₃S₂ exhibit high thermodynamic and dynamic stabilities. In order to verify the accuracy of the regression model, the magnetic and electronic properties of the three compounds were calculated with the adoption of DFT. As shown in Fig. 10, the projected density of states (DOS) and spin density of Os₂Cl₈, Fe₃GeSe₂, and Mn₄N₃S₂ were obtained by the GGA+U approach. The band gap of Os₂Cl₈ is 1.54 eV, which is larger than that of C2DB (0.21 eV). The net magnetic moment of the Os₂Cl₈ unit cell is 8.075μ_B. The predicted magnetic moment (8.435μ_B) is closer to the calculated value based on DFT than the value in the C2DB database. The magnetic moment is mainly originated from Os ions, and a small portion originates from Cl ions. For Fe₃GeSe₂ and Mn₄N₃S₂, the DOS indicates metallic characteristic; the calculated magnetic moments are 7.933μ_B and 12.844μ_B, respectively, which agree well with the predicted ones, where the magnetic moment mainly originates from the unpaired electrons in the 3D orbitals of TMs Fe and Mn.

The high-throughput DFT calculations for the 615 candidates took approximately 57,600 core-hours on an high-throughput calculations (HPC) cluster running on Intel Xeon Gold 5218 processors. The training and screening of the ML classification and regression models took approximately 560 core-hours on an Intel (R) Xeon (R) Platinum 8259CL. The computational time was reduced by more than one order of magnitude during this process; thus, ML is used to accelerate the discovery of new 2DFM materials.

C2DB and 2DMP were adopted as the primary sources of the dataset. It should be noted that C2DB includes 4035 entries up to 2021; however, up to 2023.03, the C2DB includes 15,733 entries. The 2DMP dataset contains a total of 6351 data points. According to the C2DB report, there were 1617 2D materials with thermodynamic and dynamic stability in 4035 samples. Among the 1617 candidates, 1292 were non-magnetic, and 325 magnetic. To avoid an imbalance between the magnetic and nonmagnetic datasets, 800 magnetic 2D materials were randomly selected with thermodynamic and dynamic stability in the 2DMP dataset. Thus, a dataset of 2417 samples was constructed, of which 1292 were nonmagnetic and 1125 magnetic. This dataset was labeled as data-N-M and used to train the magnetic/nonmagnetic classification model. It is also found that 615 magnetic samples are with magnetic moments > 1μ_B in C2DB with 4035 entries. However, they did not include antiferromagnetic 2D samples. Even a large database of 15,733 2D materials contained only 17 samples marked as AFM. In this context, high-throughput DFT calculations for 615 FM samples were performed with the adoption of the workflow. The retrieved datasets were used to train the FM and AFM classification models. This dataset is labeled as data-F-A. An accurate FM dataset with 383 samples calculated using DFT+U was used to train the regression model of the magnetic moment, and the FM dataset was defined as data-R. Compared with previous studies, the datasets used contained a more extensive and precise variety of 2DFM materials.

The properties of a functional material can be determined either from experimental measurements or through simulations via calculations with the adoption of an ab initio approach. ML eliminates the need for expensive programs by predicting the target properties of novel 2D materials using suitable feature descriptors based on prior data. The selection of a suitable set of feature descriptors that determine the target property is of vital importance in all ML computations. Since all the properties of a material are ultimately functions of its structure, composition, and elemental properties, we adopt feature descriptors that encode this information³¹. For the magnetic-nonmagnetic classification model, the feature set includes descriptors of two types: the first one is the Materials Agnostic Platform for Informatics and Exploration (MAGPIE)³², which was proposed by Ward et al. and used to obtain the elemental property labeled as EP. The second feature is the composition, which includes the number of space groups. For the FM/AFM classification model, the feature-descriptor space was expanded by including the number of d orbital electrons on the TM atoms (n_d) which are strongly correlated for magnetization, and the atom-centered symmetry function (ACSF), which depends on the local coordination environment around a specific atom³³. In addition, the ionic radius, electronegativity, and dipole polarization of the TM and nearest-neighbor coordination atoms were taken from the Python Mendeleev package³⁴, which are largely responsible for FM/AFM classification. For the regression model of magnetization, we added Smooth Overlap of Atomic Positions (SOAP) features and sine Coulomb matrix features associated with the structure, which were labeled as STRUCT³⁵. The unpaired d orbital electrons of TM was considered as an important feature descriptor.

Our high-throughput first-principles calculations for data F-A were performed in the framework of DFT as implemented in the Vienna ab-initio simulation package code^36,37. The generalized gradient approximation (GGA)³⁸ with DFT‒D3 was used to describe the exchange correlation, which was a semiempirical dispersion‒correction method to correct the van der Waals (vdw) interactions³⁹. GGA+U correction was applied to strongly correlate the TM d orbitals⁴⁰. The U value was determined through linear response theory which ensures the reliability of the qualitatively calculated results in the current work. The corresponding U values of each TM are presented in Table SI. The plane-wave energy cutoff was set to 500 eV. The Brillouin zone integration was sampled with a Monkhorst-Pack⁴¹esh of 8 × 8 × 1. All the structure parameters were sequentially relaxed so that the Hellmann-Feynman forces were < 0.01 eV/Å, and the total energy changes converged to < 10⁻⁵ eV. The conjugate gradient algorithm⁴²as used to optimize the structure. A vacuum distance of > 15 Å was set between adjacent slabs to eliminate spurious interactions. Structure and spin density visualization and analysis were performed adopting the VESTA code⁴³.