Accurate and scalable exchange-correlation with deep learning

ArXiv: 2506.14665

🎯 Pitch

A deep learning functional matches the accuracy of expensive hybrid DFT while costing as little as semi-local methods—breaking a decades-old trade-off. Skala learns non-local interactions directly from data using only cheap meta-GGA inputs, slashing GMTKN55 error to 2.8 kcal/mol without computing exact exchange. Its performance also improves predictably as training data expands, opening a path toward systematically more predictive DFT without increasing computational complexity.

---

1. Executive Summary

This paper introduces Skala, a deep learning–based exchange-correlation (XC) functional for density functional theory that departs from the traditional Jacob's ladder paradigm by learning non-local electronic representations directly from data—rather than relying on increasingly costly hand-engineered features—while accepting only cheap semi-local (meta-GGA) inputs. Trained on ~400,000 high-accuracy CCSD(T)/CBS energy differences from the Microsoft Research Accurate Chemistry Collection and evaluated on the GMTKN55 benchmark, Skala achieves a WTMAD-2 error of 2.8 kcal/mol, surpassing state-of-the-art range-separated hybrid functionals (e.g., ωB97M-V at 3.23 kcal/mol) at the computational cost characteristic of semi-local DFT. The architecture's non-local interaction module—which passes information between integration grid points via coarse points at atomic centers—reduces average error by ~50% compared to purely local models while preserving asymptotic O(N³) scaling, establishing that learned non-locality can substitute for exact Hartree–Fock exchange without incurring its O(N⁴) cost. Systematic data ablations reveal that Skala's accuracy improves monotonically as training data expands in both size and chemical diversity, and that physical constraints—such as the positivity of the kinetic correlation energy—emerge spontaneously from data only when the training set is sufficiently broad, establishing a framework for continuously improvable DFT whose accuracy scales with data rather than algorithmic complexity.

2. Context and Motivation

The Core Problem: The Unresolved Accuracy–Cost Trade-off in Density Functional Theory

The central problem this paper tackles is the fundamental limitation that has constrained density functional theory for decades: all practical exchange-correlation functionals sit on a rigid accuracy–cost efficiency curve, and improving one dimension inevitably sacrifices the other. This trade-off is not incidental—it is baked into the very framework, known as Jacob's ladder, that the field has used to develop XC approximations since the 1990s.

To understand why this is so consequential, we need to appreciate DFT's role. DFT is the workhorse of computational chemistry and materials science because it offers a tractable way to compute the electronic energy of molecules and materials. The exact reformulation by Kohn and Sham (1965) decomposes the total energy into several well-understood components (kinetic energy, electron-nuclear attraction, classical electrostatic repulsion) plus one remaining unknown term—the exchange-correlation (XC) functional—that accounts for the quantum many-body effects: exchange (from the Pauli exclusion principle) and correlation (electrons adjusting their motion to avoid each other). The XC functional is formally a map from the electron density ρ(r), a function of three spatial coordinates, to a scalar energy. It is universal, meaning it has the same functional form for all molecules and materials regardless of composition.

If we knew the exact XC functional, DFT would give exact answers at a computational cost that scales as O(N³) with system size—cheap enough for molecules with hundreds of atoms. But we don't know it, and the central theoretical challenge of DFT since its inception has been to approximate it well.

Jacob's Ladder: The Accurate–Efficient Trade-off Embedded in Practice

The prevailing paradigm for building XC approximations is Jacob's ladder, formalized by Perdew and Schmidt in 2001. The ladder has five rungs, each adding a new type of ingredient to construct the functional:

• Rung 1 (LDA): Uses only the local electron density ρ(r). The simplest and cheapest approximation, but typically too inaccurate for chemical applications (errors ~30–40 kcal/mol for atomization energies).

• Rung 2 (GGA): Adds the gradient of the density, ∇ρ(r). This captures some non-local effects—how the density changes in space—while remaining computationally efficient. Functionals like revPBE and B97M-V sit here. Typical errors fall to ~5–10 kcal/mol for main-group thermochemistry.

• Rung 3 (meta-GGA): Further adds the Kohn-Sham kinetic energy density τ(r). This provides additional semi-local information about the electronic structure without fundamentally changing the scaling. Functionals like r2SCAN reside at this level.

• Rung 4 (hybrid): Introduces a fraction of exact Hartree–Fock exchange, computed from the Kohn-Sham orbitals rather than the density alone. This is the most significant accuracy gain, since exact exchange inherently captures non-local quantum effects that semi-local ingredients miss. Examples include B3LYP (the most cited functional in history), M06-2X, and range-separated hybrids like ωB97M-V. The cost, however, jumps: evaluating exact exchange scales as O(N⁴) in standard implementations, making hybrid functionals substantially more expensive than meta-GGAs for large systems.

• Rung 5 (double-hybrid): Adds perturbative second-order correlation (MP2) on top of hybrid DFT. This pushes accuracy further but at O(N⁵) cost, restricting double-hybrids to relatively small molecules.

The crucial point is that all widely used functionals live somewhere on this ladder, and moving up the ladder improves accuracy but increases computational cost. This trade-off is the organizing principle of functional development: if you need better accuracy, you climb a rung and pay the computational price. The paper's Figure 1(c) makes this visible—the WTMAD-2 error on the GMTKN55 benchmark drops from ~8.4 kcal/mol (revPBE, GGA) to ~5.6 (B97M-V, meta-GGA) to ~3.2 (ωB97M-V, range-separated hybrid), but each step comes with asymptotic cost increases.

Why This Trade-off Matters: The "Chemical Accuracy" Gap

The practical consequence is that no existing functional simultaneously achieves chemical accuracy (~1 kcal/mol) across the breadth of main-group chemistry while remaining computationally affordable for routine use on medium-to-large molecules. Hybrid functionals, which come closest, still leave errors of 3–4 kcal/mol averaged over broad benchmark sets, and individual errors can be much larger for specific reaction types. This gap has real-world consequences:

• Drug discovery: Screening candidate molecules computationally before synthesis requires reliable predictions of binding energies. Errors of a few kcal/mol correspond to orders of magnitude in binding affinity, meaning that computational screens pass too many false positives to experimental validation, wasting laboratory resources.

• Catalysis design: Predicting whether a catalyst will lower a reaction barrier sufficiently requires barrier height accuracies better than ~1–2 kcal/mol. Current functionals have typical barrier errors of 3–7 kcal/mol, making reliable in silico catalyst design impossible.

• Materials discovery: Identifying materials with specific properties (carbon capture, battery electrodes) depends on accurate relative energies between competing crystal polymorphs or between reactants and products. Errors at the hybrid level can reverse the predicted stability ordering.

• Force field training: Many machine-learned interatomic potentials are trained on DFT data. Errors in the underlying DFT functional propagate directly into these downstream models, limiting their predictive power regardless of how sophisticated the force field architecture is.

Beyond practical applications, there is a deeper theoretical frustration: the exact XC functional is known to be non-local in the density. This is not speculation—it is a proven property of the exact functional (Lieb, 1983). Semi-local approximations (GGA, meta-GGA) are therefore fundamentally limited in what they can capture; they rely on error cancellation to partially compensate for missing non-locality. Hybrid functionals recover some non-locality via exact exchange, but in a rigid, non-adaptable form: the fraction of exact exchange is a fixed global parameter that cannot vary across different chemical environments. This means that even the best hybrids are over-stabilizing in some bonding situations and under-stabilizing in others—a problem known for decades as the "fractional electron problem" or delocalization error.

Prior Attempts: Machine Learning and the Feature Engineering Trap

The idea of using machine learning to improve XC functionals is not new. The earliest efforts date to the 1990s (Tozer et al., 1996), and there has been a steady stream of work since. However, as the paper notes in its introduction, these efforts have not led to a functional that sees widespread adoption, and the fundamental accuracy–cost trade-off has remained intact. There are specific, interlinked reasons for this:

Scarcity of high-accuracy training data. The XC functional is not a standalone quantity—it is a component of a self-consistent computational procedure. To train an XC functional, you need reference data that comes from methods more accurate than DFT—specifically, wavefunction-based methods like coupled-cluster theory with singles, doubles, and perturbative triples, CCSD(T), extrapolated to the complete basis set limit (CCSD(T)/CBS). CCSD(T)/CBS is considered the "gold standard" of quantum chemistry, achieving ~1 kcal/mol accuracy relative to experiment for systems where a single electronic configuration dominates (single-reference systems). But CCSD(T) calculations are extremely expensive: their cost scales as O(N⁷), limiting them to molecules with roughly 10–20 non-hydrogen atoms. Generating enough high-accuracy data to train a general-purpose deep learning model across diverse chemical space has historically been infeasible—it requires both massive computational resources and specialized expertise to run these calculations at scale. Prior ML efforts were therefore confined to small datasets (hundreds to low thousands of data points), which is orders of magnitude insufficient for learning complex non-local mappings from scratch.

The feature engineering paradigm. Confined to this low-data regime, the vast majority of prior ML work adopted the classic machine learning strategy that dominated before deep learning: hand-design informative features and train a relatively simple model on them. Specifically, researchers fed standard Jacob's ladder ingredients (density, gradient, kinetic energy density) or newly designed semi-local descriptors into models such as neural networks, Gaussian processes, or kernel methods, but the models themselves operated locally—they mapped features at a single grid point to the energy contribution at that point, without learning interactions across grid points. Examples include:

• DM21 (Kirkpatrick et al., 2021): Trained a neural network as a local functional of semi-local features, incorporating exact constraints and training on fractional electron systems. It achieved a WTMAD-2 of ~3.97 kcal/mol—competitive with hybrids—but required adding exact-exchange-like features to the input, effectively building non-locality in by hand rather than learning it.

• CIDER (Bystrom & Kozinsky, 2022, 2023): Designed an expressive set of non-local feature descriptors ("convolved" density-derived quantities) and fed them into a linear model, showing that non-local features improve accuracy, but still requiring expert specification of which non-local functions matter.

• Various efforts to learn from exact densities or potentials (Dick & Fernandez-Serra, 2020; Kanungo et al., 2024) showed promise on small systems but did not scale to broad chemical space.

This feature engineering approach mirrors the state of computer vision before AlexNet (2012)—models were limited by the quality of hand-designed features, and the key transition to deep learning happened only when large datasets (ImageNet) became available to power learned feature hierarchies. The analogy is explicit in the paper (Section 2.3), which states that their work "mirrors the shift from classic machine learning with hand-designed features to deep representation learning."

Non-local architectures remain nascent. A handful of prior studies explored learning non-local interactions on the DFT grid rather than engineering non-local features. Schmidt et al. (2019) learned a non-local XC functional for 1D model systems; Margraf & Reuter (2021) used convolutional neural networks on 3D grids for small molecules; Li et al. (2021) and Kalita et al. (2022) used the Kohn-Sham equations as a regularizer; Gao et al. (2024) introduced equivariant message-passing between grid points. These proof-of-concept studies demonstrated that learned non-locality is viable in principle but were limited to model systems, narrowly defined problem classes, or very small molecules. None scaled to the diversity and size of chemical space needed for a general-purpose functional.

The self-consistency challenge. The XC functional is not evaluated in isolation but inside a self-consistent field loop: one guesses an initial density, computes the XC potential from that density, solves the Kohn-Sham equations to get a new density, and iterates until consistency. This creates a chicken-and-egg problem for training. If you train the model on fixed densities from an existing functional (say B3LYP), the model learns to be accurate on those densities but may produce poor self-consistent densities when deployed, leading to density-driven errors that degrade total energies. If you train self-consistently from the start—running the SCF loop during training and backpropagating through it—the cost is prohibitive at scale, and the training signal is poor early on when the model produces unphysical densities that cause SCF convergence failures. Prior work has explored various compromises (e.g., training on exact densities where available, using the Kohn-Sham equations as a regularizer), but no solution has been demonstrated to scale to the dataset sizes needed for broad chemical generality.

How Skala Positions Itself

The paper positions Skala as a departure from all three legs of the prior paradigm simultaneously:

1. Data scale: By undertaking an unprecedented data generation effort—~400,000 CCSD(T)/CBS energy differences spanning atomization energies, reaction barriers, non-covalent interactions, proton affinities, ionization potentials, conformers, distorted geometries, and atomic properties—the authors escape the low-data regime that forced prior work into feature engineering. They explicitly position the MSR-ACC dataset as analogous to ImageNet for computer vision: a large-scale, diverse, high-quality dataset that enables the transition from hand-designed features to learned representations.

2. Architecture: Rather than engineering non-local features, Skala adopts a neural operator approach: the model takes only standard, cheap meta-GGA inputs (density, density gradient, kinetic energy density) and learns to represent non-local interactions through its architecture. The key innovation is an information-passing mechanism between integration grid points via coarse points at atomic centers—analogous to message-passing in graph neural networks but adapted for the continuous, irregular point clouds that DFT integration grids represent. The architecture is designed to scale: the number of coarse points equals the number of atoms, the communication is local (each grid point talks only to its own atom's coarse point), and the asymptotic computational cost remains O(N³). The paper explicitly contrasts this with traditional hybrids: rather than computing exact exchange (O(N⁴)) to get non-locality, Skala learns whatever non-local effects the data demand, at lower cost.

3. Training strategy: The two-stage protocol—pre-training on fixed B3LYP densities at scale, followed by short SCF fine-tuning—decouples the regression problem from the self-consistency problem. By starting with approximate densities that are available at scale and then closing the self-consistency gap with a small number of on-the-fly SCF steps, the authors make large-scale training computationally feasible without sacrificing self-consistent accuracy.

The paper's central claim is that this combination—large data, learned non-locality, and scalable training—breaks the Jacob's ladder trade-off: Skala achieves hybrid-level or better accuracy across main-group chemistry at semi-local (meta-GGA) cost. More fundamentally, the paper positions Skala not as a static functional but as a framework for continuously improvable DFT—its accuracy scales with data, not algorithmic complexity. This reframes the challenge of "finding the divine functional" from a problem of designing better analytic forms to a problem of generating better training data, paralleling the paradigm shift that deep learning brought to computer vision, natural language processing, and many other fields.

Why Prior Hybrid Functionals, Despite Their Success, Are Not the End State

To motivate why learned non-locality is preferable to simply using exact exchange (hybrid functionals), the paper makes several points:

• Hybids contain a fixed fraction of exact exchange (typically 15–54%, or range-separated). This fraction is globally constant, meaning it cannot adapt to different chemical environments. Certain problems are over-stabilized (e.g., electron delocalization errors in charge-transfer complexes), while others remain under-stabilized.

• Computing exact exchange on a grid—as local hybrids do to make the fraction position-dependent—introduces a much larger prefactor, making local hybrid implementations often more expensive than global hybrids despite having the same formal scaling. The paper explicitly draws this analogy in Section 5 to explain why the architecture's prefactor matters as much as asymptotic scaling.

• Hybrids still require additional corrections (e.g., D3 dispersion) to handle long-range van der Waals interactions, which exact exchange does not capture. Skala currently also uses D3, but the architecture's non-local module could in principle learn dispersion if trained on appropriate data—a direction marked as future work.

• The accuracy of hybrid functionals has plateaued: the best range-separated hybrids (ωB97M-V) achieve WTMAD-2 of ~3.2 kcal/mol, and more recent efforts (COACH, 2026) report only incremental improvements. Breaking below ~2.5 kcal/mol across broad chemical space appears to require new mechanisms, not just refinements of the hybrid paradigm.

Reconciling Conflicting Findings in the Literature

One implicit motivation the paper addresses is the mixed track record of ML functionals. Prior work reached apparently contradictory conclusions: some found that ML functionals could approach hybrid accuracy (Kirkpatrick et al., 2021), while others found that learned functionals overfit, generalize poorly to unseen chemistry, or fail in self-consistent calculations. The paper's position is that these failures largely stem from insufficient data diversity, not from fundamental limitations of the approach. By training on ~400,000 data points spanning diverse reaction types, elements, and chemical environments—and by demonstrating systematic improvement as data diversity increases (Figure 2, right)—the paper makes the case that the primary bottleneck has been data, not methodology.

Summary of the Motivational Landscape

The remainder of the paper is devoted to substantiating these claims: showing that the architecture does learn non-local effects (Section 3.1), that accuracy improves monotonically with data (Section 3.2), that physical constraints emerge from data (Section 3.3), that the functional performs well self-consistently (Section 4), and that the computational cost indeed remains in the semi-local regime (Section 5).

3. Technical Approach

3.1 Reader Orientation

Skala is a deep neural network that takes the place of the exchange-correlation functional inside a standard density functional theory calculation — it receives the electron density of a molecule, represented on a grid of points in 3D space, and outputs the exchange-correlation energy, which is the missing piece needed to compute the total electronic energy. The system solves the problem of designing an XC functional that is simultaneously accurate across diverse chemistry and computationally efficient by learning non-local quantum mechanical effects directly from data rather than hand-engineering them from expensive wavefunction-based ingredients, using only cheap semi-local features as input and a scalable message-passing architecture to approximate the multi-body interactions between electrons at different locations in space.

3.2 Big-Picture Architecture (Diagram in Words)

The Skala system has five major components that work in sequence:

1. Input Feature Extractor — reads the electron density and its derived quantities (density itself, gradient norm, kinetic energy density, all spin-resolved) on each point of a large numerical integration grid (~10⁴–10⁶ points), converts them to log-space, applies a small multilayer perceptron to learn a spin-symmetrized hidden representation per grid point.

2. Non-Local Interaction Module — the core innovation: a stack of three layers that enable grid points to communicate indirectly through coarse points placed at atomic centers. Each layer (i) projects local features onto radial and spherical harmonic basis functions, (ii) aggregates ("downsamples") these projections across each atom's grid onto its coarse point, (iii) applies an equivariant linear transformation and a symmetric tensor product to capture higher-body-order correlations on the coarse points, and (iv) projects the processed coarse features back ("upsamples") to the fine grid via the same basis functions, combining with the original features through a skip connection.

3. Output Model — a purely local MLP that takes the non-locally enriched features at each grid point (after three layers of message passing) and maps them to a scalar enhancement factor, which is passed through a scaled sigmoid to enforce bounds.

4. Integration — multiplies the enhancement factor at each grid point by the LDA exchange energy density and the integration weight, then sums over all grid points to produce the total exchange-correlation energy $E_{\text{xc}}$.

5. Training Pipeline — a two-stage protocol: pre-training on ~400,000 fixed B3LYP densities (bypassing the self-consistent field loop) followed by short fine-tuning using the model's own SCF densities, with a reaction-energy regression loss.

Information flows as follows: a molecule's density matrix (from B3LYP during pre-training, or from Skala's own SCF during fine-tuning) generates density-derived features on the integration grid → the features are log-transformed and spin-symmetrized → the non-local module enriches each grid point's representation with information from other grid points on the same atom → the output MLP produces an enhancement factor per grid point → numerical integration yields $E_{\text{xc}}$ → this is combined with the other Kohn-Sham energy components to produce the total energy → the total energies of reactants and products are subtracted to produce a reaction energy, which is compared against CCSD(T)/CBS reference values via a weighted MSE loss.

3.3 Roadmap for the Deep Dive

• First, the exchange-correlation energy expression (Equation 1) that defines the functional form Skala learns — this anchors all subsequent architectural decisions in the physics of DFT.
• Second, the input features on the integration grid — what semi-local quantities are extracted from the density matrix, how they are preprocessed, and why log-transformation matters.
• Third, the spin-symmetrization and initial local feature extraction — the first learned component of the architecture, which produces a rotation-invariant, spin-channel-invariant hidden representation.
• Fourth, the non-local interaction module in full detail — the downsampling-upsampling mechanism through coarse points, the radial and spherical harmonic basis, the equivariant linear mixing, the symmetric contraction, and the skip connections, including the theoretical justification for why this approximates multi-body interactions.
• Fifth, the output model and the final energy integration — how the per-point enhancement factor becomes $E_{\text{xc}}$ and how exact constraints (Lieb-Oxford bound) are enforced.
• Sixth, the two-stage training protocol — why pre-training on fixed B3LYP densities followed by SCF fine-tuning is necessary, the loss function, the adaptive dataset sampling mechanism, and the optimizer configuration.

3.4 Detailed, Sentence-Based Technical Breakdown

This is primarily an empirical architecture and training paper whose core idea is that non-local exchange-correlation effects can be learned from semi-local inputs using a scalable neural-operator-inspired message-passing architecture, provided the training data is sufficiently large and chemically diverse.

---

The Exchange-Correlation Energy Form

The fundamental quantity Skala must learn is the exchange-correlation energy $E_{\text{xc}}[\rho]$, a functional that maps the electron density $\rho(\mathbf{r})$ (a function of three spatial coordinates) to a scalar energy. The paper models this as an integral of an energy density weighted by an enhancement factor that multiplies the Local Density Approximation (LDA) exchange energy density:

$E_{\text{xc}}^\theta[\rho] = -\frac{3}{4}\left(\frac{6}{\pi}\right)^{1/3} \int \left[\rho^{(\uparrow)}(\mathbf{r})^{4/3} + \rho^{(\downarrow)}(\mathbf{r})^{4/3}\right] f_\theta[\mathbf{x}[\rho]](\mathbf{r})\,d\mathbf{r}$

where $\rho^{(\uparrow)}(\mathbf{r})$ and $\rho^{(\downarrow)}(\mathbf{r})$ are the electron densities of the two spin channels (spin-up and spin-down electrons) at position $\mathbf{r}$, $\mathbf{x}[\rho](\mathbf{r})$ is a vector of semi-local features derived from the density at that point (to be specified below), and $f_\theta$ is the enhancement factor — a scalar-valued functional parameterized by a deep neural network with parameters $\theta$.

What it computes: for each point $\mathbf{r}$ in space, this expression takes the local spin densities $\rho^{(\uparrow)}$ and $\rho^{(\downarrow)}$, raises them to the power 4/3 (the LDA exchange form), multiplies by the enhancement factor $f_\theta$ (which is a learned multiplier that can range from 0 to 2, as we will see), and integrates over all space weighted by quadrature weights (discretized on the integration grid). The prefactor $-\frac{3}{4}(6/\pi)^{1/3}$ is the LDA exchange constant, chosen so that when $f_\theta = 1$ the expression reduces to LDA exchange — the simplest possible XC approximation that treats the density as a uniform electron gas. The enhancement factor is therefore the correction to LDA exchange that captures all the complex physics: gradient corrections (GGA), kinetic energy corrections (meta-GGA), and non-local effects that would traditionally require exact exchange (hybrid functionals).

Why this form: there are several interlocking design decisions embedded in this equation:

• Separation from LDA: by factoring out the LDA exchange energy density $-\frac{3}{4}(6/\pi)^{1/3}(\rho^{(\uparrow)4/3} + \rho^{(\downarrow)4/3})$, the enhancement factor $f_\theta$ becomes a dimensionless scalar close to 1. This makes the learning problem well-conditioned — the neural network predicts a small multiplicative correction rather than the full energy magnitude, which makes initialization and training more stable. If the network outputs 1, the functional defaults to LDA, a reasonable baseline.

• High-density uniform coordinate scaling: a fundamental exact constraint on the XC functional is that under uniform scaling of the density $\rho_\gamma(\mathbf{r}) = \gamma^3\rho(\gamma\mathbf{r})$ (compressing the density by factor $\gamma$), the exchange energy scales as $E_x[\rho_\gamma] = \gamma E_x[\rho]$ in the high-density limit $\gamma \to \infty$. The LDA exchange form $\rho^{4/3}$ correctly captures this scaling, so by factoring the LDA exchange prefactor, Skala inherits this constraint automatically regardless of what the network learns. Building in the exact high-density scaling is energetically critical — it governs the core electron behavior in atoms and molecules, where densities are highest.

• Size-consistency: the energy is written as an integral over space. If two molecules are infinitely separated, the density is the sum of the two individual densities, and the integral separates into two independent integrals because the integration grid for each molecule is independent. This means the XC energy is additive for non-interacting fragments — a property called size-consistency that is essential for predicting reaction energies correctly (you don't want the functional to artificially stabilize or destabilize the separated reactants relative to the bound product).

• Lieb-Oxford lower bound: the exact XC functional is known to satisfy a lower bound — the exchange-correlation energy cannot be arbitrarily negative relative to the density; there is a rigorous lower bound proportional to $\int \rho^{4/3}d\mathbf{r}$. By constraining $f_\theta \in [0, 2]$ (as we will see in the output model section), Skala enforces this bound, preventing unphysically large negative XC energies that would artificially stabilize systems. The upper bound of 2 is chosen following Kirkpatrick et al. (2021), who showed that most physically relevant enhancement factors fall within this range.

• Spin decomposition: treating spin-up and spin-down densities separately in the prefactor (rather than using the total density) is physically correct because exchange only occurs between electrons of the same spin. The spin channels are additive: the LDA exchange for spin-up electrons depends on $\rho^{(\uparrow)4/3}$, not on $(\rho^{(\uparrow)} + \rho^{(\downarrow)})^{4/3}$. Getting this right matters for open-shell systems (radicals, transition metals) where spin polarization is significant.

The notation $f_\theta[\mathbf{x}[\rho]](\mathbf{r})$ (with square brackets) indicates that the enhancement factor is a neural functional: it is not merely a function of local features $\mathbf{x}(\mathbf{r})$ at point $\mathbf{r}$, but can depend on the features at other points as well, through the non-local message-passing architecture. This is the crucial distinction from prior ML functionals: $f_\theta$ is not $f_\theta(\mathbf{x}(\mathbf{r}))$ (local) but $f_\theta[\mathbf{x}[\rho]](\mathbf{r})$ (non-local functional). The architecture that enables this non-local dependence while maintaining scalability is the core technical contribution.

In practice, the integral is discretized on the numerical integration grid with $G$ points at positions $\mathbf{r}_i$ and associated quadrature weights $w_i$:

$E_{\text{xc}}^\theta[\rho] \approx -\frac{3}{4}\left(\frac{6}{\pi}\right)^{1/3} \sum_{i=1}^{G} h_{\text{enh},i}\left[\rho^{(\uparrow)}(\mathbf{r}_i)^{4/3} + \rho^{(\downarrow)}(\mathbf{r}_i)^{4/3}\right] w_i$

where $h_{\text{enh},i} = f_\theta[\mathbf{x}[\rho]](\mathbf{r}_i)$ is the enhancement factor predicted by the network at grid point $i$. The sum runs over all grid points of all atoms in the molecule. The grid weights $w_i$ are determined by the quadrature scheme (Treutler-Ahlrichs, Mura-Knowles, Gauss-Chebyshev, or Delley radial grids combined with Lebedev angular grids), and the network is trained with grid augmentation — four different radial schemes at different grid levels — to regularize against overfitting to any specific quadrature.

---

Input Features: Semi-Local Density Descriptors on the Integration Grid

The input to Skala at each grid point is a 7-dimensional vector of meta-GGA features — the standard semi-local ingredients available at the third rung of Jacob's ladder:

$\mathbf{x}[\rho](\mathbf{r}_i) = \left(\rho^{(\uparrow)}(\mathbf{r}_i), \rho^{(\downarrow)}(\mathbf{r}_i), \|\nabla\rho^{(\uparrow)}(\mathbf{r}_i)\|_2, \|\nabla\rho^{(\downarrow)}(\mathbf{r}_i)\|_2, \tau^{(\uparrow)}(\mathbf{r}_i), \tau^{(\downarrow)}(\mathbf{r}_i), \|\nabla\rho^{(\uparrow)}(\mathbf{r}_i) + \nabla\rho^{(\downarrow)}(\mathbf{r}_i)\|_2\right)$

where the components are:

• $\rho^{(\uparrow)}(\mathbf{r}_i)$ and $\rho^{(\downarrow)}(\mathbf{r}_i)$: the electron density of spin-up and spin-down electrons at grid point $i$. These are the most fundamental quantities — they tell you how many electrons of each spin are present at each point in space. When the system is closed-shell (same number of up and down electrons), these are equal.

• $\|\nabla\rho^{(\uparrow)}(\mathbf{r}_i)\|_2$ and $\|\nabla\rho^{(\downarrow)}(\mathbf{r}_i)\|_2$: the Euclidean norm of the gradient of the spin densities. The gradient norm measures how rapidly the density changes in space — it is large near nuclei (where the density drops sharply) and in bonding regions, and small in the slowly-varying tails of the electron distribution. This is the key ingredient that distinguishes GGA from LDA: density-gradient information captures the inhomogeneity of the electron gas, which is essential for describing molecular bonding.

• $\tau^{(\uparrow)}(\mathbf{r}_i)$ and $\tau^{(\downarrow)}(\mathbf{r}_i)$: the Kohn-Sham kinetic energy density for each spin channel, defined as $\tau^{(\sigma)}(\mathbf{r}) = \frac{1}{2}\sum_{j}^{\text{occ}} |\nabla\phi_j^{(\sigma)}(\mathbf{r})|^2$, where $\phi_j^{(\sigma)}$ are the occupied Kohn-Sham orbitals of spin $\sigma$. The kinetic energy density provides information beyond the density and its gradient — it can distinguish between bonding regions (where multiple orbitals contribute), lone-pair regions (where one localized orbital dominates), and empty space (where all orbitals have small gradients). This is the additional ingredient that defines the meta-GGA rung.

• $\|\nabla\rho^{(\uparrow)}(\mathbf{r}_i) + \nabla\rho^{(\downarrow)}(\mathbf{r}_i)\|_2$: the norm of the total density gradient (the gradient of the sum, which equals the sum of the gradients). This is included alongside the spin-resolved gradient norms to provide the model with both individual spin-channel information and combined information, facilitating learning of spin-dependent effects.

Why these seven quantities specifically: they represent the maximum amount of semi-local information available at the meta-GGA rung without computing non-local quantities like the exact exchange energy density or the Laplacian of the density (which is a meta-GGA ingredient used in some traditional functionals but not included here — the authors mention in Supplementary Information that they omit the Laplacian to keep the feature set focused). Critically, these features can all be computed from the density matrix at $O(N^3)$ cost — they do not require evaluating exact exchange (O(N⁴)) or MP2 correlation (O(N⁵)). The entire philosophy of Skala is that non-local physics should be learned by the architecture, not hard-coded into more expensive features.

Log transformation and numerical stability: before feeding these features into the neural network, the paper applies a log transform with a small constant offset to avoid numerical issues with near-zero density values:

$\mathbf{x}_i^{(\uparrow,\downarrow)} = \log\left(\left[\rho_i^{(\uparrow)}, \rho_i^{(\downarrow)}, \|\nabla\rho_i^{(\uparrow)}\|_2, \|\nabla\rho_i^{(\downarrow)}\|_2, \tau_i^{(\uparrow)}, \tau_i^{(\downarrow)}, \|\nabla\rho_i^{(\uparrow)} + \nabla\rho_i^{(\downarrow)}\|_2\right] + \epsilon\right)$

where $\epsilon = 10^{-5}$ prevents taking the logarithm of zero (which would occur for $\rho_i^{(\downarrow)}=0$ in fully spin-polarized regions). The log transformation serves two purposes:

• Dynamic range compression: electron densities vary by many orders of magnitude — near a nucleus, the density can be $\sim 10^5$ atomic units, while in the valence region it drops to $\sim 10^{-3}$ and in the far tail to $\sim 10^{-8}$. Without log transformation, the network would need to learn to handle inputs spanning eight orders of magnitude, which is numerically challenging. Log-transformed inputs have a much more manageable range, making the optimization landscape smoother.

• Multiplicative relationships become additive: in electronic structure, many physical relationships are multiplicative in the density (e.g., exchange energy scales as $\rho^{4/3}$, correlation energy has logarithmic density dependence in the high-density limit). By taking logs, the network can learn multiplicative effects via additive weight combinations, which is easier for gradient-based optimization than learning multiplicative functions directly.

The superscript $(\uparrow,\downarrow)$ on $\mathbf{x}_i$ emphasizes the spin ordering in the input vector: the features are arranged with spin-up quantities first, then spin-down. As we will see next, the model processes both spin orderings and averages them to enforce spin-channel permutation invariance.

Grid structure and per-atom partitioning: the $G$ grid points are not an unstructured point cloud — they are organized into atom-centered sub-grids. Each atom $j$ has its own set of grid points $\{\mathbf{r}_k^{(j)}\}$ with associated per-atom weights $\tilde{w}_k^{(j)}$, and a position-dependent partition function $\pi(\mathbf{r}, \mathbf{R}_j; \{\mathbf{R}_{j'}\})$ (Becke partition) that smoothly assigns every point in space to nearby atoms with weights that sum to 1. The total energy is then a sum over atoms of per-atom energy contributions:

$E_{\text{xc}}[\rho] = \sum_j \sum_{i \in \mathcal{G}_j} \text{[energy density at grid point }i\text{ of atom }j\text{]} \cdot \tilde{w}_i^{(j)} \cdot \pi(\mathbf{r}_i^{(j)}, \mathbf{R}_j; \{\mathbf{R}_{j'}\})$

where $\mathcal{G}_j$ denotes the set of grid points belonging to atom $j$. This per-atom decomposition is crucial for the non-local module: each grid point communicates only with the coarse point of its own associated atom, which provides a natural locality constraint that keeps the architecture scalable (linear in the number of atoms).

Grid augmentation during training: to regularize against numerical sensitivity to grid placement (a practical issue with ML functionals that can be "brittle" to grid choice), the authors compute the density features on four distinct radial integration schemes at two grid levels: Treutler-Ahlrichs, Mura-Knowles, Gauss-Chebyshev, and Delley radial grids at PySCF grid level 1 during pre-training and level 3 during fine-tuning. Additionally, at each training step, the space-partitioning scheme (the Becke partition function $\pi$) is randomized among three options: standard Becke partition with Treutler-Ahlrichs radii adjustment, the Stratmann-Scuseria-Frisch (SSF) scheme, or the Laqua-Kussmann-Ochsenfeld (LKO) linear-scaling scheme. This augmentation exposes the model to different spatial weight distributions during training, forcing it to learn a functional that is robust to grid and partition choices rather than overfitting to the idiosyncrasies of any single scheme.

---

Part One: Spin Symmetrization and Initial Local MLP

The first learned component of the Skala architecture processes the log-transformed 7-dimensional input features independently at each grid point to produce a spin-symmetrized hidden representation. This component consists of a small multi-layer perceptron (MLP) applied twice with different spin orderings:

$h_i = \frac{f_{\text{repr}}(\mathbf{x}_i^{(\uparrow,\downarrow)}) + f_{\text{repr}}(\mathbf{x}_i^{(\downarrow,\uparrow)})}{2}$

where $f_{\text{repr}}$ is a 2-layer MLP (two linear layers with Swish activation functions in between), and $\mathbf{x}_i^{(\uparrow,\downarrow)}$ and $\mathbf{x}_i^{(\downarrow,\uparrow)}$ denote the same 7 features but with the spin channels swapped (spin-up features come first in one case, spin-down first in the other).

The explicit form of $f_{\text{repr}}$ as given in the supplementary information is:

$f_{\text{repr}}(\mathbf{x}) = \sigma\left(W_2\,\sigma(W_1\mathbf{x} + \mathbf{b}_1) + \mathbf{b}_2\right)$

where $\sigma$ is the Swish activation function (SiLU): $\sigma(z) = z \cdot \text{sigmoid}(z)$, $W_1$ is a $D_{\text{hid}} \times 7$ weight matrix mapping from 7 input features to the hidden dimension $D_{\text{hid}} = 256$, $\mathbf{b}_1$ is a bias vector of size $D_{\text{hid}}$, $W_2$ is a $D_{\text{hid}} \times D_{\text{hid}}$ weight matrix (same-dimensional mapping), and $\mathbf{b}_2$ is a $D_{\text{hid}}$-dimensional bias. The output $h_i \in \mathbb{R}^{256}$ is a learned hidden representation at grid point $i$.

What it computes: for each grid point independently, this MLP maps the 7 raw semi-local features to a 256-dimensional hidden vector that captures local patterns in the density, gradient, and kinetic energy that are predictive of the exchange-correlation energy. By running the same MLP on both spin orderings and averaging, the representation becomes invariant to which spin channel is called "up" versus "down" — a physical symmetry (spin labels are arbitrary; the physics depends only on the two spin densities, not on which is labelled which).

Why this design:

• Local processing before non-local communication: the initial MLP learns a rich local representation before any information is exchanged between grid points. This is analogous to the initial convolution layers in a CNN: they extract elementary features (edges, textures) before the deeper layers combine them into more complex patterns. Here, the local MLP learns to encode combinations of density, gradient, and kinetic energy that are locally diagnostic of exchange-correlation effects (e.g., recognizing a covalent bond from the density depletion and kinetic energy enhancement in the bonding region).

• Swish activation: the choice of Swish over ReLU is likely motivated by its smoothness (Swish has continuous derivatives everywhere, unlike ReLU which has a discontinuity at zero). Since the XC energy must be a smooth functional of the density (the energy surface should not have kinks as a function of nuclear coordinates), smooth activation functions are preferred to avoid introducing artificial non-smoothness in the energy landscape. Swish also avoids the "dying ReLU" problem (neurons that output zero for all inputs and stop learning), which matters for stable training.

• Spin symmetrization before non-local processing: by making the representation spin-symmetric before it enters the non-local module, the expensive message-passing operations (which dominate the computational cost) are performed once on a single 256-dimensional vector per grid point, rather than twice on spin-specific features. Since the non-local layers involve projecting onto spherical harmonics, tensor products, and contractions, this cuts the non-local computation cost roughly in half. The spin-specific information is not lost — it is encoded by the local MLP into patterns in the 256-dimensional vector that the non-local layers can still read.

• Hidden dimension 256: this is a key hyperparameter that balances expressivity against cost. Each grid point's representation at this stage is a 256-dimensional vector; after non-local processing, it remains 256-dimensional; before being projected to a scalar enhancement factor. For comparison, the non-local branch that processes messages through coarse points operates in a much lower dimension ($D_{\text{nonl}} = 16$, as we will see), making the bulk of the model's capacity concentrated in the local representation learning, with the non-local module providing a lightweight but crucial mechanism for communication.

---

Part Two: Non-Local Interaction Through Coarse Points

This is the core architectural innovation of Skala — the mechanism that allows the functional to learn non-local exchange-correlation effects without computing expensive wavefunction-based ingredients like exact exchange.

Motivation and the Coarse Point Idea

In a standard (non-local) functional such as a hybrid, the non-local exchange interaction between two electrons at positions $\mathbf{r}$ and $\mathbf{r}'$ is approximated by computing the exact Hartree-Fock exchange energy, which involves a double integral over all pairs of grid points $\iint \frac{|\sum_j \phi_j(\mathbf{r})\phi_j^*(\mathbf{r}')|^2}{|\mathbf{r} - \mathbf{r}'|} d\mathbf{r}\,d\mathbf{r}'$. For $G \approx 10^4\text{–}10^6$ grid points, this double sum is $O(G^2)$ — computationally prohibitive. Hybrid functionals avoid this by evaluating the exchange in the basis-set representation (as a sum over occupied orbitals), which scales as $O(N^4)$ with system size — better than $O(G^2)$ but still more expensive than the $O(N^3)$ of semi-local DFT. The question Skala answers is: can we capture the essential non-local physics without the full $O(G^2)$ all-to-all interaction?

The answer is a coarse-grained message-passing scheme: instead of having every grid point talk to every other grid point, the model creates a small set of intermediate "coarse points" — one per atom, placed at the atomic nucleus — that serve as communication hubs. Each grid point sends information to its atom's coarse point; the coarse point aggregates information from all its grid points; after some processing, the coarse point sends information back to its grid points. This reduces the communication cost from $O(G^2)$ to $O(G \cdot N_{\text{atoms}})$ (since each grid point communicates only with one coarse point, and each coarse point aggregates only from its own atom's grid points), which is $O(N^4)$ in worst case but with a tiny prefactor because $N_{\text{atoms}} \ll G$ and the aggregation uses efficient basis-function projections (spherical harmonics and radial Gaussians) whose cost scales with the number of basis functions (16 radial, up to $\ell_{\text{max}}=3$ angular), not with the number of atoms.

The theoretical motivation — that a single downsampling-upsampling pass can approximate any two-body interaction kernel and that stacking multiple layers with symmetric tensor products can capture higher-body-order correlations — is formalized in Supplementary Information Section A.6 and discussed later in this breakdown.

Structure of One Non-Local Layer

Skala applies $L = 3$ non-local layers in sequence (the same layer structure repeated three times with independent learnable parameters). Each layer takes a $D_{\text{hid}} = 256$-dimensional hidden vector at each grid point and outputs a 256-dimensional enriched vector. The layer is composed of five sequential operations:

Step 1: Pre-downsampling projection. Before sending information to the coarse points, the high-dimensional hidden representation is projected to a much lower dimension for computational efficiency:

$h_i^{\text{pre-down}} = \sigma\left(W_{\text{pre-down}} h_i^{(l)} + \mathbf{b}_{\text{pre-down}}\right)$

where $h_i^{(l)} \in \mathbb{R}^{256}$ is the input hidden vector at grid point $i$ (for the first non-local layer, this is the spin-symmetrized output from Part One; for subsequent layers, it is the output of the previous non-local layer), $W_{\text{pre-down}}$ is a $D_{\text{nonl}} \times D_{\text{hid}}$ matrix (with $D_{\text{nonl}} = 16$), and $\sigma$ is the Swish activation. The output is a 16-dimensional vector per grid point — a 16× compression that dramatically reduces the cost of all subsequent spherical-harmonic projections.

Why 16-dimensional: the non-local branch is the computational bottleneck because it involves tensor products with spherical harmonics (which have $2\ell+1$ components per $\ell$, and the basis size grows as $C\cdot\sum_\ell (2\ell+1)$ where $C=16$ radial channels). By compressing to 16 channels, the cost of these projections is kept manageable. The local branch (256-dimensional) still captures most of the expressive power through deep MLPs; the non-local branch provides a correction that enriches the local representation with longer-range context.

Step 2: Downsampling — aggregating grid information onto coarse points. For each atom $j$, a set of coarse features is computed by projecting the 16-dimensional grid features onto a product basis of radial functions and spherical harmonics, then summing over the atom's grid points:

$H_{j\ell c} = \sum_{k \in \mathcal{G}_j} \phi_c(\|\mathbf{r}_{kj}\|) \mathbf{Y}_\ell(\hat{\mathbf{r}}_{kj}) \sum_{c'} W_{\text{down},\ell}^{cc'} h_{k c'}^{\text{pre-down}} \tilde{w}_k$

where the symbols are:

• $\mathbf{r}_{kj} = \mathbf{r}_k - \mathbf{R}_j$ is the vector from the coarse point (nucleus $j$) to grid point $k$, $\|\mathbf{r}_{kj}\|$ is its Euclidean distance, and $\hat{\mathbf{r}}_{kj} = \mathbf{r}_{kj}/\|\mathbf{r}_{kj}\|$ is the unit direction vector.
• $\phi_c(r)$ is the $c$-th radial basis function, with $c = 1,\ldots,16$ (same as $D_{\text{nonl}}$). The specific radial basis chosen is a second-moment Gaussian: $\phi_c(r) = \frac{2}{\text{dim} \cdot (2\pi s_c^2)^{\text{dim}/2}} \frac{r^2}{2s_c^2} \exp(-r^2/(2s_c^2))$ where $\text{dim}=3$ and $s_c$ are 16 scale coefficients evenly spaced between 0.3023 $a_0$ and 2.192 $a_0$ (chosen so that two standard deviations span the range of covalent radii from Pyykkö and Atsumi's compilation for elements up to 118).
• $\mathbf{Y}_\ell(\hat{\mathbf{r}})$ is the vector of all $(2\ell+1)$ spherical harmonics of order $\ell$ evaluated at direction $\hat{\mathbf{r}}$. The spherical harmonic orders used are $\ell = 0, 1, 2, 3$ ($\ell_{\text{max}} = 3$).
• $W_{\text{down},\ell}^{cc'}$ is a learnable $D_{\text{nonl}} \times D_{\text{nonl}} = 16 \times 16$ weight matrix per spherical harmonic order $\ell$ that mixes the 16 input channels before projection onto the product basis.
• $h_{kc'}^{\text{pre-down}}$ is the $c'$-th channel of the pre-downsampled feature at grid point $k$.
• $\tilde{w}_k$ is the per-atom integration weight for grid point $k$.
• $H_{j\ell c} \in \mathbb{R}^{2\ell+1}$ is the coarsened feature vector — it has $2\ell+1$ components (one per spherical harmonic $m = -\ell,\ldots,\ell$) for each channel $c$ and each $\ell$. This is an equivariant feature: if you rotate the molecule, $H_{j\ell c}$ transforms according to the Wigner D-matrix of order $\ell$.

What this operation computes physically: for each coarse point (nucleus), this computes a set of multipole moments of the local feature distribution on the atom's grid. For $\ell = 0$, this is a monopole — the isotropic average of the feature weighted by radial distance and aggregated over the atom's grid points. For $\ell = 1$, this is a dipole — it captures the first-order angular anisotropy (which direction the feature distribution is biased toward). For $\ell = 2$, a quadrupole; for $\ell = 3$, an octupole. Because each channel $c$ uses a different radial basis function $\phi_c$ with a different length scale $s_c$, these moments are radially resolved: channels with small $s_c$ capture features near the nucleus (core region), while channels with large $s_c$ capture features in the valence and tail regions. The learnable mixing matrix $W_{\text{down},\ell}^{cc'}$ allows the model to learn which linear combinations of the 16 local feature channels to project onto which radial basis functions.

Why the $r^2$ factor in the radial basis: the $r^2$ factor in $\phi_c(r)$ suppresses contributions from grid points very close to the nucleus (where $r \to 0$). Near-core densities are extremely high and dominated by the nuclear cusp (the density has a discontinuous derivative at the nucleus), and the exchange-correlation physics in this region is dominated by the simple LDA form (because core electrons behave like a uniform electron gas with slowly-varying density). The $r^2$ factor de-emphasizes this region, allowing the non-local module to focus on the bonding and valence regions (roughly 0.5–4 $a_0$ from the nucleus) where non-local exchange and correlation effects are most important.

Rotational equivariance: the downsampled features $H_{j\ell c}$ are equivariant — if you rotate the molecule, they transform according to the spherical harmonic representation of the rotation group SO(3). This is a fundamental design principle: the XC functional must be rotationally invariant (the energy of a molecule should not depend on its orientation in space), but the intermediate representations can be equivariant (they rotate with the molecule) before eventually being contracted to invariants. Equivariant message passing, pioneered in neural force fields (e.g., NequIP, MACE), has proven vastly more data-efficient than invariant-only message passing because it preserves directional information that can be used to compute anisotropic interactions.

Step 3: Coarse-point mixing and symmetric contraction. After downsampling, the coarse features undergo two transformations on each coarse point independently:

Equivariant linear mixing: $H'_{j\ell c} = \sum_{c'} W_{\text{mix},\ell}^{cc'} H_{j\ell c'}$

where $W_{\text{mix},\ell}$ is a $D_{\text{nonl}} \times D_{\text{nonl}}$ learnable weight matrix for each $\ell$, initialized to the identity. This mixes information across the 16 radial channels within each spherical harmonic order without mixing different $\ell$ (which would break equivariance — tensors of different orders cannot be added).

Symmetric contraction — capturing higher-body-order correlations: After linear mixing, a symmetric contraction of correlation order $\nu = 3$ is applied. Following the MACE architecture (Batatia et al., 2022), the symmetric contraction takes $\nu$-fold symmetric tensor products of the equivariant feature vectors $H'_{j}$ (across all $\ell$ and $c$) and contracts them to scalar invariants using the appropriate Clebsch-Gordan coefficients. The contraction of order $\nu$ generates features that depend on $\nu$-body correlations of the underlying grid-point distribution.

In plain language: the downsampled features $H'_{j}$ encode spatial distributions (multipole moments) of the local XC features around nucleus $j$. The symmetric contraction computes scalar invariants that describe how $\nu$ such distributions interact — for $\nu = 2$, this would capture pairwise interactions (e.g., the correlation between the monopole moment and the quadrupole moment); for $\nu = 3$ (as used here), this captures three-way interactions. This is directly analogous to constructing three-body descriptors in interatomic potentials, but applied to the electron density features rather than atomic positions.

Why correlation order $\nu = 3$: in the MACE literature, correlation order $\nu$ corresponds to the body order of interactions that can be represented. $\nu = 1$ captures one-body terms (local only); $\nu = 2$ captures two-body; $\nu = 3$ captures three-body, and so on. For interatomic potentials, two-body (pairwise) interactions capture much of the physics, but three-body terms are needed for angular dependence (e.g., in water, the H-O-H angle energy depends on the three-body geometry). For an XC functional, non-local exchange between two electrons is fundamentally a two-body interaction, but the effective XC functional (which integrates out the wavefunction) can require higher-body terms. Using $\nu = 3$ enables up to four-body interactions within each non-local layer (starting from the two-body downsampling step plus $\nu$-body contraction). The paper found this sufficient for main-group chemistry accuracy.

Important: shared weights across coarse points. Unlike the original MACE architecture, which uses element-dependent weights to distinguish atom types (carbon vs. oxygen vs. iron), Skala's symmetric contraction uses shared weights across all coarse points. This is because the XC functional depends only on the electron density, not on the nuclear charge directly — the density already encodes the element identity through the number of electrons and the shape of the density distribution. Using shared weights enforces this principle and reduces the parameter count.

Step 4: Upsampling — sending processed coarse information back to the grid. After processing on the coarse points, the enriched features are projected back onto the fine grid using the same radial and spherical harmonic basis:

$h'_{ic} = \sum_\ell \phi_c(\|\mathbf{r}_{ij}\|) \left\langle \mathbf{Y}_\ell(\hat{\mathbf{r}}_{ij}), \sum_{c'} W_{\text{up},\ell}^{cc'} \mathbf{H}_{j\ell c'} \right\rangle$

where $\mathbf{r}_{ij} = \mathbf{r}_i^{(j)} - \mathbf{R}_j$ is the vector from the coarse point to the grid point (opposite direction from downsampling), $\langle\cdot,\cdot\rangle$ denotes the Euclidean inner product between two $(2\ell+1)$-dimensional vectors (the spherical harmonic evaluated at the grid point direction and the processed coarse feature vector), and $W_{\text{up},\ell}^{cc'}$ is another learnable mixing matrix per $\ell$. The output $h'_{ic}$ is a 16-dimensional vector at each grid point $i$ that encodes non-local information from all other grid points on the same atom.

What this operation computes: for each grid point $i$ on atom $j$, this projects the processed coarse features (which contain aggregated information from all grid points on atom $j$) back onto the specific grid point's location, weighted by how strongly the grid point "feels" each angular and radial mode. The inner product $\langle\mathbf{Y}_\ell(\hat{\mathbf{r}}_{ij}), \mathbf{H}_{j\ell c'}\rangle$ quantifies the alignment between the pattern of the coarse feature at order $\ell$ and the direction of the grid point from the nucleus — if the coarse feature "points" in the same direction as the grid point, the alignment is strong and the upsampled contribution is large.

Interpretation as a learned non-local kernel: as shown in Supplementary Information Theorem 1, any rotationally invariant two-body kernel $\kappa(\mathbf{r}_i, \mathbf{r}_k)$ can be expanded as a sum over radial bases and spherical harmonics of the form $\sum_{c_1,\ell} \phi_{c_1}(\|r_i\|) \langle \mathbf{Y}_\ell(\hat{\mathbf{r}}_i), \sum_{c_2} C_{c_1,\ell,c_2} \phi_{c_2}(\|r_k\|) \mathbf{Y}_\ell(\hat{\mathbf{r}}_k) \rangle$. The downsampling step computes the inner integral $\sum_k \phi_{c_2}(\|r_k\|)\mathbf{Y}_\ell(\hat{\mathbf{r}}_k) \tilde{h}_k$ (with $\tilde{h}_k$ being the learnable feature at $\mathbf{r}_k$). The upsampling step applies the outer linear combination. The linear mixing weights $W_{\text{down},\ell}$ and $W_{\text{up},\ell}$ together with the channel mixing learn the effective coefficients $C_{c_1,\ell,c_2}$. In other words, the non-local layer learns an effective non-local interaction kernel directly from data, without specifying its functional form in advance. The radial basis with 16 Gaussian scales provides enough flexibility to represent a wide range of kernel shapes, and the spherical harmonic expansion to $\ell_{\text{max}}=3$ captures angular anisotropy up to octupolar order.

Step 5: Post-upsampling processing and skip connection. The 16-dimensional upsampled feature is post-processed and combined with the original 256-dimensional input via a skip connection:

$h_i^{\text{post-up}} = \sigma(W_{\text{post-up}} h'_i + \mathbf{b}_{\text{post-up}})$

$h_i^{(l+1)} = \sigma\left(W_{\text{cat}} \left[h_i^{(l)}, \exp(-\rho_i) \cdot h_i^{\text{post-up}}\right] + \mathbf{b}_{\text{cat}}\right)$

where $W_{\text{post-up}}$ is a $D_{\text{nonl}} \times D_{\text{nonl}} = 16 \times 16$ weight matrix, $h_i^{(l)}$ is the original 256-dimensional input (before this non-local layer), $h_i^{\text{post-up}} \in \mathbb{R}^{16}$ is the upsampled non-local feature, $\rho_i = \rho_i^{(\uparrow)} + \rho_i^{(\downarrow)}$ is the total electron density at the grid point, $\exp(-\rho_i)$ is an exponential density-dependent gating factor, $[\cdot,\cdot]$ denotes concatenation of the 256-dimensional local feature and the gated 16-dimensional non-local feature into a 272-dimensional vector, and $W_{\text{cat}}$ is a $D_{\text{hid}} \times (D_{\text{hid}} + D_{\text{nonl}}) = 256 \times 272$ weight matrix that projects back to 256 dimensions.

What this operation computes: the skip connection preserves the local information from the input, while the non-local branch provides a correction term. The density gating $\exp(-\rho_i)$ suppresses non-local corrections in regions of high density (near the nucleus), where LDA exchange is known to be a good approximation and non-local effects are small. The weight $W_{\text{cat}}$ learns how to combine the local and non-local features to produce the enriched representation.

Why stacking $L=3$ layers: a single downsampling-upsampling pass captures two-body interactions (Theorem 1 in the supplement). With the symmetric contraction of order $\nu=3$, each layer captures up to $\nu+1 = 4$-body correlations. Stacking $L=3$ such layers with independent parameters extends the effective receptive field: in the first layer, a grid point receives information from all other grid points on its atom; in the second layer, after another round of downsampling and upsampling, the coarse features have been influenced by the first layer's non-local features, so the upsampled information at layer 2 encodes multi-step interactions. However, since each atom's grid points only ever communicate with that atom's coarse point (there is no cross-atom message passing between coarse points), the non-locality is restricted to within-atom exchanges of density features. Cross-atom non-locality enters indirectly: the density features themselves depend on the entire molecular density (through the SCF solution), so information from other atoms is already encoded in the semi-local features at a given grid point. The non-local layers then allow grid points on a given atom to "pool" this cross-atom information through their shared coarse point.

Parameter count of the non-local module: each non-local layer has approximately:

• $W_{\text{pre-down}}$: $16 \times 256 = 4,096$ weights
• $W_{\text{down},\ell}$ (4 matrices, one per $\ell$): $4 \times 16 \times 16 = 1,024$ weights
• $W_{\text{mix},\ell}$ (4 matrices): $4 \times 16 \times 16 = 1,024$ weights
• $W_{\text{up},\ell}$ (4 matrices): $4 \times 16 \times 16 = 1,024$ weights
• $W_{\text{post-up}}$: $16 \times 16 = 256$ weights
• $W_{\text{cat}}$: $256 \times 272 = 69,632$ weights
• Plus biases, symmetric contraction parameters, etc. Across three layers, the non-local module has roughly $3 \times (69,632 + \text{small terms}) \approx$ ~200,000 parameters. The paper reports that Skala has 385,217 total parameters, with 331,265 in the local branch and the remainder (~54,000) in the non-local branch — consistent with our rough estimate given the recurrent structure (the three layers share structure but not weights).

---

Part Three: Output Model and Energy Integration

After three non-local layers, each grid point has a 256-dimensional enriched hidden vector $h_i^{(L)}$. The final component of the architecture is a purely local MLP that maps this to a scalar enhancement factor:

$h_{\text{enh},i} = \sigma_{\text{out}}\left(W_4\,\sigma(W_3 h_i^{(L)} + \mathbf{b}_3) + b_4\right)$

where $W_3$ is a $D_{\text{hid}} \times D_{\text{hid}} = 256 \times 256$ weight matrix, $\mathbf{b}_3$ is a 256-dimensional bias, $\sigma$ is the Swish activation, $W_4$ is a $1 \times D_{\text{hid}} = 1 \times 256$ row vector (projecting to a scalar), $b_4$ is a scalar bias, and $\sigma_{\text{out}}$ is a scaled sigmoid activation function:

$\sigma_{\text{out}}(x) = \frac{2}{1 + \exp(-x/2)}$

Why the scaled sigmoid: the standard sigmoid $1/(1+\exp(-x))$ maps to $(0, 1)$. The scaled version maps to $(0, 2)$, centered around 1 (since $\sigma_{\text{out}}(0) = 1$). This enforces the Lieb-Oxford lower bound: the enhancement factor is non-negative, so the exchange energy density (which has a negative sign in front — see Equation 1) is bounded below. The upper bound of 2 is a recommended value from Kirkpatrick et al. (2021) that covers physically plausible enhancement factors. When the logit is close to zero (the network outputs small values), $h_{\text{enh}} \approx 1$, and the functional defaults to LDA exchange — a sensible baseline that the network can deviate from as needed.

The final energy integration: $E_{\text{xc}}^\theta[\rho] = -\frac{3}{4}\left(\frac{6}{\pi}\right)^{1/3} \sum_{i=1}^{G} h_{\text{enh},i}\left[\rho_i^{(\uparrow)4/3} + \rho_i^{(\downarrow)4/3}\right] w_i$

where the sum runs over all grid points (all atoms), $w_i$ are the merged integration grid weights (incorporating the space-partitioning function), and $h_{\text{enh},i}$ is the scalar enhancement factor at grid point $i$. This gives the total exchange-correlation energy in Hartree atomic units.

Computational scaling: the cost is dominated by the per-grid-point operations. The initial local MLP, the three non-local layers, and the output MLP all scale linearly in the number of grid points $G$. Since $G$ is proportional to the number of atoms for a given grid quality, the overall asymptotic scaling is $O(N_{\text{atoms}})$ for the XC energy evaluation — i.e., $O(N^3)$ when combined with the other Kohn-Sham steps (Coulomb matrix construction in $O(N^3)$, diagonalization in $O(N^3)$). The non-local communication cost is $O(G \cdot \text{[basis size]})$ rather than $O(G^2)$ because each grid point talks only to one coarse point. The constant factor is determined by the network architecture parameters ($D_{\text{hid}}=256$, $D_{\text{nonl}}=16$, $L=3$, $\ell_{\text{max}}=3$, 16 radial basis functions), which are fixed and independent of molecule size.

---

The Two-Stage Training Protocol

The training of Skala is not a standard supervised learning problem because the XC functional is used inside a self-consistent field (SCF) loop at inference time. The densities on which the functional is evaluated during training are not the same as the densities it will produce when deployed self-consistently. The paper addresses this with a two-stage protocol.

Stage 1: Pre-training on Fixed B3LYP Densities

What happens: for each molecule in the training set, the density $\rho_{\text{B3LYP}}$ and the Kohn-Sham orbitals are computed once using the B3LYP functional (hybrid GGA, 20% exact exchange) with the def2-QZVP basis set (or ma-def2-QZVP for anions). These densities and orbitals are frozen — they are not recomputed during training. The XC energy predicted by Skala is evaluated on this fixed density, and the total energy is:

$E_{\text{tot}}^\theta[M] = E_{\text{tot}-\text{xc}}[M] + E_{\text{xc}}^\theta[\rho_{\text{B3LYP}}]$

where $E_{\text{tot}-\text{xc}}[M]$ is the total B3LYP energy minus the B3LYP XC energy (precomputed and fixed). This term includes the kinetic energy, electron-nuclear attraction, Hartree electrostatic energy, and B3LYP's D3(BJ) dispersion correction. By adding $E_{\text{xc}}^\theta$ instead of the B3LYP XC energy, the total energy reflects Skala's XC functional while keeping all other energy components at the B3LYP level.

A reaction energy is then computed as the stoichiometric combination of total energies:

$\Delta E = \sum_M c_M E_{\text{tot}}^\theta[M]$

where $c_M$ are the stoichiometric coefficients (positive for products, negative for reactants). For a simple reaction $A \to B$, $\Delta E = E_{\text{tot}}^\theta[B] - E_{\text{tot}}^\theta[A]$.

Loss function — weighted mean squared error on reaction energies: $\mathcal{L} = \mathbb{E}\left[\frac{\left|\Delta E - \Delta E_{\text{ref}}\right|^2}{10^{-4}E_h + |\Delta E_{\text{ref}}|}\right]$

where $\Delta E_{\text{ref}}$ is the reference reaction energy from CCSD(T)/CBS calculations, $E_h$ (Hartree) is the atomic unit of energy (1 Hartree = 627.5 kcal/mol), and $10^{-4}E_h$ (approximately 0.06 kcal/mol) is a small constant that prevents division by zero when the reference energy is near zero (as happens for some conformational energy differences).

What it computes: each reaction contributes a squared error normalized by its reference energy magnitude. For a reaction with $|\Delta E_{\text{ref}}| = 100$ kcal/mol, the denominator is $0.06 + 100/\text{[conversion]} \approx 100$, so the effective loss weight is $1/100$. For a reaction with $|\Delta E_{\text{ref}}| = 1$ kcal/mol, the denominator is $0.06 + 1 \approx 1$, giving a weight of ~1. This inverse-energy weighting ensures that large-energy reactions (like total atomization energies of hundreds of kcal/mol) do not dominate the loss at the expense of small-energy reactions (like conformational changes of ~1 kcal/mol) that are equally important for chemical accuracy. Without this weighting, the optimizer would focus on fitting the large energies and ignore the subtle energy differences that matter for conformational preferences and non-covalent interactions.

Why train on reaction energies, not total energies: reference total energies from wavefunction methods are generally not as accurate as reaction energy differences, because systematic errors in the wavefunction treatment (basis set incompleteness, missing higher-order excitations) largely cancel when subtracting energies of reactants and products. By training on reaction energies directly, the model learns to reproduce the energy differences that matter for chemistry without needing to fit absolute total energies. This is the standard approach in XC functional development — functionals are judged by their reaction energy errors, not total energy errors.

Why use B3LYP densities rather than a simpler functional: B3LYP densities are of reasonable quality — better than LDA or GGA densities — and are available at scale for the entire training dataset. The cost of generating B3LYP densities for ~400,000 reactions is amortized over training and re-training, making it a one-time computational investment. The choice of B3LYP over, say, PBE or r2SCAN is pragmatic: B3LYP is the most widely available and well-tested hybrid functional, and its densities are known to be reasonably accurate for main-group chemistry (dipoles, geometries). However, as the paper discusses in Section 4, training on B3LYP densities creates a gap when the functional is deployed self-consistently — Skala's own SCF densities will differ from B3LYP densities, and the accuracy on B3LYP densities is an upper bound on self-consistent accuracy. This motivates the fine-tuning stage.

Grid augmentation during pre-training: to prevent the model from overfitting to the specific grid used in training and to make the functional robust to grid choice (an important practical consideration for DFT users who may choose different grid settings), the pre-training computes the density features on four different radial integration schemes at grid level 1 (coarser grids, ~fewer points per atom) and randomizes the space-partitioning scheme among three options. The stochasticity regularizes the model against grid-dependent artifacts.

Stage 2: SCF Fine-Tuning

After pre-training, the model produces accurate reaction energies when evaluated on B3LYP densities, but may perform worse when evaluated self-consistently — a phenomenon known as density-driven error: even if the functional is perfect when given the exact density, the density it produces self-consistently may be degraded, leading to an overall energy error that combines functional error and density error.

The fine-tuning procedure (Supplementary Information B.5):

• The pre-trained model weights are loaded, and the optimizer is reinitialized (Adam + Muon, as detailed below).
• For each training step, a reaction is sampled from the dataset, and a self-consistent field calculation is run for each molecule in the reaction using Skala's own XC functional. The SCF loop starts from a B3LYP initial guess and iterates to convergence using DIIS (Direct Inversion in the Iterative Subspace) with a window of 8 iterations, a density convergence threshold of $10^{-8}$, an energy change threshold of $5 \times 10^{-6}E_h$, and an orbital gradient norm threshold of $10^{-3}E_h$, capped at 40 iterations. The SCF uses Treutler-Ahlrichs radial grids at level 3 (finer grids than pre-training) with the def2-QZVP basis set and density fitting.
• The converged self-consistent density $\rho_{\text{Skala}}$ is used to compute the total energy $E_{\text{tot}}^\theta$ and the reaction energy $\Delta E$.
• The same weighted MSE loss (Equation 31) is minimized, but now the model's own SCF densities are the input, closing the gap between training conditions and deployment conditions.

Crucial property — the Envelope Theorem: the paper invokes the envelope theorem to justify that no backpropagation through the SCF loop is needed. Let $C^* = \arg\min_C E_{\text{tot}}^\theta[C]$ be the minimizing orbital coefficients (density matrix). At the minimum, the gradient of $E_{\text{tot}}$ with respect to $C$ is zero (first-order stationarity condition). Therefore, the gradient of the minimum total energy with respect to the model parameters $\theta$ is simply the gradient of the energy evaluated at the optimal orbitals:

$\nabla_\theta E_{\text{tot}} = \nabla_\theta E_{\text{tot}}^\theta[C^*]$

That is, you can take the gradient of the energy expression with respect to $\theta$ as if the density were fixed, even though the density depends on $\theta$ through the SCF procedure. This avoids the expensive and numerically unstable process of backpropagating through dozens of SCF iterations, making SCF fine-tuning computationally feasible. The training signal comes from the fact that $C^*$ changes as $\theta$ changes across training steps, so the model sees evolving SCF densities over the course of fine-tuning.

Fine-tuning hyperparameters: constant learning rate $10^{-5}$ for both Muon and Adam, per-parameter RMS gradient clipping with threshold $10^{-4}$, exponential moving average (EMA) on model weights with decay 0.999, batch size of 1 reaction per GPU on a single A100 GPU, run for 20,000 steps. The adaptive sampling probabilities from pre-training are frozen at their final pre-training values. The SCF fine-tuning is monitored by tracking (i) the gap between the accuracy on B3LYP densities vs. SCF densities and (ii) the dipole moment error on a validation set (a sensitive probe of density quality).

What happens during fine-tuning (Figure 3a): the paper reports a striking empirical observation. Over the first ~2,000 fine-tuning steps, SCF convergence is erratic (some molecules fail to converge, and the converged subset shows fluctuating dipole errors). After ~2,000 steps, all systems converge, and from this point onward, both the reaction energy error and the dipole moment error improve simultaneously. This is important because it suggests that the functional is not engaging in error compensation — it is not improving energies by producing worse densities that happen to cancel functional errors. Instead, the functional appears to be approaching a property of the exact functional: improving the energy necessarily improves the density, because the exact functional has the density as its functional derivative. The final Skala dipole errors (Table 3b) fall well below B3LYP's dipole errors, confirming that the SCF densities are genuinely improved, not just the energies.

Dataset Sampling and the Relative Excess Loss (REL) Mechanism

The training dataset is heterogeneous (spanning atomization energies, barrier heights, non-covalent interactions, atomic properties, etc., with wildly varying dataset sizes) and imbalanced. Uniform sampling would cause training to be dominated by the largest datasets (e.g., ~120,000 TAEs) while underrepresenting smaller but chemically important sets (e.g., proton affinities, transition metals). The paper introduces a two-level sampling strategy with adaptive reweighting:

Initial sampling probabilities: each dataset $i$ is assigned a probability $p_i \propto \alpha_{c(i)} \cdot |D_i|$, where $|D_i|$ is the number of reactions in the dataset, $c(i)$ is the category (atomization energies, barrier heights, non-covalent interactions, etc.), and $\alpha_c$ is a per-category weight. The category weights $\{\alpha_c\}$ are solved such that the total sampling probability allocated to each of the nine categories matches prescribed target proportions.

Relative Excess Loss (REL): every 25,000 training steps, the current model is evaluated on each dataset's holdout set, producing a mean absolute error (MAE) $\ell_i$. The REL for dataset $i$ is:

$R_i = \frac{\max(\ell_i - t_i, 0)}{n_i}$

where $t_i$ is the minimum MAE among baseline functionals (revPBE, r2SCAN, B97M-V, B3LYP, M06-2X, ωB97X-V, ωB97M-V) evaluated on the same holdout set, and $n_i$ is the difference between the median baseline MAE and the minimum baseline MAE ($n_i = \text{median}_i - t_i$).

What REL measures: $R_i = 0$ means Skala is at least as accurate as the best baseline functional on dataset $i$. $R_i = 1$ means Skala performs at the median baseline level. $R_i > 1$ means Skala underperforms even the median baseline. The total REL $R = \sum_i R_i$ counts, in a soft sense, how many datasets have not yet reached best-baseline accuracy.

Adaptive reweighting: after computing REL, the sampling probabilities are updated via multiplicative exponential weighting:

$p_i \leftarrow p_i \cdot \exp(\eta R_i)$

with $\eta = 0.025$, followed by renormalization $p_i \leftarrow p_i / \sum_j p_j$. Datasets where Skala is already better than the best baseline ($R_i = 0$) receive no additional weight; datasets where Skala is still behind ($R_i > 0$) receive increased sampling probability. This provides a principled stopping criterion: when REL approaches zero on all datasets, training has saturated the available signal, and further improvement would require more data or a better architecture.

Why this is superior to uniform sampling or loss-based reweighting: uniform sampling would overfit to the largest datasets while underfitting the smaller ones. Loss-based reweighting (up-weighting high-loss examples) can be counterproductive because some datasets have inherently higher irreducible error due to reference label noise (even CCSD(T)/CBS has finite accuracy). The REL-based approach uses baseline functional performance as a soft target: it stops trying to reduce error on a dataset once Skala surpasses the best traditional functional, since further reduction may fit noise rather than signal.

Optimizer Configuration and Training Details

Optimizer split — Muon for hidden weights, Adam for biases and output layer:

The paper uses a hybrid Muon-Adam strategy. Model parameters are split into two groups:

• Muon group: all hidden weight matrices (parameters with ≥2 dimensions, excluding the final output layer). Muon runs standard SGD with momentum ($\text{momentum} = 0.95$) and then forces the update matrix to be orthogonal (via Newton-Schulz iteration), which has been shown to improve training stability and convergence in large-scale deep learning (Keller Jordan, 2024; Liu et al., 2025).
• Adam group: all low-dimensional parameters (biases, 1D parameters) and the final output layer's weights. Adam uses $\beta = (0.9, 0.95)$ and $\epsilon = 10^{-10}$.

Both groups share a cosine learning rate schedule with linear warmup:

$\text{lr}(t) = \begin{cases} \text{lr}_{\text{max}} \times \frac{t}{t_{\text{warmup}}}, & t \leq t_{\text{warmup}} \\ \text{lr}_{\text{min}} + \frac{\text{lr}_{\text{max}} - \text{lr}_{\text{min}}}{2}\left(1 + \cos\left(\frac{t - t_{\text{warmup}}}{T - t_{\text{warmup}}}\pi\right)\right), & t > t_{\text{warmup}} \end{cases}$

with peak learning rates $\text{lr}_{\text{max}}^{\text{Muon}} = 7 \times 10^{-4}$, $\text{lr}_{\text{max}}^{\text{Adam}} = 1.5 \times 10^{-4}$, warmup period $t_{\text{warmup}} = 50,000$ steps, and total training steps $T = 1,000,000$. Gradient clipping is applied with threshold 0.0001. An exponential moving average of model weights is maintained with decay 0.9999 for pre-training and 0.999 for fine-tuning. Training uses 8 A100 GPUs with a minibatch size of 1 reaction per GPU (effective batch size 8) in a distributed data parallel configuration.

Why Muon: the orthogonalization step in Muon has been empirically shown to reduce the number of training steps needed for convergence and to improve final accuracy in large transformers. For Skala, the hidden weight matrices in the local MLPs and non-local layers have the right structure (2D matrices), while biases and the final projection layer are better handled by standard Adam (which adapts per-parameter learning rates based on gradient history).

Parameter initialization: the local MLP weights are initialized using Xavier (Glorot) uniform initialization, and biases are set to zero. The non-local layer weights are initialized with a modified Xavier scheme that accounts for weight sharing across tensor components of the same spherical harmonic order (from the e3nn library, Geiger & Smidt, 2022).

Training Data Summary

The training set (Extended Data Table 6) comprises approximately 409,424 reactions after removing overlap with test sets (GMTKN55, W4-17) and reserving holdout validation sets. The data spans:

• Set A (thermochemistry, basic properties, distortions, transition metals): ~284,000 reactions, including ~118,000 total atomization energies (MSR-ACC/TAE), ~16,000 ionization potentials, ~13,000 electron affinities, ~20,000 proton affinities, ~19,000 distorted geometries, plus atomic properties and GDB9-W1-F12 atomization energies.
• Set B (reactions and barrier heights): ~114,000 reactions, mainly from MSR-ACC/Reactions (elementary reaction paths with up to 13 non-hydrogen atoms) plus BH9 (public barrier heights).
• Set C (non-covalent interactions and conformations): ~41,000 reactions, including ~41,000 conformers (MSR-ACC/Conf), ~2,200 water dimer structures, ~3,700 NCI clusters, and ~43,000 NCIAtlas + DES370K interactions.
• Other (atomics, transition metals): ~16,000 atomic properties and small transition-metal datasets.

4. Key Insights and Innovations

Innovation 1: Learned non-locality as a substitute for Hartree–Fock exchange — breaking the Jacob's ladder paradigm

The organizing principle of exchange-correlation functional development for three decades has been Jacob's ladder: to get better accuracy, you must climb to a higher rung, which adds a more expensive ingredient — and that ingredient has almost always been exact Hartree–Fock exchange in some form (global, range-separated, or local). This paper's most fundamental conceptual move is to reject the premise that non-locality must be provided by exact exchange — and to demonstrate that a neural network can learn whatever non-local interactions the data demand directly from cheap semi-local inputs, at a cost far below that of computing the Hartree–Fock integrals.

This is not an incremental refinement of the hybrid paradigm. It is an alternative to it. Hybrid functionals improve accuracy by adding a fixed fraction of exact exchange — a physically motivated but rigid non-local ingredient that is identical for all systems and all chemical environments. Skala instead learns an adaptive, data-driven non-local interaction kernel whose effective range and functional form emerge from the training data rather than from an analytic expression derived from the Coulomb operator. The evidence is in the architecture's design: the non-local layers implement a learnable two-body kernel via the downsampling–upsampling mechanism (Section 3.3.4), and the symmetric contraction extends this to higher-body-order correlations — but the kernel's shape is determined entirely by the learned weights, not by $1/r_{12}$. This is a qualitatively different source of non-locality: it is whatever the data say is needed, not what physicists in the 1990s argued should be added.

The significance of this move extends beyond the accuracy numbers. For decades, the DFT community has debated the "optimal" fraction of exact exchange — 20% (B3LYP), 25% (PBE0), 54% (M06-2X), or range-separated with system-dependent parameters. Skala sidesteps this debate entirely. The effective non-local exchange measured by the Grimme molecular "thermometer" (Extended Data Figure 9a) shows Skala achieving ~54% effective non-local exchange — close to the CCSD(T) reference of 60%, and competitive with the best hybrids — without a single line of code that computes a Hartree–Fock exchange integral. This is evidence that the learned kernel captures physically meaningful non-local exchange effects, not just an arbitrary correction that happens to work on the test set.

The paper explicitly frames this as disrupting the accuracy–cost trade-off that defines Jacob's ladder (Figure 1a vs. Figure 1c). The ladder says: accuracy comes from climbing rungs, and each rung costs more. Skala says: accuracy comes from data and learned representations, and the computational cost is determined by the architecture, not by the physics of the ingredients. This is a reframing of the problem from "what physics should we add?" to "what architecture can learn from data?" — the same shift that deep learning brought to computer vision (from hand-designed filters to learned convolutional kernels) and natural language processing (from feature engineering to learned token representations).

The ~50% error reduction from adding the non-local branch to the local baseline (Figure 2a, left) quantifies the value of this learned non-locality: a purely local model with the same training data achieves GGA-level accuracy on Diet GMTKN55, while the non-local model reaches hybrid-level. This is a clean ablation that isolates the contribution of the architectural innovation from the contribution of the large dataset — the data are the same, but the model's ability to use them depends critically on having a mechanism to learn non-local interactions.

A subtle but important distinction: this is not the first paper to propose learning non-local XC effects. Schmidt et al. (2019), Margraf & Reuter (2021), and Gao et al. (2024) all explored non-local architectures for XC. What distinguishes Skala is the combination of scalability (the coarse-point mechanism keeps the cost at O(N³) with a manageable prefactor, tested on molecules with up to 910 atoms / 19,800 basis functions in Figure 4) and generality (tested on the full GMTKN55 benchmark, not model systems or narrow problem classes). Prior non-local ML functionals either could not scale to production molecular sizes or were only tested on restricted chemical spaces. Skala demonstrates that learned non-locality can work at scale, on real benchmarks, and against the best traditional functionals.

Finally, the paper positions this as a deliberate departure from feature engineering. The input features are exactly the standard meta-GGA ingredients — density, gradient, kinetic energy density, nothing more. Any non-locality in the functional comes from the architecture's ability to mix information across grid points, not from cleverly designed non-local input descriptors (as in CIDER or DM21). This is the "deep representation learning" bet: given enough data, the model will learn better representations of non-locality than a human can design. The results validate this bet for main-group chemistry.

Innovation 2: Data-driven emergence of physical constraints — an alternative to hard-coding exact conditions

Traditional XC functional development has placed enormous emphasis on exact constraints — known mathematical properties that the exact functional must satisfy, such as the uniform electron gas limit, coordinate scaling relations, the Lieb-Oxford bound, and positivity conditions on correlation energy components. The standard methodology (epitomized by the SCAN functional, Sun et al., 2015) is to carefully design functional forms that satisfy as many exact constraints as possible, and then fit any remaining free parameters to data. This approach has been enormously successful — SCAN satisfies 17 known exact constraints and is widely regarded as the best non-empirical meta-GGA. The assumption embedded in this methodology is that satisfying constraints is necessary for robustness and transferability, and that data alone cannot be trusted to learn them.

Skala takes the opposite approach: it builds in only the most energetically critical constraints (the Lieb-Oxford lower bound via the bounded sigmoid, high-density uniform coordinate scaling via the LDA prefactor, and size-consistency via the integral form) and leaves everything else to be learned from data. The model makes no attempt to enforce, for example, the positivity of the kinetic correlation energy $T_c$ — a constraint derived from the fact that electron correlation increases the kinetic energy (electrons move faster to avoid each other). The implicit bet is that with enough diverse training data, the model will learn these constraints as emergent properties because violating them produces physically incorrect predictions that the data penalize.

Figure 2b provides a striking test of this bet. When trained only on the MSR-ACC/TAE25 dataset (~730 atomization energies of small molecules), Skala's predicted $T_c$ for atoms violates the positivity constraint: the kinetic correlation energy is negative for several elements, meaning the model has not learned that correlation should increase kinetic energy. But when trained on the full dataset (~400,000 energy differences spanning diverse chemistry), the constraint emerges: $T_c$ becomes positive for all atoms, and the values are quantitatively close to reference data obtained by reverse-engineering CCSD densities in a pure Kohn-Sham framework (Figure 14 in the supplement shows this is consistent across five random seeds). This is not a small effect — the constraint goes from being largely violated to consistently satisfied purely as a function of training data diversity.

The intellectual significance of this finding is that it suggests a different relationship between data and constraints than the traditional view. The traditional view is: constraints are prior knowledge that compensate for limited data; they ensure the functional behaves physically in regions of chemical space where training data are sparse. The Skala result suggests a more nuanced picture: some constraints may be emergent properties of sufficient data diversity, not separate requirements that must be imposed from outside. This does not mean all constraints are learnable — the Lieb-Oxford bound and coordinate scaling are enforced because they govern the high-density limit where training data are effectively absent (core electrons in heavy atoms). But for constraints that govern chemically accessible regions — like $T_c$ positivity for valence electrons in atoms — the data may speak loudly enough that explicit enforcement is unnecessary.

This has practical implications for future functional development. If constraints can be learned, then the design space for ML functionals is much larger than previously thought: you do not need to restrict the functional form to guarantee constraint satisfaction; you need only enough diverse data and a sufficiently expressive architecture, and the constraints will emerge. This is a fundamentally different design philosophy from the Perdew tradition of constraint-based functional construction, and it opens the door to architectures that would be impossible to constrain analytically.

A key nuance the paper acknowledges: $T_c$ in Skala is evaluated in the generalized Kohn-Sham (GKS) framework, where meta-GGA functionals have more freedom to move energy between kinetic and exchange-correlation components. In a GKS framework, $T_c$ is not required to be positive — a negative $T_c^{\text{GKS}}$ could be compensated by a more negative $E_{\text{xc}}$ without violating any physics. The remarkable observation is that Skala nevertheless learns $T_c^{\text{GKS}} \approx T_c^{\text{KS}} > 0$ — it chooses the physically more natural decomposition, even though it has the freedom to do otherwise. This suggests that the data favor the pure Kohn-Sham partitioning, providing an empirical argument for the physical meaningfulness of the pure-KS decomposition that is independent of theoretical arguments.

This insight is one of the most conceptually novel in the paper because it inverts the traditional relationship between theory and data in DFT: rather than theory constraining what the data can teach, the data are used to test which theoretical constraints are truly necessary and which are automatically satisfied by nature's own patterns. This is a rare example of machine learning not just matching human-designed approaches but revealing something about the underlying physics — that some exact constraints may be redundant with sufficient chemical information.

Innovation 3: The two-stage training protocol as a scalable solution to the self-consistency gap

Training an XC functional is not a standard regression problem because the functional is evaluated inside a self-consistent field loop at inference time. The density on which the functional is trained is not the density it will produce when deployed — this is the self-consistency gap, and managing it has been one of the central practical challenges in ML functional development. Prior work has explored various solutions: training self-consistently from the start (prohibitively expensive at scale; SCF convergence is unreliable with untrained models), using exact densities from wavefunction methods (only available for very small systems), regularizing with the Kohn-Sham equations (adds complexity and cost), or simply training on fixed approximate densities and accepting the self-consistency gap (which can lead to density-driven errors that degrade total energies).

The paper's solution — a long pre-training phase on fixed B3LYP densities followed by a short fine-tuning phase on self-consistent densities — is not conceptually complex, but the empirical observation about what happens during fine-tuning is what makes it an innovation rather than just an engineering trick. Figure 3a shows the dynamics: during the first ~2,000 fine-tuning steps, SCF convergence is erratic, and the converged subset shows unstable dipole errors. After ~2,000 steps, all systems converge, and from that point onward, both the reaction energy error and the dipole moment error improve simultaneously. The reaction energy MAE drops from ~1.5 kcal/mol to ~1.0 kcal/mol, while the dipole RMSE drops from ~6% to ~5%.

This simultaneous improvement is significant because it suggests the functional is not resorting to error compensation — the well-known pathology where an approximate functional improves its energy predictions by producing worse densities that happen to cancel the functional's intrinsic error. Error compensation is the dark side of empirical fitting in DFT: it can produce excellent benchmark numbers while degrading the underlying physical quantities, leading to failures on properties that were not directly trained on. The classic cautionary tale is Medvedev et al. (2017), who showed that as functionals have become more empirically fitted over time, their self-consistent densities have actually gotten worse for some systems, even as energy errors improved. The paper's observation that energies and densities improve together during Skala's fine-tuning is evidence that the model is approaching a genuine physical improvement, not exploiting cancellation.

The theoretical justification for why this should happen is the envelope theorem: at the SCF minimum, the gradient of the total energy with respect to model parameters does not require backpropagating through the SCF loop because the density is stationary with respect to orbital variations. This means the fine-tuning can treat the SCF densities as fixed for gradient computation, even though they evolve across training steps. The paper leverages this to make SCF fine-tuning computationally feasible — no expensive unrolling of the SCF iterations — while still closing the gap between training and deployment densities.

What makes this a genuine innovation rather than routine transfer learning is the scale at which it works. Prior work (e.g., Li et al., 2021; Kalita et al., 2022) explored training with self-consistency but was limited to small datasets and small molecules. Skala demonstrates that a simple two-stage protocol, with pre-training on approximate densities at scale and brief fine-tuning on self-consistent densities, is sufficient to produce a functional that is both accurate in energies and faithful in densities across broad chemical space. The final dipole errors (Table 3b: RMSE 4.43% for Skala vs. 5.84% for ωB97M-V, 7.09% for B3LYP) confirm that the self-consistent densities are genuinely improved, not just the energies.

This protocol is likely to become a template for future ML functional training because it decouples two problems that are difficult to solve simultaneously: learning the regression mapping (which needs scale) and achieving self-consistency (which needs quality). By sequencing them, the paper shows you can have both without either the computational cost of full self-consistent training or the accuracy degradation of ignoring the gap.

Innovation 4: Verifier-free accuracy — establishing that learned functionals generalize without overfitting to benchmark chemistry

A persistent concern with empirically fitted functionals — whether traditional (M06-2X, with ~30 parameters fit to databases) or ML-based (DM21, with thousands of neural network parameters) — is that their strong benchmark performance reflects overfitting to the training distribution rather than learning the underlying physics. This concern is magnified for deep learning models, which have the capacity to memorize patterns that do not generalize. The DFT community has developed specific diagnostics for this: testing on systems with different chemical bonding from the training set, evaluating on properties not directly trained on (geometries, dipoles), and checking that improvements in energy errors are not accompanied by degradation in density errors.

Skala passes these diagnostics at a scale that reset expectations for what a learned functional can achieve. The paper's test-set overlap removal is aggressive: all systems with more than two atoms whose molecular graphs match those in GMTKN55 or W4-17 are excluded from training (Section 7.2). This means that for the majority of the test molecules, Skala has never seen any molecule with the same connectivity — the test evaluates compositional generalization to new compounds, not configurational generalization to new geometries of the same compounds. This is a stricter standard than many prior ML functional evaluations, which sometimes train and test on different conformations of the same molecules.

The evidence for genuine generalization comes from multiple angles:

• GMTKN55 breakdown (Table 1): Skala is the best functional on 32 of 55 subsets, more than all other functionals combined. The subsets span qualitatively different chemistry: atomization energies (W4-11), proton affinities (PA26), non-covalent interactions (S66, WATER27), barrier heights (BH76, BHDIV10), self-interaction-error-sensitive systems (SIE4x4), and artificial "mindless" molecules (MB16-43). Performing well on all of these simultaneously is not something a model can achieve by memorizing — the chemical patterns are too diverse.

• W4-17 atomization energies (Extended Data Figure 8a): on the 183 single-reference systems, Skala achieves 0.92 kcal/mol MAE — within chemical accuracy. The multi-reference systems (17 molecules) show larger errors (4.61 kcal/mol), which the paper attributes to lack of multi-reference training data, not to overfitting — the model degrades gracefully rather than failing catastrophically.

• Equilibrium geometries (Table 3c): Skala's geometry optimization accuracy matches or exceeds the best hybrid functionals across four independent benchmark sets (LMGB35, HMGB11, CCse21, W4-11-GEOM). Geometry optimization involves forces — gradients of the energy with respect to nuclear positions — which were not explicitly trained on. Good geometry performance means the potential energy surface is not just accurate at energy minima but has the right curvature in all directions, a much stronger test than single-point energies alone.

• The GEO metric (Supplementary Figure 15a): the additional error from using Skala to both optimize the geometry and compute the energy — instead of using reference geometries for energy evaluation — is ~0.02 kcal/mol on W4-11-GEOM, comparable to ωB97M-V. This quantifies what chemists care about in practice: when you use the functional for a full structure prediction + energy calculation, how much additional error do you incur? Skala's GEO being small and comparable to the best hybrids means its potential energy surfaces are globally accurate, not just at the equilibrium geometries seen in training.

• SN2 reaction barriers (Extended Data Figure 9b): these are a classic stress test for XC functionals because they involve partially broken bonds in the transition state, where non-local exchange effects are critical. Skala achieves errors ~2.5× smaller than the best meta-GGAs, comparable to the best hybrids, without computing exact exchange — direct evidence that the learned non-local kernel captures the exchange stabilization of the transition state.

The underlying innovation here is not a specific technique but a standard of evidence. By testing on held-out chemical graphs, held-out reaction types, properties not in the training objective (geometries, dipoles, GEO), and by showing systematic improvement trends with data rather than sharp cliffs at the training distribution boundary, the paper makes a case that the functional has learned something about exchange and correlation that generalizes, rather than memorizing benchmark answers. This shifts the burden of proof for future ML functionals: demonstrating compositional generalization at this level of rigor becomes the expected standard, not an aspirational goal.

Furthermore, the data ablation study (Figure 2a, right, and Figure 2b) shows monotonic improvement with data diversity: adding more types of chemistry to the training set never hurts and often substantially helps. This is precisely the pattern one expects from a model that is learning generalizable features — it is in the underfitting regime, where more data improves representations — rather than in the overfitting regime, where more data would not help because the model is already memorizing. This positions Skala as a framework for systematically improvable DFT: the path to even better accuracy is simply more and more diverse training data, not a better functional form.

Coupled with the observation that physical constraints emerge from data (Innovation 2), this suggests a synthesis: with sufficient data diversity, a sufficiently expressive architecture with minimal built-in constraints can learn both the empirical patterns of chemical accuracy and the underlying physical constraints that ensure robustness — without the explicit constraint engineering that has defined functional development for decades. Whether this synthesis holds for constraints that govern asymptotic limits far from the training data (e.g., the uniform electron gas at metallic densities, or the high-density limit of correlation) remains an open question, but for main-group chemistry — where training data are abundant — the evidence is compelling.

5. Experimental Analysis

Evaluation Methodology

Dataset. The primary evaluation uses the GMTKN55 benchmark (Goerigk et al., 2017), comprising 55 subsets covering five categories: basic properties, thermochemistry, kinetics, intermolecular non-covalent interactions, and conformational energies. A representative 100-reaction subset, Diet GMTKN55 (Gould, 2018), is used for faster iteration in ablation studies. For atomization energies specifically, the W4-17 benchmark (Karton et al., 2017) provides 200 diverse total atomization energies of small molecules labeled at CCSDTQ5/CBS level (the W4 protocol), including 183 single-reference and 17 multi-reference cases. Additional evaluation sets test density quality (dipole moments: 151 structures from Hait & Head-Gordon, 2018), equilibrium geometries (LMGB35, HMGB11, CCse21, W4-11-GEOM), and strained conformer energies (Wiggle150). Training data overlap with GMTKN55 and W4-17 is removed based on molecular graphs for all systems with more than two atoms (Section 7.2).

Base model. All experiments use Skala-1.1, the best model from a data-ablation-driven development process, with ~385,000 parameters total (331,265 in the local branch). Prior versions (Skala-1.0, trained on ~150,000 data points with a less refined architecture) are superseded. The ablation studies comparing local vs. non-local architectures train model instances with different random seeds (affecting initialization and data ordering) from the same training data and evaluate on Diet GMTKN55.

Metrics. The primary metric is WTMAD-2 (weighted total mean absolute deviation) in kcal/mol on GMTKN55, which combines weighted MAEs across all 55 subsets to produce a single number capturing overall accuracy for general main-group chemistry. The weights follow the scheme in Goerigk et al. (2017). For individual subsets and benchmarks, unweighted MAE is reported. For dipole moments, regularized RMSE (following Hait & Head-Gordon, 2018) is reported as a percentage. For geometries, average absolute errors in bond lengths (Å) and bond angles (degrees) against semiexperimental or high-level reference structures. For the GEO metric, MAE of the energy difference between evaluating a functional on its own optimized vs. reference geometry. SCF convergence statistics (percentage converged per retry level, average iterations) are reported in Tables 3 and 4.

Baselines. Three GGA and meta-GGA functionals: revPBE (Zhang & Yang, 1998) with D3(BJ), r2SCAN (Furness et al., 2020) with D3(BJ), B97M-V (Mardirossian & Head-Gordon, 2015) with VV10. Four hybrid functionals: B3LYP (Becke, 1993; Stephens et al., 1994) with D3(BJ), M06-2X (Zhao & Truhlar, 2008) with D3(0), ωB97X-V (Mardirossian & Head-Gordon, 2014) with VV10, ωB97M-V (Mardirossian & Head-Gordon, 2016) with VV10. For geometry optimization benchmarks, GFN2-xTB (Bannwarth et al., 2019) is also included. Baselines are evaluated with def2-QZVP basis sets (ma-def2-QZVP for anions) using the PySCF package.

Compute accounting. Computational cost is measured empirically as wall-clock time per SCF step (XC integration only, excluding Coulomb matrix construction, diagonalization, and SCF convergence) as a function of molecular size, measured using the GauXC library on GPU (NVIDIA A100) and CPU (16 OpenMP threads, AMD EPYC 7V13). Molecules range from 907 to 19,816 atomic orbitals (47–910 atoms). Fitted power-law exponents $n$ in $a N_{\text{orbitals}}^n$ characterize asymptotic scaling. All computations use def2-TZVP basis set and GauXC's GM3 integration grid.

Cross-validation / statistical protocol. The data ablation study in Figure 2a trains five models per condition (different random seeds affecting initialization and data ordering) and reports both individual model instances (dots) and averages (bars), providing a measure of inter-seed variance. The model with the lowest total REL (Section 2.4 of the Further Details) on the validation sets is selected as Skala-1.1. For the GMTKN55 benchmark, standard practice is used: evaluate on all subsets, report WTMAD-2. No cross-validation over test-set splits is performed for the main benchmark results. The Diet GMTKN55 evaluation used for ablations contains 100 representative reactions designed to correlate well with the full WTMAD-2.

---

Main Quantitative Results

GMTKN55: Overall Accuracy Across Main-Group Chemistry

Headline result (Figure 1c, Table 1): Skala-1.1 achieves a WTMAD-2 of 2.80 kcal/mol on the full GMTKN55 benchmark, surpassing all tested functionals including the best hybrid, ωB97M-V (3.23 kcal/mol), and the best meta-GGA, B97M-V (5.56 kcal/mol). This represents a ~13% relative improvement over the best hybrid and approximately a 2× reduction in error relative to the best meta-GGA.

The breakdown by individual subset (Table 1) reveals the breadth of Skala's advantage: it is the best-performing functional on 32 of 55 subsets, more than all other functionals combined (ωB97M-V: 11, ωB97X-V: 6, M06-2X: 3, B3LYP: 1, B97M-V: 0, r2SCAN: 0, revPBE: 0). Skala is the best or second-best on 44 of 55 subsets.

Performance by category:

• Basic properties (18 subsets): Skala achieves the lowest MAE on 13 of 18. On SIE4x4 (self-interaction-error-related problems), a long-standing challenge for semi-local functionals, Skala achieves 12.68 kcal/mol — worse than M06-2X (8.67), ωB97X-V (11.48), and ωB97M-V (10.75), but dramatically better than revPBE (23.42) and r2SCAN (18.17). On DC13 (13 difficult cases for DFT), Skala's 2.20 kcal/mol is far below the next best (B97M-V: 5.20; ωB97M-V: 5.40).

• Reaction energies (9 subsets): Skala is best or second-best on 7 of 9. On MB16-43 (decomposition energies of artificial "mindless" molecules — a test of robustness to unusual bonding), Skala achieves 2.57 kcal/mol, compared to 14.52 for ωB97M-V, 14.69 for r2SCAN, and 35.97 for B97M-V. This subset is a stress test for non-local exchange: the mindless molecules contain unusual bonding patterns not seen in typical training sets, and the large improvement (>5× over ωB97M-V) suggests Skala's learned non-locality generalizes to out-of-distribution bonding. On ISO34 (isomerization energies of medium-sized organics), Skala achieves 0.22 kcal/mol vs. 0.62 for ωB97M-V. On C60ISO (relative energies of C60 isomers), Skala's 8.16 kcal/mol is substantially worse than B3LYP (2.40) and B97M-V (4.85), indicating a failure mode — possibly related to the large system size and extensive conjugation, which may lie outside the training distribution.

• Barrier heights (6 subsets): Skala is best on 2 subsets (BHDIV10: 0.53 kcal/mol, WCPT18: 0.41 kcal/mol) and competitive on the others. The overall barrier height MAE is substantially lower than any meta-GGA and close to the best hybrids. On BH76 (barrier heights of various reaction types), the progression is clear: revPBE 7.84 → r2SCAN 7.01 → B97M-V 4.15 → ωB97M-V 1.33 vs. Skala 1.75 kcal/mol. Skala's error of 1.75 is higher than ωB97M-V (1.33) but far below any meta-GGA.

• Intermolecular NCIs (12 subsets): Skala is best on 3 subsets (HAL59: 0.23, S66: 0.07, PNICO23: 0.19), competitive on most others, but shows notably higher errors on CHB6 (cation-neutral dimers: 1.29 vs. ωB97M-V 0.96) and WATER27 (1.25 vs. B97M-V 0.77, ωB97M-V 1.10).

• Intramolecular NCIs (conformational energies, 10 subsets): Skala is best on 4 subsets (BUT14DIOL: 0.06, SCONF: 0.09, IDISP: 0.60, ACONF: 0.07). On IDISP (intramolecular dispersion interactions), Skala's 0.60 kcal/mol substantially outperforms all other functionals (next best: ωB97M-V: 1.63; B97M-V: 3.18), despite Skala using a fixed D3 dispersion correction rather than learning dispersion — this suggests the learned non-local kernel is capturing some short-to-medium-range correlation effects that complement the D3 long-range correction.

The significance of the 2.80 kcal/mol WTMAD-2 is not just that it beats prior functionals, but that it does so at the computational cost of a meta-GGA (see Computational Cost section below). The traditional accuracy–cost trade-off would place a functional with this accuracy at the hybrid or double-hybrid level — Skala achieves it at the third rung.

W4-17: Atomization Energies

Headline result (Extended Data Figure 8a): On the full W4-17 benchmark (200 molecules), Skala achieves an MAE of 1.23 kcal/mol, compared to 2.06 for ωB97M-V, 2.90 for B97M-V, 4.51 for r2SCAN, and 8.09 for revPBE. On the 183 single-reference systems (where the reference W4 protocol is most reliable), Skala reaches 0.92 kcal/mol — within the ~1 kcal/mol threshold for chemical accuracy. On the 17 multi-reference systems, Skala's error of 4.61 kcal/mol is higher than B3LYP (5.25 is actually worse — B3LYP is 5.25, M06-2X is 9.70, ωB97M-V is 6.36), but note the ordering: Skala (4.61) < B3LYP (5.25) < ωB97M-V (6.36) < B97M-V (6.95) < M06-2X (9.70). The multi-reference error of 4.61 is consistent with the training data containing only single-reference-labeled molecules (%TAE[(T)] < 6%); Skala does not learn to handle strong static correlation because it has never seen it.

The singlereference subset accuracy of 0.92 kcal/mol is a strong result for a semi-local-cost functional: it exceeds the accuracy of ωB97M-V (1.66) by ~45% on this chemically important class. Since atomization energies form the backbone of thermochemical predictions (reaction energies are linear combinations of atomization energies), this accuracy on TAEs directly translates to accuracy in predicting heats of formation and reaction enthalpies.

Holdout Set Generalization

Headline result (Extended Data Figure 7): On holdout splits of the training data spanning diverse chemical properties, Skala achieves MAEs that are consistently lower than all baseline functionals. On MSR-ACC/TAE25 holdout (730 reactions), Skala's MAE is 0.36 kcal/mol vs. 3.76 (M06-2X, the best baseline) — a 10× reduction. On MB2061 holdout (103 reactions), 3.56 vs. 13.58 (B3LYP). On Water2510 holdout (115 reactions), 0.07 vs. 0.18 (ωB97X-V and ωB97M-V tied). On DES370K holdout (217 reactions), 0.14 vs. 0.23 (both ωB97X-V and ωB97M-V). The error distributions (histograms in Figure 7a) show that Skala's errors are tightly concentrated around zero, with few large outliers, compared to the broader and more skewed distributions of the baselines.

The holdout results confirm two things: (1) Skala is not overfit to its training data — it generalizes to held-out molecules from the same data distribution with errors well below those of traditional functionals, and (2) the improvements on GMTKN55 and W4-17 are not simply a consequence of those benchmarks being "easier" — Skala genuinely reduces error across diverse reaction types.

Computational Cost: Empirical Scaling

Headline result (Figure 4): On GPU (NVIDIA A100), Skala's wall-clock time for XC integration per SCF step scales as $O(N_{\text{orbitals}}^{1.42})$, compared to 1.55 for LDA (SVWN5), 1.54 for revPBE (GGA), 1.34 for r2SCAN (meta-GGA), and 1.80 for B3LYP (hybrid), 1.70 for M06-2X (meta-hybrid). On CPU, Skala scales as $O(N_{\text{orbitals}}^{1.60})$ vs. 2.06 for LDA, 2.06 for revPBE, 1.91 for r2SCAN, 1.90 for B3LYP, 1.88 for M06-2X.

At 2,228 orbitals (113 atoms), GPU timings per SCF step are approximately: Skala ~0.15 s, r2SCAN ~0.12 s, revPBE ~0.08 s, B3LYP ~0.50 s (extrapolated — hybrids could not be run beyond ~2,000 basis functions with GauXC's sn-LinK). At 8,448 orbitals (380 atoms) on GPU: Skala ~0.6 s, r2SCAN ~0.5 s, revPBE ~0.3 s. At 19,816 orbitals (910 atoms) on GPU: Skala ~2.4 s, r2SCAN ~2.1 s, revPBE ~1.2 s. Hybrids could not be run at these system sizes with the available sn-LinK implementation.

The key empirical findings:

• GPU: Skala is within ~25–30% of r2SCAN's cost, and 3–6× cheaper than hybrid functionals at the largest comparable system sizes. The fitted exponent of 1.42 is slightly higher than r2SCAN's 1.34, indicating Skala's cost grows marginally faster with system size, perhaps due to the non-local module's overhead — but both are far below the 1.80 of B3LYP.

• CPU: Skala is approximately 2× slower than r2SCAN for small systems (512 orbitals: ~0.75 s vs. ~0.38 s), but the ratio decreases toward unity as system size increases (the lines converge), because at large sizes the cost is dominated by atomic orbital evaluation (common to all functionals) rather than XC-specific operations. At 4,096 orbitals, Skala is ~1.4× r2SCAN; at 8,448 orbitals, ~1.3×.

• The asymptotic scaling exponents confirm the design goal: the non-local architecture does not introduce higher-order scaling — the exponents (1.42 GPU, 1.60 CPU) are below quadratic and comparable to other meta-GGAs. The prefactor penalty (~30% on GPU, ~30–100% on CPU depending on system size) is the price of learned non-locality. For practical purposes, the 3–6× speedup over hybrid functionals (which would need at least 4× more time at 8,000+ orbitals, if they could even run) more than compensates.

The paper explicitly notes (Section 5, main text, and Section D.7) that hybrid functionals could not be run for molecules with over 2,000 basis functions using GauXC's sn-LinK exact exchange, so the cost advantage over hybrids at large sizes is effectively infinite within that implementation. Production codes with optimized exact exchange (e.g., RI-JK, COSX) would reduce but not eliminate this gap.

Dipole Moments: Self-Consistent Density Quality

Headline result (Table 3b): Skala-1.1 achieves a regularized RMSE of 4.43% on the dipole moment benchmark of Hait & Head-Gordon (2018), outperforming all tested functionals including the best range-separated hybrids: ωB97X-V (5.18%), ωB97M-V (5.84%), B3LYP (7.09%), M06-2X (7.73%). The breakdown by spin state shows Skala performs consistently: non-spin-polarized systems (NSP: 4.23%) and spin-polarized systems (SP: 4.80%).

This result is noteworthy for two reasons: (1) dipole moments are a direct probe of electron density quality — a functional that predicts accurate energies but poor densities will show large dipole errors. Skala's dipole errors being lower than all tested hybrids indicates that its self-consistent densities are genuinely improved, not degraded by error compensation. (2) The dipole benchmark was not directly trained on — Skala was trained on reaction energies, not dipole moments. The improvement in dipoles during SCF fine-tuning (Figure 3a) is an emergent property of the simultaneous energy-and-density improvement the paper documents.

The dipole RMSE improvement over B3LYP (7.09% → 4.43%) is particularly significant because Skala was pre-trained on B3LYP densities — if SCF fine-tuning merely reproduced B3LYP densities, the dipole error would be close to B3LYP's. The fact that Skala substantially outperforms B3LYP on dipoles is direct evidence that SCF fine-tuning moves the densities away from the B3LYP starting point toward a more accurate distribution.

Equilibrium Geometries

Headline result (Table 3c): Skala-1.1 achieves geometry optimization accuracy on par with or exceeding the best hybrid functionals across four independent benchmark sets:

• LMGB35 (light main-group bond lengths): Skala 0.008 Å, vs. ωB97M-V 0.008, r2SCAN 0.006, B3LYP 0.007.
• HMGB11 (heavy main-group bond lengths): Skala 0.009 Å, vs. ωB97M-V 0.010, ωB97X-V 0.040, r2SCAN 0.012.
• CCse21 bond lengths: Skala 0.003 Å, vs. ωB97M-V 0.005, B97M-V 0.005, ωB97X-V 0.005.
• CCse21 bond angles: Skala 0.18°, vs. ωB97M-V 0.18°, ωB97X-V 0.24°, B3LYP 0.38°.
• W4-11-GEOM bond lengths: Skala 0.005 Å, vs. ωB97M-V 0.005, ωB97X-V 0.006, r2SCAN 0.005.
• W4-11-GEOM bond angles: Skala 0.18°, vs. ωB97M-V 0.25°, ωB97X-V 0.32°, B3LYP 0.37°.

The box plots in Table 3c show the quartiles of the error distribution, revealing that Skala's error distribution is not just shifted to lower mean values but also has smaller spread (tighter quartiles) on several sets, particularly CCse21 bond lengths and W4-11-GEOM bond angles.

The GEO metric (Supplementary Figure 15a) quantifies the energy penalty from using self-optimized rather than reference geometries: Skala's GEO on W4-11-GEOM is ~0.02 kcal/mol on average, comparable to ωB97M-V (~0.01) and ωB97X-V (~0.02), and far lower than revPBE (~0.35) or r2SCAN (~0.10). This means that when Skala is used both to optimize the structure and compute the energy — the standard use case in computational chemistry — it introduces negligible additional error relative to evaluating on reference geometries.

The distortion training data ablation (Extended Data Figure 8c and Supplementary Figure 15b) is critical context: without the MSR-ACC/Distortion dataset (~19,000 distorted geometries), Skala's geometry optimization accuracy degrades to approximately revPBE levels (bond length errors roughly double). Training on distorted structures along vibrational modes effectively teaches the model the shape of the potential energy surface near equilibrium, which it needs to correctly locate the gradient zero during geometry relaxation. The fact that this targeted data augmentation is the difference between GGA-level and hybrid-level geometry performance illustrates both the power and the limitation of the data-driven approach: for properties not well-represented in the base training data, accuracy degrades gracefully (defaulting to semi-local quality) rather than catastrophically.

SCF Convergence Behavior

Headline result (Tables 3 and 4): Skala-1.1 converges in the first SCF attempt for 98.87% of the 2,734 evaluated structures (W4-17 + GMTKN55 + dipoles), compared to 99.93–100.00% for traditional functionals. The additional convergence attempts (damping, level shifting, Newton solver, manual intervention) resolve all but 0.18% of cases. The average number of SCF iterations in the first attempt is 9.2 for Skala vs. 8.0–8.9 for traditional functionals — a small but measurable increase (~0.3–1.0 extra iterations on average) consistent with the slightly less smooth behavior expected from a neural functional.

The slightly lower first-attempt convergence rate and slightly higher iteration count are the practical cost of the neural functional's non-smoothness relative to analytic forms. However, the paper's retry protocol (six fallback strategies, detailed in Section D.4) handles these cases robustly, and the overall convergence behavior is within acceptable ranges for production use. Using a def2-TZVP basis set (smaller than the default QZVP) improves first-attempt convergence to 99.38% and reduces average iterations to 8.9 (Table 4), suggesting that basis-set size affects convergence more for Skala than for traditional functionals.

Transition Metals: Preliminary Generalization

Headline result (Extended Data Figure 10): On transition-metal benchmarks where Skala has seen minimal or no training data:

• 3d4dIPSS (ionization potentials and spin splittings of 3d/4d atoms, seen in training): Skala 7.19 kcal/mol vs. ωB97M-V 6.37, B3LYP 7.69, r2SCAN 12.52.
• TMD10 (dimer bond energies, seen in training): Skala 9.34 vs. B97M-V 4.31, ωB97X-V 5.18, ωB97M-V 7.50.
• DAPD (Pd-containing diatomics, seen in training): Skala 1.45 vs. B3LYP 3.54, B97M-V 6.67, ωB97M-V 7.85 — a strong result suggesting good transfer to 4d metals.
• MOR41 (closed-shell metal-organic reactions, unseen): Skala 3.29 vs. ωB97M-V 2.13, ωB97X-V 2.18, r2SCAN 3.22.
• MOBH35 (metal-organic barrier heights, unseen): Skala 2.08 vs. ωB97M-V 1.84, B97M-V 1.82, ωB97X-V 2.29.
• 3dTMV (3d transition metal electrocatalysts, unseen): Skala 8.01 vs. B97M-V 8.00, B3LYP 7.34, ωB97M-V 13.93.
• ROST59 (open-shell transition metal reactions, unseen): Skala 3.38 vs. ωB97X-V 2.67, r2SCAN 3.20, ωB97M-V 2.72.

The pattern is clear: on benchmarks where transition metals were included in training (even minimally — atomic properties and dimer bonds), Skala performs competitively; on unseen transition-metal benchmarks, its performance is mid-pack, roughly at the r2SCAN/B3LYP level. This is consistent with the paper's thesis that Skala's accuracy is data-driven: without sufficient transition-metal training data, it defaults to semi-local quality. As a preview of future scaling, these results indicate that expanding the transition-metal training data is the primary bottleneck to extending Skala beyond main-group chemistry.

---

Ablation Studies and Robustness Checks

Non-local branch contribution (Figure 2a, left): Purely local models (non-local branch ablated, identical architecture otherwise, trained on the full dataset in Extended Data Table 6) achieve an average WTMAD-2 of ~5.6 kcal/mol on Diet GMTKN55 (averaged over 5 seeds), compared to ~2.8 kcal/mol for the full non-local Skala architecture — a ~50% error reduction. The local-only models show greater inter-seed variance (two seeds overlap nearly exactly, but the spread across seeds is wider than for the non-local models), indicating that non-local representations not only improve accuracy but also make training more stable across random initializations.

This ablation is the key architectural justification for the paper's central claim that learned non-locality substitutes for exact exchange. A local model with the same input features, same training data, and comparable parameter count achieves only meta-GGA-level accuracy. The non-local module is what pushes it into hybrid territory.

Data composition ablation (Figure 2a, right, and Extended Data Table 6): Training on progressively more diverse chemistry:

• Set A only (thermochemistry, basic properties, distorted geometries, transition metals — ~284,000 reactions): average WTMAD-2 on Diet GMTKN55 is ~8.4 kcal/mol (GGA-level; large seed variance). One seed reaches ~6.5, indicating potential for better results with more tuning, but the inter-seed spread is wide.
• Sets A + B (adds reactions and barrier heights — ~398,000 reactions): average WTMAD-2 drops to ~6.6 kcal/mol (approximately B97M-V level) with reduced seed variance.
• Sets A + B + C (adds NCIs and conformers — full ~409,000 reactions, i.e., Skala-1.1): average WTMAD-2 drops to ~2.8 kcal/mol with the lowest seed variance.

The monotonic improvement with data diversity — both in average error and in inter-seed consistency — is the paper's strongest evidence for the "systematically improvable" claim. Each new category of chemistry (barrier heights, non-covalent interactions) provides the model with information it cannot extract from thermochemistry alone, and the error keeps dropping. This resembles the scaling behavior of deep learning in other domains (vision, language) where performance improves with dataset size and diversity, motivating the versioned-release model (Skala-1.0, -1.1, etc.).

Distortion data for geometry optimization (Extended Data Figure 8c and Supplementary Figure 15): Without the MSR-ACC/Distortion dataset (~19,000 distorted geometries), Skala's Wiggle150 MAE for strained conformer energies degrades from ~0.73 to ~1.0–1.6 kcal/mol (depending on seed). The GEO metric on W4-11-GEOM degrades from ~0.02 to ~0.2–1.0 kcal/mol, and the bond length/angle errors roughly double. This ablation confirms that predicting accurate geometries requires explicit training on the shape of the potential energy surface near minima — equilibrium structures alone are insufficient.

Physical constraint emergence (Figure 2b and Supplementary Figure 14): Training on MSR-ACC/TAE25 only (~730 atomization energies of small molecules) produces a model with consistently negative $T_c$ values (violating the positivity constraint) for atoms H, B, C, N, O, F, Ne, Na, Mg, Al, Si, P, S, Cl, Ar — essentially every element tested. Training on the full dataset produces positive $T_c$ for all atoms, with values close to reference data from reverse-engineered CCSD densities. This pattern is consistent across five random seeds (Figure 14): the full-training-set models all produce positive $T_c$, while the TAE25-only models all violate the constraint. The paper interprets this as emergent constraint satisfaction from data diversity — an intriguing finding that motivates the minimal-constraint design philosophy.

Grid convergence (Supplementary Figure 13): Skala's MAE on Diet GMTKN55 and W4-17 changes by ≤0.035 kcal/mol when moving from grid level 3 to level 6 (a ~2× increase in radial grid points and ~4× in angular points), comparable to the grid sensitivity of traditional functionals (revPBE: ~0.02, ωB97M-V: ~0.01). This confirms that the architecture has a well-defined continuum limit and does not overfit to the specific grid resolution used in training.

Spin-symmetry breaking for multi-reference systems (Supplementary Table 5h and Figure 12): For the 17 multi-reference molecules in W4-17, allowing spin-symmetry breaking by adding random noise to the initial density matrix reduces Skala's TAE errors for 1BN, B2, and C2 — the three molecules where spin-symmetry breaking has the largest effect across all functionals. For most other multi-reference molecules, spin-symmetry breaking has minimal effect on Skala's energies. This indicates that Skala behaves similarly to traditional functionals in its response to symmetry breaking, and that the remaining multi-reference errors (~4.6 kcal/mol) are due to missing physics (static correlation) rather than an inability to access symmetry-broken solutions.

Specific challenging subsets (Table 6 in Supplementary): The reactions in GMTKN55 where Skala's error differs most from ωB97M-V (both directions) are dominated by MB16-43 (mindless molecules) and SIE4x4 (self-interaction error). Skala dramatically outperforms ωB97M-V on many MB16-43 reactions (e.g., 4.30 vs. 56.33 kcal/mol, 0.42 vs. 45.23 kcal/mol), suggesting its learned non-locality handles unusual bonding much better than fixed-fraction exact exchange. Conversely, ωB97M-V outperforms Skala most on SIE4x4 reactions (e.g., -67.61 vs. -7.25 kcal/mol for the notorious $2\text{H}_6\text{N}_2^+ \to \text{H}_3\text{N} + 2\text{H}_3\text{N}^+$ reaction), indicating that self-interaction correction remains a relative weakness for Skala — consistent with its lack of explicit self-interaction error mitigation and its reliance on D3 dispersion rather than wavefunction-based correlation.

---

Critical Assessment

Claim 1: Skala surpasses state-of-the-art hybrid functionals in accuracy across GMTKN55 while retaining semi-local computational cost. This claim is supported with notable caveats. The WTMAD-2 of 2.80 kcal/mol (Skala) vs. 3.23 (ωB97M-V) represents a genuine improvement, and the per-subset breakdown (Table 1) shows Skala dominating across thermochemistry and kinetics. However, the "semi-local cost" claim relies primarily on the empirical timings in Figure 4, which show Skala to be 25–30% slower than r2SCAN on GPU and up to 2× slower on CPU for small systems. This is indeed substantially cheaper than hybrids (3–6×), but it is not at parity with the fastest meta-GGAs — there is a real prefactor cost to the non-local architecture. Whether Skala is considered "semi-local cost" depends on whether one groups it with r2SCAN (1.34–1.91 exponent, 0.12–0.38 s at 2,228 orbitals) or with B3LYP (1.80–1.90 exponent, ~0.50 s). The paper's framing as "semi-local cost" is defensible if the criterion is asymptotic scaling (exponents of 1.42 and 1.60 are solidly in the O(N³) regime), but the prefactor penalty is real and should be acknowledged. Furthermore, the hybrid timings are based on GauXC's sn-LinK, which is not the fastest exact-exchange implementation available — production codes such as ORCA with RIJCOSX or Q-Chem with RI-J would show smaller cost differences. The paper would be strengthened by reporting timings from at least one optimized production code to bound the practical cost advantage.

Claim 2: The non-local architecture is responsible for the accuracy gains, not just the large dataset. This claim is strongly supported by the model ablation in Figure 2a (left). The purely local variant trained on the full dataset achieves ~5.6 kcal/mol WTMAD-2 on Diet GMTKN55 vs. ~2.8 for the non-local model — a 2× difference. This is a clean ablation: same data, same local architecture capacity (331,265 parameters in the local branch), same training protocol, and the only difference is the ~54,000 non-local parameters. The improvement is large and consistent across seeds. The local model's performance (~5.6 kcal/mol) is consistent with prior ML functionals (DM21 achieved 3.97 on full GMTKN55, but with exact-exchange-like features) — confirming that without learned non-locality, one cannot escape the semi-local accuracy ceiling regardless of dataset size.

Claim 3: Skala's accuracy improves systematically with training data. This claim is supported with qualifications. The data composition ablation (Figure 2a, right) shows monotonic improvement as more chemistry is added to training. However, the ablation uses three coarse data groupings (A, A+B, A+B+C), not a truly systematic scaling study (e.g., 10%, 25%, 50%, 100% of data). The claim that Skala is "systematically improvable" — that more data will continue to improve accuracy — rests on extrapolation from three data points plus the general deep-learning scaling evidence from other domains. The paper would be substantially strengthened by showing learning curves (accuracy vs. dataset size) on Diet GMTKN55 for multiple dataset sizes within a fixed data mixture, which would provide direct evidence for the scaling trend rather than the categorical ablation presented. The emergence of the $T_c$ positivity constraint (Figure 2b) provides supporting evidence that data diversity, not just size, matters — but this is a single constraint, and it is unclear whether other exact constraints would similarly emerge with sufficient data or would require explicit architectural enforcement.

Claim 4: Skala achieves hybrid-level accuracy at semi-local cost, breaking the Jacob's ladder trade-off. This is the paper's headline claim and it is supported but requires careful boundary-drawing. The evidence that Skala breaks the trade-off on GMTKN55 is compelling: it achieves lower WTMAD-2 than any hybrid at lower computational cost. However, the "trade-off" is defined by asymptotic scaling plus practical runtime, and Skala does introduce a new cost — training. The training cost (1 million pre-training steps on 8×A100, plus 20,000 fine-tuning steps) is a one-time investment amortized over inference use, but it is a real cost that traditional functional development does not incur in the same way (traditional functionals require human effort for design and fitting, not GPU-hours). For a practical user, the relevant cost is inference time, and on that metric Skala is cheaper than hybrids. The broader claim that Skala "breaks" the ladder should be contextualized: it breaks it for main-group chemistry with abundant training data. For transition metals (Extended Data Figure 10), Skala's accuracy on unseen benchmarks degrades to r2SCAN/B3LYP levels — the trade-off reappears where data are scarce.

Genuine weaknesses in the experimental design:

• Single model family: all results are for Skala-1.1 trained on MSR-ACC data. There is no comparison to other neural XC functionals (DM21, CIDER, etc.) on the same benchmarks, which would help disentangle architecture contributions from data contributions. The paper compares only to traditional analytic functionals.

• Dispersion correction is fixed, not learned: Skala uses B3LYP's D3(BJ) correction parameters — the same dispersion model used during pre-training on B3LYP densities. This means the non-local module is not learning long-range dispersion; it is learning everything else. The strong performance on IDISP (intramolecular dispersion, 0.60 vs. ωB97M-V's 1.63) suggests the non-local kernel captures some medium-range correlation, but the architecture likely benefits from the D3 correction for long-range van der Waals. The paper acknowledges this as future work but it complicates the claim that the architecture "learns non-locality" — it learns some non-locality, with the tail handled by a classical correction.

• Difficulty estimation / density dependence not ablated: the model is trained on B3LYP densities and fine-tuned on self-consistent densities, but there is no ablation showing how sensitive the final accuracy is to the choice of pre-training density functional. Could Skala have been trained on PBE densities instead and achieved the same final accuracy after SCF fine-tuning? This matters for transferability: if Skala depends on B3LYP-quality starting densities during pre-training, it may not generalize to systems where B3LYP densities are poor (strongly correlated systems, excited states).

• The ~400,000 training reactions are heavily dominated by a few datasets: ~120,000 TAEs, ~154,000 reaction paths, ~41,000 conformers. The adaptive sampling mitigates imbalance but the total information content is skewed toward thermochemistry. The holdout results (Figure 7) show excellent generalization on holdout splits of the same data distributions, but this does not test generalization to chemistry completely absent from training — e.g., heavier elements, excited states, periodic systems.

• SCF convergence issues are real but manageable: the 98.87% first-attempt convergence rate and the need for a six-level retry protocol (Tables 3 and 4) indicate that Skala is less robust than traditional functionals in SCF calculations. The paper is transparent about this, providing detailed convergence statistics, but a practical user switching from B3LYP to Skala will encounter more frequent convergence failures, particularly with larger basis sets. The paper would benefit from case studies showing what types of systems cause SCF difficulties and whether simple workarounds (e.g., starting from a converged B3LYP density) alleviate them.

Missing experiments that would have strengthened the paper:

• Learning curves for dataset size at fixed data mixture, to substantiate the "systematically improvable" claim beyond the coarse A/B/C ablation.
• Training data efficiency comparison: how much data does Skala need to match ωB97M-V? The ablation shows that Set A (~284,000 reactions, GGA-level) is insufficient and Set A+B+C (~409,000) surpasses hybrids. The crossover point (how many reactions needed to match the best hybrid) is not identified.
• Generalization to solids / periodic systems: all benchmarks are molecular. DFT is widely used for solids (crystals, surfaces). Can Skala, trained only on molecules, generalize to periodic boundary conditions? This requires a qualitatively different test.
• Direct comparison to other ML functionals: evaluating DM21 or CIDER on the same benchmarks with the same basis sets and protocols would clarify whether Skala's gains come from the architecture, the data, or the training protocol.
• Sensitivity to basis set: all evaluations use def2-QZVP; results with smaller basis sets (def2-TZVP, def2-SVP) would inform practical users who cannot afford QZVP for large systems.
• Calibration of error estimates: Skala produces an energy prediction but no uncertainty estimate. Traditional functionals have known error patterns; users have intuition for when to trust them. Skala's error characteristics on out-of-distribution chemistry are not characterized beyond the transition-metal preview.

Despite these limitations, the experimental evidence for Skala's main claim — that learned non-locality from large-scale data can produce a functional that competes with or surpasses hybrid functionals on main-group chemistry at lower computational cost — is robust. The GMTKN55 results (Table 1), W4-17 single-reference accuracy of 0.92 kcal/mol, dipole moments (Table 3b), geometry optimizations (Table 3c), and computational cost measurements (Figure 4) together form a consistent picture: Skala is a viable, sometimes superior, alternative to range-separated hybrids for main-group chemical applications. The paper is transparent about the boundaries of this conclusion (multi-reference systems, transition metals, the dependence on D3 dispersion), and the versioned-release model (Skala-1.0 → Skala-1.1) with documented improvements provides a framework for tracking progress as these boundaries are pushed outward by expanded training data.

6. Limitations and Trade-offs

The Difficulty Estimation Problem: Training Data Generation at Scale Requires Wavefunction Methods That Are Themselves the Bottleneck

The assumption or constraint. Skala's training depends on ~400,000 CCSD(T)/CBS energy differences spanning diverse chemistry — a dataset whose generation required an industrial-scale campaign using the W1-F12 and W1w thermochemical protocols (Section 2.3, Section 7.2). Generating these reference labels required running high-level wavefunction calculations whose cost scales as O(N⁷) with system size. The paper explicitly acknowledges that extending this approach to multi-reference systems is blocked by a lack of scalable reference methods:

"A key challenge in this endeavor will be extending coverage to multi-reference and strongly correlated systems, where generating accurate reference data at scale remains an obstacle that will require new scientific and computational advances to overcome." (Section 6)

The consequence. The very bottleneck that Skala is designed to circumvent at inference time — the high computational cost of accurate quantum chemistry — remains the bottleneck for creating Skala in the first place. This creates a fundamental asymmetry: Skala can replace expensive hybrid DFT calculations (O(N⁴) inference cost) for systems within its training distribution, but it cannot extend accuracy to chemical regimes where generating CCSD(T)/CBS training labels is infeasible. Specifically, multi-reference systems (where a single Slater determinant is a poor approximation — common in transition metals, bond-breaking, excited states) are largely absent from the training data because the W1-F12/W1w protocols themselves break down. The 17 multi-reference molecules in W4-17 show Skala's MAE of 4.61 kcal/mol (Extended Data Figure 8a), compared to 0.92 kcal/mol on the 183 single-reference molecules — a 5× degradation that the paper attributes directly to missing training data. For systems with strong static correlation (e.g., many transition-metal complexes, dissociating bonds, biradicals), Skala cannot be trained because reference labels at CCSD(T)/CBS accuracy are unreliable or unavailable, creating a hard ceiling on where the data-driven paradigm can reach.

What evidence exists in the paper. The W4-17 breakdown (Extended Data Figure 8a) quantifies the single-reference vs. multi-reference gap. The transition-metal benchmarks (Extended Data Figure 10) show that on unseen TM datasets (MOR41, MOBH35, 3dTMV, ROST59), Skala performs at roughly r2SCAN/B3LYP levels — consistent with defaulting to semi-local accuracy when training data are absent. The paper states that training data are excluded if %TAE[(T)] exceeds 6% (Section C.1), which explicitly filters out multi-reference character that CCSD(T) cannot reliably handle. The training data composition (Extended Data Table 6) shows transition metals contributing only 35 reactions total — negligible compared to ~118,000 for main-group atomization energies.

Mitigation status. The paper partially acknowledges this (Section 6), framing it as a future direction contingent on "new scientific and computational advances" in multi-reference electronic structure. No architectural mitigation is attempted — the model does not flag when it encounters multi-reference character, and there is no fallback mechanism to a more appropriate wavefunction method. The paper suggests the versioned-release model (Skala-1.0 → Skala-1.1) will incorporate progress as reference data improve, but this passes the bottleneck to the electronic structure theory community rather than solving it within the ML framework. A practitioner evaluating a molecule with suspected multi-reference character (e.g., a transition-metal catalyst with a low-lying excited state) has no way to know whether Skala's prediction is reliable or falls into the 4.61 kcal/mol error regime.

---

Difficulty Estimation Cost: the Pre-Training Density Functional Is a Frozen, Unoptimized Choice With Unquantified Sensitivity

The assumption or constraint. Skala is pre-trained on B3LYP densities (Section 2.1, Section B.1): the total energy is computed by taking the B3LYP total energy, subtracting the B3LYP XC contribution, and adding Skala's own XC energy evaluated on the B3LYP density. The choice of B3LYP as the pre-training density functional is pragmatic — "B3LYP densities are of reasonable quality" and "can be generated at scale for our entire dataset" (Section 2.1) — but there is no ablation studying alternative choices (PBE, r2SCAN, or even a Hartree–Fock density). The SCF fine-tuning phase is then tasked with closing the gap between accuracy on B3LYP densities and accuracy on self-consistent Skala densities, but the pre-training distribution fundamentally shapes what the model learns.

The consequence. If B3LYP densities are systematically poor for certain bonding situations — and they are known to be, particularly for charge-transfer complexes (due to delocalization error in B3LYP's 20% exact exchange), anions (where B3LYP over-delocalizes the extra electron), and systems with significant static correlation — then Skala may inherit biases from the pre-training density even after SCF fine-tuning. The fine-tuning protocol (20,000 steps with batch size 1, Section B.5) is brief compared to the 1,000,000-step pre-training, raising the question of whether the model has capacity to fully escape biases in the pre-training distribution. More subtly, the entire training objective (reaction energy regression on B3LYP densities) may encourage the model to learn a functional that reproduces B3LYP's density-driven error cancellation patterns — patterns that are correct for B3LYP densities but may not translate to self-consistent accuracy in a way that the short fine-tuning phase can fully correct.

What evidence exists in the paper. The dipole moments provide indirect evidence: after SCF fine-tuning, Skala's dipole RMSE (4.43%, Table 3b) is substantially better than B3LYP's (7.09%), indicating that the self-consistent densities have moved away from B3LYP densities. The simultaneous improvement of energies and dipoles during fine-tuning (Figure 3a) suggests the model is not locked into B3LYP patterns. However, there is no direct ablation comparing Skala trained on B3LYP densities vs. trained on, say, r2SCAN or Hartree–Fock densities. The paper does not report how much fine-tuning is needed to erase density-functional-dependence from pre-training — whether 20,000 steps are enough or whether 100,000 would yield further improvement. The GEO metric (~0.02 kcal/mol, Supplementary Figure 15a) shows that Skala's potential energy surface around equilibrium is accurate, but this tests the functional's behavior near B3LYP geometries (which are themselves reasonable), not in regions where B3LYP densities are known to be poor.

Mitigation status. The authors acknowledge the dependency implicitly by providing SCF fine-tuning results (Figure 3a) that monitor the energy accuracy gap between B3LYP and self-consistent densities, showing it narrows but does not close completely. They do not claim that the choice of pre-training functional is irrelevant. The paper does not explore whether the model could be fine-tuned for longer, or whether a second fine-tuning phase on PBE densities could further improve robustness. For a practitioner, the takeaway is that Skala's accuracy is most reliable for chemistry similar to what B3LYP describes well — which is broad but not universal — and there is no quantification of how performance degrades on systems where B3LYP densities are known to be problematic (e.g., charge-transfer excitation energies, although these are not directly an XC functional property, the density quality matters).

---

Difficulty Estimation Cost: Training Data Generation Dominates the Total Compute Budget and Is Not Amortized

The assumption or constraint. The paper's headline result — Skala achieves hybrid-level accuracy at semi-local cost — measures only inference-time computational cost. The one-time cost of generating the ~400,000 CCSD(T)/CBS reference energies, computing B3LYP densities for every training molecule, and running 1,000,000 gradient steps on 8 A100 GPUs is excluded from the cost comparison against traditional functionals. The training data generation, in particular, involved running the W1-F12 and W1w protocols (which require multiple CCSD(T) calculations with basis-set extrapolation) on molecules with up to 9 non-hydrogen atoms (MSR-ACC/TAE) and up to 13 non-hydrogen atoms (MSR-ACC/Reactions). The paper states this was done at "an unprecedented volume" (Section 2.3) but does not provide the total core-hours or GPU-hours consumed to generate the reference data.

The consequence. The cost of generating training data is relevant for two reasons. First, it determines the barrier to entry for extending Skala to new chemical domains: if a research group wants to train a version of Skala for, say, f-block chemistry or excited states, they must replicate the data generation campaign — which may be infeasible without access to the computing resources of a company like Microsoft. This makes Skala less "open" than traditional functionals (which require no training data) and more "open" than the specific model checkpoint (which is released under MIT license). Second, the data generation cost factors into the total cost of ownership: if Skala-1.1 required 10 million CPU-hours of CCSD(T) calculations to train, then the comparison against hybrid functionals (which required zero training GPU-hours) is an apples-to-oranges comparison for total resource consumption, even if inference is cheaper.

What evidence exists in the paper. The paper does not report the total computational cost of generating MSR-ACC. The supplementary information (Section C) provides the protocols (basis sets, extrapolation formulas, software packages) but not the total core-hours or GPU-hours. The training cost is partially reported: 1,000,000 pre-training steps on 8 A100 GPUs (Section B.4), which at ~0.5 seconds per step gives roughly 140 GPU-days. The SCF fine-tuning adds ~20,000 steps on 1 A100 GPU (~6 GPU-days). But these numbers exclude the cost of computing the reference CCSD(T) energies and the B3LYP densities for all training molecules — by far the dominant contribution. A rough order-of-magnitude estimate: W1-F12 for a molecule with 5 non-hydrogen atoms requires roughly ~100–500 CPU-hours; for 120,000 such calculations, this is 12–60 million CPU-hours, orders of magnitude more than the training cost.

Mitigation status. The paper acknowledges that data generation is a major undertaking — it is described as a "large-scale effort" (Section 2.3) and the dataset is given a name (MSR-ACC) with a dedicated publication (Ehlert et al., 2025). The versioned-release model (Skala-1.0 → Skala-1.1) amortizes the data generation cost over successive model improvements. The release of the MSR-ACC/TAE25 subset (molecules with up to 5 non-hydrogen atoms and a single molecular fragment) under a public license partially addresses the barrier to entry by providing a benchmark-quality subset for reproducible research. However, the full training data (MSR-ACC/TAE, /Reactions, /Conf, /NCI, /Distortion, etc.) are not released, meaning the broader community cannot reproduce Skala-1.1 or extend it without independently generating equivalent data. For a practitioner evaluating whether to invest in Skala for their application, the total cost should include the inference cost savings over hybrid functionals amortized over the expected number of calculations — a calculation the paper does not perform.

---

Difficulty Estimation Cost: Convergence Robustness — Skala Is More Difficult to Converge Self-Consistently Than Traditional Functionals

The assumption or constraint. The paper assumes that Skala can be deployed as a drop-in replacement for traditional XC functionals in standard SCF workflows. However, neural network functionals are inherently less smooth in their density dependence than analytic functional forms, which manifests as increased difficulty in SCF convergence. The paper's own statistics (Tables 3 and 4) show that Skala requires a six-level retry protocol (DIIS → damping → level shifts at 0.1, 0.3, 0.5 → Newton solver → manual intervention) to converge 99.82% of 2,734 evaluated structures, compared to first-attempt convergence rates of 99.93–100% for traditional functionals.

The consequence. In practical deployment, Skala's lower first-attempt convergence rate (98.87% vs. 99.93–100%) means a small but measurable fraction of calculations will require user intervention or automated fallback logic — costing both compute time and human effort. The average SCF iteration count (9.2 for Skala vs. 8.0–8.9 for traditional functionals, Table 4) adds an additional ~5–15% overhead to each converged calculation. The retry protocol itself introduces complexity: a production workflow that calls Skala through PySCF or GauXC must handle cases where the first SCF attempt fails and re-attempt with different settings, which may interact badly with automated high-throughput screening pipelines. Furthermore, the convergence statistics are measured for Skala-1.1 with the specific SCF settings described in Section D.4 (DIIS window size 8, convergence thresholds 5×10⁻⁶ Eₕ for energy and 10⁻³ Eₕ for orbital gradient). Different codes, basis sets, or convergence thresholds may produce different — potentially worse — convergence behavior. The paper notes that using a smaller basis set (def2-TZVP) improves first-attempt convergence to 99.38% and reduces average iterations to 8.9 (Table 4), suggesting that the convergence difficulty is basis-set-dependent, but this interaction is not systematically explored.

What evidence exists in the paper. Tables 3 and 4 provide detailed convergence statistics across 2,734 structures (W4-17: 211 structures, GMTKN55: 2,372 structures, dipoles: 151 structures). Table 3 breaks down which retry level was ultimately successful for the small fraction that did not converge in the first attempt: 0.18% required damping, 0.33% required level shift 0.1, 0.22% required level shift 0.3, 0.22% required Newton solver, and 0.18% required manual convergence. Table 4 reports average SCF iterations for the first attempt only. The paper does not report convergence behavior on larger molecules (the timing data in Figure 4 includes molecules up to 910 atoms, but convergence statistics for those are absent) or on transition-metal systems (where SCF convergence is notoriously more difficult even for traditional functionals due to small HOMO-LUMO gaps).

Mitigation status. The paper is transparent about convergence statistics and provides the retry protocol in full detail (Section D.4), which is helpful for reproducibility. The orbital gradient descent fallback (Section D.5) provides a robust last-resort converger, but it requires double-precision computation and is not integrated into standard DFT codes — it was used manually for the 0.18% of cases. The paper suggests that "slightly less smooth behavior than traditional functionals" is "expected for a deep-learned functional" and characterizes all variations as "within acceptable ranges." For a practitioner, the additional convergence friction is manageable but real: it adds complexity to automated workflows and may cause failures on systems not represented in the convergence statistics (large molecules, transition metals, open-shell species). The paper does not provide guidance on which types of systems are most likely to cause convergence difficulties or how to adjust SCF settings proactively.

---

Difficulty Estimation Cost: The Training Data Distribution Skews Heavily Toward Equilibrium, Single-Reference, Main-Group Chemistry — Creating a Sharp Accuracy Cliff Outside This Regime

The assumption or constraint. Skala's training data are dominated by molecules and reactions from main-group chemistry at or near equilibrium geometries, with single-reference electronic structure character. The largest datasets are total atomization energies (~118,000), reaction paths (~114,000), and conformers (~41,000) — all computed for molecules with up to 13 non-hydrogen atoms from elements H, Li–F, Na–Cl (Extended Data Table 6). The paper explicitly excludes systems with %TAE[(T)] > 6% (multi-reference character) from training. Even within main-group chemistry, the training distribution is concentrated around equilibrium structures and low-energy conformers (within 10 kcal/mol of the minimum).

The consequence. While Skala generalizes well within main-group chemistry as measured by GMTKN55 and W4-17, the accuracy can degrade sharply for chemistry that falls outside the training distribution's support. The paper provides direct evidence for three such regimes:

• Multi-reference systems: The 17 multi-reference molecules in W4-17 show a 5× error increase (4.61 vs. 0.92 kcal/mol, Extended Data Figure 8a).
• Highly strained conformations: On Wiggle150 (strained conformers of drug-sized molecules), Skala achieves 0.73 kcal/mol MAE — competitive but not dominant (ωB97M-V achieves 0.87, Extended Data Figure 8b). Importantly, the distortion data ablation (Extended Data Figure 8c) shows that without explicitly training on distorted geometries, Wiggle150 errors degrade to ~1.0–1.6 kcal/mol, indicating that even equilibrium main-group training data are insufficient for predicting energies away from minima.
• Large conjugated systems: On C60ISO (relative energies of C60 isomers), Skala's error of 8.16 kcal/mol is far higher than B3LYP (2.40) or B97M-V (4.85) — Table 1. These systems (60 carbon atoms, extensive π-conjugation, multiple low-energy isomers) are far outside the training distribution in both size and electronic structure, and Skala's accuracy degrades to well below hybrid quality.
• Transition metals (unseen): On MOR41, MOBH35, 3dTMV, and ROST59 — datasets with no training overlap — Skala performs at roughly r2SCAN/B3LYP level (Extended Data Figure 10), not at hybrid level. The transition-metal accuracy is "on par with established functionals" but not state-of-the-art, reflecting the near-absence of transition metals in training (35 reactions total).

More broadly, all training molecules contain 13 or fewer non-hydrogen atoms. For larger molecules, Skala is extrapolating — the functional should be size-consistent by construction (the integral form ensures this), but the learned non-local kernel is trained on atomic environments characteristic of small molecules and may misrepresent the electronic structure of large, delocalized systems. The performance on C60ISO is a warning sign: Skala does not gracefully handle large conjugated systems, possibly because its non-local interactions are limited to within-atom coarse-point communication (there is no cross-atom message passing), and the learned effective non-local exchange fraction may not be appropriate for extensive π-systems.

What evidence exists in the paper. The evidence is distributed across multiple results: the single-reference/multi-reference breakdown in Extended Data Figure 8a, the transition-metal benchmarks in Extended Data Figure 10, the C60ISO subset of GMTKN55 (Table 1), and the Wiggle150 results (Extended Data Figure 8b). The data composition ablation (Figure 2a, right) provides indirect evidence: adding barrier heights and non-covalent interactions to the training data improves performance on those specific categories, suggesting that the model does not extrapolate from atomization energies alone — it needs explicit examples of each chemical phenomenon. The paper does not provide a systematic study of accuracy vs. molecular size (e.g., plotting MAE as a function of number of atoms for GMTKN55 subsets, which range from diatomics to C60), which would quantify the extrapolation behavior.

Mitigation status. The paper positions this limitation explicitly as the motivation for the versioned-release model: Skala-1.0 → Skala-1.1 improved accuracy by expanding the training data, and future versions will incorporate more chemistry. The transition-metal preview (35 reactions included intentionally, Section 6) demonstrates "strong transferability" on seen benchmarks (DAPD: 1.45 kcal/mol, the best result by far) and "par with established functionals" on unseen ones — a deliberate design to show that even minimal training signal carries some transfer. For multi-reference systems, the paper is honest that this "will require new scientific and computational advances to overcome." However, for a practitioner, the current version (Skala-1.1) has a well-defined scope: main-group chemistry with up to ~13 non-hydrogen atoms, single-reference electronic structure, near-equilibrium geometries. Outside this scope, accuracy degrades — sometimes to meta-GGA levels (transition metals, large conjugated systems), sometimes worse (multi-reference). The paper does not provide a diagnostic for determining whether a given molecule is within scope, leaving users to rely on chemical intuition or benchmarking against a trusted hybrid functional.

---

Difficulty Estimation Cost: the Architecture's Non-Locality Is Constrained to Within-Atom Communication — Long-Range Correlation Beyond Dispersion Is Not Learned

The assumption or constraint. The non-local module in Skala enables grid points to communicate only with the coarse point of their own associated atom (Section A.4, Section A.5). There is no message passing between coarse points of different atoms. The paper justifies this by noting that the density features themselves depend on the full molecular density through the SCF solution, so cross-atom information is implicitly encoded in the semi-local features. Additionally, the D3 dispersion correction (using B3LYP parameters, Section 2.2) handles the longest-range van der Waals interactions. The non-local module is designed to learn the remaining non-local physics — predominantly shorter-range exchange and correlation effects that hybrid functionals traditionally capture with exact exchange.

The consequence. This design choice means the non-local module is not learning long-range correlation beyond what D3 provides, nor is it learning cross-atom non-local exchange of the type that would be captured by an all-to-all interaction kernel. The model cannot, for example, learn that the exchange interaction between two electrons on atom A and atom B depends on the specific bonding environment of both atoms in a way that goes beyond what the density features at atom A's grid points already encode about atom B. This may explain why Skala underperforms the best hybrids on certain non-covalent interaction subsets: CHB6 (cation-neutral dimers, 1.29 kcal/mol vs. 0.96 for ωB97M-V), WATER27 (1.25 vs. 0.77 for B97M-V), and AHB21 (0.41 vs. 0.28 for ωB97M-V) — see Table 1. These involve interactions where the detailed electronic structure of the monomers is perturbed by the intermolecular interaction in ways that depend non-locally on the partner's electron distribution — a cross-atom effect that the within-atom non-local architecture would struggle to capture. In contrast, the architecture excels at describing non-local effects within a single atom's environment, which is precisely what matters for covalent bond energies — and indeed Skala dominates thermochemistry and kinetics subsets.

More speculatively, the within-atom constraint may limit Skala's ability to describe charge-transfer excitations, where an electron moves from one atom to another and the exchange interaction between the donor and acceptor orbitals is inherently cross-atom. Charge-transfer energetics are a known weakness of semi-local functionals and a primary motivation for range-separated hybrids; Skala's architecture, lacking explicit cross-atom non-locality and relying on D3 for long-range correlation, may inherit some of these limitations despite its strong performance on ground-state thermochemistry.

What evidence exists in the paper. The per-subset GMTKN55 results (Table 1) indirectly reveal this limitation: Skala underperforms the best hybrids on several NCI subsets despite excelling on thermochemistry. The molecular "thermometer" for effective non-local exchange (Extended Data Figure 9a) shows Skala achieving 53.7% effective non-local exchange — close to the CCSD(T) reference of 60% — suggesting the learned within-atom non-locality does capture substantial exchange physics. But this metric measures an aggregate effect; it does not distinguish between within-atom and cross-atom contributions. The paper does not provide an ablation where cross-atom coarse-point communication is enabled vs. disabled, which would directly quantify the contribution of cross-atom non-locality. The theoretical framework in Section A.6 proves that a single downsampling-upsampling pass through a coarse point can approximate any two-body kernel — but only if that kernel operates on features that themselves encode cross-atom information. If the kernel needs to directly couple grid points on different atoms (e.g., $\kappa(\mathbf{r}_i^{\text{atom A}}, \mathbf{r}_j^{\text{atom B}})$), the current architecture cannot represent it without cross-atom coarse-point connections.

Mitigation status. The paper is explicit that it "do[es] not yet attempt to model dispersion explicitly" and uses a fixed D3 correction (Section 2.2). It frames learning dispersion within this architecture as future work. The paper does not propose an architectural extension to enable cross-atom message passing, nor does it discuss how such an extension would affect the computational cost scaling (cross-atom all-to-all communication would be O(N_atoms²), which could compromise the O(N³) scaling). The paper's framing — that the architecture suffices for "thermochemistry and kinetics" while dispersion is deferred — implicitly bounds the scope of what the current architecture is designed to capture. For a practitioner, the practical implication is that Skala is most reliable for problems dominated by local bonding (atomization energies, barrier heights, isomerizations) and may be less reliable for problems dominated by long-range correlation (large supramolecular complexes, π-stacking, extended van der Waals systems) where D3(BJ) with B3LYP parameters is the primary source of non-locality and the learned kernel contributes little.

7. Implications and Future Directions

How This Work Changes the Landscape

This paper changes the landscape of exchange-correlation functional development by reframing the problem from a search for better analytic forms to a problem of generating better training data. This is a paradigm shift in methodology, not just an incremental accuracy gain — it mirrors the transition that occurred in computer vision when AlexNet (2012) replaced hand-designed SIFT features with learned convolutional filters, and the paper draws this analogy explicitly (Section 2.3). The consequence is that progress in DFT accuracy is now gated primarily by the scale and diversity of high-accuracy reference data, not by human ingenuity in designing functional forms that balance exact constraints against empirical flexibility.

The concrete reframing is this: Jacob's ladder has been the organizing principle of XC functional development for three decades because it encodes a physical truth — that capturing non-local exchange and correlation requires ingredients that go beyond the local density and its gradients. The ladder says: if you want better accuracy, you must climb a rung, which adds a more expensive ingredient — almost always exact Hartree–Fock exchange. Skala demonstrates that this trade-off is not a law of nature; it is an artifact of the feature-engineering paradigm. When a neural network is given enough data, it can learn non-local exchange-correlation effects directly from semi-local inputs, at a cost far below that of computing the Hartree–Fock integrals. The evidence is in the numbers: a 2.80 kcal/mol WTMAD-2 on GMTKN55 (Table 1) at computational cost within ~30% of r2SCAN on GPU (Figure 4), compared to 3.23 for ωB97M-V at 3–6× higher cost. The ladder is not broken in the sense that non-locality is unnecessary — it is broken in the sense that you no longer need to pay the computational price of exact exchange to get the benefits of non-locality.

This shift has several downstream effects on the research landscape:

Data generation replaces functional design as the primary bottleneck. The paper's data ablation (Figure 2a, right) shows monotonic improvement as more categories of chemistry are added: Set A (thermochemistry only, ~284,000 reactions) yields ~8.4 kcal/mol WTMAD-2 on Diet GMTKN55 — GGA-level; adding barrier heights and reactions (Set B, ~398,000 reactions) drops this to ~6.6 — meta-GGA level; adding non-covalent interactions and conformers (Set C, full ~409,000 reactions) drops it to ~2.8 — hybrid level. This is a learning curve, not a saturation curve: there is no sign of diminishing returns. The implication is that the fastest path to better accuracy is not a cleverer architecture but more comprehensive training data — a message that redirects research effort from functional form development (which has been dominant for 60 years) to reference data generation and curation.

The versioned-release model supersedes the "DFT zoo." Traditional DFT has accumulated hundreds of functionals — Goerigk et al. (2017) assessed 200 — with no mechanism for retiring obsolete ones. Practitioners select functionals based on conflicting benchmark studies and domain-specific lore. Skala's versioned development (1.0 → 1.1, with documented improvements from expanding the training set from ~150,000 to ~400,000 data points) provides a different model: a single continuously improving functional that supersedes its predecessors through demonstrable accuracy gains. If this model is adopted by the community, it could simplify the "which functional should I use?" question that has plagued computational chemistry for decades — there would be one functional (Skala-N or its successor), and the question becomes "what version?"

Physical constraints are demoted from architectural requirements to emergent properties. This is the most philosophically significant shift in the paper. Traditional constraint-based functional development (the Perdew tradition, culminating in SCAN with 17 exact constraints) treats exact conditions as non-negotiable — a functional that violates them is considered physically unsound regardless of its empirical performance. Skala inverts this: it enforces only the most energetically critical constraints (Lieb-Oxford bound, high-density scaling, size-consistency) and lets the data teach the rest. Figure 2b shows that the positivity of the kinetic correlation energy T_c — a standard exact constraint — emerges spontaneously when training data are sufficiently diverse, but is violated when training is restricted to atomization energies only. This does not mean all constraints are learnable; the paper wisely retains the constraints that govern asymptotic limits (high-density scaling, Lieb-Oxford) where training data are effectively absent. But it does mean that the space of admissible functional architectures is much larger than previously thought — you do not need to restrict the functional form to guarantee every known constraint; you need enough data and the constraints will emerge if they are physically necessary for accurate predictions. This finding, if replicated across other constraints and other ML functionals, would fundamentally change how physicists think about the relationship between theory and data in electronic structure.

The paper reconciles contradictory findings in the ML-DFT literature. Prior ML functionals reached apparently conflicting conclusions: some showed promise (DM21: WTMAD-2 of 3.97 kcal/mol, approaching hybrid accuracy), while others found that learned functionals overfit, generalize poorly, or fail in self-consistent calculations. Skala's results suggest these contradictions were largely artifacts of data scale. DM21 was trained on a few thousand data points with carefully engineered features; it achieved respectable accuracy but required explicit non-local features (exact-exchange-like descriptors). Simpler models trained on similarly small datasets achieved worse results. Skala shows that with ~400,000 data points, even a model with only standard meta-GGA inputs can surpass the best hybrids — and that the transition from GGA-level to hybrid-level accuracy occurs somewhere between ~300,000 and ~400,000 training reactions (the gap between Set B and Set C in Figure 2a). This explains why prior work in the low-data regime could not break the accuracy–cost trade-off: they were on the flat part of the learning curve, where more data is the missing ingredient.

The accuracy–cost trade-off is reframed as a data–cost trade-off. The traditional view: accuracy costs FLOPs at inference time (climbing Jacob's ladder). The Skala view: accuracy costs reference data at training time, and inference FLOPs are decoupled from accuracy. This reframes the economics of DFT. A research group deciding whether to adopt Skala must weigh the one-time cost of generating or acquiring training data against the recurring savings from cheaper inference. For high-throughput applications (drug screening, materials discovery, training set generation for force fields), the inference savings dominate and Skala is compelling. For one-off calculations on small molecules, traditional hybrids may remain cheaper in total cost because they require no upfront data generation investment. The paper does not provide the total training data generation cost (CCSD(T)/CBS calculations are enormously expensive — the paper describes 12–60 million CPU-hours for the atomization energies alone, Section 2.3), making this calculus difficult for practitioners.

The work redirects research attention from architecture search to data curation. The non-local module in Skala is remarkably simple — 3 layers of downsampling-upsampling through coarse points, operating in only 16 dimensions. The total non-local parameter count is ~54,000 out of 385,217. The architecture is not the source of the accuracy gains; the scale and diversity of the training data are. The implication is that future effort should focus on generating reference data for currently unrepresented chemical regimes — multi-reference systems, transition metals, excited states, periodic systems — rather than on incrementally improving the message-passing architecture. The paper's transition-metal preview (Extended Data Figure 10) shows that even 35 training reactions produce competitive results on seen benchmarks and reasonable transfer to unseen ones; expanding this to hundreds or thousands of transition-metal training reactions would likely yield a Skala that handles transition metals at hybrid accuracy.

A note on what this work does not change. The paper does not claim that machine learning has "solved" DFT. The exact functional remains unknown, and Skala's accuracy — while state-of-the-art — is still 2.80 kcal/mol WTMAD-2, well above the ~1 kcal/mol target for chemical accuracy in covalent bond-making/breaking. The multi-reference problem is not addressed; indeed, the paper's training protocol explicitly filters out multi-reference character, and Skala's errors on the 17 multi-reference molecules in W4-17 are 4.61 kcal/mol (Extended Data Figure 8a). The transition-metal performance on unseen benchmarks is at the r2SCAN/B3LYP level, not at hybrid level. And the largest molecules in the training set have 13 non-hydrogen atoms; performance on C60ISO (8.16 kcal/mol, Table 1) indicates extrapolation to larger, delocalized systems is unreliable. These boundaries are clearly marked in the paper, and they define the research agenda for the next several years.

---

Follow-Up Research This Work Enables

Cheap difficulty estimation for Skala: training a model to predict Skala's error on a new molecule without running a reference calculation. The paper demonstrates that Skala's accuracy is not uniform across chemical space — it varies from ~0.2 kcal/mol on isomerizations to ~8.2 kcal/mol on C60 isomers (Table 1) to ~4.6 kcal/mol on multi-reference systems (Extended Data Figure 8a). A practical deployment of Skala needs to know when to trust its predictions and when to fall back to a more expensive method (e.g., a hybrid functional or explicitly correlated wavefunction method). A natural follow-up would train a lightweight classifier that takes as input the molecular structure (or Skala's own features, such as the per-point enhancement factor statistics or the non-local feature norms) and predicts Skala's expected MAE on that molecule, using the GMTKN55 and W4-17 benchmarks as labeled data. The classifier could flag molecules likely to have multi-reference character (using, e.g., the %TAE[(T)] diagnostic that the training pipeline already computes) or large delocalized π-systems (where the within-atom non-local module may be insufficient). A strong result would show that the classifier can identify 90% of molecules where Skala's error exceeds 2 kcal/mol while only flagging 20% of molecules where Skala is actually reliable, enabling a two-tier computation strategy where most molecules use cheap Skala and only flagged molecules incur the cost of a hybrid or wavefunction calculation.

Cross-atom coarse-point message passing: an architectural extension to enable learning long-range exchange and correlation beyond D3. The current architecture restricts non-local communication to within each atom's sub-grid (Section A.4, A.5). Coarse points at atomic centers never exchange messages with each other. This design choice preserves O(N³) scaling but may limit Skala's ability to capture non-local effects that depend directly on the relative positions and electronic structures of different atoms — phenomena like charge-transfer exchange, through-space conjugation, and many-body dispersion beyond the pairwise D3 model. An architectural extension would add cross-atom message passing between coarse points, perhaps using a similar equivariant framework (spherical harmonic projections of the distance vectors between coarse points, with learnable radial filters). The key question is whether the accuracy gains justify the computational cost increase, which could push scaling toward O(N_atoms²) in the non-local module. A strong experiment would compare Skala with and without cross-atom coarse-point communication on benchmarks that stress long-range exchange and correlation: the S66 non-covalent interaction set (where Skala already performs well at 0.07 kcal/mol, Table 1, so gains may be small), charge-transfer complexes (where semi-local functionals are notoriously poor), and the C60ISO set (where Skala's 8.16 kcal/mol error suggests missing physics). If cross-atom communication reduces the C60ISO error to below 4 kcal/mol without increasing the computational cost beyond, say, 2× r2SCAN, it would demonstrate that the within-atom restriction is the primary bottleneck for large delocalized systems. If it does not help, that would suggest the missing physics is something else — perhaps extensive static correlation in the π-system that no non-local exchange kernel can capture.

Scaling laws for XC functional accuracy: a systematic study of how Skala's error varies with training dataset size and composition. The paper's data ablation (Figure 2a, right) uses three coarse groupings (Sets A, A+B, A+B+C) and shows monotonic improvement, but this is not a proper scaling law. A systematic study would train Skala on random subsets of the full training data at sizes ranging from 10,000 to 400,000 reactions (e.g., 1%, 2.5%, 5%, 10%, 25%, 50%, 100%), with fixed data mixture proportions, and measure WTMAD-2 on Diet GMTKN55 and W4-17 for each size. This would produce a learning curve — error as a function of training set size — that could be fit to a power law Decay ∝ Nᵅ. The exponent α would characterize how efficiently Skala converts data into accuracy. If α is close to −0.5 (the typical deep learning scaling regime), then halving the error would require 4× more data — expensive but feasible. If α is shallower (e.g., −0.2), then reaching chemical accuracy (1 kcal/mol) would require orders of magnitude more data than is currently available, suggesting that architectural improvements (cross-atom non-locality, explicit dispersion learning, multi-reference handling) are needed in addition to more data. The study could also vary the data mixture independently — e.g., holding total dataset size constant while varying the fraction of barrier heights vs. NCIs — to determine whether some categories yield more "learning per reaction" than others. A finding that, say, 10,000 barrier heights improve overall accuracy more than 10,000 additional atomization energies would guide future data generation campaigns. This is analogous to the Chinchilla scaling laws for language models (Hoffmann et al., 2022) and would provide the first principled guidance for allocating computational resources between generating more data and improving the model architecture in computational chemistry.

Self-consistent training from scratch with better initialization: can we bypass the B3LYP pre-training dependency? The two-stage training protocol (pre-training on fixed B3LYP densities, then SCF fine-tuning) works but introduces a dependency on B3LYP that is not fully characterized. The paper does not ablate the choice of pre-training density functional, and there is no evidence that 20,000 fine-tuning steps are sufficient to erase B3LYP-specific biases. An alternative approach, enabled by the Skala architecture and training framework, would be to train from scratch using self-consistent densities produced by a sequence of increasingly accurate functionals — a curriculum learning strategy. Start with LDA (trivial, always converges), train Skala to reproduce CCSD(T) reaction energies on LDA densities, then use the trained Skala to generate slightly better self-consistent densities, retrain on those, and iterate. The key question is whether this self-consistent curriculum can converge to the same accuracy as the B3LYP pre-training protocol while eliminating the dependence on a specific traditional functional. A negative result — the curriculum fails to converge or converges to lower accuracy — would indicate that B3LYP densities contain essential physical information (correct spin states, reasonable bond lengths, proper charge distributions) that Skala cannot learn from scratch on LDA densities, suggesting a fundamental limitation of the data-driven approach: you need a "good enough" starting density to bootstrap. A positive result — the curriculum matches or exceeds B3LYP-pre-trained accuracy — would demonstrate that the framework is truly independent of traditional functionals and could be applied to regimes where B3LYP is unreliable (e.g., strong correlation).

Can Skala fine-tune on experimental data? The current training uses exclusively wavefunction reference energies (CCSD(T)/CBS), which limits the training to molecules small enough for coupled-cluster calculations (~10–20 non-hydrogen atoms). For larger molecules, experimental thermochemical data (heats of formation, reaction enthalpies) are available but carry measurement uncertainties and may include vibrational, thermal, and solvent effects not captured by the electronic energy alone. A natural extension would fine-tune a pre-trained Skala on experimental reaction energies for larger molecules, using the electronic energy as computed by Skala plus standard thermochemical corrections (zero-point energy, thermal enthalpy) to compare against experiment. This tests whether the learned functional representations transfer from the CCSD(T)/CBS regime (small molecules, gas-phase, vibrationless) to the experimental regime (larger molecules, condensed phase, finite temperature). A strong positive result would show that fine-tuning on a few hundred experimental reaction enthalpies reduces Skala's error on those reactions without degrading performance on the GMTKN55 benchmarks, indicating that the representations capture physically meaningful electronic effects that generalize across data sources. A negative result — fine-tuning helps on the experimental reactions but substantially degrades GMTKN55 performance — would indicate that experimental data introduce systematic biases (solvent effects, measurement errors) that conflict with the electronic structure information learned from wavefunction data, implying that training on experiment and training on theory are incompatible objectives for a single functional.

Skala as a drop-in replacement for B3LYP in high-throughput screening: a comprehensive benchmark on drug-discovery-relevant chemistry. The paper's evaluation focuses on the GMTKN55 and W4-17 benchmarks, which are designed to stress-test general main-group chemistry but do not directly represent the molecular types and property predictions relevant to pharmaceutical drug discovery — large, flexible molecules with multiple functional groups, conformational ensembles, protonation states, and tautomers. A practical evaluation would run Skala (and baseline functionals: B3LYP-D3(BJ), ωB97M-V, r2SCAN) on standard drug-discovery benchmarks: the PDBbind refined set (protein-ligand binding affinities, where the electronic energy contributes to the interaction), the TorsionNet500 set (conformational energies of drug-like fragments), and the SAMPL host-guest binding challenges. The key metrics would be ranking correlation (Spearman ρ) rather than absolute MAE, since drug discovery pipelines typically rank candidates rather than requiring absolute binding energies. Given Skala's strong performance on conformational energies in GMTKN55 (0.09 kcal/mol on SCONF for sugar conformers, Table 1), it may excel at the TorsionNet500 benchmark, where semi-local functionals are known to struggle with intramolecular dispersion. A finding that Skala matches ωB97M-V's ranking accuracy on protein-ligand binding at r2SCAN-level cost would make it immediately attractive for pharmaceutical virtual screening, where the cost of hybrid functionals prohibits their routine use on large compound libraries.

---

Practical Applications and Downstream Use Cases

High-throughput virtual screening for drug discovery. Modern drug discovery pipelines routinely screen libraries of 10⁶–10⁹ candidate molecules against a protein target, typically using docking scores or force-field-based methods because DFT with hybrid functionals is prohibitively expensive at this scale. The result is that computational screens pass too many false positives to experimental validation — a problem the paper explicitly identifies in its introduction. Skala offers a direct path to upgrading the energy evaluation in these pipelines from force-field accuracy (~2–5 kcal/mol errors for non-covalent interactions) to hybrid-DFT accuracy (~0.5–1.0 kcal/mol on the S66 benchmark, Table 1) at a computational cost increase of only ~30% over a meta-GGA like r2SCAN (Figure 4, GPU: Skala ~0.15 s/SCF step vs. r2SCAN ~0.12 s for 113-atom systems). For a library of 10⁶ molecules, evaluating each with a single-point Skala calculation at the def2-TZVP level (~2,000–3,000 basis functions for drug-sized molecules, roughly the 2,228-orbital data point in Figure 4) would cost ~150,000 GPU-hours — feasible on a moderate cluster. The same library evaluated with ωB97M-V at ~0.50 s/step would cost ~500,000 GPU-hours and may exhaust memory for the largest molecules. The practical benefit is a reduction in false-positive rate in computational screening, which directly translates to fewer wasted synthesis and assay cycles in the laboratory. The paper's strong performance on non-covalent interactions — particularly S66 (0.07 kcal/mol vs. 0.14 for ωB97M-V) and HAL59 (0.23 vs. 0.28 for ωB97M-V) — suggests Skala is well-suited to the types of intermolecular interactions that dominate protein-ligand binding.

Training set generation for machine-learned interatomic potentials. Many state-of-the-art force fields (ANI, MACE, NequIP, Allegro) are trained on DFT data — typically at the ωB97M-V or ωB97X-D3 level for accuracy, or at the PBE/r2SCAN level for cost reasons. The accuracy of the resulting force field is bounded by the accuracy of the DFT functional used to generate its training data. Skala offers the opportunity to generate training data at hybrid-level accuracy (or better) at semi-local cost, which could improve force field quality without increasing the computational budget for data generation. For example, generating 1 million energy and force evaluations for a MACE training set using Skala at def2-TZVP would cost roughly the same as generating the same data with r2SCAN but would encode ωB97M-V-quality potential energy surfaces. The paper's GEO metric (~0.02 kcal/mol, Supplementary Figure 15a) and geometry optimization accuracy (Table 3c) confirm that Skala's potential energy surfaces are globally accurate, not just at energy minima — essential for force field training, which requires accurate forces across the entire configuration space. The improvement in force field accuracy would cascade downstream to molecular dynamics simulations of protein folding, materials failure, and reaction dynamics — all of which currently inherit the errors of the underlying DFT functional. Given that DFT data generation is often the computational bottleneck in force field development, a 3–6× cost reduction (Skala vs. hybrids on GPU, Figure 4) with equal or better accuracy is immediately actionable.

Replacing B3LYP in routine computational chemistry workflows. B3LYP is the most cited functional in history — not because it is the most accurate, but because it represents an acceptable compromise between cost and accuracy for a wide range of chemical problems. Skala's accuracy (2.80 WTMAD-2 vs. B3LYP's 6.38, Table 1) and computational cost (within ~30% of r2SCAN, which is typically ~2–3× cheaper than B3LYP) position it as a candidate to replace B3LYP as the "default functional" in electronic structure codes. For a typical computational chemistry project — say, optimizing the geometry of a 50-atom organic molecule and computing its reaction barrier with a small substrate — B3LYP/def2-TZVP might require ~100 SCF steps for geometry optimization plus ~30 steps for the transition state search, totaling ~130 SCF steps. On GPU, B3LYP with sn-LinK would take ~65 seconds (at ~0.5 s/step for a molecule of this size, extrapolating Figure 4). Skala would take ~20 seconds (at ~0.15 s/step) — a 3× speedup — while providing accuracy that is substantially better than B3LYP across thermochemistry, kinetics, and non-covalent interactions. The practical benefit to a bench chemist is faster turnaround for routine calculations, enabling more extensive conformational sampling, larger basis sets, or explicit solvent modeling within the same computational budget. The barrier to adoption is integration into production codes — the paper notes that Psi4 and CP2K are already integrating Skala through GauXC, and PySCF integration is available through Python/C/Fortran interfaces (Section 5). As these integrations mature, the transition from B3LYP to Skala for routine main-group calculations becomes technically straightforward.

Accelerating catalyst and materials discovery. Computational catalyst design requires screening thousands of candidate materials for their ability to lower reaction barriers — a task that demands accurate barrier heights (typically 1–2 kcal/mol errors are needed for reliable predictions). The best hybrid functionals (ωB97M-V: 1.33 kcal/mol on BH76, Table 1) achieve this for main-group reactions, but their cost prohibits large-scale screening. Skala's barrier height accuracy (1.75 kcal/mol on BH76, 0.53 kcal/mol on BHDIV10) is competitive with ωB97M-V, and its computational cost is 3–6× lower. For a screening campaign evaluating 5,000 catalyst candidates, each requiring a transition-state search (~50 SCF steps for the bare reaction plus ~100 steps for each of the 5,000 catalyst variants), the total SCF step count is ~500,000. At Skala's cost (~0.15 s/step on GPU), this is ~21 GPU-hours. At ωB97M-V's cost (~0.5 s/step), it is ~69 GPU-hours. The difference of ~48 GPU-hours per screening campaign is significant but not transformative — the larger advantage is that Skala can handle larger catalyst models (more atoms) before hitting memory or time limits. Since transition-metal catalysts are the most industrially relevant but also the most computationally demanding, extending this application to transition metals requires first solving the data scarcity problem identified in Extended Data Figure 10. However, even for main-group organocatalysis (which is an active area of green chemistry), Skala's combination of hybrid-level barrier heights and semi-local cost is immediately useful.

---

When to Prefer This Method

The paper positions Skala explicitly against the traditional Jacob's ladder paradigm — it is an alternative to climbing the ladder, not an additional rung. The decision between Skala and traditional functionals depends on the specific chemical problem, the available computational resources, and the tolerance for the slightly less robust SCF convergence behavior documented in Tables 3 and 4. The paper's results support the following decision framework:

• Prefer Skala for main-group thermochemistry and kinetics on molecules with up to ~50–100 atoms. The GMTKN55 results (Table 1) show Skala dominating these categories, and the computational cost advantage over hybrids (3–6× on GPU) grows with system size. The def2-TZVP basis set provides good cost/accuracy balance and shows better SCF convergence (99.38% vs. 98.87% for def2-QZVP, Tables 3 and 4).

• Prefer Skala when the computational budget is dominated by the number of SCF calculations rather than individual calculation cost — high-throughput screening, training set generation, conformational sampling. The cost per SCF step is within 30% of r2SCAN, enabling orders of magnitude more calculations than hybrid functionals for the same budget. The accuracy is hybrid-level or better, so there is no accuracy penalty for scaling up the number of calculations.

• Prefer traditional hybrid functionals (ωB97M-V, ωB97X-V) for systems with significant multi-reference character, extensive π-conjugation beyond ~20 atoms, or transition metals with little training data coverage. The W4-17 multi-reference MAE (4.61 kcal/mol, Extended Data Figure 8a), C60ISO error (8.16 kcal/mol, Table 1), and unseen transition-metal benchmark performance (Extended Data Figure 10, at r2SCAN/B3LYP level) all indicate that Skala degrades to semi-local quality outside its training distribution. Hybrid functionals, with their physically motivated exact-exchange fraction, may extrapolate more robustly to these regimes even if their average accuracy on main-group benchmarks is lower.

• Prefer Skala when accurate electron densities (dipole moments) matter as much as accurate energies. The dipole RMSE of 4.43% (Table 3b) is lower than all tested hybrids (ωB97M-V: 5.84%, ωB97X-V: 5.18%), indicating that Skala's self-consistent densities are genuinely improved — not just its energies. For applications where dipole moments, polarizabilities, or electrostatic potentials are the quantities of interest, Skala offers state-of-the-art accuracy at lower cost than hybrids.

• Prefer traditional meta-GGAs (r2SCAN) when robustness and code maturity are paramount — educational settings, exploratory calculations on unknown systems, integration with codes that do not support GauXC. Skala's SCF convergence issues (98.87% first-attempt, requiring a six-level retry protocol) and dependence on the GauXC library for efficient grid operations mean it is not yet as "plug-and-play" as r2SCAN, which is implemented in virtually every electronic structure code via libxc and converges robustly. For a researcher exploring a new chemical system where convergence problems would be difficult to diagnose, the slightly worse accuracy of r2SCAN may be preferable to the slightly less robust convergence of Skala.

• Prefer Skala when the goal is to generate reference-quality data for downstream machine learning. The versioned-release model and systematic improvability with data mean that Skala's accuracy will improve over time, and using it to generate training data for force fields or property prediction models ensures that those downstream models inherit the best available accuracy. This "cascade of accuracy transfer" is explicitly highlighted in the paper's conclusion (Section 6).