Paper deep dive

Coupled-cluster molecular properties across the main group that extrapolate beyond training size

Wenhao He, Xu Chen, Noah Song, Haowei Xu, Tim S. Hindges, Bohan Li, Zihan Lin, Yu Yao, Avetik R. Harutyunyan, Fang Liu, Yao Wang, Hao Tang, Ju Li

Year: 2026Venue: arXiv preprintArea: physics.chem-phType: PreprintEmbeddings: 68

Intelligence

Status: succeeded | Model: Gemma-4-26B-A4B | Prompt: intel-v1 | Confidence: 94%

Last extracted: 8/20/2026, 4:25:10 AM

Summary

The paper introduces MEHnet-MG, an equivariant neural network that predicts an effective one-electron Hamiltonian to derive multiple molecular properties (energy, optical gap, dipole, quadrupole, polarizability, Mulliken charges, Mayer bond orders) with coupled-cluster accuracy. Trained on a dataset of 41,939 molecules across nine main-group elements (H, C, N, O, F, Si, P, S, Cl) labeled at the CCSD(T) level, the model uses a B3LYP/def2-SVP baseline and adds learned corrections. It significantly outperforms DFT methods (BP86, B3LYP, DSD-PBEP86) in accuracy while maintaining low computational cost (~25 ms per molecule). Crucially, by deriving properties from a predicted Hamiltonian rather than pooling atomic features, the model exhibits correct size-scaling and can extrapolate to larger systems (up to 58 atoms) where traditional pooling-based architectures fail.

Entities (10)

MEHnet-MG · model · 98%B3LYP · method · 95%CCSD(T) · method · 95%Chlorine · element · 95%Optical Gap · property · 95%Phosphorus · element · 95%Polarizability · property · 95%Sulfur · element · 95%EquiformerV2 · model · 90%Oligothiophenes · chemical-system · 90%

Relation Signals (10)

MEHnet-MG → coverselements → Sulfur

confidence 95% · across nine main-group elements, including the under-served... sulfur

MEHnet-MG → coverselements → Chlorine

confidence 95% · across nine main-group elements, including the under-served... chlorine chemistries.

MEHnet-MG → coverselements → Phosphorus

confidence 95% · across nine main-group elements, including the under-served phosphorus

MEHnet-MG → predictsproperty → Optical Gap

confidence 95% · derives a broad suite of properties from it (energy, optical gap...)

MEHnet-MG → predictsproperty → Polarizability

confidence 95% · derives a broad suite of properties from it (...polarizability...)

MEHnet-MG → trainedondatalevel → CCSD(T)

confidence 95% · The model is trained on a new in-house dataset of multi-property labels computed at the CCSD(T) level

MEHnet-MG → usesbaselinemethod → B3LYP

confidence 95% · predicts an effective one-electron Hamiltonian from one inexpensive B3LYP/def2-SVP calculation

MEHnet-MG → backbonearchitecture → EquiformerV2

confidence 90% · The backbone is a spherical-harmonic graph-attention network (EquiformerV2...)

Cypher Suggestions (0)

No Cypher suggestions yet.

Abstract

Abstract:Coupled-cluster theory defines the accuracy standard for molecular electronic-structure properties but scales too steeply for routine application, whereas density-functional theory is affordable yet systematically biased. We resolve this trade-off with a single equivariant network, MEHnet-MG, that predicts an effective one-electron Hamiltonian from one inexpensive B3LYP/def2-SVP calculation and derives a broad suite of properties from it (energy, optical gap, dipole, quadrupole, polarizability, Mulliken atomic charges, and Mayer bond orders) at coupled-cluster accuracy across nine main-group elements, including the under-served phosphorus, sulfur, and chlorine chemistries. The model is trained on a new in-house dataset of multi-property labels computed at the CCSD(T) level for all nine elements. On a held-out test set, it reduces the error of every property by a factor of 3.8 to 230 relative to semi-local, hybrid, and double-hybrid DFT (referenced to composite CCSD(T)/cc-pVTZ; Methods), while adding only ~25 ms wall time per molecule, delivering coupled-cluster-quality predictions at the cost of a single DFT calculation. Critically, deriving every property from a predicted Hamiltonian rather than pooling per-atom features builds the correct size-scaling into the model architecture: on pi-conjugated oligothiophenes it matches finite-field CCSD polarizability and the EOM-CCSD optical gap to ~2% at the largest sizes where those references remain affordable (44 and 37 atoms, where a single CCSD field point already costs ~500x the model's entire inference) and extrapolates the corrected trends to 58-atom chains, a regime where pooling-based architectures fail by construction. Accurate extrapolation is therefore set by the model's inductive bias rather than by the training data.

PDF

Open source PDF →Open local PDF →

Trouble viewing inline? Open PDF directly →

Full Text

64,546 characters extracted from source content.

Expand or collapse full text

Coupled-cluster molecular properties across the main group that extrapolate beyond training size Wenhao He Thanks: These authors contributed equally. Affiliation: Center for Computational Science and Engineering, Massachusetts Institute of Technology, Cambridge, MA 02139, USA Xu Chen†footnotemark: Affiliation: Department of Chemistry, Emory University, Atlanta, GA 30322, USA Noah Song†footnotemark: Affiliation: Department of Chemistry, Emory University, Atlanta, GA 30322, USA Haowei Xu Affiliation: Department of Nuclear Science and Engineering, Massachusetts Institute of Technology, Cambridge, MA 02139, USA Tim S. Hindges Affiliation: Department of Nuclear Science and Engineering, Massachusetts Institute of Technology, Cambridge, MA 02139, USA Bohan Li Affiliation: Center for Computational Science and Engineering, Massachusetts Institute of Technology, Cambridge, MA 02139, USA Zihan Lin Affiliation: Department of Materials Science and Engineering, Massachusetts Institute of Technology, Cambridge, MA 02139, USA Yu Yao Affiliation: Department of Nuclear Science and Engineering, Massachusetts Institute of Technology, Cambridge, MA 02139, USA Avetik R. Harutyunyan Affiliation: Honda Research Institute USA, San Jose, CA 95134, USA Fang Liu Affiliation: Department of Chemistry, Emory University, Atlanta, GA 30322, USA Yao Wang Affiliation: Department of Chemistry, Emory University, Atlanta, GA 30322, USA Hao Tang Thanks: Correspondence: haot@mit.edu Affiliation: Department of Materials Science and Engineering, Massachusetts Institute of Technology, Cambridge, MA 02139, USA Ju Li Thanks: Correspondence: liju@mit.edu Affiliation: Center for Computational Science and Engineering, Massachusetts Institute of Technology, Cambridge, MA 02139, USA Affiliation: Department of Materials Science and Engineering, Massachusetts Institute of Technology, Cambridge, MA 02139, USA Affiliation: Department of Nuclear Science and Engineering, Massachusetts Institute of Technology, Cambridge, MA 02139, USA Abstract Coupled-cluster theory defines the accuracy standard for molecular electronic-structure properties but scales too steeply for routine application, whereas density-functional theory is affordable yet systematically biased. We resolve this trade-off with a single equivariant network, MEHnet-MG, that predicts an effective one-electron Hamiltonian from one inexpensive B3LYP/def2-SVP calculation and derives a broad suite of properties from it (energy, optical gap, dipole, quadrupole, polarizability, Mulliken atomic charges, and Mayer bond orders) at coupled-cluster accuracy across nine main-group elements, including the under-served phosphorus, sulfur, and chlorine chemistries. The model is trained on a new in-house dataset of multi-property labels computed at the CCSD(T) level for all nine elements. On a held-out test set, it reduces the error of every property by a factor of 3.83.8 to 230230 relative to semi-local, hybrid, and double-hybrid DFT (referenced to composite CCSD(T)/c-pVTZ; Methods), while adding only ∼ 25 ms wall time per molecule, delivering coupled-cluster-quality predictions at the cost of a single DFT calculation. Critically, deriving every property from a predicted Hamiltonian rather than pooling per-atom features builds the correct size-scaling into the model architecture: on π-conjugated oligothiophenes it matches finite-field CCSD polarizability and the EOM-CCSD optical gap to ∼ 2% at the largest sizes where those references remain affordable (4444 and 3737 atoms, where a single CCSD field point already costs ∼ 500× the model’s entire inference) and extrapolates the corrected trends to 58-atom chains, a regime where pooling-based architectures fail by construction. Accurate extrapolation is therefore set by the model’s inductive bias rather than by the training data. 11footnotetext: These authors contributed equally.22footnotetext: Correspondence: haot@mit.edu33footnotetext: Correspondence: liju@mit.edu Predicting the electronic properties of molecules from their atomic structure is a foundational capability of the chemical sciences, underpinning the rational design of drugs, catalysts, electrolytes, and functional materials [1, 2]. As computational screening increasingly precedes and guides experiments, the rate of discovery is set by how accurately and inexpensively electronic-structure properties can be evaluated across the vast space of candidate compounds [3]. A method that combines the accuracy of high-level wavefunction theory with a cost low enough for routine application to large molecular libraries would transform this process. Developing such approaches has therefore long been a central goal of computational chemistry. The trade-off between accuracy and efficiency is the clearest between high-level wavefunction-based quantum chemistry method and density-functional theory. Coupled-cluster theory with perturbative triples, CCSD(T), is the de facto gold standard for single-reference molecules [4], but its ⁡(N7)O(N^7) scaling confines it to small systems and precludes high-throughput use. Density-functional theory [5] is far more affordable and is the most widely applied electronic-structure method, yet its accuracy is functional-dependent and subject to systematic errors that no single functional eliminates. Generalized-gradient approximations, in particular, exhibit a delocalization (self-interaction) error that degrades response properties and worsens with extended π-conjugation [6, 7]. Ascending the hierarchy of functionals reduces some errors while introducing others, at steadily increasing cost. Decades of methodological development have narrowed, but not closed, the gap between affordable and accurate electronic-structure methods. Machine learning offers a route around this trade-off, because a trained model’s inference cost is set by its input, not by its supervision: a network trained on coupled-cluster labels whose input is a single inexpensive DFT calculation (here, a Δ -learning correction to the DFT Fock matrix) reproduces correlated-level accuracy at DFT cost. Two broad strategies have been pursued, and each, we argue, addresses only one half of the problem. The first learns the electronic structure itself, regressing the Kohn–Sham Hamiltonian from geometry and then diagonalizing it, so that properties follow from quantum mechanics with the correct dependence on system size. This route has been developed for periodic solids by DeepH and its equivariant successors [8, 9]. For molecules, related approaches include SchNOrb [10], PhiSNet [11], QHNet [12], and, more recently, HELM, which spans much of the periodic table [13]. Because the supervision is the reference Hamiltonian matrix, however, the accuracy is capped at the functional that produced the labels, and the model accelerates DFT rather than surpassing it. The second strategy learns the properties directly. Equivariant property networks such as SpookyNet [14], PaiNN [15], and AIMNet2 [16], and Δ -learning models more broadly [17], can be trained against correlated or experimental references and so are not capped in accuracy. Their limitation lies instead in how each property is produced: it is read out of pooled atomic features, a sum or mean over the graph. This pooling fixes the size-scaling of every property to be strictly extensive or strictly intensive, so a response that grows super-linearly with system size cannot be represented, and accuracy degrades once the model is applied beyond its training sizes. A separate output head must also be trained for each property, in contrast to the single electronic-structure object from which an entire suite of observables follows by construction [10, 11]. Most of these models are, moreover, confined to first- and second-row organic chemistry (C, H, N, O). Indeed, Si, P, S, and Cl on the same footing are rare, and many widely used potentials omit P and S entirely. A third strategy improves the exchange–correlation functional itself: deep-learning functionals such as DM21 [18] and Skala [19] attain near-coupled-cluster accuracy for energies at semi-local cost and span the main group. Their supervision and validated benchmarks are, however, energetic (DM21 additionally constrains fractional-charge and fractional-spin behaviour): every other observable still follows from the resulting Kohn–Sham solution, and accuracy at the correlated level for densities, response properties, or excitations is neither trained for nor yet demonstrated. Our model occupies the intersection not addressed by these strategies. Like MEHnet, the Hamiltonian-learning approach introduced for hydrocarbons by Tang et al. [20], our model—MEHnet-MG, for the main group—predicts an effective one-electron Hamiltonian and derives every property from it using exact quantum-mechanical operators, inheriting the correct size-scaling by construction. Critically, it is supervised on the observables rather than on a reference Hamiltonian matrix: a single learned correction to one cheap B3LYP/def2-SVP single point is trained so that the derived properties match a correlated reference; the model is therefore not capped at DFT accuracy and could, in principle, target any level of theory. Our contributions are fourfold: (i) deriving every property from the predicted Hamiltonian by exact operators builds the correct size-scaling into the architecture itself, and the model demonstrably extrapolates out of distribution: it tracks finite-field CCSD polarizability on π-conjugated oligothiophenes to ∼ 2% out to the 44-atom CCSD limit and continues the corrected trend to 58-atom chains, all from sub-25-atom training molecules, where pooling-based architectures diverge (Fig. 5), because the inductive bias, not the data, sets the size-scaling; (i) it reaches CCSD(T)/CCSD/EOM-CCSD-level accuracy on seven ground- and response-properties at the cost of one DFT single point, all from one baseline rather than a per-property construction; (i) it spans nine main-group elements, treating the under-served P/S/Cl chemistries on the same footing as the organic elements; and (iv) it does so with an architecture built for the main group: a higher-angular-resolution equivariant backbone (EquiformerV2, lmax=4l_ =4, as required by the d-orbital Fock blocks of the third-row elements), Δ -learning extended to every predicted object—an explicit energy-correction head and a screening matrix anchored to the baseline response, so that an untrained correction reproduces the baseline polarizability exactly—and a hybrid baseline whose near-unit-slope Kohn–Sham gap lets the optical gap be read directly from the frontier eigenvalues, with no learned rescaling (SI). Two principles organize these results: the model’s architecture determines whether a property extrapolates beyond the training distribution, and the level of supervision sets the accuracy it can reach. Together they place a broad suite of molecular properties at the favourable corner of the accuracy–cost plane across the nine main-group elements covered here. Results Model architecture and dataset Figure 1: Model overview. From a molecule graph, an equivariant network (input embedding followed by K E(3)-convolution layers) produces node and edge features that feed three correction heads: ΔEθ E_θ, Δθ _θ, and Δθ _θ. Each correction is added to the corresponding baseline (EDFTE_DFT, DFTH_DFT, DFTT_DFT) from the cheap (TD)DFT (B3LYP/def2-SVP) reference to give the total energy E, the effective Hamiltonian H, and the screening matrix T. Properties are then read off by the defining quantum-mechanical operator rather than pooled from atomic features: the energy from E; the dipole, quadrupole, Mulliken atomic charges, Mayer bond orders, and optical gap from H; and the polarizability from the screening matrix T. The dashed path denotes gradients (response properties obtained by differentiation). The full forward pass adds only ∼ 25 ms (GPU) on top of the DFT baseline, so coupled-cluster-quality properties are obtained at the cost of one DFT single point. Rather than correcting each target property directly, our model applies a learned correction to a cheap physical baseline at the level of a few underlying electronic-structure objects, and derives every property from them (Fig. 1). A single B3LYP/def2-SVP calculation supplies the baseline: a total energy EDFTE_DFT, a Fock matrix DFTH_DFT, and a screening matrix DFTT_DFT. The ground-state single point provides EDFTE_DFT and DFTH_DFT, and hence every property except the polarizability; only DFTT_DFT requires an additional linear-response (TDDFT) step on the same baseline, so the TDDFT cost is incurred only when α is requested. An E(3)-equivariant graph network reads the molecule graph and produces, through separate heads, three corrections to this baseline, ΔEθ E_θ, Δθ _θ, and Δθ _θ, which are added back to give a corrected energy E, an effective Hamiltonian H, and a screening matrix T. Its backbone is a spherical-harmonic graph-attention network (EquiformerV2 [21]; ∼ 5.6 M parameters, four equivariant interaction layers, spherical-harmonic degree lmax=4l_ =4, a 66 Å cutoff); on top of this backbone we add equivariant heads, built with e3n [22], that output the three corrections. Full settings are given in Methods. The essential design choice is that the network outputs corrections to these physical objects, not the target properties themselves. Because the baseline is inexpensive and the network adds only milliseconds (Fig. 3a), the total cost is essentially that of the hybrid-DFT single point, while the accuracy approaches the coupled-cluster labels the network is trained on. A single such model spans all nine main-group elements (H, C, N, O, F, Si, P, S, Cl) and yields seven properties from one forward pass: the total electronic energy E, the gap EgE_g, and Mulliken atomic charges [23] qiq_i (scalars); the dipole μ (a vector); the quadrupole Q and polarizability α (rank-2 tensors, predicted equivariantly); and pairwise Mayer bond orders [24] bijb_ij. Each corrected object yields its properties by the operator that defines it, rather than through a separate property head. The corrected energy E=EDFT+ΔEθE=E_DFT+ E_θ is the total electronic energy. Diagonalizing the effective Hamiltonian =DFT+ΔθH=H_DFT+ _θ gives molecular-orbital energies and coefficients, and hence the one-particle density matrix, from which follow the optical gap (the frontier HOMO–LUMO eigenvalue difference), the dipole and quadrupole (multipole operators traced against the density), and the atomic charges and bond orders. The polarizability uses the corrected screening matrix =DFT+ΔθT=T_DFT+ _θ: a bare sum-over-states response built from the eigenstates of H is screened by T to give the reported tensor (Methods, Eqs. (2)–(3)), and because DFTT_DFT comes from the baseline DFT response, an uncorrected model reproduces the baseline polarizability. Because each property is the exact functional of objects assembled from local corrections, its dependence on system size is fixed by quantum mechanics rather than by a pooling rule, which is what allows the model to extrapolate collective responses beyond the training sizes (Fig. 5). MEHnet-MG keeps this effective-Hamiltonian read-out from its hydrocarbon predecessor [20] but rebuilds the machinery around it. The two-layer EGNN backbone of MEHnet (lmax=2l_ =2) becomes the four-layer EquiformerV2 above, carrying angular momenta to lmax=4l_ =4 as the d×d× d Fock blocks of Si–Cl require; the correction set gains an explicit energy head ΔEθ E_θ (MEHnet read the energy as the occupied-eigenvalue sum); the screening matrix becomes a Δ -learning quantity anchored to the baseline response, =DFT+ΔθT=T_DFT+ _θ, where MEHnet learned T from scratch; and the hybrid baseline allows the gap to be the bare frontier eigenvalue difference, retiring MEHnet’s learned gap rescaling, Eg=(1+G1)(εL−εH)+G2E_g=(1+G_1)( _L- _H)+G_2. Retiring the rescaling is more than a simplification. The G1,G2G_1,G_2 correctors were attention-pooled learned scalars—precisely the pooling-style read-out this paper argues has no guaranteed size behaviour—and they were the last property that bypassed the Hamiltonian. With them removed, every property in MEHnet-MG flows through the operator read-out, so the size transferability above applies to the full suite; the out-of-distribution gap of Fig. 5b is its direct test. A row-by-row comparison with MEHnet, including accuracy, is given in the SI. The model is trained on a new in-house dataset built from 44,41244,412 three-dimensional structures from PubChem and passed through a three-stage quality filter (Methods). Designed to broaden chemical coverage, the dataset includes compounds containing silicone, phosphorus, sulfur, chlorine, and fluorine in addition to the organic elements H, C, N, O. (Fig. 2a). The retained molecules are small, neutral, closed-shell singlets of 22–2424 atoms (≤9≤ 9 heavy, median 1414; Fig. 2b). In terms of element coverage, hydrogen (99.3%99.3\% of molecules) and carbon (98.1%98.1\%) are near-ubiquitous. The set is also rich in the heteroatoms that most omit: sulfur in 58.1%58.1\% of molecules, nitrogen in 56.8%56.8\%, oxygen in 47.6%47.6\%, chlorine in 21.9%21.9\%, phosphorus in 15.2%15.2\%, silicon in 9.5%9.5\%, and fluorine in 6.5%6.5\% (Fig. 2a). Every molecule is labeled at property-appropriate coupled-cluster levels (energies at composite CCSD(T)/c-pVTZ—canonical CCSD(T)/c-pVDZ plus a DLPNO-CCSD(T) basis-set correction to c-pVTZ, Eq. (1)—polarizabilities at finite-field CCSD/c-pVDZ, gaps at EOM-CCSD/c-pVDZ, and dipoles, quadrupoles, Mulliken charges and Mayer bond orders from the coupled-cluster density at the same composite level; Methods), and the resulting ground-truth properties span a broad chemical range (Fig. 2c,d). After filtering, 41,93941,939 molecules remain. These are partitioned into training, validation, and a held-out test set of 959959 molecules, which we report on throughout this work. Figure 2: Dataset composition and coverage. a Elemental composition of the dataset across the nine main-group elements included in this work (periods 1–3). Shaded as a heatmap to show the number of molecules in the dataset containing a given element. b Distribution of molecular size, measured as atoms per structure, with median size 1414. c UMAP embedding of the CCSD(T) property space, colored by heavy-atom count. d Distributions of the seven CCSD(T) ground-truth properties across the dataset: total ground state energy E, the optical gap EgapE_gap, the dipole ||| μ|, the quadrupole ‖\|Q\|, the isotropic polarizability αiso _iso, the per-atom Mulliken charge q, and the per-bond Mayer bond order (dashed lines mark medians). Note that all panels cover the 42,91242,912 quality-controlled structures, but the final train/val/test set contains only the 41,93941,939 that additionally pass baseline and outlier filtering (Methods). Accuracy benchmark across all properties Table 1: Global mean absolute error (MAE) versus the composite CCSD(T)/c-pVTZ ground truth (Methods) on the 959-molecule test set, per property and method. Lower is better; best in bold. Units: E in kcal mol-1, EgapE_gap in eV, |μ||μ| in Debye; ‖Q‖\|Q\| and αiso _iso in a.u.; charge in e. Each method is evaluated at the basis set of the property’s coupled-cluster reference: c-pVTZ for E (compared as an atomization energy), |μ||μ|, ‖Q‖\|Q\|, q and bond order, and c-pVDZ for αiso _iso and EgapE_gap (finite-field CCSD and EOM-CCSD, respectively). DSD-PBEP86 is the D3(BJ)-corrected double hybrid. property BP86 B3LYP DSD-PBEP86 MEHnet-MG (this work) E (kcal mol-1) 63.0 13.2 15.1 0.27 EgapE_gap (eV) 1.32 0.61 4.84 0.067 |μ||μ| (Debye) 0.160 0.159 0.109 0.022 ‖Q‖\|Q\| (a.u.) 0.389 0.347 0.233 0.045 αiso _iso (a.u.) 4.71 2.79 1.18 0.121 q (e) 0.0205 0.0256 0.0220 0.0054 bond order 0.0556 0.0588 0.0251 0.0066 Figure 3: Bulk accuracy and cost. a Per-molecule wall time per single-point evaluation versus atom count (log–log) with power-law fits; the model shares its DFT baseline’s shallow scaling, while CCSD(T) and the double-hybrid scale steeply. b Per-property error radar on a shared log scale (inner → smaller error vs. CCSD(T)). c Per-element MAE across the nine supported elements, per-atom averaged; units: Mulliken charge in e, E in kcal mol-1/atom, EgapE_gap in eV, μ in D/atom, Q and αiso _iso in a.u./atom. As in Table 1, each functional is evaluated at the basis of the property’s coupled-cluster reference (c-pVTZ for E, |μ||μ|, ‖Q‖\|Q\|, q and bond order; c-pVDZ for αiso _iso and EgapE_gap), all referenced to composite CCSD(T)/c-pVTZ (Methods), finite-field CCSD/c-pVDZ and EOM-CCSD/c-pVDZ; panel a times the c-pVTZ single point of each functional against DLPNO-CCSD(T)/c-pVTZ, the c-pVTZ component of the composite reference. On the 959-molecule test set, we compare the model against three widely used density functionals (GGA BP86, hybrid B3LYP, and double-hybrid DSD-PBEP86). Each functional is evaluated at the basis set of the property’s coupled-cluster reference (c-pVTZ for energies, dipoles, quadrupoles, Mulliken charges and Mayer bond orders; c-pVDZ for polarizabilities and optical gaps), referenced to the composite CCSD(T)/c-pVTZ, finite-field CCSD/c-pVDZ, and EOM-CCSD/c-pVDZ ground truth (Table 1, Fig. 3). The model reduces the mean absolute error of every property relative to the DFT methods by factors ranging from ∼ 2.52.5 (per-atom charges and bond orders) to over two orders of magnitude (energies). The total energy reaches sub-chemical accuracy (∼ 10−310^-3 Ha, well below 11 kcal mol-1), more than two orders of magnitude below BP86, and the tensorial and response properties (quadrupole, polarizability) improve by roughly an order of magnitude, the regime where DFT performs worst. Simply ascending the functional hierarchy does not close this gap, whereas the learned correction does: on this set B3LYP is not consistently better than BP86, while the double-hybrid DSD-PBEP86 — each functional now evaluated at a basis matched to its coupled-cluster reference — is the most accurate of the three for the tensorial and response properties (dipole, quadrupole, polarizability, and bond order), consistent with its design, whereas B3LYP is best for the energy. The double hybrid’s Kohn–Sham orbital gap, however, remains a poor proxy for the EOM-CCSD optical gap. MEHnet-MG nonetheless leads every property in the table. Deriving every property from one predicted Hamiltonian by exact operators is what lets MEHnet-MG report the complete suite at one accuracy level. The correction also generalizes uniformly across chemistry, with no degradation for the rarer elements F and Si relative to C/H/N/O (Fig. 3c), evidence that it has learned transferable chemistry rather than memorizing the common organic motifs. The hydrocarbon-only predecessor MEHnet [20] is compared per-property in the Supplementary Information (Tables S12 and S13); because MEHnet reports in-distribution RMSE on hydrocarbons against a c-pVDZ reference, it is not directly commensurable with the columns here and is kept out of Table 1. This accuracy comes at the cost of using the B3LYP baseline. Figure 3a reports the per-molecule wall time of a single-point energy at each method’s reference level (serial ORCA; the additional cost of computing the full response-property suite is shown in Supplementary Fig. S8). Fitted as t=ANτt=A\,N^τ, the DLPNO-CCSD(T)/c-pVTZ component of the composite reference is one to two orders of magnitude more expensive than the DFT baselines at every size in the set (prefactor A≈13A≈ 13 s versus 0.30.3–0.50.5 s) and is the steepest-scaling method shown (τ≈1.6τ≈ 1.6), reaching ∼ 10310^3 s on the largest molecules. These fitted exponents are effective slopes over the 55–2424-atom test range and sit well below the familiar formal asymptotics (⁡(N3–4)O(N^3--4) for hybrid DFT, ⁡(N7)O(N^7) for canonical CCSD(T)): at these sizes the wall time is dominated by size-independent overheads (SCF setup, integral generation, I/O) and reduced by integral screening and density fitting (RIJCOSX), and DLPNO-CCSD(T) in particular trades the canonical ⁡(N7)O(N^7) for near-linear asymptotic scaling at a ∼ 30×30× prefactor, appearing as a vertically offset, nearly parallel line rather than a steeper curve. The textbook exponents only emerge at larger sizes: on the 99–4444-atom oligothiophenes the measured canonical-CCSD slope steepens to τ≈4.6τ≈ 4.6 (Supplementary Fig. S7). The model rides its B3LYP/def2-SVP baseline: the neural correction adds only 2424–2828 ms per molecule on a single NVIDIA A100 GPU, essentially flat in size, so its cost tracks the baseline DFT curve, roughly two orders below CCSD(T); the deep-learning functional Skala [19] is comparably cheap (a meta-GGA-cost single point). We report wall-times because they are what a practitioner experiences, but they conflate hardware (ORCA on CPU, the network and Skala on GPU); expressed as accelerator-time the model’s cost is dominated by the classical baseline either way (SI). The model therefore sits at the favourable corner of the accuracy–cost plane (Fig. 1d): coupled-cluster accuracy at hybrid-DFT cost, ∼ orders below CCSD(T) at these sizes and a separation that widens with system size, without bound once the reference turns intractable, as in the larger oligothiophenes below. Accuracy aside, no released machine-learning model even delivers this combination of chemistry and properties. Table 2 maps representative released models onto the benchmark’s nine-element chemistry and property suite: across the three families of ML approaches to molecular properties, none simultaneously spans the elements, outputs the full property suite, and targets a correlated reference. A complementary scope map—which charge states, spin states, and geometry regimes each of these models is trained or validated on—is given in the Supplementary Information. Table 2: What released machine-learning models can deliver on this benchmark’s chemistry and property suite. A capability/coverage map, not a head-to-head error table (for accuracy see Table 1): across the three families of ML approaches to molecular properties, no released model simultaneously spans the nine-element chemistry, outputs the full property suite, and targets a correlated reference. Elements states each released model’s own coverage—several exceed this work’s scope—with, in parentheses, how many of this benchmark’s nine elements (H, C, N, O, F, Si, P, S, Cl) it spans: QHNet lacks Si, P, S, Cl; the released SpookyNet (QM7-X) lacks F, Si, P; MACE-OFF23 lacks Si. Symbols: ✓ validated output of the cited model; ✓ obtainable from the model’s predicted electronic-structure objects (Hamiltonian, density, or charges) but not demonstrated in the cited work; – not accessible. Properties: total/atomization energy E, HOMO–LUMO gap EgE_g, dipole μ, quadrupole Q, polarizability α, atomic charge q, bond order (BO); this work reports Mulliken charges and Mayer bond orders, whereas the population convention differs among the other entries (SpookyNet, e.g., predicts its own learned partial charges). HELM’s weights are announced but, as of this writing, not public (its dataset is); MACE-OFF23 dipoles are validated only in its separate −μ-μ variant; SpookyNet dipoles follow from its predicted partial charges but are not benchmarked. The Reference column is the correlated level the tabulated properties are fit to: for the machine-learned functionals (DM21, Skala) only the energy targets a correlated reference (E), whereas MEHnet v2 targets one for every property (all). Output properties Method Elements E EgE_g μ Q α q BO Reference Public Electronic-structure (Hamiltonian) read-out MEHnet v2 (this work) 9 (9/9) ✓ ✓ ✓ ✓ ✓ ✓ ✓ CCSD(T) (all) ✓ QHNet [12] 5 (5/9) ✓ ✓ ✓ ✓ ✓ ✓ ✓ B3LYP ✓ HELM [13] 58 (9/9) ✓ ✓ ✓ ✓ ✓ ✓ ✓ ω 97M-V × Machine-learned density functional DM21 [18] H–Kr (9/9) ✓ ✓ ✓ ✓ ✓ ✓ ✓ CCSD(T)/exp. (E) ✓ Skala [19] H–Ar (9/9) ✓ ✓ ✓ ✓ ✓ ✓ ✓ CCSD(T) (E) ✓ Direct property (pooling) read-out SpookyNet [14] 6 (6/9) ✓ – ✓ – – ✓ – PBE0+MBD ✓ MACE-OFF23 [25] 10 (8/9) ✓ – ✓ – – – – ω 97M-D3(BJ) ✓ UMA / OMol25 [26] 83 (9/9) ✓ – – – – – – ω 97M-V ✓ Agreement with experiment The benchmarks so far are referenced to a computational ground truth. As an independent check, we compare predicted gas-phase dipole magnitudes against experimentally measured values from the NIST Computational Chemistry Comparison and Benchmark Database (CCCBDB), for twelve molecules spanning the full element scope (HF, HCl, H2S, PH3, SiH4, CH3F, CH3Cl, CHF3, CHCl3, PF3, PCl3, SO2; experimental dipoles 00–1.891.89 Debye). We use the dipole moment because it is sensitive to electronic-structure quality yet insensitive to basis-set incompleteness in the c-pVXXZ series, so the experimental value is a clean target (unlike atomization energies, where c-pVTZ incompleteness dominates; discussed in the SI). Figure 4 shows the parity plot against experiment. Over the twelve molecules the model predicts the measured dipole magnitudes to a mean absolute error of 0.0780.078 Debye—close to the 0.0460.046 Debye that a canonical CCSD(T)/def2-TZVPP calculation achieves against the same experiments, and an improvement over its B3LYP/def2-SVP baseline, at the baseline’s cost. The model thus reproduces a real, independently measured observable across H, C, N, O, F, Si, P, S, and Cl chemistries at DFT cost. Figure 4: Agreement with experiment. Predicted versus measured (NIST CCCBDB) gas-phase dipole magnitudes for twelve molecules spanning the model’s element scope; the line is y=xy=x. Mean absolute error 0.0780.078 Debye (RMSE 0.0980.098 Debye), compared with the B3LYP baseline and 0.0460.046 Debye for a canonical CCSD(T)/def2-TZVPP reference. SiH4 is at the origin by TdT_d symmetry. Out-of-distribution size scaling and its mechanism Figure 5: Out-of-distribution size scaling on oligothiophenes TnT_n (n=1n=1–88, 99–5858 atoms). Both panels share the same equivariant backbone, training data, and local correction; only the read-out differs. a Isotropic polarizability αiso _iso versus size: the Hamiltonian read-out of this work (MEHnet-MG, red) agrees with finite-field CCSD/c-pVDZ (black, dashed) to ∼ 2% everywhere CCSD is computable (n≤6n≤ 6, up to 44 atoms) and continues the trend to 58 atoms, whereas the direct property read-out (blue; the pooling-style Δ -prediction ablation) over-polarizes, drifting toward the BP86 baseline (grey). b Optical gap versus size: MEHnet-MG (red) reproduces the EOM-CCSD/c-pVDZ gap where that reference is affordable (T1T_1–T5T_5, black; within 0.040.04 eV beyond the monomer) and saturates with length; a Kohn–Sham (DFT) gap is not an optical gap and is not shown. The mechanism (frontier-orbital delocalization) and the longitudinal αxx _x are analysed in the Supplementary Information. The benchmarks above probe molecules from the training distribution. A more demanding test is transfer to systems larger than every training molecule and outside the training size distribution: α,α′α,α -linked oligothiophenes TnT_n (n=1n=1–88; thiophene to octithiophene, 99–5858 atoms), the sulfur-rich π-conjugated motif of organic electronics. These are doubly out-of-distribution: larger than the training molecules (≤9≤ 9 heavy and 2424 total atoms, already exceeded by bithiophene T2T_2) and realising extended π-conjugation absent from training, with coupled-cluster references that turn intractable as the chain grows. The isotropic polarizability grows super-linearly along the backbone (Fig. 5a; a ∼ 12-fold rise from T1T_1 to T8T_8, ∼ 21-fold for the longitudinal αxx _x, SI). MEHnet-MG tracks finite-field CCSD/c-pVDZ across the affordable range (signed αiso _iso error within ∼ 2% out to T6T_6; the longitudinal αxx _x stays within ∼ 3%, SI) and continues the corrected trend to T8T_8, while the uncorrected DFT baseline over-polarizes by a length-growing margin (+43%+43\% at T4T_4 for a GGA). The decisive control is the read-out ablation (blue): the same backbone, training data, and local correction, with a direct property read-out in place of the Hamiltonian read-out, reproduces CCSD on short chains but then drifts back toward the BP86 baseline, the deviation widening monotonically to T8T_8. Only the read-out differs, so the divergence is attributable to it alone. The optical gap is the complementary observable (Fig. 5b): where the polarizability grows super-linearly, the gap closes monotonically with conjugation and saturates toward a finite polymer-limit value, an intensive-like quantity that the same read-out also gets right. Our accuracy claim here is deliberately narrow. The gap label is, by construction, the lowest EOM-CCSD/c-pVDZ excitation (Methods), so the like-for-like reference is EOM-CCSD/c-pVDZ on the same geometries; against it the model transfers well out of distribution, reproducing the EOM-CCSD gap within 0.040.04 eV at T2T_2–T5T_5. The T1T_1 monomer deviates by −0.37-0.37 eV; we verified this is not a state-assignment artefact—the lowest EOM root at T1T_1 is cleanly the HOMO→ single excitation (92%92\% singles character), exactly the state the frontier-gap read-out targets—but an ordinary prediction error on the smallest chain, lying in the tail of the in-distribution test error (gap RMSE 0.150.15 eV). Both the model and this reference lie ∼ 1 eV above experimental UV–Vis maxima. That offset originates in the limitations of the EOM-CCSD/c-pVDZ labels, not in the network: c-pVDZ lacks diffuse functions, which alone account for 0.410.41 eV of the offset at T1T_1 (SI), and the comparison further pits vertical gas-phase excitations against solution band maxima. The model’s gap accuracy is therefore bounded by the level of its labels, a point we return to in the Discussion. A network with a finite cutoff RcR_c can nevertheless capture a response that grows over lengths far exceeding RcR_c, because each property is read from the eigenstates of the predicted Hamiltonian, not from local features: although the learned correction Δθ _θ is local, diagonalizing =DFT+ΔθH=H_DFT+ _θ yields frontier orbitals that delocalize over the conjugation length ξ≫Rcξ R_c, and the polarizability is the linear response of these delocalized states (Methods, Eqs. (2)–(3)). This is measured, not assumed: the frontier quantities of the predicted Hamiltonian keep evolving with chain length far beyond the 66 Å cutoff, the HOMO→ transition dipole growing as ∼N0.65 N^0.65 and the gap shrinking as ∼N−0.33 N^-0.33 out to T8T_8 (SI); their combination, αxx∼|⟨ψH|x^|ψL⟩|2/(εL−εH) _x | _H| x| _L |^2/( _L- _H), makes αxx _x grow as ∼N1.6 N^1.6, super-linearly (τ>1τ>1). The response length is set by the global eigenproblem (ξ), not by the network cutoff; the full derivation is in the SI, and a distributed-polarizability decomposition there confirms the real-space picture: 6868–97%97\% of αxx _x is collective inter-atomic charge flow rather than local atomic dipoles. A pooled read-out cannot reach this regime. Summing local contributions frozen beyond RcR_c locks the size-scaling to strictly extensive (τ=1τ=1) or intensive (τ=0τ=0), so a super-linear collective response is structurally out of reach; even charge-weighted (charge × position) dipole read-outs remain capped at linear, and only genuinely non-local schemes (global attention, charge equilibration) escape the bound, but must then learn the size dependence from data they lack at large N (SI). This is exactly the ablation of Fig. 5a: no amount of training data lifts it, because the missing ingredient is the global eigenproblem, not a richer local representation. Discussion The results show that a single equivariant model can place a broad suite of molecular properties at the favourable corner of the accuracy–cost plane, coupled-cluster accuracy, at the cost of a single DFT single point, across the main group. These results rest on the two principles set out in the Introduction. First, the level of supervision sets the accuracy ceiling: MEHnet-MG is trained on the observables, matched to a correlated (CCSD(T)) reference, with the cheap B3LYP/def2-SVP single point serving only as a baseline the network corrects (so it need only represent the smooth, relatively small difference between hybrid DFT and coupled cluster). Because it never fits a reference Hamiltonian matrix, its accuracy is bounded by the reference level, not by DFT, which is what separates it from Hamiltonian-regression surrogates. Second, the architecture sets the size-scaling: the network outputs its correction as a local change to the Hamiltonian, and every property is read off by the exact quantum-mechanical operator acting on the eigenstates of that Hamiltonian. The learnable part (the Hamiltonian correction Δθ _θ) is strictly local and depends only on chemical environments that recur across sizes and chemistries (which is why per-element accuracy does not degrade, and why the model transfers from ≲ 25-atom training molecules to 58-atom chains), but the property is computed globally, by diagonalization and linear response, so its dependence on system size is fixed by quantum mechanics rather than by a pooling choice. An architecture that instead pools per-atom contributions can only produce strictly extensive (sum) or strictly intensive (mean) scaling; it cannot represent a collective response that grows super-linearly with length, and, because that growth lies outside the small-molecule training distribution, no amount of data can supply it. Correct size-scaling is, in this sense, an architectural property, not a learnable one. The oligothiophene case makes this concrete. The network does not merely echo its baseline but actively removes the DFT delocalization error, tracking finite-field CCSD polarizability to ∼ 2% out to the 44-atom CCSD limit, as accurately at T6T_6 as in-distribution, and continuing the corrected super-linear trend beyond it, where even one CCSD field point costs 5.75.7 h (∼ 500× the model’s entire inference; Supplementary Fig. S7). It can do so because the longitudinal response of these chains is dominated by collective, delocalized charge transfer (Supplementary Fig. S3; 68–97% of αxx _x), which the global linear response of the predicted Hamiltonian captures even though the learned correction is local, exactly the contribution a pooling or induced-dipole surrogate would miss. A same-backbone read-out ablation makes this concrete: an otherwise identical model with a direct property read-out drifts back toward the over-polarized baseline as the chain grows, while the Hamiltonian read-out does not (Fig. 5). The model is strongest where its targets are clean, single-reference coupled-cluster properties: total energies, polarizabilities, atomic charges, and bond orders. The clear weak link is the optical gap, whose label is an EOM-CCSD/c-pVDZ excitation; the model reproduces that target faithfully even out of distribution, but the target itself is not a physical gold standard for conjugated systems, and the model inherits its ∼ 1 eV offset from experiment. This is a property of the labels, not the surrogate, and it points directly to how the dataset should evolve (see below). Practically, the method enables coupled-cluster-quality prediction of the six ground-state and response properties at high throughput for the nine main-group elements covered, including the P/S/Cl space relevant to ligands, agrochemicals, flame retardants, and sulfur- and phosphorus-based materials, at a cost dominated by a single cheap DFT call. The optical gap is the exception: it tracks its EOM-CCSD/c-pVDZ label faithfully but inherits that label’s offset from experiment, so the gap should be treated as a relative, label-consistent quantity rather than a substitute for a higher-level spectroscopic prediction. A GPU-accelerated DFT baseline would reduce the residual cost further, but the limiting accuracy is set by the labels, not the surrogate. The model’s accuracy is bounded by the level of theory of its training labels, and our analysis identifies the optical gap as the binding constraint: its EOM-CCSD/c-pVDZ label is not a gold standard for π-conjugated systems, because the small basis lacks diffuse functions, the method omits triples, and the relevant excited states acquire double-excitation character that grows with chain length [27, 28]. The model reproduces this label faithfully, so improving the gap requires improving the labels, not the network. This motivates a concrete dataset roadmap. (i) Add diffuse functions (aug-c-pVDZ or def2-TZVPD) for the response and excited-state targets, the single largest and cheapest accuracy gain; on thiophene this alone recovers 0.410.41 eV of the gap offset. (i) Make the fundamental gap (IP−-EA), obtainable from the EOM-IP/EA-CCSD data we already compute, the headline gap: it is a single-electron process, free of the double-excitation pathology, and robust with system size. (i) Build out-of-distribution holdouts and a graded series of conjugated molecules by design, so transfer is measured rather than discovered post hoc. (iv) Calibrate against a higher-level reference (C3/aug-c-pVTZ) on a few hundred molecules to pin the residual label error, rather than recomputing the entire set at that level. Full coupled-cluster triples for excitations (C3/CCSDT, ⁡(NOPEN7–8))O(N^7--8))) are intractable for systems of this size, which is why we bound the gap ceiling with a basis-set study and literature benchmarks (SI) rather than with explicit triples. Two further limitations are worth naming explicitly. First, the training/validation/test split is by shuffled dataset-index range (Methods) rather than by scaffold or size; while every molecule appears in exactly one set, this does not guarantee absence of close chemical near-neighbours between train and test, and scaffold-split benchmarks would more conservatively measure transfer to genuinely unseen chemistry. Second, our comparison against other machine-learning models is necessarily heterogeneous: we map released checkpoints onto the benchmark’s chemistry and property suite (Table 2) and report accuracy for the machine-learned functional Skala alongside conventional functionals (Table 1), but the released models differ in training set, reference level, and target chemistry, so these are not like-for-like error comparisons. We cleanly isolate the role of the read-out itself through the same-backbone ablation (direct vs. Hamiltonian read-out; Fig. 5); a controlled retraining of competing architectures on our common dataset and CCSD(T) reference would extend this to a full head-to-head accuracy benchmark and is left to future work. Methods Dataset and reference levels The dataset comprises 41,93941,939 closed-shell, neutral main-group molecules (2–24 atoms, ≤9≤ 9 heavy) built from H, C, N, O, F, Si, P, S, and Cl, retained after a three-stage quality filter from 44,41244,412 molecular structures collected from PubChem (one 3D structure per compound; SI). The molecules are partitioned by dataset-index range into training (40,98540,985; indices 20012001–44,41244,412), validation (968968; 10011001–20002000) and test (959959; 11–10001000). Each molecule appears in exactly one set. Because the structures were pseudorandomly shuffled and relabeled 11–44,41244,412 before any reference calculations (SI), this index-range split is random with respect to chemistry and molecule size; it is, however, not scaffold-split or size-stratified, and we therefore do not guarantee absence of close chemical near-neighbours between train and test, a limitation that the out-of-distribution oligothiophene benchmark (Fig. 5) is partly designed to address. Each property is labeled at a property-appropriate level. The electronic energy—and likewise the dipole, quadrupole, Mulliken atomic charges, and Mayer bond orders, taken from the relaxed coupled-cluster density—is a composite estimate of canonical CCSD(T)/c-pVTZ, assembled from three coupled-cluster runs per molecule as Xlabel=XCCSD⁡(T)c-pVDZ+XDLPNOcc-pVTZ−XDLPNOcc-pVDZ,X_label\;=\;X^c-pVDZ_CCSD(T)\;+\;X^c-pVTZ_DLPNO\;-\;X^c-pVDZ_DLPNO, (1) i.e., a canonical CCSD(T)/c-pVDZ value plus a DLPNO-CCSD(T) [29] basis-set correction from c-pVDZ to c-pVTZ; we write “composite CCSD(T)/c-pVTZ” for this level throughout. The isotropic and full polarizability is labeled from finite-field CCSD/c-pVDZ (numerical second derivative of the energy with respect to a static field, by 7-point central differences); the lowest singlet vertical excitation (the gap EgE_g) from EOM-CCSD/c-pVDZ [30] (lowest root); and ionization potentials and electron affinities are from EOM-IP/EA-CCSD/c-pVDZ. Two basis-set-dependence caveats apply: (i) the polarizability reference is c-pVDZ without diffuse functions and is therefore itself biased relative to a complete-basis CCSD limit by a few percent (SI), so “∼ 2% of CCSD” should be read as “∼ 2% of this systematically biased reference”; (i) Mulliken charges and Mayer bond orders are basis-dependent partitioning schemes: the B3LYP/def2-SVP baseline and the composite CCSD(T)/c-pVTZ target are strictly different quantities, and the model learns the correction between them rather than predicting an absolute basis-independent “physical” partial charge. The Δ -learning baseline (and the baseline Fock matrix DFTH_DFT that the network corrects) is a single B3LYP/def2-SVP calculation. All reference calculations were performed with ORCA 6.0 [31]; protocols and input templates are in the SI. Model architecture Rather than reading each target property out of pooled atomic features, the model predicts an effective one-electron Hamiltonian and derives every property from it by the same quantum-mechanical operator that defines it in the reference calculation. This extends the molecular-Hamiltonian-learning formulation of Tang et al. [20] (demonstrated there for hydrocarbons) to the nine-element main group, and is the structural source of the model’s size transferability (Fig. 5, Discussion). The backbone is an SO(3)-equivariant, spherical-harmonic graph-attention network (EquiformerV2 [21], in the eSCN lineage [32]): four interaction layers, a radial cutoff of 6.06.0 Å with up to 2020 neighbours per atom, 6464 spherical channels, 88 attention heads, spherical-harmonic degree lmax=mmax=4l_ =m_ =4—the degree required to represent the d×d× d on-site Fock blocks that the third-row elements (Si, P, S, Cl) introduce—and RMS spherical-harmonic normalization with gate activations (∼ 5.6 M parameters). Its equivariant node features are mapped by e3n [22] heads, per element for on-site blocks and pairwise (combined with edge spherical harmonics) for off-site blocks, onto a correction Δθ _θ to the baseline Fock matrix in the atomic-orbital basis. The effective Hamiltonian is =DFT+ΔθH=H_DFT+ _θ, where DFTH_DFT is the B3LYP/def2-SVP Kohn–Sham matrix of the cheap baseline single point; solving the eigenvalue equations of H yields molecular-orbital energies εn\ _n\ and coefficients, hence the one-particle density matrix. Every property is then evaluated from H and its eigenstates by the defining operator, not by a separate learned head: the total electronic energy as the sum over occupied levels; the gap as the frontier (HOMO–LUMO) eigenvalue difference; the dipole and quadrupole as the multipole operators traced against the density; the polarizability by linear response (detailed below); and Mulliken atomic charges and Mayer bond orders from the density matrix. The polarizability is built in two steps from the eigenpairs εn,ψn\ _n, _n\ of H. A bare independent-particle sum-over-states response, αxy0=∑i∈occ∑a∈virt⟨ψi|x^|ψa⟩⟨ψi|y^|ψa⟩εa−εi,α^0_xy=4\!\! _i \; _a _i|\, x\,| _a \, _i|\, y\,| _a _a- _i, (2) is screened (depolarization) by the screening matrix =DFT+ΔθT=T_DFT+ _θ to give the reported tensor, =(+0)−10, α= (I+ α^0T )^-1 α^0, (3) where the baseline screening matrix DFTT_DFT is obtained from a linear-response (TDDFT) calculation on the baseline (so an uncorrected model reproduces the baseline polarizability; this TDDFT step is needed only when the polarizability is requested) and Δθ _θ is the learned correction. Because the read-out is the exact quantum-mechanical functional of a Hamiltonian assembled from a local, finite-cutoff correction, intensive and extensive properties inherit their correct system-size dependence by construction, rather than from a pooling choice (sum, which forces strict extensivity, or mean, which forces strict intensivity). This is what allows the model to extrapolate the size-dependence of collective properties beyond the training set (Fig. 5). With the B3LYP baseline the Kohn–Sham gap already tracks the correlated optical gap with near-unit slope, so, unlike a GGA baseline, no multiplicative gap rescaling is needed. Training The model was trained on the full 40,98540,985-molecule training partition by distributed data parallelism across 16 NVIDIA A100 GPUs (4 nodes), minimizing a weighted sum of per-property mean-squared-error losses (the per-property weights balance the disparate property scales; values in the released config), using the Adam optimizer [33] with a step-decay schedule taking the learning rate from 3×10−33× 10^-3 to 10−410^-4 over the 2,0002,000-epoch schedule. Regularization comprised attention dropout (0.10.1) and stochastic depth (drop-path 0.050.05). Each epoch is one pass over the training partition in 100100 minibatches (∼ 410 molecules per global batch); training ran for 2,0002,000 epochs (2×1052× 10^5 gradient steps) in ∼ 72 h of wall-clock time, and the checkpoint with the lowest validation loss is used throughout. Full hyperparameters are in the released configuration file. Training stability and level tracking. A Hamiltonian read-out introduces a failure mode absent from property-head models: as the learned correction Δθ _θ shifts the spectrum during training, two molecular orbitals can cross, and a strict energy ordering then reassigns which orbital is the HOMO discontinuously, producing a jump in the density, and in every response property derived from it, that destabilizes the gradient. We addressed this with a diabatic level-tracking scheme: at each step the occupied subspace is identified not by the nen_e lowest eigenvalues but by maximal overlap of each eigenstate with the baseline occupied density, so the occupied/virtual assignment varies continuously through a crossing. This is distinct from, and complementary to, the perturbation-theory treatment of near-degenerate eigenvalue gradients [20]: it removes the discontinuity in the forward assignment, not the singularity in the backward pass. Level tracking was essential for stable training on the unfiltered crawl; once the dataset was quality-filtered (Methods, SI), surviving crossings were rare enough that plain energy ordering sufficed, and the released model dispenses with it. We document the scheme because it is a cheap, general remedy for training Hamiltonian-prediction models on noisier data. Benchmark systems and timing Oligothiophene, chloroalkane, and chloropolyene geometries were generated programmatically (planar, idealized; SI). Reference properties on these geometries used ORCA 6.0 [31]: analytic BP86 polarizabilities via the ELPROP module, finite-field CCSD/c-pVDZ polarizabilities (7-point central differences), EOM-CCSD/c-pVDZ excitations (lowest three roots), and BP86 and B3LYP single points (the latter is the model’s baseline). The distributed-polarizability decomposition (Supplementary Fig. S3) partitions the analytic finite-field DFT response density via the ELPROP module. Neural-network forward-pass timings were measured on a single NVIDIA A100 (40 GB) with batch size 1, five warmup passes, and the minimum of three timed passes per molecule; DFT-baseline wall-times are the ORCA TOTAL RUN TIME on 7 CPU cores at dataset-generation time. Supplementary Information Supplementary Information is provided as a separate document (SI.tex), covering: the chloroalkane/chloropolyene companion size-extrapolation; the full inference-timing analysis; the basis-set ceiling of the optical-gap reference (c-pVDZ vs. aug-c-pVDZ); the atomization-energy/basis- incompleteness decomposition; per-property and per-element error tables; and full computational protocols. Data availability The 959959-molecule held-out test set—the three-dimensional structures together with every coupled-cluster label reported here (composite CCSD(T)/c-pVTZ, finite-field CCSD/c-pVDZ and EOM-CCSD/c-pVDZ)—underlies all quantitative results in this work and is openly available in a Zenodo archive whose DOI will be minted upon publication. The same record provides the train/validation/test and scaffold split definitions and the PubChem identifiers of all 44,41244,412 source structures, so the full corpus can be reconstructed from public inputs using the released generation pipeline (see Code availability). The coupled-cluster labels for the 41,93941,939-molecule training and validation partitions are not redistributed with that record; they can be regenerated from the released identifiers and pipeline, and are available from the corresponding author on reasonable request. Code availability The model implementation, the training and inference pipeline, the trained MEHnet-MG weights, the ORCA input templates and data-generation scripts, and the three-stage quality-control and filtering code are openly available at https://github.com/He-Wenhao/ML_electronic_Wenhao under an OSI-approved open-source license; a versioned release will be archived on Zenodo upon publication. Acknowledgements This work was supported by the Honda Research Institute. N.S. and Y.W. acknowledge support from U.S. Department of Energy, Office of Science, Basic Energy Sciences, under Early Career Award No. DE-SC0024524. This research used resources of the National Energy Research Scientific Computing Center (NERSC), a U.S. Department of Energy Office of Science User Facility at Lawrence Berkeley National Laboratory, operated under Contract No. DE-AC02-05CH11231 (the Perlmutter system), which provided the primary computational resources for this work. Additional computing resources were provided by the Frontera system at the Texas Advanced Computing Center (TACC) at The University of Texas at Austin. Author contributions W.H. designed and implemented the model, performed the experiments and analysis, and wrote the manuscript. H.T. contributed to the methodology. X.C. implemented the distributed-training parallelization and optimized GPU utilization. N.S. generated the training dataset. T.S.H. optimized observable back-propagation pipelines. Z.L. and B.L. studied the alignment of the model with experiment. Y.Y. contributed the neural-network training methodology. F.L., Y.W. and J.L. provided computational resources, supervised the project, acquired funding, and edited the manuscript. A.R.H. co-initiated the research theme, co-formulated the research goals, and acquired funding. All authors discussed the results and commented on the manuscript. Competing interests The authors declare no competing interests. References [1] Stefano Curtarolo, Gus L. W. Hart, Marco Buongiorno Nardelli, Natalio Mingo, Stefano Sanvito, and Ohad Levy. The high-throughput highway to computational materials design. Nature Materials, 12:191–201, 2013. [2] Keith T. Butler, Daniel W. Davies, Hugh Cartwright, Olexandr Isayev, and Aron Walsh. Machine learning for molecular and materials science. Nature, 559:547–555, 2018. [3] O. Anatole von Lilienfeld, Klaus-Robert Müller, and Alexandre Tkatchenko. Exploring chemical compound space with quantum-based machine learning. Nature Reviews Chemistry, 4:347–358, 2020. [4] Krishnan Raghavachari, Gary W. Trucks, John A. Pople, and Martin Head-Gordon. A fifth-order perturbation comparison of electron correlation theories. Chemical Physics Letters, 157:479–483, 1989. [5] W. Kohn and L. J. Sham. Self-consistent equations including exchange and correlation effects. Physical Review, 140:A1133–A1138, 1965. [6] Aron J. Cohen, Paula Mori-Sánchez, and Weitao Yang. Insights into current limitations of density functional theory. Science, 321:792–794, 2008. [7] Paula Mori-Sánchez, Aron J. Cohen, and Weitao Yang. Many-electron self-interaction error in approximate density functionals. The Journal of Chemical Physics, 125:201102, 2006. [8] He Li, Zun Wang, Nianlong Zou, Meng Ye, Runzhang Xu, Xiaoxun Gong, Wenhui Duan, and Yong Xu. Deep-learning density functional theory Hamiltonian for efficient ab initio electronic-structure calculation. Nature Computational Science, 2:367–377, 2022. [9] Xiaoxun Gong, He Li, Nianlong Zou, Runzhang Xu, Wenhui Duan, and Yong Xu. General framework for E(3)-equivariant neural network representation of density functional theory Hamiltonian. Nature Communications, 14:2848, 2023. [10] K. T. Schüt, M. Gastegger, A. Tkatchenko, K.-R. Müller, and R. J. Maurer. Unifying machine learning and quantum chemistry with a deep neural network for molecular wavefunctions. Nature Communications, 10:5024, 2019. [11] Oliver T. Unke, Mihail Bogojeski, Michael Gastegger, Mario Geiger, Tess Smidt, and Klaus-Robert Müller. SE(3)-equivariant prediction of molecular wavefunctions and electronic densities. In Advances in Neural Information Processing Systems (NeurIPS), volume 34, pages 14434–14447, 2021. [12] Haiyang Yu, Zhao Xu, Xiaofeng Qian, Xiaoning Qian, and Shuiwang Ji. Efficient and equivariant graph networks for predicting quantum Hamiltonian. In International Conference on Machine Learning (ICML), volume 202 of PMLR, 2023. [13] Manasa Kaniselvan, Benjamin Kurt Miller, Meng Gao, Juno Nam, and Daniel S. Levine. Learning from the electronic structure of molecules across the periodic table, 2025. [14] Oliver T. Unke, Stefan Chmiela, Michael Gastegger, Kristof T. Schüt, Huziel E. Sauceda, and Klaus-Robert Müller. SpookyNet: Learning force fields with electronic degrees of freedom and nonlocal effects. Nature Communications, 12:7273, 2021. [15] Kristof T. Schüt, Oliver T. Unke, and Michael Gastegger. Equivariant message passing for the prediction of tensorial properties and molecular spectra. In International Conference on Machine Learning (ICML), 2021. [16] Dylan M. Anstine, Roman Zubatyuk, and Olexandr Isayev. AIMNet2: a neural network potential to meet your neutral, charged, organic, and elemental-organic needs. Chemical Science, 2024. [17] Raghunathan Ramakrishnan, Pavlo O. Dral, Matthias Rupp, and O. Anatole von Lilienfeld. Big data meets quantum chemistry approximations: The δ-machine learning approach. Journal of Chemical Theory and Computation, 11:2087–2096, 2015. [18] James Kirkpatrick, Brendan McMorrow, David H. P. Turban, Alexander L. Gaunt, James S. Spencer, Alexander G. D. G. Matthews, Annette Obika, Louis Thiry, Meire Fortunato, David Pfau, Lara Román Castellanos, Stig Petersen, Alexander W. R. Nelson, Pushmeet Kohli, Paula Mori-Sánchez, Demis Hassabis, and Aron J. Cohen. Pushing the frontiers of density functionals by solving the fractional electron problem. Science, 374:1385–1389, 2021. [19] Microsoft Research AI for Science. Accurate and scalable exchange-correlation with deep learning, 2025. Skala exchange–correlation functional; https://github.com/microsoft/skala. [20] Hao Tang, Brian Xiao, Wenhao He, Pero Subasic, Avetik R. Harutyunyan, Yao Wang, Fang Liu, Haowei Xu, and Ju Li. Approaching coupled-cluster accuracy for molecular electronic structures with multi-task learning. Nature Computational Science, 5(2):144–154, 2025. [21] Yi-Lun Liao, Brandon Wood, Abhishek Das, and Tess Smidt. EquiformerV2: Improved equivariant transformer for scaling to higher-degree representations. In International Conference on Learning Representations (ICLR), 2024. [22] Mario Geiger and Tess Smidt. e3n: Euclidean neural networks. arXiv:2207.09453, 2022. [23] R. S. Mulliken. Electronic population analysis on LCAO–MO molecular wave functions. I. The Journal of Chemical Physics, 23:1833–1840, 1955. [24] I. Mayer. Charge, bond order and valence in the ab initio SCF theory. Chemical Physics Letters, 97:270–274, 1983. [25] Dávid P. Kovács, J. Harry Moore, others, and Gábor Csányi. MACE-OFF: Transferable short-range machine learning force fields for organic molecules. Journal of the American Chemical Society, 2025. [26] Brandon M. Wood, Muhammed Shuaibi, Adeesh Kolluru, et al. UMA: A family of universal models for atoms, 2025. [27] Pierre-François Loos, Martial Boggio-Pasqua, Anthony Scemama, Michel Caffarel, and Denis Jacquemin. Reference energies for double excitations. Journal of Chemical Theory and Computation, 15:1939–1956, 2019. [28] Mark A. Watson and Garnet Kin-Lic Chan. Excited states of butadiene to chemical accuracy: Reconciling theory and experiment. Journal of Chemical Theory and Computation, 8:4013–4018, 2012. [29] Christoph Riplinger and Frank Neese. An efficient and near linear scaling pair natural orbital based local coupled cluster method. The Journal of Chemical Physics, 138:034106, 2013. [30] John F. Stanton and Rodney J. Bartlett. The equation of motion coupled-cluster method. a systematic biorthogonal approach to molecular excitation energies, transition probabilities, and excited state properties. The Journal of Chemical Physics, 98:7029–7039, 1993. [31] Frank Neese. Software update: The ORCA program system—version 5.0. WIREs Computational Molecular Science, 12:e1606, 2022. [32] Saro Passaro and C. Lawrence Zitnick. Reducing SO(3) convolutions to SO(2) for efficient equivariant GNNs. In International Conference on Machine Learning (ICML), 2023. [33] Diederik P. Kingma and Jimmy Ba. Adam: A method for stochastic optimization. In International Conference on Learning Representations (ICLR), 2015.

Instant research discovery