Machine learning augmented design of 2D magnet with planar cyclo-tetranitrogen: ambient thermal stability from quantum to mesoscale

We first optimized the structure of N₄ clusters with the molecular point groups of T_d, D₂h, C₂h and C_s at B3LYP-D3/def2TZVP level and the results are shown in Fig. 1a. N₄-D₂h is identified as the metastable cluster and N₄-C_2h is found to be the most stable species (Supplementary Table 1). As illustrated in Fig. 1b, the D₂h cyclo-N₄ cluster is a singlet ground state owing to its highest occupied molecular orbital (MO14), which exhibits significant N–N bonding character. Furthermore, one- and two-electron reductions populate the next orbital (MO15), thereby stabilizing the anionic species in with the elongated N–N bonds. However, once we embed the such D₂h topology into the 2D periodic coordinated system with Fe, cyclo-N₄ exhibits a perfect square motif with an optimized N − N bond length of 1.37 Å (PBE + U level), as shown in Fig. 1c and Supplementary Fig. 1. Additionally, calculations with several exchange–correlation functionals (Supplementary Table 2) consistently predict an N–N bond length of ~1.37 Å, which lies between the typical values for N = N double bonds (≈1.23 Å) and N–N single bonds (≈1.43 Å). For instance, at ambient pressure, the N = N bond in [N₂]^2- ions in BaN₂ is 1.23 Å and the calculated N − N bond lengths in [N₂]^4- in PtN₂ and OsN₂ are 1.41 Å and 1.43 Å respectively^19,20. Therefore, although the polymeric cyclo-N₄ framework bonds have been reported in the cluster model and bulk compounds previously, the t-FeN₄ monolayer possesses a unique single N − N bond with perfectly planar square N₄ unit rarely observed/predicted in 2D confirmed space.

a The optimized structure of N₄-T_d, N₄-D₂h, N₄-C₂h, and N₄-C_s. Bond lengths are given in angstrom. Values were obtained from B3LYP-D3/def2TZVP computations. b Frontier molecular orbitals of the neutral cyclo-N₄ unit. c Crystal structure of t-FeN₄ monolayer viewed along c-axis. The bond lengths are given in angstrom and bond angle in degree. d Critical points in the unit cell of the t-FeN₄ monolayer (shown as smaller balls: bonds, red; rings, yellow). e Plot of the ELF map viewed from the vertical direction. ELF ranges from 0.0 (blue) to 1.0 (red). f Deformation density for the t-FeN₄ monolayer. Red and green refer to electron accumulation and depletion regions. The isovalue is 0.005 e/ų. g Phonon dispersion spectrum and (h) phonon density of states of the t-FeN₄ monolayer. i Variation of the \({\angle }_{N-N-N}\) angle in t-FeN₄ supercell as a function of the AIMD simulations. The green trace shows the \({\angle }_{N-N-N}\) angle values as recorded from the AIMD simulation, while the pink trace represents the fitted trend through data points.

An intriguing inquiry arises regarding the bonding character of such unique 2D (4.6²)₄(6⁴) cpa topology²¹. It is of particular interest to determine if the cyclo-N₄ specie form ionic bonds with the cations, or possibly coordination or covalent bonds. To address this, we conducted a Quantum Theory of Atoms in Molecules (QTAIM) analysis²². To study the bonding characteristics, we computed the Laplacian ∇²_ρ(r) of the N − N and Fe−N bonds in the t-FeN₄ monolayer structure, as depicted in Table 1. The corresponding bond critical points (BCPs) are shown in Fig. 1d. The type and strength of chemical bond can be characterized by the charge density the ρ(r) and the curvature ∇²_ρ(r) at the BCP. The calculated values of ∇²ρ(r) for the N − N and Fe−N interactions are −0.897 a.u. and 0.505 a.u., respectively. According to established QTAIM criteria, A negative ∇²_ρ(r) at the BCP, together with a relatively large ρ(r), is taken as a hallmark of shared‑electron (covalent) bonding, whereas a positive ∇²_ρ(r) coupled with a smaller ρ(r) indicates a closed‑shell (ionic or coordination) interaction^22,23. The results indicated that this value is between bulk AlN₄ and Li₂N₄ (Supplementary Table 3), closely match the strongly negative value we observe, confirming robust N − N covalency of the structure. Next, the bonding character of t-FeN₄ monolayer is further understand by electron localization function (ELF) mapping on the cyclo-N₄ plane. Figure 1e and Supplementary Fig. 2 display a prominent accumulation of electron density between each nitrogen pair, confirming the presence of four N − N bonds in cyclo-N₄. Meanwhile, the N lone pairs are assigned to regions around the four N atoms and bonded to adjacent Fe atoms, which exhibit moderate ELF values (<0.5), consistent with partial electron delocalization typical the bonds. Deformation-density plots show pronounced electron accumulation in the N − N bonding region and comparatively diffuse accumulation around Fe−N, corroborating the stronger covalent nature of N − N bonding and the mixed ionic/covalent nature of Fe−N bonding, as shown in Fig. 1f. Moreover, Crystal Orbital Hamilton Population (COHP) analyses (Supplementary Fig. 3) reveal large negative integrated COHP (–ICOHP) for N − N bonds (5.63 eV), indicating strong covalent character²⁴. While the ICOHP value of 2.32 eV for Fe–N bond suggests moderate covalent character compared to N − N bonds. Bader charge analysis indicates that each Fe atom donates approximately 1.14 e to the nearby N atoms. This observation is also consistent with the deformation density plot, providing further confirmation of the charge redistribution within the t-FeN₄ monolayer.

The stability of the t-FeN₄ monolayer has been thoroughly validated through comprehensive analyses, including the calculated formation energy, phonon dispersion spectrum (Figs. 1g, h), ab-initio molecular dynamics (AIMD) simulations (Supplementary Fig. 4), and elastic constants (Supplementary Note 1 and Supplementary Table 4). Interestingly, the N₄ ring retains a nearly perfect square configuration even at finite temperature, as shown in Fig. 1i and Supplementary Fig. 5, which indicates the high stability of the cyclo-tetranitrogen in the t-FeN monolayer. Detailed discussion of the structural energy stability, mechanical stability/properties, and the stability of the t-FeN₄ structure under aqueous solution can be found in Supplementary Note 1 and Supplementary Figs. 6–9. Moreover, based on the recent empirical correlation proposed by Bondarchuk et al. for N₅-containing systems, the decomposition temperature (\({T}_{{dec}}\)) can be estimated from the longest N–N bond length (\({l}_{N-N}\)) using the equation \({T}_{{dec}}(^\circ {\rm{C}})=-4073\times {[l}_{N-N}({\text{\AA }})]+5594.7\)²⁵. For the t-FeN₄ monolayer, we obtain \({T}_{{dec}}\) ≈ 15 °C. This estimate further supports the ambient-temperature stability of the cyclo-N₄ unit in the t-FeN₄ monolayer. Therefore, even energetically the MN₄ structure with cyclo‑N₄ unit represents a metastable configuration¹⁸, our multiscale computational analyses demonstrate that the proposed t‑FeN₄ structure exhibits robust stability in both dynamical and thermodynamic dimensions.

Beyond unique structural and chemical bonding properties, the 2D confinement induces emergent magnetic phenomena through spin-polarized t-FeN₄ monolayer. We then investigated in detail the electronic structure and magnetic properties of the t-FeN₄ monolayer structure. Using the supercell model, we evaluated three potential magnetic configurations (Fig. 2a) to identify the magnetic ground state of the t-FeN₄ monolayer. These configurations included one ferromagnetic (FM) state and two antiferromagnetic (AFM) states (stripy-AFM and checkerboard-AFM, respectively). The relative energies of the different magnetic states, calculated using various methods, are summarized in Supplementary Table 5. Among them, the stripy-AFM state was determined to be the ground state. This energy difference is attributed to the competition between the ferromagnetic coupling of nearest neighbours and the antiferromagnetic coupling of next-nearest neighbours, emphasizing the delicate balance of exchange interactions in low-dimensional systems. The magnetism primarily originates from the Fe atom, as demonstrated by the spin charge density (Supplementary Fig. 10), with each Fe atom displaying a magnetic moment of 2.0 μB. Based on spin-polarized HSE06 functional, our band structure calculations reveal that, despite all configurations exhibiting semiconductor behavior, the FM state uniquely displays half-metallicity, as shown in Figs. 2b–f (cf. density of states using HSE06 functional and band structure using PBE + U functional in Supplementary Figs. 11–16). In the FM configuration, while the spin-down channel remains metallic, the spin-up channel develops an insulating gap of 3.39 eV. This pronounced asymmetry is associated with crystal field-induced splitting of the Fe_3d orbitals, where the d-orbital plays a dominant role in the spin-down conduction near the Fermi level (Supplementary Fig. 17). Such behaviour is highly advantageous for spintronic applications, as it facilitates 100% spin-polarized transport, which is essential for efficient spin injection and detection^26,27.

a Three different magnetic states: FM, checkerboard-AFM, and stripy-AFM, respectively. b Brillouin zone path with labelled high-symmetry k-points of the t-FeN₄ monolayer. c Spin-orbit coupling (SOC) induced magnetic moment orientations, illustrating in-plane (φ) and out-of-plane (θ) anisotropy. d–f Spin-polarized band structures for the three magnetic states based on HSE06 functional. g, h MAE analysis reveals strong out-of-plane anisotropy compared to in-plane driven by SOC-mediated orbital rehybridization. i The lattice constant and magnetic moment of the t-FeN₄ monolayer at different doping density. j The relative energy of three different magnetic arrangements and (k) MAE of t-FeN₄ monolayer as a function of carrier density. l The specific heat as a function of temperature for the t-FeN₄ monolayer.

Magnetic anisotropy energy (MAE) is another key factor in achieving long-range magnetic ordering in 2D materials, with higher MAE values indicating greater thermal stability of the magnetic state^28,29,30. The MAE is quantified as the energy differential between in-plane (\({E}_{x}\)) and out-of-plane (\({E}_{z}\)) ferromagnetic orientations, expressed mathematically as \({MAE}={E}_{x}-{E}_{z}\). A positive MAE value signifies preferential alignment of the magnetization axis perpendicular to the plane (out-of-plane easy axis), while a negative MAE indicates in-plane magnetization dominance. This directional dependence arises from spin-orbit coupling effects that modulate orbital hybridization in the 2D confined geometry. The MAE values corresponding to rotations in the x–y and x–z planes are illustrated in Fig. 2g, h, respectively. The out-of-plane easy axis exhibits a remarkably high MAE of 393 μeV/f.u. (cf. Supplementary Fig. 18 for 3D schematic of MAE), surpassing values reported for prototypical 2D magnets such as CrI₃ (24.68 μeV/f.u.)³¹, CrS₂ (255.57 μeV/f.u.)³², and CoGa₂S₄ (47 μeV/f.u.)³³.

We further investigate doping as a strategic approach to modulate its magnetic behaviour. As illustrated in Fig. 2i, electron doping induces lattice expansion, whereas hole doping drives contraction, a phenomenon attributable to charge-dependent Coulomb screening effects: Electron doping enhances interatomic repulsion via occupied antibonding states, whereas hole doping strengthens the Fe–N bonds by depleting electrons from bonding orbitals near the Fermi level. Notably, magnetic moments exhibit an inverse trend—electron doping suppresses the Fe-centered moment (from 2.0 μB to 1.75 μB), while hole doping enhances it (up to 2.3 μB at +0.3 e/f.u.). This dichotomy arises from doping-driven shifts in 3d-orbital occupancy, where electron filling reduces unpaired spins, whereas hole creation amplifies spin polarization. Remarkably, the stripy-AFM ground state persists across a broad doping range (−0.2 to +0.3 e/f.u.), as shown in Fig. 2j. Only under extreme electron doping (−0.3 e/f.u.) does the system transition to FM ordering, driven by weakened superexchange interactions and enhanced direct Fe–Fe ferromagnetic coupling. Crucially, AIMD simulations confirm structural integrity even at these doping extremes (Supplementary Figs. 19, 20), with no bond dissociation observed over 3 ps at 300 K. The doping strategy also enables precise control over magnetic anisotropy. For example, heavy electron doping reduces the MAE magnitude and eventually inverts its sign (−754 μeV/f.u. at −0.3 e/f.u.), switching the easy axis from out-of-plane to in-plane. This tunability stems from doping-induced redistribution of d-orbital populations, which modifies the relative contributions of second-order SOC terms. Such programmable magnetic anisotropy, coupled with ambient stability, positions t-FeN₄ as a versatile platform for spin-orbit torque devices and non-volatile memory architectures, where electric-field-controlled doping could enable dynamic switching between distinct magnetic states.

For spintronic devices to function effectively under ambient conditions, their critical temperatures must exceed room temperature. To evaluate this critical parameter for the t-FeN₄ monolayer, we performed large-scale Monte Carlo simulations based on the anisotropic Heisenberg model, explicitly accounting for both nearest-neighbor (J₁) and next-nearest-neighbor (J₂) exchange interactions (Supplementary Fig. 21). The interaction parameters, derived from spin-polarized PBE + U calculations, revealed strong antiferromagnetic coupling with J₁ = − 0.85 meV and J₂ = − 19.7 meV, where the dominant J₂ term governs the stripy-AFM ground state through competing super-exchange pathways mediated by the N₄ rings. The temperature-dependent magnetization curve (cf. Fig. 2l), obtained through finite-size scaling analysis of the 30 × 30 supercell, demonstrates a Néel temperature of 600 K—a value that not only exceeds room temperature but also surpasses some prototypical 2D magnets such as CrCl₃ (10–17 K)³⁴, FeCl₂ (47 ± 8 K)³⁵, FePS₃ (120 K)³⁶, CrMoC₂S₆ (556 K)³⁷, etc. Therefore, the combination of ultrahigh Néel temperature and strong magnetic anisotropy suggests t-FeN₄ could serve as a promising material for room-temperature spin logic gates, non-volatile magnetic memory, and other quantum-enhanced technologies requiring persistent magnetic ordering without cryogenic constraints^28,38,39,40.

The demonstrated magnetic tunability motivates development of accurate computational tools for simulating realistic material configurations. However, traditional DFT-based approach becomes computationally prohibitive for studying temperature-dependent behaviour in extended systems containing more than thousands of atoms^41,42. The workflow for constructing the MLP for the t-FeN₄ monolayer and its subsequent application is comprehensively depicted in Fig. 3a. Detailed computational protocols for both the MLP construction and the AIMD simulations are provided in the Methods section and Supplementary Figs. 22–25. This dual approach is designed to rigorously assess the model’s capability to replicate critical properties—namely, the energies, forces, and structural configurations—which are essential for accurately reproducing the complex characteristics of the materials⁴³.

a Flowchart of the training procedure of the MLP model. Comparison of (b) energy and (c) averaged atomic forces for the configurations from DFT and MLP. d Phonon dispersion and (e) phonon density of states of t-FeN₄ monolayer computed with MLP. f The radial distribution function (RDF) for the Fe–Fe, Fe–N, and N–N distance in t-FeN₄ monolayer at 300 K. Solid lines indicate AIMD simulations, whereas dashed lines show the results of MLP trajectories. g Same for the angular distribution functions (ADF) at 300 and 800 K, computed using a radial cut-off of 2.0 Å.

Figures 3b, c present a direct comparison between the MLP model predictions and density functional theory (DFT) calculations for the t-FeN₄ monolayer. The results demonstrate excellent agreement, with the MLP model achieving a root mean square error (RMSE) of 14.1 meV per atom for energy and 0.016 eV Å⁻¹ for forces. These low error margins underscore the high fidelity of the MLP in capturing the intricate interatomic interactions within the system⁴³. Further, the dynamic stability of the FeN₄ monolayer is substantiated by the phonon analysis presented in Fig. 3d, e. The phonon dispersion curves obtained via the MLP show no imaginary frequencies, a clear indication of the dynamic stability of the structure. The phonon density of states analysis reveals that the low-frequency vibrational modes are primarily attributed to the Fe–N bonds, whereas the high-frequency modes predominantly arise from the vibrations within the N₄ ring. This nuanced understanding of the vibrational characteristics is critical for predicting the material’s thermal and mechanical behaviour.

Figure 3f compares the radial distribution functions (RDFs) computed from both AIMD and MLP-MD simulations for the t-FeN₄ monolayer supercell. The RDF analysis provides valuable insights into the spatial correlations and local bonding environments, reflecting the accuracy with which the MLP model replicates the DFT-derived structural features. Complementing this, Fig. 3g illustrates the angular distribution profiles at various temperatures. Notably, at 800 K, the increased thermal perturbations result in a broader angular distribution—a phenomenon that is effectively captured by the MLP-based molecular dynamics simulations. This consistency across temperature regimes further validates the robustness of the MLP model.

Overall, the integration of MLP with both DFT and AIMD simulations not only validates the model’s accuracy in describing the fundamental properties of the t-FeN₄ monolayer but also provides a powerful framework for exploring the material’s dynamic behaviour under varied conditions. This comprehensive approach is instrumental in advancing our understanding of complex material systems and paves the way for the rational design of novel materials with tailored properties.

The generalisation capability of the MLP is critical to ensure prediction accuracy, transferability, robustness against overfitting, and uncertainty quantification, thereby facilitating further utilisation in materials discovery. Therefore, we first evaluated the performance of the well-trained MLP for periodic t-FeN₄ structures of varying sizes, as illustrated in Figs. 4a–c and Supplementary Fig. 26. As the number of atoms increases, we observe a corresponding increase in the effective out-of-plane extension along the z-axis. This behaviour is primarily attributed to the cumulative effects of thermal fluctuations and minor surface reconstructions that become more pronounced in larger systems. Remarkably, despite these effects, the overall structural integrity is maintained even for systems as large as 57.0 nm², indicating a robust stability of the t-FeN₄ framework, as shown in Supplementary Figs. 27–29. To rigorously compare the MLP predictions with DFT results, we performed a 5-ps AIMD simulation under room temperature. The AIMD results, as shown in Fig. 4d, demonstrate that both the N–N bond lengths and bond angles remain almost perfectly square throughout the simulation. In our subsequent analysis, we applied the MLP to systems with different lattice sizes and directly compared the results to the AIMD benchmarks. As depicted in Figs. 4e, f, the MLP accurately reproduces the AIMD outcomes, yielding average N–N bond lengths of ~1.37 Å and bond angles of 90°. This high degree of correspondence confirms that the MLP possesses exceptional generalization capability; it accurately captures the thermodynamic and dynamic behaviour of the t-FeN₄ monolayer across a range of system sizes, even when the extent of out-of-plane fluctuations increases with the number of atoms. To further test the transferability of the MLP, we systematically simulated various non-periodic t-FeN₄ structures, including nanosheets, nanomatrices, and nanoribbons, as shown in Figs. 4g–i and Supplementary Figs. 30–32. Notably, these simulations encompass structures with edges and potential defect sites—configurations that were not included in the original training dataset of the MLP. Despite this, the MLP was able to reproduce high-precision finite-temperature behaviour across these diverse morphologies. The predicted structural parameters and dynamical evolution under finite temperature conditions closely match those obtained from AIMD, thereby reinforcing the reliability of the MLP across different dimensional scales and defect configurations.

Periodic t-FeN₄ configurations with the system size of (a) 3.2 × 3.2 nm², (b) 7.9 × 7.9 nm², and (c) 11.9 × 11.9 nm², respectively. Colour coding indicates the atomic positions in the direction normal to the layer. d Statistics of the N–N bond lengths and bond angles obtained from AIMD simulation at room temperature. e, f Statistical averages of the N–N bond length and bond angle obtained from the MLP-MD simulations at room temperature. g–i Systematic tests on t-FeN₄ nanosheets, nanomatrices, and nanoribbons. Colour coding indicates the atomic positions in the direction normal to the layer.

These systematic studies provide compelling evidence for the exceptional generalization ability of the MLP. The model not only reproduces the essential structural and dynamical properties of periodic t-FeN₄ monolayers, but also accurately simulates complex nanostructured systems that include edges, defects, and variations in atomic arrangement. Moreover, the MLP’s ability to efficiently simulate systems containing tens of thousands of atoms—while maintaining a high degree of accuracy—suggests that it can be effectively applied to large-scale atomistic simulations that are computationally prohibitive for conventional first-principles methods. The fact that even defected t-FeN₄ structures and various nanostructures remain stable at finite temperatures suggests that these materials are intrinsically robust. This robustness indicates that t-FeN₄ monolayers, with their unique structural and magnetic features, have a high likelihood of being successfully synthesized in the laboratory. Furthermore, the capability of the MLP to accurately predict the thermal, structural, and dynamic behaviour of these materials under realistic conditions can serve as a valuable guide for experimentalists aiming to design and fabricate novel low-dimensional magnetic systems for spintronic and magnetic storage applications.

In summary, we demonstrate a fundamental advance in polynitrogen stabilization strategies: 2D confinement in t-FeN₄ monolayer ambiently stabilizes planar cyclo-N₄ units through the combined action of geometric constraints and charge transfer, thus potentially achievable without high-pressure synthesis. The developed MLP exhibits excellent generalization across structural motifs—including extended monolayers, nanoribbons, and defective configurations—achieving quantum-mechanical accuracy while enabling large-scale simulations of temperature-dependent behaviour previously inaccessible to first-principles method. The emergent giant out-of-plane magnetic anisotropy and high Néel temperature, establish t-FeN₄ as a tunable platform for designing spin-filtering devices and non-volatile memory elements. By elucidating how dimensional reduction suppresses N ≡ N bond reformation through orbital hybridization control, this work provides a generalizable blueprint for stabilizing metastable nitrogen topologies, while the validated MLP framework bridges critical gaps between quantum-level precision and mesoscale materials engineering, collectively advancing both polynitrogen chemistry and magnetic material design.

Considering that the positive example of the synthesis of 2D polynitrogen layers under high pressure⁹, it is expected that such 2D system can be prepared under ambient conditions in future studies. Our findings unveil an alternative design principle for the stabilisation of polynitrogen, whereby 2D confinement suppresses the bond dissociation pathways that are inherent to bulk systems. This approach may guide the development of ambient‑stable polynitrogen architectures without relying solely on high‑pressure approach, and could inform future efforts to synthesize functional polynitrogen materials.

Ab initio calculations

For different N₄ species, the structure optimizations were carried out by using the B3LYP exchange-correlation functional^44,45. All cluster calculations were performed in the Gaussian 09 program suite⁴⁶. The def2-TZVP basis set⁴⁷ was employed for nitrogen during structure optimizations. The DFT-D3 method⁴⁸ was used to include the dispersion correction. All reported energies are obtained by single point energy calculations with G4 method⁴⁹.

All the periodic DFT calculations were performed using the Vienna Ab initio Simulation (VASP, 5.3.2 version) package⁵⁰. In these simulations, the projected augmented wave (PAW) method⁵¹ was adopted, and the generalized gradient approximation (GGA) in the form of the Perdew–Burke–Ernzerhof (PBE) functional⁵² was used for the exchange-correlation potential. We systematically evaluated the influence of different functionals (PBE, PBE + U, PBEsol, and B3LYP) on the lattice parameters, including both spin-polarized and non-spin-polarized calculations, which yielded comparable lattice constants (Supplementary Table 2). Hence, unless explicitly stated otherwise, all calculations presented here are based on DFT + U results. The on-site Coulomb interactions for localized d electrons in iron were treated using the Hubbard corrections following Dudarev’s approach⁵³, with U = 4.0 eV and J = 1.0 eV chosen to reproduce experimental bulk FeN₄ configuration¹¹. An energy cutoff of 850 eV and a 12 × 12 × 1 Monkhorst–Pack k-point mesh (as determined by convergence tests, Supplementary Fig. S33) were used for Brillouin zone integrations. For an accurate band structure analysis, the hybrid Heyd–Scuseria–Ernzerhof (HSE06) functional⁵⁴ was additionally employed, and spin–orbit coupling (SOC) effects were included in MAE calculations.

Electron charge density analysis

The electron charge density was analyzed topologically for t-FeN₄ monolayer using the Quantum Theory of Atoms in Molecules (QTAIM) approach²². This theory defines a solid using a zero-flux surface of the electron density gradient ∇ρ(r). The charge density distribution ρ(r) and its principal curvatures (eigenvalues of the Hessian matrix) at the bond critical point (BCP) provide insight into bonding type and properties. The second derivative of the electron density (Laplacian) ∇²ρ(r) reveals electron concentration at the BCP. Results showed a direct correlation between ∇ρ(r) and ∇²ρ(r) at the BCP and bond order and strength. Electron localization function (ELF) ranges from 0 to a maximum of 1, with values close to 1 indicating localized electrons. High ELF values are seen in bonds, cores, and lone pairs, with maxima appearing in the middle of bonds (e.g., the N-N bond in N₂ molecule).

Training data preparation

The finite temperature ab initio molecular dynamics (AIMD) simulations were performed to generate datasets for machine learning training by using time steps of 1 fs in 4 × 4 × 1 supercell (containing 80 atoms) of t-FeN₄ monolayer. Then the supercell structure was simulated by the Nosé–Hoover thermostat⁵⁵ in an NVT (N-number of particles, V-Volume, and T-temperature) ensemble. The Brillouin zone was sampled with a 1 × 1 × 1 Monkhorst−Pack k-points grid⁵⁶. We performed three different AIMD simulations by considering the different temperatures of 300, 500, and 800 K, to assess the impact of the temperature. The total dataset consists of around 30,000 configurations.

MLP training and validation

The smooth version of the DeePMD model with an inner cutoff of 1 Å and an outer cutoff of 6.9 Å was adopted in this work, and the DeePMD-kit package⁵⁷ was used for training and validation. The embedding network has the ResNet-like architecture⁵⁸ with a size of (25, 50, 100).

MLP-based MD simulations

For MD simulations, an optimized interface of DeePMD potential has been implemented in the Large-scale Atomic/Molecular Massively Parallel Simulator (Lammps) molecular dynamics package⁵⁹. All subsequent simulation through the DeepMD protentional was carried out by Lammps code with Nosé-Hoover style thermostat. The timestep in all MD simulations was 1 fs. We performed the ns-level MD simulation at 300 K temperature for the deferent 2D systems with the different atomic numbers (from 320 to 103,680 atoms) and different structural types, covering ideal periodic configurations, nanoribbons, nanosheets, and nanosheets with defects. (cf. Figure 4 and Supplementary Figs. 26–32).

Phonon calculation

The phonon properties were computed using two python-based tools: phonopy for VASP and phonoLAMMPS for LAMMPS^50,59,60. The phonon bands for the 19 × 19 supercell of t-FeN₄ monolayer were first calculated with a 12 × 12 × 1 mesh point sampling of the Brillouin-zone. The partial phonon density of states (Ph DOS) was then obtained from the calculations.

Visualization

MD trajectories and model structures were visualized using VESTA and OVITO software^61,62.

The central codes used in this paper are VASP and DeePMD-kit. Detailed information related to the license and user guide are available at https://www.vasp.at and https://github.com/deepmodeling/deepmd-kit.

Authors and Affiliations

1. School of Materials Science and Engineering, Chongqing Jiaotong University, Chongqing, PR China

   Dong Fan, Ke Zheng & Hongfang Li

2. KFMCH, Faculty of Science, Charles University, Prague, Czech Republic

   Junjie He

3. Key Laboratory of Low Dimensional Materials and Application Technology of Ministry of Education, School of Materials Science and Engineering, Xiangtan University, Xiangtan, China

   Pengbo Lyu

Contributions

D.F. conceived the study. D.F., P.L. and H.L, performed simulations. D. F., K. Z., H. Li, P. L., and J. H. discussed the results and contributed to the manuscript preparation. D.F. led the writing and editing of the manuscript. D.F. and P.L. supervised the study.

Corresponding authors

Correspondence to Dong Fan or Pengbo Lyu.

Competing interests

The authors declare no competing interests.

Open Access This article is licensed under a Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 International License, which permits any non-commercial use, sharing, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if you modified the licensed material. You do not have permission under this licence to share adapted material derived from this article or parts of it. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by-nc-nd/4.0/.

Reprints and permissions