Large Supercells Research Articles

Atomistic first-principles modeling of single donor spin-qubit

Using an impurity atom in crystal silicon as a spin-1/2 qubit has been made experimentally possible recently where the impurity atom acts as a quantum dot (QD). Quantum transport in and out of such a donor QD occurs in the sequential tunneling regime where a physical quantity of importance is the charging (addition) energy, which measures the energy necessary for adding an electron into the donor QD. In this work, we present a first-principles method to quantitatively predict the addition energy of the donor QD. Using density functional theory (DFT), we determine the impurity states that serve as the basis set for subsequent exact diagonalization calculation of the many-body states and energies of the donor QD. Due to the large effective Bohr radius of the conduction electrons in Si, very large supercells containing more than 10 000 atoms must be used to obtain accurate results. For the donor QD of a phosphorus impurity in bulk Si, the combined DFT and exact diagonalization predicts the first addition energy to be 53 meV, in good agreement with the corresponding experimental value. For the donor QD of an arsenic impurity in Si, the first addition energy is predicted to be 44.2 meV. The calculated many-body wave functions provide a vivid electronic picture of the donor QD.

Domain Dynamics Response to Polarization Switching in Relaxor Ferroelectrics.

Nanoscale polar regions, or nanodomains (NDs), are crucial for understanding the domain structure and high susceptibility of relaxors. However, unveiling the evolution and function of NDs during polarization switching at the microscopic level is of great challenge. The experimental in situ characterization of NDs under electric-field perturbations, and computational accurate prediction of the dipole switching within a sufficiently large supercell, are notoriously tricky and tedious. These difficulties hinder a full understanding of the link between micro domain dynamics and macro polarization switching. Herein, the real-time evolution of NDs at the nanoscale is observed and visualized during polarization switching in an exemplary relaxor system of Bi5- xLaxMg0.5Ti3.5O15. Two fundamentally different domain switching pathways and dynamic characteristics are revealed: one steep, bipolar-like switching between two degenerate polarization states; and another flat, multi-step switching process with a thermodynamically stable non-polar mesophase mediating the degenerate polarization states. The two are determined by the distinct Landau energy landscapes that are strongly dependent on the intrinsic domain configurations and interdomain interactions. This work bridges the gap between micro domain dynamics and macro polarization switching, providing a guiding principle for the strategic design and optimization of relaxors.

Comparative study of atomic models of plastic deformation mechanism in ductile inorganic semiconductors: A theoretical discussion

Recently developed ductile inorganic semiconductors have opened a new avenue toward potential applications such as flexible electronics. Significant studies have been conducted based on different atomic models to understand the deformation mechanism using first-principle calculations; however, the comparative study on these atomic models remains unexplored. In this paper, taking Ag2S as an example, we simulated the slipping process using bulk, slab, and tilt-cell models. Systematic supercell tests were performed to investigate the convergence of the slip/cleavage energy based on the three atomic models. It was found that a reasonably large supercell is required to converge the slip/cleavage energy, and the ratio between minimum cleavage and maximum slip energy converges to a comparable value among different atomic models, where the tilt-cell and slab models show slightly higher energy ratio values compared with the bulk model. However, the comparison of different exchange-correlation functionals and van der Waals corrections indicates that the calculations of slip and cleavage energies are sensitive to the choice of calculation methods, highlighting the importance of using the same method for comparing slip and cleavage energies of different materials. This work provides insights into understanding different atomic models of ductility mechanisms in ductile inorganic semiconductors.

Computational insights into the structural, thermodynamic and transport properties of CaF2-MgF2 binary fluoride system at high temperatures

The structural, thermodynamic and transport properties of the CaF2-MgF2 molten salt system were investigated with ab initio molecular dynamics (AIMD), system-specific neural network interatomic potentials (NNIPs) and universal PreFerred Potentials (PFP). We trained an NNIP model using AIMD data as input and used this potential to efficiently simulate the interactions within a large supercell in a temperature range of 1273–1773 K. The Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) code was employed to validate our trained NNIP model. The Matlantis software with universal PFP is also presented to prove its feasibility for MD calculations and can be considered as a useful alternative simulation tool for higher-order systems where existing potentials are not readily available. We calculated structural and thermodynamic properties including radial distribution function (RDF), angular distribution function (ADF), specific heat capacity, ionic self-diffusivity, and viscosity. Our results indicate that the system exhibited a high degree of structural disorder, with the Ca, Mg, and F ions forming a liquid solution. Using PFP, the positions of the first peak in RDFs for Ca-F and Mg-F pairs are only slightly left-shifted (<0.05 Å), and the estimated viscosity of the melt decreases from 4.613 mPa·s to 1.846 mPa·s with an increase in temperature from 1273 K to 1773 K, in agreement with the NNIP trained specifically for CaF2-MgF2. Our results provide valuable insights into the properties of the CaF2-MgF2 system at high temperatures and serve as predictive models for the development of new electrolytes that could be used for silicon epitaxy by adding silica.

Development of All-Atom Empirical Potentials for Phloroglucinol (Phg) and Hexachlorocyclotriphosphazene (HCCP) Based on their Crystal Phase Structures.

Two existing generic force fields have been augmented with partial charges and tuned in order to give intercompatible all-atom empirical potentials that can satisfactorily represent the known crystal phase structures of the organic phloroglucinol (Phg) (C6H6O3) and inorganic hexachlorocyclotriphosphazene (HCCP) (N3P3Cl6) molecules at several temperatures. It has been proposed that HCCP-Phg network polymers could act as efficient H2 barrier layers in hydrogen storage tanks for cars. However, essential requirements for modeling such networks are adequate representations of both monomers in their pure dense phases using a common form of force field. Tests of their ability to maintain stable crystal structures have been made using classical molecular dynamics (MD) simulations on large 800-molecule supercells. The force fields have been optimized to match the densities calculated from the experimental unit cell dimensions at ambient conditions as well as the intermolecular potential energy, as estimated from experimental enthalpies of sublimation. For Phg, the crystal structure is stabilized by a network of hydrogen bonds and the Coulombic interactions contribute to over 55% of the total intermolecular potential energy. In contrast, the crystal structure of HCCP is intrinsically stabilized by the van der Waals terms. Both optimized force fields reproduce very well the orthorhombic symmetry of their respective crystals under constant-pressure NPT conditions. The model parameters tuned at ambient temperature also give reasonable agreement with crystallographic data at lower temperatures.

Revisiting Aspirin Polymorphic Stability Using a Machine Learning Potential.

In this study, we present a systematic computational investigation to analyze the long-debated free energy stability of two well-known aspirin polymorphs, denoted as Form I and Form II. Specifically, we developed a strategy to collect training configurations covering diverse interatomic interactions between representative functional groups in aspirin crystals. Utilizing a state-of-the-art neural network interatomic potential (NNIP) model, we trained an accurate machine learning potential to simulate aspirin crystal dynamics under finite temperature conditions with ∼0.46 kJ/mol/molecule accuracy. Employing the trained NNIP model, we performed thermodynamic integration to assess the free energy difference between aspirins I and II, accounting for the anharmonic effects in a large supercell consisting of 512 molecules. For the first time, our results convincingly demonstrated that Form I is more stable than Form II at 300 K, ranging from 0.74 to 1.83 kJ/mol/molecule, aligning with experimental observations. Unlike the majority of previous simulations based on (quasi)harmonic approximations in a small super cell, which often found degenerate energies between aspirins I and II, our findings underscore the importance of anharmonic effects in determining polymorphic stability ranking. Furthermore, we proposed the use of the rotational degrees of freedom of methyl and ester/phenyl groups in aspirin crystals as characteristic motions to highlight rotational entropic contribution that favors the stability of Form I. From the structural perspective, we also found that the subtle free energy difference can be used to explain the distinct thermal expansion responses as observed in both experimental and simulation data. Beyond the aspirin polymorphism, we anticipate that such entropy-driven stabilization can be broadly applicable to many other organic systems, suggesting that our approach holds great promise for stability studies in small-molecule drug design.

Insights into core-mantle differentiation from bulk Earth melt simulations

The earth is thought to have gone through complex physicochemical changes during the accretion and magma ocean stages. To better understand this evolution process at the fundamental level, we investigate the behavior of a bulk earth melt system by simulating the composition Fe35.7Mg19.0Si15.2O30.2 (in wt%) at high pressure. A deep neural network potential trained by first-principles data can enable accurate molecular dynamics simulation of large supercells that greatly enhances sampling for reliable evaluation of elemental partitioning. Our simulated system undergoes a phase separation in which the four elements clump together to different extents into two major domains. Based on the coordination and space-decomposition analyses, the inferred composition at 3000 K and 29.1 GPa contains 96.2, 0.1, 1.9 and 1.7 wt% of Fe, Mg, Si, and O, respectively, for the one domain and the corresponding elemental proportions are 3.0, 29.7, 22.0, and 45.3 wt% for the other domain. The predicted segregation thus leads to the formation of an iron-rich phase which corresponds to the metallic core and a magma ocean phase which corresponds to the silicate mantle. The metallic domain incorporates more silicon and more oxygen whereas the magma ocean domain gains more iron oxides at higher temperatures. Our predicted compositions compare favorably with those derived from experimental work for the equilibrium state metal and silicate reacting under high-pressure conditions.

Nonadiabatic Molecular Dynamics in Momentum Space Beyond Harmonic Approximation: Hot Electron Relaxation in Photoexcited Black Phosphorus.

We simulated hot-electron relaxation in black phosphorus using the nonadiabatic molecular dynamics (NA-MD) approach with a non-Condon effect in momentum space beyond the harmonic approximation. By comparing simulations at the Γ point in a large supercell with those using a few k-points in a smaller supercell─while maintaining the same number of electronic states within the same energy range, we demonstrate that both setups yield remarkably consistent energy relaxation times, regardless of the initial state energy. This consistency arises from the complementary effects of supercell size in real space and the number of k-points in the reciprocal space. This finding confirms that simulations at a single k-point in large size supercells are an effective approximation for NA-MD with a non-Condon effect. This approach offers significant advantages for complex photophysics, such as intervalley scattering and indirect bandgap charge recombination, and is particularly suitable for large systems without the need for a harmonic approximation.

ScaleLat: A chemical structure matching algorithm for mapping atomic structure of multi-phase system and high entropy alloys

ScaleLat (Scale Lattice) is a computer program written in C for performing the atomic structure analysis of multi-phase system or high entropy alloys (HEAs). The program implements an atomic cluster cell extraction algorithm to obtain all symmetry independent characteristic atomic cluster cells for the complex atomic configurations which are usually obtained from molecular dynamics or kinetic Monte-Carlo simulations at nanoscale or mesoscopic scale. ScaleLat implements an efficient and unique chemical structure matching algorithm to match all extracted atomic clusters from a large supercell (>104 atoms) to a representative small one (∼ 103 or less), providing the possibility to directly use the highly accurate quantum mechanical methods to study the electronic, magnetic, and mechanical properties of multi-component alloys for complex microstructures. We demonstrate the capability of ScaleLat code by conducting both the atomic structure matching analysis for Fe-12.8 at.% Cr binary alloy and equiatomic CrFeCoNiCu high entropy alloy, successfully obtaining the representative supercells containing 102∼103 atoms for two systems. The reliability of the proposed chemical structure matching scheme is tested and confirmed by calculating the electronic structures of both examples using trial supercells with various sizes. Overall, ScaleLat program provides a universal platform to efficiently map all essential chemical structures of large complex atomic structures to a relatively easy-handling small supercell for quantum mechanical calculations of various user interested properties. Program SummaryProgram Title: ScaleLatCPC Library link to program files: https://doi.org/10.17632/cv6wsxy938.1Developer's repository link: https://github.com/NanLi-xjtu/ScaleLat.gitLicensing provisions: MITProgramming language: CNature of Problem: Very large supercells containing more than 104 atoms must be used to describe the atomic structure of multi-phase materials or high entropy alloys, and their electronic and mechanical properties are almost not possible to be investigated using quantum mechanical calculations within conventional computational hardware. Decreasing the total number of atoms in a smaller representative supercell which preserves all essential chemical structures of the original large supercell is the key to tackling the problem.Solution to the Problem: An atomic cluster cell extraction algorithm is realized to thoroughly obtain all symmetry independent characteristic chemical structures of multi-phase system or high entropy alloys. Utilizing the direct atom swapping method, the efficient chemical structure matching procedure is developed to map all extracted characteristic atomic cluster cells of benchmark structure to the small supercell by minimizing their differences in both the types and relative fractions of atomic cluster cell sets.

Orbital ordering patterns in KBF (B = Sc, Ti, Fe, Co) perovskites.

The orbital ordering (OO) resulting from the partial occupancy of the subshell of the transition metals in KBF (B = Sc, Ti, Ffe, Co) perovskites, and the many possible patterns arising from the coupling between the B sites, have been investigated at the quantum mechanical level ( Gaussian type basis set, B3LYP hybrid functional) in a 40 atoms supercell. The numerous patterns are distributed into 162 classes of equivalent configurations. For each fluoroperovskite, one representative per class has been calculated. The four compounds behave similarly: an identical dependence of the energy and volume (or cell parameters) on the OO pattern; the spanned energy interval is small (1 to 2 mE per formula unit), suggesting that most of the configurations are occupied at room and even at low temperature. A linear model, taking into account the relative orbital order in contiguous sites, reproduces the energy order in the full set for each compound, suggesting that it could be used for studying OO in larger supercells.

Accelerating defect predictions in semiconductors using graph neural networks

First-principles computations reliably predict the energetics of point defects in semiconductors but are constrained by the expense of using large supercells and advanced levels of theory. Machine learning models trained on computational data, especially ones that sufficiently encode defect coordination environments, can be used to accelerate defect predictions. Here, we develop a framework for the prediction and screening of native defects and functional impurities in a chemical space of group IV, III–V, and II–VI zinc blende semiconductors, powered by crystal Graph-based Neural Networks (GNNs) trained on high-throughput density functional theory (DFT) data. Using an innovative approach of sampling partially optimized defect configurations from DFT calculations, we generate one of the largest computational defect datasets to date, containing many types of vacancies, self-interstitials, anti-site substitutions, impurity interstitials and substitutions, as well as some defect complexes. We applied three types of established GNN techniques, namely crystal graph convolutional neural network, materials graph network, and Atomistic Line Graph Neural Network (ALIGNN), to rigorously train models for predicting defect formation energy (DFE) in multiple charge states and chemical potential conditions. We find that ALIGNN yields the best DFE predictions with root mean square errors around 0.3 eV, which represents a prediction accuracy of 98% given the range of values within the dataset, improving significantly on the state-of-the-art. We further show that GNN-based defective structure optimization can take us close to DFT-optimized geometries at a fraction of the cost of full DFT. The current models are based on the semi-local generalized gradient approximation-Perdew–Burke–Ernzerhof (PBE) functional but are highly promising because of the correlation of computed energetics and defect levels with higher levels of theory and experimental data, the accuracy and necessity of discovering novel metastable and low energy defect structures at the PBE level of theory before advanced methods could be applied, and the ability to train multi-fidelity models in the future with new data from non-local functionals. The DFT-GNN models enable prediction and screening across thousands of hypothetical defects based on both unoptimized and partially optimized defective structures, helping identify electronically active defects in technologically important semiconductors.

First-principles calculation of lattice distortion, electronic structure, and bonding properties of GeTe-based and PbSe-based high-entropy chalcogenides

The massive amount of wasted heat energy from industry has pushed the development of thermoelectric (TE) materials that directly convert heat into electricity to a new level of concern. Recently, multicomponent alloys such as GeTe-based and PbSe-based high-entropy (HE) chalcogenides have attracted a great deal of attention due to their potential application as TE materials. The nature of the interatomic bonding, lattice distortion (LD), and the electronic structure in this class of materials is not fully understood. Herein, we report a comprehensive computational investigation of nine GeTe-based HE alloys with eight metallic elements (Ag, Pb, Sb, Bi, Cu, Cd, Mn, and Sn) with large supercells of 1080 atoms each; seven PbSe-based HE solid solutions: Pb0.99−ySb0.012SnySe1−2xTexSx (x = 0.1, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45, with y = 0) with supercells of 1000 atoms each; and five Pb0.99−ySb0.012SnySe1−2xTexSx (y = 0.05, 0.1, 0.15, 0.2, 0.25 with x = 0.25) solid solutions with supercells of 1000 atoms each. All these HE models are theoretically investigated for the first time. The electronic structure, interatomic bonding, charge transfer, and lattice distortion (LD) are investigated by first-principles calculations based on density functional theory. Multicomponent HE alloys can cause a significant LD, which affects their mechanical, thermal, and TE properties. The calculations for the GeTe-based HE chalcogenides showed that they are semiconductors with a narrow bandgap, except for m8, which has a semi-metallic characteristic, and this makes them good candidates for TE applications. For most of these models, the Fermi level shifts upward and locates deeply in the conduction bands, resulting in the enhancement of the electrical conductivity (σ). The bonding properties showed that most bonds in m5 are more dispersed, indicating highest LD and lower lattice thermal conductivity. For PbSe-based HE solid solutions, the LD calculations showed that the models Pb0.99Sb0.012Se0.5Te0.25S0.25 and Pb0.89Sb0.012Sn0.1Se0.5Te0.25S0.25 have the higher LD, and thus a lower lattice thermal conductivity. Such investigations are in high demand since it enables us to design new HE chalcogenides for TE applications. We use the novel concept of total bond order density as a single quantum mechanical metric to characterize the internal cohesion of these HE alloys and correlate with calculated properties, especially the mechanical properties. This work provides a solid database for HE chalcogenides and a road map for many potential applications. Moreover, the computational procedure we developed can be used to design new HE chalcogenides for specific TE applications.

Role of phonon scattering and bonding in resolving lattice thermal conductivity ambiguities of β-Ga2O3.

Understanding the lattice thermal conductivity (κl) of β-Ga2O3 is very intriguing owing to its advantages in high-voltage and high-temperature applications. Despite several attempts, the underlying mechanism and causes of the notable discrepancies found in the κl values of β-Ga2O3 along [100] and [001] directions calculated using first principles remained unresolved. We demonstrate that the understanding of the nature of chemical bonding is crucial to overcome the inconsistency in theoretically reported κl values. In low-symmetry structures such as β-Ga2O3, the nature of the interactions is primarily long-range; therefore, a sufficiently large supercell inclusive of various bonding characteristics is required to capture relevant phonon wavelengths. Bonding nature-aware structure modeling allows precise estimation of acoustic and optical mode contributions towards κl. Additionally, phonon mean free path analysis confirms that considering only third-order interaction terms is adequate to determine the κl of β-Ga2O3. The calculated κl values are in excellent agreement with experimentally reported values in all three directions. Our results establish that the bonding nature and its influence on phonon scattering are essential to consider in calculating κl accurately.

Accelerating GW Calculations of Point Defects with the Defect-Patched Screening Approximation.

The GW approximation has been widely accepted as an ab initio tool for calculating defect levels with the many-electron effect included. However, the GW simulation cost increases dramatically with the system size, and unfortunately, large supercells are often required to model low-density defects that are experimentally relevant. In this work, we propose to accelerate GW calculations of point defects by reducing the simulation cost of many-electron screening, which is the primary computational bottleneck. The random-phase approximation of many-electron screening is divided into two parts: one is the intrinsic screening, calculated using a unit cell of pristine structures, and the other is the defect-induced screening, calculated using the supercell within a small energy window. Depending on specific defects, one may only need to consider the intrinsic screening or include the defect contribution. This approach avoids the summation of many conduction states of supercells and significantly reduces the simulation cost. We have applied it to calculate various point defects, including neutral and charged defects in two-dimensional and bulk systems with small or large bandgaps. The results are consistent with those from the direct GW simulations. This defect-patched screening approach not only clarifies the roles of defects in many-electron screening but also paves the way to fast screen defect structures/materials for novel applications, including single-photon sources, quantum qubits, and quantum sensors.

Optical properties of charged defects in monolayer MoS2

We present theoretical calculations of the optical spectrum of monolayer MoS2 with a charged defect. In particular, we solve the Bethe–Salpeter equation based on an atomistic tight-binding model of the MoS2 electronic structure which allows calculations for large supercells. The defect is modelled as a point charge whose potential is screened by the MoS2 electrons. We find that the defect gives rise to new peaks in the optical spectrum approximately 100–200 meV below the first free exciton peak. These peaks arise from transitions involving in-gap bound states induced by the charged defect. Our findings are in good agreement with experimental measurements.

Trajectory sampling and finite-size effects in first-principles stopping power calculations

Real-time time-dependent density functional theory (TDDFT) is presently the most accurate available method for computing electronic stopping powers from first principles. However, obtaining application-relevant results often involves either costly averages over multiple calculations or ad hoc selection of a representative ion trajectory. We consider a broadly applicable, quantitative metric for evaluating and optimizing trajectories in this context. This methodology enables rigorous analysis of the failure modes of various common trajectory choices in crystalline materials. Although randomly selecting trajectories is common practice in stopping power calculations in solids, we show that nearly 30% of random trajectories in an FCC aluminum crystal will not representatively sample the material over the time and length scales feasibly simulated with TDDFT, and unrepresentative choices incur errors of up to 60%. We also show that finite-size effects depend on ion trajectory via “ouroboros” effects beyond the prevailing plasmon-based interpretation, and we propose a cost-reducing scheme to obtain converged results even when expensive core-electron contributions preclude large supercells. This work helps to mitigate poorly controlled approximations in first-principles stopping power calculations, allowing 1–2 order of magnitude cost reductions for obtaining representatively averaged and converged results.

Special glass structures for first principles studies of bulk metallic glasses

The atomic-level structure of bulk metallic glasses is a key determinant of their properties. An accurate representation of amorphous systems in computational studies has traditionally required large supercells that are unfortunately computationally demanding to handle using the most accurate ab initio calculations. To address this, we propose to specifically design small-cell structures that best reproduce the local geometric descriptors (e.g., pairwise distances or bond angle distributions) of a large-cell simulation. We rely on molecular dynamics (MD) driven by empirical potentials to generate the target descriptors, while we use reverse Monte Carlo (RMC) methods to optimize the small-cell structure. The latter can then be used to determine mechanical and electronic properties using more accurate electronic structure calculations. The method is implemented in the Metallic Amorphous Structures Toolkit (MAST) software package.

Effect of Framework Composition and NH3 on the Diffusion of Cu+ in Cu-CHA Catalysts Predicted by Machine-Learning Accelerated Molecular Dynamics.

Cu-exchanged zeolites rely on mobile solvated Cu+ cations for their catalytic activity, but the role of the framework composition in transport is not fully understood. Ab initio molecular dynamics simulations can provide quantitative atomistic insight but are too computationally expensive to explore large length and time scales or diverse compositions. We report a machine-learning interatomic potential that accurately reproduces ab initio results and effectively generalizes to allow multinanosecond simulations of large supercells and diverse chemical compositions. Biased and unbiased simulations of [Cu(NH3)2]+ mobility show that aluminum pairing in eight-membered rings accelerates local hopping and demonstrate that increased NH3 concentration enhances long-range diffusion. The probability of finding two [Cu(NH3)2]+ complexes in the same cage, which is key for SCR-NOx reaction, increases with Cu content and Al content but does not correlate with the long-range mobility of Cu+. Supporting experimental evidence was obtained from reactivity tests of Cu-CHA catalysts with a controlled chemical composition.

Multiple k-Point Nonadiabatic Molecular Dynamics for Ultrafast Excitations in Periodic Systems: The Example of Photoexcited Silicon.

With the rapid development of ultrafast experimental techniques for the research of carrier dynamics in solid-state systems, a microscopic understanding of the related phenomena, particularly a first-principle calculation, is highly desirable. Nonadiabatic molecular dynamics (NAMD) offers a real-time direct simulation of the carrier transfer or carrier thermalization. However, when applied to a periodic supercell, there is no cross-k-point transitions during the NAMD simulation. This often leads to a significant underestimation of the transition rate with the single-k-point band structure in a supercell. In this work, based on the surface hopping scheme used for NAMD, we propose a practical method to enable the cross-k transitions for a periodic system. We demonstrate our formalism by showing that the hot electron thermalization process by the multi-k-point NAMD in a small silicon supercell is equivalent to such simulation in a large supercell with a single Γ point. The simulated hot carrier thermalization process of the bulk silicon is compared with the recent ultrafast experiments, which shows excellent agreements. We have also demonstrated our method for the hot carrier coolings in the amorphous silicons and the GaAlAs_{2} solid solutions with the various cation distributions.

Efficient ab initio stacking fault energy mapping for dilute interstitial alloys

Density Functional Theory (DFT) is the prevalent first principles computational method for determining the stacking fault energy (SFE) of face centered cubic (fcc) metals and alloys. Due to several theoretical and computational challenges, SFE determination for interstitial alloys with alloying elements such as carbon, nitrogen, and hydrogen, has so far been limited to few studies at relatively high interstitial content. We propose a new method, rooted in the axial interaction model, that allows rapid and robust mapping of SFE for virtually arbitrary interstitial contents. Instead of computing the total energy of a very large supercell to represent dilute interstitial solutions, representative interstitial-affected and bulk regions are treated separately at the equivalent volume. The SFE is obtained by balancing the SFE values of the regions with a lever rule approach. The method matches SFE values from the axial interaction model within ≤4 mJ.m−2 error, as validated for non-magnetic fcc Fe-N and paramagnetic fcc Fe-N and AISI 304 alloys. The significantly reduced computational workload and equidistant SFE mapping vs. interstitial content down to extremely low values allows accurate fitting of the SFE vs. interstitial content with only few datapoints. This further improves the computational efficiency. So far DFT-based SFE mapping was limited to purely substitutional alloys; we demonstrate the first-time DFT-based SFE mapping in fcc AISI 304 vs. N and Ni, revealing a non-additive contribution of N and Ni to the SFE. Finally, the remaining challenges and future application for high-throughput DFT SFE computation in interstitial alloys is discussed.