Paper deep dive
Learning the Kohn-Sham map with neural operators for quasi-linear scaling density functional theory
Danish Khan, Maurice D. Hanisch, Nikolai Argatoff, Evan Xie, Sandeep Sharma, Anima Anandkumar
Intelligence
Status: succeeded | Model: Gemma-4-26B-A4B | Prompt: intel-v1 | Confidence: 92%
Last extracted: 8/26/2026, 4:58:50 AM
Summary
This paper introduces the Kohn-Sham Fourier Neural Operator (Kohn-Sham FNO), a machine learning method for orbital-free Density Functional Theory (DFT). It addresses the cubic scaling of traditional Kohn-Sham DFT by learning the Kohn-Sham map, which predicts electron density from the Kohn-Sham potential, using an SE(3)-equivariant Fourier Neural Operator. Trained on 8,504 molecules and solids, the model enables stable, quasi-linear scaling self-consistent field (SCF) calculations without explicit orbital diagonalization, achieving accuracy comparable to standard DFT and scaling to systems with up to 82,500 valence electrons.
Entities (8)
Relation Signals (7)
Kohn-Sham Fourier Neural Operator → isa → Fourier Neural Operator
confidence 95% · In this work, we introduce a new domain-invariant, SE(3)-equivariant FNO architecture... Kohn-Sham FNO learns the KS solution operator
Kohn-Sham Fourier Neural Operator → maps → Kohn-Sham potential
confidence 95% · It maps a Kohn-Sham potential directly to the corresponding density
Kohn-Sham Fourier Neural Operator → maps → electron density
confidence 95% · It maps a Kohn-Sham potential directly to the corresponding density
Kohn-Sham Fourier Neural Operator → replaces → orbital diagonalization
confidence 92% · The model replaces the orbital-solution step... enabling stable quasi-linear scaling SCFs
Kohn-Sham Fourier Neural Operator → achievesaccuracyof → Density Functional Theory
confidence 90% · reproducing densities, electronic spectra, and structural observables at Kohn–Sham DFT accuracy
Kohn-Sham Fourier Neural Operator → enables → quasi-linear scaling
confidence 90% · enabling stable quasi-linear scaling SCFs... Linear-scaling SCFs additionally allow converging magnesium dislocation densities
Kohn-Sham Fourier Neural Operator → trainedon → QM9
confidence 85% · Trained jointly on 8,504 molecules and solids... small-organic QM9
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:Kohn--Sham density functional theory (DFT) underpins electronic-structure simulations, but repeated orbital diagonalizations lead to cubic scaling, restricting quantum calculations to modest scales only. Eliminating these auxiliary orbitals while retaining Kohn--Sham accuracy is the central goal of orbital-free DFT, but both analytical and machine-learning methods have so far fallen short. Prior learning approaches either try to learn the variational kinetic-energy functionals, which are ill-conditioned, or directly predict the ground state, which extrapolate poorly to larger systems. Instead, we identify the Kohn--Sham map as the right learning target for orbital-free DFT. It maps a Kohn--Sham potential directly to the corresponding density and noninteracting kinetic energy, quantities otherwise obtained through an orbital diagonalization. Focusing on the density component in this work, a domain-invariant $\mathrm{SE}(3)$-equivariant Fourier neural operator learns to predict it from the potential as input on real-space grids, enabling stable quasi-linear scaling SCFs. Trained jointly on 8,504 molecules and solids, a single model generalizes to out-of-distribution organic molecules, insulators, and metals. For the first time, the same method converges SCFs across these systems without explicitly constructing Kohn--Sham orbitals, while reproducing densities, electronic spectra, and structural observables at Kohn--Sham DFT accuracy. Linear-scaling SCFs additionally allow converging magnesium dislocation densities containing up to 82,500 valence electrons on a single GPU.
Tags
Links
- Source: https://arxiv.org/abs/2608.23895v1
- Canonical: https://arxiv.org/abs/2608.23895v1
Trouble viewing inline? Open PDF directly →
Full Text
59,671 characters extracted from source content.
Expand or collapse full text
Learning the Kohn–Sham map with neural operators for quasi-linear scaling density functional theory Danish Khan 1† , Maurice D. Hanisch 1† , Nikolai Argatoff 1,3 , Evan Xie 1 , Sandeep Sharma 2,4 , and Anima Anandkumar 1* 1 Department of Computing and Mathematical Sciences, California Institute of Technology 2 Division of Chemistry and Chemical Engineering, California Institute of Technology 3 Department of Mathematics, ETH Z ̈urich 4 Marcus Center for Theoretical Chemistry, Pasadena, CA 91125, USA * Corresponding author. Email: anima@caltech.edu † Equal contribution. August 26, 2026 Abstract Kohn–Sham density functional theory (DFT) underpins electronic-structure simulations, but repeated orbital diagonalizations lead to cubic scaling, restricting quantum calculations to modest scales only. Eliminating these auxiliary orbitals while retaining Kohn–Sham ac- curacy is the central goal of orbital-free DFT, but both analytical and machine-learning methods have so far fallen short. Prior learning approaches either try to learn the varia- tional kinetic-energy functionals, which are ill-conditioned, or directly predict the ground state, which extrapolate poorly to larger systems. Instead, we identify the Kohn–Sham map as the right learning target for orbital-free DFT. It maps a Kohn–Sham potential directly to the corresponding density and noninteracting kinetic energy, quantities otherwise ob- tained through an orbital diagonalization. Focusing on the density component in this work, a domain-invariant SE(3)-equivariant Fourier neural operator learns to predict it from the potential as input on real-space grids, enabling stable quasi-linear scaling SCFs. Trained jointly on 8,504 molecules and solids, a single model generalizes to out-of-distribution or- ganic molecules, insulators, and metals. For the first time, the same method converges SCFs across these systems without explicitly constructing Kohn–Sham orbitals, while reproduc- ing densities, electronic spectra, and structural observables at Kohn–Sham DFT accuracy. Linear-scaling SCFs additionally allow converging magnesium dislocation densities contain- ing up to 82,500 valence electrons on a single GPU. 1 Introduction More than six decades after Hohenberg and Kohn established the ground-state density as the fundamental variable of electronic structure, extending predictive first-principles calculations to mesoscopic scales remains a central challenge [1]. Kohn and Sham made density functional the- ory (DFT) computationally practical, leading to its ubiquitous use across chemistry, materials science, condensed-matter physics, and biophysics [2–5]. Their formulation introduces an aux- iliary system of non-interacting electrons whose orbitals must be solved for and orthogonalized at every self-consistent field (SCF) iteration. These operations scale nominally as O(N 3 ) with N electrons, restricting the length scales accessible to first-principles simulations of extended defects, interfaces, electrochemical environments, and biological systems [6, 7]. The diagonalization of Kohn–Sham Hamiltonians is among the most frequently repeated and computationally consequential eigenvalue problems in modern science. As an example, nearly 1 arXiv:2608.23895v1 [physics.chem-ph] 24 Aug 2026 30% of the workload at the US National Energy Research Scientific Computing Center in 2018 was attributed to DFT calculations [8]. The search for linear-scaling electronic structure methods has a long history. Several efforts leverage Kohn’s nearsightedness principle to truncate the density matrix (DM) through localized orbitals, leading to sparse-matrix operations, while others use DM purification methods to impose idempotency [9, 10]. Other approaches include Fermi-operator expansions, selected inversion, and stochastic trace estimation [11, 12]. These methods extend the reach of KS-DFT, but their performance depends on density-matrix decay, sparsity, dimensionality, or stochastic error. There is therefore no general framework for near-linear scaling DFT with KS accuracy across molecules, metals, insulators, and chemically heterogeneous systems. Hohenberg and Kohn originally formulated DFT as a variational theory of the density alone. At fixed spatial resolution, the size of the density representation grows linearly with system size, offering a route to linear or quasi-linear scaling assuming a similar cost of applying the energy functionals. Orbital-free DFT (OF-DFT) seeks to recover this density-only formulation by minimizing the energy directly over the density [13]. Its central unknown is the non-interacting kinetic energy density functional (KEDF). KS-DFT evaluates this quantity from the auxiliary orbitals, whereas OF-DFT requires an explicit density functional T s [n]. Existing analytical KEDFs are efficient but lack the accuracy and transferability needed to describe shell structure, chemical bonding, and nonlocal response across molecules and materials [13]. Machine-learned KEDFs have improved this description [14, 15], but stable self-consistent use requires accurate functional derivatives, which are considerably harder to learn than energies evaluated on pre- scribed densities [16]. OF-DFT thus defines the goal, namely a density-only solver, but leaves open the best object to learn. In contrast to OF-DFT, models that infer a converged density from atomic structure [17–19] take the direct approach of learning the composite ground-state (GS) map G GS : v ext 7→ n ⋆ . Here, v ext and n ⋆ denote the ionic external Coulomb potential and corresponding ground-state density, respectively. In exact theory, this map involves the interacting many-electron ground- state problem. KS-DFT determines n ⋆ through a sequence of simpler noninteracting problems: each SCF iteration solves for the density associated with a KS potential and uses that den- sity to construct the next potential, allowing later iterations to refine earlier estimates. A direct GS prediction model must instead learn the endpoint of an arbitrarily long, XC-specific trajectory. A similar direct endpoint-learning strategy underlies machine-learned interatomic potentials, which learn the energy map v ext 7→ E[n ⋆ ] [20]. Alternatives that use near-converged quantum-mechanical (from e.g. semi-empirical calculations) features as input rather than being learned directly like Orbitall [21] greatly simplify this map and have been shown to outper- form state-of-the-art MLIPs. The benefit of retaining intermediate computation is also evident through inference-time reasoning in large language models, where difficult answers are con- structed through intermediate steps rather than in one prediction [22]. Models that predict the converged DFT Hamiltonian from atomic structure are also direct ground-state models [23–25]. In KS-DFT, the Hamiltonian can be directly constructed from the density alone. Predicting it separately therefore introduces an unnecessary basis-dependent matrix whose number of entries grows quadratically with the system size. NeuralSCF [26] learns the composite update from the density at a KS-DFT SCF iteration and atomic structure to the next density, and iterates this map to self-consistency, outperforming a matched one- shot ground-state predictor. Because the KS potential is not supplied explicitly, however, the model must learn both its XC-dependent construction, non-local Hartree screening, and then the subsequent non-interacting solution. Additionally, it depends on the XC approximation used to generate the learned SCF trajectory, as well as the basis representation used for the density, limiting its transferability between molecules and periodic systems. 2 ASelf-consistent field (SCF) cycle Molecules Build Kohn-Sham potential O(N logN) Solve O(N 3 ) DFT Kohn-Sham FNO O(N logN) Ours Density solution Self- consistent? Solids Density guess No Ground- state Yes BModel outline Kohn-Sham potentialDensity solution 푡= 1, 2, . . . , 퐿 · · · · · · · · · · · · F 푅 휃 푡 F −1 푧 푡 (r) Fourier layer 푊 휃 푡 + 휎 CAccuracy ΓXYΣΓZΣ 1 NPY 1 Z|X P −3 0 3 6 Energy (eV) CaTiO 3 (I4/mcm) E g : 2.573 eVE g : 2.469 eV Kohn-Sham FNODFT DExtrapolation ≤9202530354045 Heavy atoms 0 15 30 45 Density error (%) QM9 Training QMugs Kohn-Sham FNODirect FNO EScaling 020k40k60k Number of electrons 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 Time (h) Magnesium dislocations DFT Kohn-Sham FNO Figure 1: Summary of the Kohn–Sham Fourier Neural Operator (Kohn-Sham FNO) based self-consistent field (SCF) cycles for density optimization. (A) Conventional DFT and Kohn-Sham FNO workflows. The model replaces the orbital-solution step (see SI) while retaining the rest of a standard DFT SCF cycle. Our SCF implementations for both molecules and solids closely follow Quantum ESPRESSO (QE) 7.5 [27, 28] with minor modifi- cations for retaining quasi-linear scaling. All DFT calculations are also performed using QE. See SI. (B) Schematic of the model mapping an input potential to the density produced by solving the corresponding Kohn-Sham equations. See SI. (C) Band structures for tetragonal perovskite CaTiO 3 (MC3D [29] ID mc3d-69107; I4/mcm, No. 140) obtained using the Kohn-Sham FNO (blue) converged densities and regular DFT (black). E g values denote the corresponding band gaps. See SI. (D) Density error versus molecular size for self-consistent and direct ground-state prediction on small-organic QM9 [30] and larger drug-like QMugs molecules [31]. See Fig. 3. (E) Measured total SCF cycle wall time versus electron count for QE and Kohn-Sham FNO. See SI. 3 Our Approach: A natural and universal learning target already exists in every Kohn–Sham calculation. At each SCF iteration, the density n defines an effective potential v KS [n]. Solving its non-interacting Hamiltonian and occupying the eigenstates produces the output density, − 1 2 ∇ 2 + v KS [n](r) φ p (r) = ε p φ p (r)=⇒ n out (r) = X p f p |φ p (r)| 2 (1) The SCF converges when n out = n, but every preceding iteration also supplies an exact sample of the same forward universal operator, G (n,T) KS : v KS [n](r)7−→ n out (r),T s [n out ] ,(2) where the electron number, boundary conditions, and occupations are implicit. Both the density and kinetic energy are therefore one non-interacting solve away from a prescribed KS potential, and together form a complete orbital-free DFT target. Stable density optimization has historically been the harder test for OF-DFT [13, 16] and hence in this work we learn only the density component, i.e. v KS [n](r)7−→ n out (r), denotedG KS for performing stable orbital-free SCFs at inference. Consequently, we obtain total energies and electronic spectra via a fixed-density post-SCF diagonalization in the following results. This forward formulation is the inverse of the relation used in traditional OF-DFT. For a local KS potential, minimization of the total energy at fixed electron number gives δT s [ρ] δρ(r) ρ=n out = μ− v KS [n](r),(3) where μ enforces the electron number. A KEDF therefore maps the density to the kinetic energy, and its functional derivative recovers the generating effective potential up to the additive constant μ : n out (r) 7→ μ− v KS [n](r), which is the inverse of the eq. 2 operator. By contrast, the two outputs in Eq. 2 map the generating potential directly to its density and kinetic energy. Compared to direct prediction models, v KS [n] 7→ n out is a more elementary map defined via a single non-interacting Hamiltonian solution. Compared to MLIPs, a simpler forward counterpart is v KS [n ⋆ ]7→ T s [n ⋆ ], after which the remaining energy terms are explicit. Section 2.1 discusses this energy analogy further, while Fig. 3 directly compares forward KS and one-shot density prediction. Equation 2 gives the local form of the KS operator while our calculations use nonlocal pseudopotentials [32, 33] to smooth the core region and allow us to use uniform-grid fast Fourier transforms (FFT). For each fixed set of nonlocal projectors, we combine their operator with the kinetic energy in a generalized non-interacting kinetic functional defined by constrained search [34, 35]. This preserves the joint forward map and its inverse Euler relation for the pseudopotential problem. See SI. This defines one operator across SCF iterations, chemical composition, geometry, and exchange-correlation (XC) approximations. These inputs modify the KS potential but not the non-interacting solution operator being learned. Potential-density pairs from different XC approximations can therefore be combined in training. We test this separation by applying a PBE-trained model [36] to PBEsol SCFs [37] without retraining [Fig. 4]. The KS map is an operator between spatial fields, making neural operators a natural model class for this framework. Neural operators learn maps between function spaces rather than between fixed-dimensional vectors. We specifically utilize Fourier neural operator (FNO) since the Fourier layers of an FNO capture global nonlocality with quasi-linear O(N g logN g ) scaling [38, 39] on the same grid used to construct the Kohn-Sham potential. The complete update therefore scales asO(N g logN g ) for N g grid points. Alternative architectures, e.g., graph neural networks rely either on finite spatial cutoffs or on globally coupled operations whose cost grows quadratically with system size, e.g., transformers. 4 In this work, we introduce a new domain-invariant, SE(3)-equivariant FNO architecture. Unlike the standard FNO, whose spectral weights are tied to a fixed domain size, the domain- invariant FNO evaluates radially factorized filters as functions of physical reciprocal-space ra- dius. The same learned filters can therefore be sampled on the reciprocal-space grids of differ- ently sized domains, allowing one architecture to treat molecules and materials across system sizes. The radial filters and spherical mode truncation make the spectral convolution rotation equivariant; together with its translation equivariance, this yields full SE(3) equivariance. We first show that learning the KS-map leads to significantly improved stability, and chemi- cal extrapolation than the inverse OF-DFT and direct ground-state learning approaches respec- tively. Following this, we show, for the first time, a single method converging SCF calculations without invoking the auxiliary Kohn-Sham orbitals across organic molecules of varying sizes as well as insulating and metallic systems spanning the first 5 rows of the periodic table. The corresponding converged densities are shown to produce observables nearly indistinguishable in quality from a regular KS-DFT calculation. This is verified using electronic and structural observables, including band structures, band gaps, density of states, equation of state curves, equilibrium volumes, and bulk moduli. The computational advantage over conventional KS- DFT is particularly large for periodic systems since the Kohn–Sham FNO inference loop evolves only the k-independent, lattice-periodic potential and density on the unit cell and requires no explicit k-point sampling [40] skipping N iterations × N k diagonalizations. Finally, after training on cells with at most 364 atoms, we obtain stable SCFs for magnesium dislocations containing up to 8,250 atoms on one GPU, empirically verifying the expected near-linear scaling. 2 Results & Discussion Unless otherwise stated, results are obtained using a single domain-invariant, radially factorized FNO trained on density-potential pairs from 2,004 molecules and 6,500 solids at their equilib- rium geometries (Fig. 3A). When the domain-invariant FNO learns the KS solution operator and is embedded in the SCF cycle, we refer to the resulting model & method interchangeably as Kohn-Sham FNO for brevity. At inference, the model replaces the KS eigensolver within a fixed-point SCF iteration to optimize the densities. Fig. 1 summarizes the SCF workflow and model architecture, with details of the training data, architecture, optimization procedure, and quasi-linear scaling SCF implementation. See SI. Algorithm 1 Self-consistent field cycle with Kohn–Sham FNO. Input: Atomic structure, convergence tolerance ε, and maximum iteration count I max . 1 Construct v loc ext and n SAD . Initialize n 0 = n SAD and the mixing factor α 0 . 2 for i = 0,...,I max − 1 do 3Construct v KS = v loc ext + v H [n i ] + v xc [n i ]. 4Predict en FNO i = n SAD + F θ [v KS ,v loc ext ] and project n FNO i =P N e [en FNO i ]. 5Evaluate the relative fixed-point residual R rel i = 1 N e R |n FNO i (r)− n i (r)| dr. 6if R rel i < ε then return n ⋆ = n FNO i . 7Mix to obtain a new density n i+1 =P N e [(1− α i )n i + α i n FNO i ]. 8 end for 9 return unconverged Output: A self-consistent density n ⋆ , potential v KS [n ⋆ ], and Hamiltonian h KS [n ⋆ ] or an explicit conver- gence failure. 2.1 Stable optimization & improved extrapolation through the forward map The learning target determines both the numerical problem presented to the model and how its errors enter the final density. We compare three targets that seek to eliminate the orbital 5 A 020406080 SCF iteration 0 5 10 15 Density error, | ∆ n | (%) Kohn–Sham FNO F θ :v KS [n](r)7→n out (r) 020406080 SCF iteration Kinetic Potential FNO F θ :n out (r)7→v T s [n](r) B 0 4 8 12 0 0.5 1 Eigenmode quantile −8 −4 0 4 8 Forward response χ s =δn/δv KS 0 4 8 12 0 0.5 1 Eigenmode quantile Inverse response −χ −1 s ≈δv T s /δn SCF iteration 0 50 100 150 200 250 300 350 Data points per bin Jacobian eigenspectrum (log 10 ) Figure 2: Learning the forward orbital-free KS map avoids inverse-response insta- bility. (A) Density error relative to PBE [36] ground-states for 100 held-out small-organic molecules from the QM9 dataset [30]. See SI. Kohn–Sham FNO (blue) maps the KS potential to its output density and is iterated to a fixed-point. The kinetic-potential FNO (red) uses the same domain-invariant FNO model as the backbone and learns the density to kinetic-potential relation and is used to solve the Euler equation based on the implementation in DFTpy [40]. See SI. Both models use the same training data (Fig. 3A) and FNO backbone. Thin curves show individual molecules, bold curves show their mean, and circles and crosses mark converged and failed calculations, respectively. (B) Binned eigenspectra of the projected forward response χ s (left) and inverse response −χ −1 s (right) over 1,100 PBE SCF states for the same molecules. Modes are ordered from strongest to weakest forward response; points mark bin medians and surface color gives the number of eigenvalues per bin. The responses were computed using all- electron calculations in Gaussian basis sets with the PySCF [41] and KS-PIES [42] packages. See SI. 6 bottleneck: the forward KS map used herein, its inverse kinetic-potential map used in traditional orbital-free density optimizations, and direct prediction of the ground-state density from the atomic environment. The forward KS map lies between the other two. It neither inverts the potential-to-density response nor compresses the complete nonlinear ground-state problem into a single prediction. We first test the consequences of this reformulation for density-optimization stability and then for chemical generalization. For a non-interacting v-representable density, the Euler equation identifies the kinetic po- tential as v T s [n out ] = μ−v KS [n]. See SI. Learning a kinetic potential therefore reverses the same relation learned by the Kohn–Sham FNO. The distinction is consequential. The forward map is defined by a one-particle ground-state problem for every KS potential. Its inverse is defined only for non-interacting v-representable densities. Recovering the corresponding potential gen- erally requires an iterative inversion such as the Wu-Yang method [43], and a complete set of constraints characterizing the v-representable domain is unknown [44, 45]. Trial densities gen- erated during orbital-free optimization can therefore leave the v-representable domain sampled by exact inverse labels. A scalar kinetic-energy functional is subject to the same issue when used variationally. Its forward counterpart is the explicit map v KS [n](r)7−→ T s [n out ],(4) where n out is the density produced by the input potential. The density map v KS [n]7→ n out and the energy map in Eq. 4 therefore follow from the same single one-particle diagonalization. The pseudopotential calculations below obey the same relations with T s replaced by the generalized functional e T s . See SI. We isolate the effect of mapping direction by training a kinetic-potential FNO with the same backbone architecture and the same training data as the Kohn–Sham FNO. The forward model learns v KS [n]7→ n out , whereas the inverse model learns n out 7→ μ−v KS [n], using the generalized kinetic potential of the nonlocal-pseudopotential formulation. See SI. Starting from the same superposition-of-atomic-densities (SAD) guess, we apply the forward model in a fixed-point SCF procedure and the inverse model in a potential-only optimization. See SI. The resulting trajectories, measured against the PBE reference density using the density-error metric, are qualitatively different [Fig. 2A]. See SI. All forward-model calculations reach a self-consistent fixed point, progressively reducing the density error with respect to the reference PBE density from approximately 10–18% for the initial guess to less than 1%. By contrast, the inverse- model optimizations fail after only a few iterations: the learned Euler residual does not provide a descent direction for which the backtracking procedure can continue, and the density remains close to its initial error. This comparison uses the same training data, model architecture, as well as equivalent training and inference strategies; its primary difference is whether the KS relation is learned in the forward or inverse direction. The same difficulty arises when a learned energy functional is differentiated to optimize the density. Remme et al. obtained stable QM9 SCFs by training the combined kinetic and XC energy E TXC [n] on about 107,000 molecules and 21 perturbed SCF states per molecule, giving roughly 2.25 million training labels [16]. In their ablation, ordinary SCF training data led to a 28% convergence-failure rate, compared with no failures using the perturbed data. Kohn–Sham FNO is only trained on 59,500 total ordinary SCF training labels from 8,504 molecular and periodic structures, while converging SCFs in both domains. Although the two architectures are not directly comparable, the large contrast shows how much off-equilibrium data may be needed to learn the inverse density-to-potential map stably. The stability difference between the two models stems from how the two learned maps are conditioned. A linear perturbation of the KS potential produces δn = χ s δv KS , where χ s = δn/δv KS is the non-interacting density response. See SI. On the particle-number conserving subspace, linearization of the Euler equation instead gives δv T s /δn =−χ −1 s . We evaluated the spectrum of χ s independently for every recorded SCF iteration of the same 100 QM9 molecules 7 using all-electron KS response calculations. See SI. The forward spectrum in Fig. 2B contains many weak-response modes with small |λ j (χ s )|. Such modes suppress potential errors in the forward map, but inversion turns them into sensitivities proportional to|λ j (χ s )| −1 and amplifies density errors by several orders of magnitude. The mirrored spectra persist throughout the SCF trajectories, so the inverse conditioning problem is not confined to the converged density. It applies whether the kinetic potential is learned directly, as in our controlled inverse model, or obtained by differentiating a learned kinetic-energy functional. We next compare the forward operator with the opposite extreme: predicting the converged density directly from the ionic environment. Direct prediction models can be accurate and transferable when symmetry, locality, and suitable density representations are built into the architecture [17–19, 46, 47]. We again hold the FNO backbone, grid representation, and training structures fixed and change the target. The direct model learns v loc ext 7→ n ⋆ Ref in one evaluation, whereas the Kohn–Sham FNO learns one non-interacting solution and recovers n ⋆ Ref through self-consistency. This distinction separates what must be learned from what can be evaluated explicitly. Direct prediction must encode the entire PBE fixed point, which involves constructing the ground-state v KS [n ⋆ ] first via an arbitrarily long non-linear SCF trajectory involving Hartree and XC feedback and long-range charge redistribution. The Kohn–Sham FNO instead receives the reconstructed Hartree and XC potentials at every iteration and repeatedly applies the same one-particle solution operator. An error in one update changes the next input potential and can be corrected by subsequent updates, while density mixing damps unstable steps. Moreover, each conventional SCF trajectory supplies multiple exact potential-density training pairs for the forward operator at no additional electronic-structure cost, rather than only its final density. On held-out QM9 molecules drawn from the same small-organic chemical space as the molec- ular training set, both approaches are accurate [Fig. 3B, left]. The Kohn–Sham FNO and direct FNO obtain density errors of 0.625% and 0.662%, respectively, with dipole errors of 0.010 and 0.011 D per electron. Their quadrupole, electrostatic-potential, Hartree-energy, and XC-energy errors are similarly close. The direct model therefore has sufficient capacity to represent the ground-state map within its training distribution. The difference emerges for QMugs drug-like molecules [31], which are larger than QM9 and introduce S, Cl, and P into molecular environ- ments absent from the molecular training set [Fig. 3B, right]. Although these elements occur in the solid-state portion of the shared training data, their bonding environments and molecular domain sizes are out of distribution. Under this combined size, composition, and chemical- environment shift, the direct-model density error rises to 9.97%, compared with 2.23% for the self-consistent Kohn–Sham FNO. The self-consistent cycle also reduces the dipole error from 0.237 to 0.026 D per electron, the quadrupole error from 0.428 to 0.031 e ̊ A 2 per electron, and the electrostatic-potential error from 40.8 to 3.28 mHa. The Hartree- and XC-energy errors fall from 5.91 to 0.267 mHa per electron and from 9.92 to 2.32 mHa per electron, respectively. Figure 3C shows that this extrapolation gap widens systematically with molecular size: from 20 to 45 heavy atoms, the direct-model density error increases from about 4% to 41%, whereas the self-consistent Kohn–Sham FNO grows only from about 1.5% to 4%. This contrasting size dependence is consistent with SCF feedback correcting intermediate errors rather than requiring a one-shot model to extrapolate the complete density of an increasingly large system. The same idea applies to energy prediction. MLIPs learn the complete ground-state map from the external potential represented by the atomic graph to the ground-state energy [20], v ext 7−→ E[n ⋆ ].(5) Once an SCF has constructed the ground-state KS potential, the forward alternative is v KS [n ⋆ ]7−→ (n ⋆ ,T s [n ⋆ ])7−→ E[n ⋆ ],(6) where the external-potential, Hartree, XC, and nuclear contributions are evaluated explicitly. The first map asks the model to infer the entire electronic ground-state solution directly from the 8 ATraining set composition 2,004 QM9 molecules + 6,500 MC3D crystals BChemical extrapolationC Kohn-Sham FNODirect FNO H 28.1 He Li 7.2 Be 3.7 B 6.0 C 31.2 N 24.1 O 49.0 F 3.3 Ne Na 7.1 Mg 10.1 Al 12.9 Si 10.3 P 7.1 S 7.9 Cl 4.4 Ar K 5.9 Ca 14.2 Sc 2.4 Ti 4.4 V 0.7 Cr 0.2 Mn 0.1 Fe 0.1 Co 0.5 Ni 1.7 Cu 1.6 Zn 10.1 Ga 3.6 Ge 5.1 As 3.8 Se 3.5 Br 2.1 Kr Rb 0.8 Sr 5.2 Y 4.1 Zr 3.5 Nb 0.6 Mo 0.6 Tc 0.1 Ru 0.4 Rh 0.7 Pd 1.4 Ag 0.8 Cd 0.7 In 0.9 Sn 1.4 Sb 1.3 Te 1.0 I 0.4 Xe 0.1 Cs 1.0 Ba 2.0 Hf 0.7 Ta 0.7 W 0.6 Re 0.2 Os 0.1 Ir 0.7 Pt 1.3 Au 1.0 Hg 0.4 Tl 0.5 Pb 0.8 Bi 0.7 PoAtRn 0 10 20 30 40 Structures containing element (%) 0 4 8 |∆n| (%) 0.0 0.1 0.2 |∆μ| (D/N e ) 0.0 0.2 0.4 ‖∆Q‖ F (e ̊ A 2 /N e ) QM9QMugs 0 20 40 |∆V ESP | (mHa) QM9QMugs 0.0 2.5 5.0 E H [∆n] (mHa/N e ) QM9QMugs 0 4 8 |∆E xc | (mHa/N e ) ≤920253035 Heavy atoms 0 8 16 24 Density error, | ∆ n | (%) QM9QMugs Figure 3: Improved extrapolation via the self-consistent Kohn–Sham map. (A) El- emental composition of the shared training set of 2,004 QM9 molecules [30] and 6,500 MC3D crystals [29]. Each value is the percentage of structures containing that element; gray denotes absence from training. (B) Mean errors of the self-consistent Kohn–Sham FNO (blue) densities and direct ground-state predicted densities on 100 held-out QM9 molecules containing the ele- ments C, H, N, O, F and 200 larger drug-like molecules from the QMugs dataset [31] containing the elements C, H, N, O, F, S, Cl, P. From left to right, the metrics are the normalized L 1 density error, Euclidean norm of the dipole-vector error, Frobenius norm of the quadrupole- tensor error, electrostatic-potential error, Coulomb-weighted density-error metric, and absolute XC-energy error; quantities labeled /N e are normalized per valence electron. (C) Density error versus number of heavy atoms. The QM9 bin at ≤9 heavy atoms consists of the same 100 test molecules from Fig. 2.1. The QMugs bins from 20− 35 heavy atoms each consist of 50 randomly sampled ground-state conformers. All errors are relative to converged Quantum ESPRESSO PBE calculations. Data construction, the direct model, and the density metric: see SI. 9 atoms; the second asks it to reproduce one non-interacting KS solution after the self-consistent potential is already known. The extrapolation advantage in Fig. 3 therefore motivates the same hypothesis for energy prediction. We do not train a forward energy model here, however, and leave this hypothesis for a future study involving a variationally consistent energy, density Kohn–Sham model. Together, Figs. 2 and 3 identify the forward KS map as the favorable intermediate learning target among the three formulations tested. Relative to the inverse kinetic-potential map, it avoids inverse-response amplification and is queried only on potentials for which the exact forward target is defined. Relative to direct ground-state prediction, it factorizes a composite, XC-dependent fixed point into repeated applications of an XC-independent one-particle solve with explicit physical feedback. In this sense, the KS map is the optimal target considered here: it removes the repeated eigensolution that dominates the cost while preserving the parts of KS-DFT that stabilize and transfer the calculation. This is not a claim of formal optimality over every possible representation, but it is the only target in these controlled comparisons that combines stable self-consistent optimization with robust extrapolation. 2.2 Periodic self-consistency without k-points Machine-learned orbital-free DFT has previously been applied to periodic systems, but demon- strations have been limited to selected material classes and have relied on local pseudopotentials constructed specifically for orbital-free calculations [48, 49]. Here, we use the same Kohn–Sham FNO that drives the molecular SCFs in the preceding section, trained jointly on molecular data and chemically diverse metallic and semiconducting systems from MC3D [29]. See SI. Its parameters, architecture, and inference procedure are unchanged across molecules, metals, and semiconductors. Moreover, it reproduces densities from the transferable nonlocal norm- conserving pseudopotentials routinely used in KS-DFT, whose orbital-dependent projectors cannot be treated directly within conventional density-only orbital-free DFT. The computational advantage over conventional KS-DFT is particularly large for periodic systems. As in conventional periodic orbital-free DFT, the Kohn–Sham FNO inference loop evolves only the k-independent, lattice-periodic potential and density on the unit cell and re- quires no explicit k-point sampling [40]. The distinction lies in how the density update is obtained. Conventional orbital-free DFT derives it from an approximate kinetic functional, whereas the Kohn–Sham FNO learns the composite KS operation that solves the Hamiltonian at every sampled k-point and contracts the occupied states into the Brillouin-zone-summed density. A single model evaluation therefore replaces all k-resolved diagonalizations in each conventional KS-DFT iteration while retaining the KS potential-density feedback. See SI. On 200 held-out MC3D crystals, the main model drives the fixed-point residual smoothly toward zero, with all but one calculation converging within 80 iterations [Fig. 4A]. The resulting mean density errors are 0.75% for semiconductors and 1.46% for metals. To test whether the converged FNO density also preserves observables that are sensitive to the ground-state Hamiltonian, we consider the held-out metallic Sr 3 SnO cubic antiperovskite (MC3D ID mc3d-54098, space group Pm 3m) [29]. Its multiband electronic structure, with several dispersive states near the Fermi level, provides a nontrivial test of both spectral and energetic fidelity. One fixed-density post-SCF diagonalization using the FNO density reproduces the self-consistent PBE band structure and density of states [Fig. 4B]. At the reference volume, its Harris total energy differs from self-consistent PBE by only 11.4 meV/cell (2.29 meV/atom). The equation of state tests transfer beyond the equilibrium geometries used for training. See SI. Across seven clamped-ion volumes spanning V/V ref = 0.94–1.06, the PBE FNO and DFT energy curves are nearly coincident: their fitted equilibrium volumes differ by 0.23% and their bulk moduli by 1.1%. Changing the volume presents a new KS potential, but the learned task remains the same one-particle solution map rather than direct prediction of a geometry-specific ground-state property. 10 AMC3D test set Kohn–Sham FNO BSr 3 SnO Kohn–Sham FNODFT 020406080 SCF iteration 10 −3 10 −2 10 −1 10 0 Fixed-point residual converged failed SemiconductorsMetals 0.0 0.5 1.0 1.5 2.0 Density error | ∆ n | (%) 0.75 1.46 ΓXMΓR X|M R −6 −3 0 3 6 Energy (eV) 05 DOS 125130135140145 Volume ( ̊ A 3 /cell) 0 20 40 60 80 100 E − E 0 (meV/cell) XCSCF V 0 ( ̊ A 3 /cell) B 0 (GPa) FNO138.7243.47 DFT138.3943.96 FNO133.0045.15 DFT132.2048.92 PBE PBEsol Figure 4: Generalization across periodic systems, SCF trajectories, geometries and XC approximation. (A) Fixed-point trajectories and final density errors for a held-out MC3D test set composed of 100 semiconductors and 100 metals as characterized by converged PBE calculations. Chemical composition of these crystals spans the same subspace of the periodic table as in Fig. 3A. Thin curves denote individual crystals, dots and crosses mark converged and failed calculations, and bars summarize the mean errors for semiconductors and metals. (B) Results for metallic cubic Sr 3 SnO (MC3D ID mc3d-54098; Pm 3m, No. 221). Left, band structure and density of states from one fixed-density post-SCF diagonalization using the FNO density (blue) and from self-consistent PBE (black dashed), each referenced to its Fermi energy. Right, clamped-ion PBE (squares) and PBEsol [37] (circles) equations of state; solid blue and dashed black curves denote FNO and DFT, respectively, and the table reports fitted equilibrium volumes V 0 and bulk moduli B 0 . The PBEsol test uses the PBE-trained model without fine- tuning and replaces only the XC potential in the KS potential construction during inference. See SI. Training, post-SCF observables, and EOS fitting: see SI. 11 We then change the XC approximation itself. PBEsol is a generalized gradient approxima- tion (GGA)-XC designed to improve the equilibrium properties of densely packed solids and their surfaces, making it a natural geometry-focused transfer test [37]. Without retraining or fine-tuning, we replace only the PBE XC term in the Kohn–Sham potential with the PBEsol potential at each FNO iteration; PBEsol is also used for energy evaluation, while all other model inputs and inference settings remain unchanged. See SI. The FNO based SCFs follow the resulting shift of the entire EOS: it predicts V 0 = 133.00 ̊ A 3 /cell and B 0 = 45.15 GPa, compared with 132.20 ̊ A 3 /cell and 48.92 GPa from self-consistent PBEsol. The corresponding errors of 0.61% and 7.7% are larger than for PBE but remain small despite the absence of any PBEsol data. This functional transfer follows from the operator definition. The noninteracting output density is determined by the complete local KS potential, not by the XC approximation used to construct it. Potential-density pairs generated with different density-only XC approximations are therefore samples of the same KS operator and can be mixed in training without a functional label. Similarly, at inference, the input KS potential to the model can be constructed using any density-dependent XC-potential approximation. Fine-tuning would improve coverage, but the learning target need not be redefined. Thus, even before learning the kinetic-energy output of the joint KS operator, the density model can replace the full orbital-based periodic SCF loop involving accurate norm-conserving pseudopotentials. Density and density-derived observables require no subsequent orbital calculation. A conventional fixed-density solve is needed only when orbital-resolved spectra are needed, rather than at every k point of every SCF iteration. 2.3 Magnesium Dislocation Finally, we move from ideal periodic crystals to extended defects at realistic length scales. A dislocation breaks primitive translational symmetry and produces a long-ranged elastic field, so first-principles calculations require supercells containing thousands of atoms. We consider ⟨c + a⟩ screw dislocations in Mg, the lightest structural metal. Its limited ductility is linked to the relative stability and cross-slip of competing pyramidal-I and pyramidal-I cores, which can be altered through alloying [50–52]. Accurately resolving their small energy difference therefore requires cell sizes at which repeated orbital diagonalization becomes prohibitive. This problem was used by the 2019 ACM Gordon Bell Prize finalist study of Das et al. to demonstrate large-scale metallic DFT with finite elements (FE) based DFT-FE package [7, 53]. A full ground-state calculation for 6,164 Mg atoms required 56 SCF iterations on 1,300 nodes of the Summit supercomputer (7,800 NVIDIA V100 GPUs), whereas the largest system, containing 10,508 atoms and 105,080 electrons, used 3,800 nodes (22,800 GPUs) for one SCF iteration. To test the practical reach of the quasi-linear Kohn–Sham FNO, we therefore attempt to converge the same class of metallic-defect SCFs on a single, modern high-memory GPU. As a test, we began with the same pretrained Kohn–Sham FNO used for the molecular and solid calculations above, without any defect-specific fine-tuning. The SCF trajectories immediately revealed that additional training was needed as the fixed-point residual rapidly increased for every tested cell and none of the calculations converged [Fig. 5B]. This warning required no reference density. The predicted density was inconsistent with the KS potential reconstructed from it, causing the next model update to move farther from a fixed point. Iterative prediction thus provides a built-in reliability diagnostic that a direct ground-state prediction model would not. The residual is not a formal error bound, and convergence alone does not certify agreement with PBE ground-state, but SCF divergence immediately identifies a model that should not be trusted. We therefore fine-tuned the same pretrained model on 1,203 Mg structures containing 20-364 atoms. See SI. The data combine bulk and strained cells, surfaces, generalized stacking faults, and local environments cut from pyramidal-I and pyramidal-I dislocation cores; the architec- ture and inference procedure were unchanged. The fine-tuned model converges every tested 12 3×10 6 10 7 3×10 7 10 8 Grid points 10 0 10 1 10 2 10 3 Time per SCF iteration (s) Quantum ESPRESSO p= 3.37 Kohn–Sham FNO p= 1.03 0255075100 SCF iteration 10 −3 10 −2 10 −1 10 0 Fixed-point residual Kohn–Sham FNO pre-trained fine-tuned 203 ̊ A 191 ̊A 8,250 Mg atoms PyrI−PyrII 600 atoms 10 ̊ A 20 40 60 80 Valence electrons ( × 10 3 ) AComputational scalingBSCF convergence CElectrostatic-potential differenceDDislocation core, magnified −0.050+0.05 (Ha) Figure 5: Near linear-scaling density optimization for large metallic dislocation cells. (A) Time per complete SCF iteration t versus the number of real-space grid points N g on logarithmic axes. Symbols show Quantum ESPRESSO (black squares) and Kohn–Sham FNO (blue circles) timings; lines are empirical power-law fits t ∝ N p g , where p is the fitted scaling exponent. Quantum ESPRESSO calculations were run on 192-core CPU nodes, whereas the FNO used one NVIDIA B300 GPU. See SI. (B) Fixed-point residuals for pyramidal-I and pyramidal-I dislocation cells using the base model (dashed), trained on the QM9 and MC3D training set in Fig. 3A, and the model fine-tuned on Mg-defect data (solid; see SI); color denotes valence-electron count. (C) Difference in the electrostatic potential v loc ext + v H [n] between the 8,250-atom pyramidal-I and pyramidal-I cells using Kohn–Sham FNO converged densities. Grey dots denote the Mg atom locations from the pyramidal-I cell. The heat map shows the 99.5th percentile difference. (D) The dislocation core at full resolution. The region boxed in (C) contains 600 atoms and is magnified roughly 3.9x. Training, fine-tuning, and inference data generation, convergence, and timing calculations specific to these results: see SI. 13 dislocation cell to the same relative-L 1 fixed-point threshold of 10 −3 used throughout this work [Fig. 5B]. See SI. The largest calculation contains 8,250 Mg atoms and 82,500 valence electrons and runs on one NVIDIA B300 GPU without reducing the real-space grid resolution. Their trajectories also contain occasional residual spikes, which we attribute to numerical instabilities when applying the FNO on such large grids. Adaptive density mixing nevertheless returns each trajectory to convergence. See SI. Timing the complete SCF update confirms the expected near-linear scaling of the approach empirically [Fig. 5A]. Power-law fits against the number of grid points give p = 1.03 for the Kohn–Sham FNO, including the largest system, compared with p = 3.37 for Quantum ESPRESSO. The FNO result is consistent with the expected O(N g logN g ) complexity of the complete FFT-based SCF iteration. FNO timings used one B300 GPU, whereas Quantum ESPRESSO timings were linearly rescaled by core count to estimate execution on a 192-core AMD EPYC 9655 node (see SI); their absolute times are therefore not a hardware-matched comparison. The upward deviation of the largest FNO point likely reflects memory pressure near the capacity of the B300, which set the maximum system size tested. Because converging PBE reference densities for the full dislocation cells with Quantum ESPRESSO was computationally prohibitive, density errors were evaluated on core-centered crops containing up to 528 Mg atoms, the largest tractable size on nodes with 750 GB of memory. See SI. Across the held-out pyramidal-I and pyramidal-I cells, the mean FNO density error remains 0.33–0.35% with no systematic growth with size [see SI]. At 8,250 atoms, where no reference is available, the difference in the FNO-derived electrostatic potential (ESP) v loc ext (r) + v H [n](r) between the two core structures remains localized around the dislocation [Fig. 5C]. Although not an independent accuracy test, this field shows that the large converged densities resolve distinct core environments. Figure 5 tests only the density component of the joint KS operator. The converged densities can be passed to a single fixed-density orbital calculation, such as a DFT-FE calculation, or paired with a future forward kinetic-energy model to obtain total energies without repeating the full orbital SCF trajectory. Here, we establish stable self-consistency for an extended metallic defect with more than 8× 10 4 electrons on one GPU. 3 Conclusion In this work, we identify the Kohn–Sham map as a complete orbital-free DFT target and construct a domain-invariant, radially factorized Fourier neural operator to learn its density component. The model replaces one non-interacting orbital solution at an arbitrary SCF itera- tion, rather than learning the inverse density-to-potential relation or compressing the complete ground-state calculation into one prediction. This factorization removes the repeated orbital diagonalization bottleneck from regular KS-DFT SCFs while retaining explicit Hartree and XC construction, density mixing, and self-consistent feedback. Controlled comparisons using the same data and model backbone show that this choice stabilizes density optimization relative to the inverse kinetic- potential map and extrapolates more reliably than direct ground-state prediction. The same Kohn–Sham FNO drives molec- ular, semiconducting, and metallic SCFs with transferable norm-conserving pseudopotentials as training reference. It can learn from truncated, unconverged reference trajectories, transfer from PBE to PBEsol without retraining, and produce densities that recover KS-DFT quality electronic observables, spectra and equations of state. After fine-tuning only on Mg cells con- taining at most 364 atoms, the model converges dislocation-cell densities containing up to 8,250 atoms and 82,500 valence electrons on one GPU, with an empirical scaling of O(N 1.03 ). The fixed-point residual additionally provides a direct diagnostic when the learned operator is being applied outside its reliable domain. The present model learns only the density component of the joint Kohn–Sham operator. 14 Orbital-resolved observables and total energies can already be obtained at KS-DFT accuracy with one fixed-density post-SCF calculation, but eliminating that final orbital solve requires learning the corresponding forward kinetic-energy (or finite-smearing free-energy) output. Be- cause both quantities are produced by the same one-particle problem, we will pursue a variation- ally consistent framework for the complete map in future work. More broadly, these results show that machine learning can extend electronic-structure calculations most effectively by replacing their dominant repeated operation while preserving the physical iteration that constructs and validates the ground state. 4 Acknowledgements A. Anandkumar is supported in part by the Bren endowed chair, ONR (MURI grant N00014- 18-12624), DARPA ExpMath HR0011, and by the AI2050 senior fellow program at Schmidt Sciences. S. Sharma was supported by the U.S. Department of Energy, Office of Science, Office of Basic Energy Sciences, Fuels from Sunlight Hub under Award Number DE-SC0021266. D. Khan acknowledges support from the Pritzker AI+Science fund and DARPA Biological Tech- nologies HR0011. M. Hanisch is supported by the Kortschak Scholars Program. E. Xie is supported through the Caltech Summer Undergraduate Research Fellowship. D. Khan and M. Hanisch acknowledge discussions with Tommaso Chiarotti, Valentin Duruisseaux, Chuwei Wang, Vignesh Bhethanabotla, Chenghan Li, and Garnet K. L. Chan. We acknowledge TACC, Schmidt Sciences and Caltech HPC for computing resources. References [1] Pierre Hohenberg and Walter Kohn. Inhomogeneous electron gas. Physical Review, 136 (3B):B864, 1964. [2] Walter Kohn and Lu Jeu Sham. Self-consistent equations including exchange and correla- tion effects. Physical Review, 140(4A):A1133, 1965. [3] Axel D Becke. Perspective: Fifty years of density-functional theory in chemical physics. The Journal of Chemical Physics, 140(18), 2014. [4] Robert O Jones. Density functional theory: Its origins, rise to prominence, and future. Reviews of Modern Physics, 87(3):897–923, 2015. [5] Bing Huang, Guido Falk von Rudorff, and O Anatole von Lilienfeld. The central role of density functional theory in the ai age. Science, 381(6654):170–175, 2023. [6] Daniel J Cole and Nicholas D M Hine.Applications of large-scale den- sity functional theory in biology.JournalofPhysics:CondensedMat- ter, 28(39):393001, aug 2016.doi:10.1088/0953-8984/28/39/393001.URL https://doi.org/10.1088/0953-8984/28/39/393001. [7] Sambit Das, Phani Motamarri, Vikram Gavini, Bruno Turcksin, Ying Wai Li, and Brent Leback.Fast, scalable and accurate finite-element based ab initio calcula- tions using mixed precision computing: 46 PFLOPS simulation of a metallic dislo- cation system. In Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis, pages 1–11, Denver Colorado, Novem- ber 2019. ACM.ISBN 978-1-4503-6229-0.doi:10.1145/3295500.3357157.URL https://dl.acm.org/doi/10.1145/3295500.3357157. 15 [8] Ryan Pederson, John Kozlowski, Ruyi Song, Jackson Beall, Martin Ganahl, Markus Hauru, Adam GM Lewis, Yi Yao, Shrestha Basu Mallick, Volker Blum, et al. Large scale quantum chemistry with tensor processing units. Journal of Chemical Theory and Computation, 19 (1):25–32, 2022. [9] Walter Kohn. Density functional and density matrix method scaling linearly with the number of atoms. Physical Review Letters, 76(17):3168, 1996. [10] David R Bowler and Tsuyoshi Miyazaki. O(n) methods in electronic structure calculations. Reports on Progress in Physics, 75(3):036503, 2012. [11] Lin Lin, Mohan Chen, Chao Yang, and Lixin He. Accelerating atomic orbital-based elec- tronic structure calculation via pole expansion and selected inversion. Journal of Physics: Condensed Matter, 25(29):295501, 2013. [12] Roi Baer, Daniel Neuhauser, and Eran Rabani. Self-averaging stochastic kohn-sham density-functional theory. Physical Review Letters, 111(10):106402, 2013. [13] Wenhui Mi, Kai Luo, SB Trickey, and Michele Pavanello. Orbital-free density functional theory: An attractive electronic structure method for large-scale first-principles simula- tions. Chemical Reviews, 123(21):12039–12104, 2023. [14] He Zhang, Siyuan Liu, Jiacheng You, Chang Liu, Shuxin Zheng, Ziheng Lu, Tong Wang, Nanning Zheng, and Bin Shao. Overcoming the barrier of orbital-free density functional theory for molecular systems using deep learning. Nature Computational Science, 4(3): 210–223, 2024. [15] Min Chen, Michele Pavanello, Wenhui Mi, Manabu Ihara, and Sergei Manzhos. Machine learning-enhanced orbital-free density functional theory. Journal of Chemical Theory and Computation, 22(7):3127–3143, 2026. [16] Roman Remme, Tobias Kaczun, Tim Ebert, Christof A Gehrig, Dominik Geng, Gerrit Gerhartz, Marc K Ickler, Manuel V Klockow, Peter Lippmann, Johannes S Schmidt, et al. Stable and accurate orbital-free density functional theory powered by machine learning. Journal of the American Chemical Society, 147(32):28851–28859, 2025. [17] Felix Brockherde, Leslie Vogt, Li Li, Mark E Tuckerman, Kieron Burke, and Klaus-Robert M ̈uller. Bypassing the kohn-sham equations with machine learning. Nature Communica- tions, 8(1):872, 2017. [18] Peter Bjørn Jørgensen and Arghya Bhowmik. Equivariant graph neural networks for fast electron density estimation of molecules, liquids, and solids. npj Computational Materials, 8(1):183, 2022. [19] Chenghan Li, Or Sharir, Shunyue Yuan, and Garnet Kin-Lic Chan. Image super-resolution inspired electron density prediction. Nature Communications, 16(1):4811, 2025. [20] Oliver T Unke, Stefan Chmiela, Huziel E Sauceda, Michael Gastegger, Igor Poltavsky, Kristof T Schutt, Alexandre Tkatchenko, and Klaus-Robert Muller. Machine learning force fields. Chemical Reviews, 121(16):10142–10186, 2021. [21] Beom Seok Kang, Vignesh C Bhethanabotla, Amin Tavakoli, Maurice D Hanisch, William A Goddard I, and Anima Anandkumar. Orbitall: a unified quantum me- chanical representation deep learning framework for all molecular systems. arXiv preprint arXiv:2507.03853, 2025. 16 [22] Jason Wei, Xuezhi Wang, Dale Schuurmans, Maarten Bosma, Fei Xia, Ed Chi, Quoc V Le, Denny Zhou, et al. Chain-of-thought prompting elicits reasoning in large language models. Advances in Neural Information Processing Systems, 35:24824–24837, 2022. [23] He Li, Zun Wang, Nianlong Zou, Meng Ye, Runzhang Xu, Xiaoxun Gong, Wenhui Duan, and Yong Xu. Deep-learning density functional theory hamiltonian for efficient ab initio electronic-structure calculation. Nature Computational Science, 2(6):367–377, 2022. [24] Xiaoxun Gong, He Li, Nianlong Zou, Runzhang Xu, Wenhui Duan, and Yong Xu. General framework for e(3)-equivariant neural network representation of density functional theory hamiltonian. Nature Communications, 14(1):2848, 2023. [25] Manasa Kaniselvan, Benjamin Kurt Miller, Meng Gao, Juno Nam, and Daniel S Levine. Learning from the electronic structure of molecules across the periodic table. arXiv preprint arXiv:2510.00224, 2025. [26] Feitong Song and Ji Feng. Neural network self-consistent fields for density functional theory. npj Computational Materials, 2026. [27] Paolo Giannozzi, Stefano Baroni, Nicola Bonini, Matteo Calandra, Roberto Car, Carlo Cavazzoni, Davide Ceresoli, Guido L Chiarotti, Matteo Cococcioni, Ismaila Dabo, et al. Quantum espresso: a modular and open-source software project for quantum simulations of materials. Journal of Physics: Condensed Matter, 21(39):395502, 2009. [28] Paolo Giannozzi, Oliviero Andreussi, Thomas Brumme, Oana Bunau, Marco Buon- giorno Nardelli, Matteo Calandra, Roberto Car, Carlo Cavazzoni, Davide Ceresoli, Matteo Cococcioni, et al. Advanced capabilities for materials modelling with quantum espresso. Journal of Physics: Condensed Matter, 29(46):465901, 2017. [29] Sebastiaan P Huber, Michail Minotakis, Marnik Bercx, Timo Reents, Kristjan Eimre, Nataliya Paulish, Nicolas H ̈ormann, Martin Uhrin, Nicola Marzari, and Giovanni Pizzi. Mc3d: The materials cloud computational database of experimentally known stoichiometric inorganics. Digital Discovery, 5(3):1114–1131, 2026. [30] Raghunathan Ramakrishnan, Pavlo O Dral, Matthias Rupp, and O Anatole Von Lilienfeld. Quantum chemistry structures and properties of 134 kilo molecules. Scientific Data, 1(1): 1–7, 2014. [31] Clemens Isert, Kenneth Atz, Jos ́e Jim ́enez-Luna, and Gisbert Schneider. Qmugs, quantum mechanical properties of drug-like molecules. Scientific Data, 9(1):273, 2022. [32] Mike C Payne, Michael P Teter, Douglas C Allan, TA Arias, and ad JD Joannopoulos. It- erative minimization techniques for ab initio total-energy calculations: molecular dynamics and conjugate gradients. Reviews of Modern Physics, 64(4):1045, 1992. [33] Donald R Hamann. Optimized norm-conserving vanderbilt pseudopotentials. Physical Review B, 88(8):085117, 2013. [34] Mel Levy. Universal variational functionals of electron densities, first-order density matri- ces, and natural spin-orbitals and solution of the v-representability problem. Proceedings of the National Academy of Sciences, 76(12):6062–6065, 1979. [35] Elliott H. Lieb. Density functionals for coulomb systems. International Journal of Quan- tum Chemistry, 24(3):243–277, 1983. doi: https://doi.org/10.1002/qua.560240302. URL https://onlinelibrary.wiley.com/doi/abs/10.1002/qua.560240302. 17 [36] John P Perdew, Kieron Burke, and Matthias Ernzerhof. Generalized gradient approxima- tion made simple. Physical Review Letters, 77(18):3865, 1996. [37] John P Perdew, Adrienn Ruzsinszky, G ́abor I Csonka, Oleg A Vydrov, Gustavo E Scuseria, Lucian A Constantin, Xiaolan Zhou, and Kieron Burke. Restoring the density-gradient expansion for exchange in solids and surfaces. Physical Review Letters, 100(13):136406, 2008. [38] Nikola Kovachki, Zongyi Li, Burigede Liu, Kamyar Azizzadenesheli, Kaushik Bhattacharya, Andrew Stuart, and Anima Anandkumar. Neural operator: Learning maps between func- tion spaces with applications to pdes. Journal of Machine Learning Research, 24(89):1–97, 2023. [39] Kamyar Azizzadenesheli, Nikola Kovachki, Zongyi Li, Miguel Liu-Schiaffini, Jean Kossaifi, and Anima Anandkumar. Neural operators for accelerating scientific simulations and de- sign. Nature Reviews Physics, 6(5):320–328, 2024. [40] Xuecheng Shao, Kaili Jiang, Wenhui Mi, Alessandro Genova, and Michele Pavanello. Dftpy: An efficient and object-oriented platform for orbital-free dft simulations. Wiley Interdisci- plinary Reviews: Computational Molecular Science, 11(1):e1482, 2021. [41] Qiming Sun, Timothy C Berkelbach, Nick S Blunt, George H Booth, Sheng Guo, Zhendong Li, Junzi Liu, James D McClain, Elvira R Sayfutyarova, Sandeep Sharma, et al. Pyscf: the python-based simulations of chemistry framework. Wiley Interdisciplinary Reviews: Computational Molecular Science, 8(1):e1340, 2018. [42] Seungsoo Nam, Ryan J McCarty, Hansol Park, and Eunji Sim. Ks-pies: Kohn–sham inversion toolkit. The Journal of Chemical Physics, 154(12), 2021. [43] Qin Wu and Weitao Yang. A direct optimization method for calculating density function- als and exchange–correlation potentials from electron densities. The Journal of Chemical Physics, 118(6):2498–2509, 2003. [44] Walter Kohn. v-representability and density functional theory. Physical Review Letters, 51 (17):1596, 1983. [45] Egor Trushin, Jannis Erhard, and Andreas G ̈orling. Violations of the v-representability condition underlying kohn-sham density-functional theory. Physical Review A, 110(2): L020802, 2024. [46] Andrea Grisafi, Alberto Fabrizio, Benjamin Meyer, David M Wilkins, Clemence Cormin- boeuf, and Michele Ceriotti. Transferable machine-learning model of the electron density. ACS Central Science, 5(1):57–64, 2019. [47] Alan M Lewis, Andrea Grisafi, Michele Ceriotti, and Mariana Rossi. Learning electron densities in the condensed phase. Journal of Chemical Theory and Computation, 17(11): 7203, 2021. [48] Fumihiro Imoto, Masatoshi Imada, and Atsushi Oshiyama. Order-n orbital-free density- functional calculations with machine learning of functional derivatives for semiconductors and metals. Physical Review Research, 3(3):033198, 2021. [49] Liang Sun and Mohan Chen. Machine learning based nonlocal kinetic energy density functional for simple metals and alloys. Physical Review B, 109(11):115135, 2024. 18 [50] Zhaoxuan Wu and W. A. Curtin.Mechanism and energetics of〈c + a〉 dis- location cross-slip in hcp metals.Proceedings of the National Academy of Sci- ences, 113(40):11137–11142, October 2016.doi:10.1073/pnas.1603966113.URL https://w.pnas.org/doi/10.1073/pnas.1603966113. [51] Zhaoxuan Wu, Rasool Ahmad, Binglun Yin, Stefanie Sandl ̈obes, and W. A. Curtin. Mechanistic origin and prediction of enhanced ductility in magnesium alloys.Sci- ence, 359(6374):447–452, January 2018.doi:10.1126/science.aap8716.URL https://w.science.org/doi/10.1126/science.aap8716. [52] Sambit Das and Vikram Gavini.Intrinsic ductility enhancement in Mg al- loys elucidated via large-scale ab-initio calculations,January 2026.URL http://arxiv.org/abs/2601.12202. arXiv:2601.12202 [cond-mat.mtrl-sci]. [53] Phani Motamarri, Sambit Das, Shiva Rudraraju, Krishnendu Ghosh, Denis Davydov, and Vikram Gavini. Dft-fe–a massively parallel adaptive finite-element code for large-scale density functional theory calculations. Computer Physics Communications, 246:106853, 2020. 19