Is physics necessary for building foundation models for chemistry or data is all you need?

Our Orbitall foundation model uses 35× less molecular training data and is a 50× smaller model than frontier UMA model but outperforms and is 100x faster and accurate in reactions with solvents.

https://arxiv.org/abs/2507.03853

Unlike language and image models, generating training data of larger chemical systems through density functional theory (DFT) is incredibly expensive.

Instead Orbitall trains on smaller systems and incorporates physics as: (1) input orbital features using cheaper semi-empirical calculations (2) symmetries.

Orbital features natively incorporate spin and charge allowing us to model open shell systems cleanly and use cheaper implicit solvation methods that UMA is unable to.

This allows it to beat Meta’s UMA model in chemically important solvent-based reactions and be 100x faster. OrbitAll accurately reproduces solvent-dependent reaction energetics and transition-state structures.

Orbitall is a molecular foundation model with physics-grounded representations that handle charges and spin states and environmental effects such as solvation and external electric fields in an elegant manner.

Foundation models for science are not built blindly on data and certainly not just with LLMs. Physics is key!

A Unified Quantum Mechanical Representation Deep Learning Framework for All Molecular Systems

Source: https://arxiv.org/html/2507.03853 \equalcont

These authors contributed equally to this work.

\equalcont

These authors contributed equally to this work.

[1]\fnmWilliam A.\surGoddard III

[2]\fnmAnima\surAnandkumar

1]\orgdivDivision of Chemistry and Chemical Engineering,\orgnameCaltech,\orgaddress\street1200 E California Blvd,\cityPasadena,\postcode91125,\stateCA,\countryUSA

2]\orgdivDivision of Engineering and Applied Science,\orgnameCaltech,\orgaddress\street1200 E California Blvd,\cityPasadena,\postcode91125,\stateCA,\countryUSA

\fnmVignesh C.\surBhethanabotla\fnmAmin\surTavakoli\fnmMaurice D.\surHanisch\fnmArimitsu\surHorikawa-Strakovsky\fnmMiguel\surNouman\fnmDanish\surKhan[email protected][email protected][[

Abstract

We introduce OrbitAll, a geometry- and physics-informed deep learning framework that encodes any molecular system with arbitrary charges, spins, and environmental effects using electronic structure information. It utilizes spin-polarized orbital features from the underlying quantum mechanical method and combines them with SE(3)-equivariant graph neural networks. OrbitAll demonstrates superior performance and generalization in predicting charged, open-shell, and solvated molecules, and robustly extrapolates to molecules significantly larger than the training data. OrbitAll achieves chemical accuracy using 10 times fewer training data than competing AI models, with approximately 103– 104speedup compared to density functional theory. Trained on a chemically diverse dataset, OrbitAll performs robustly on challenging molecular systems, and outperforms the foundational machine-learned interatomic potential, UMA, for highly charged species, despite using 35 times less molecular data and a 50-times-smaller model. After learning solvent effects, it accurately predicts solvent-dependent reaction pathways at about 100 times lower cost than explicit-solvation simulations using UMA.

keywords:

Quantum Chemistry, Quantum Mechanics, Deep Learning, DFT, Graph Neural Network, Molecular Representation, Electronic Structure

1Introduction

Quantum mechanical (QM) methods, such as density functional theory (DFT), can accurately simulate atomic and molecular systems[cohen2012challenges,pereira2017machine,mood2020methyl], but require an enormous computational cost[frison_compare_dft&semi_2008,engel2011dft,friesner2005abinitio]. Machine learning (ML)-based methods enable high-throughput atomic simulations at a fraction of the computational cost, while maintaining relative quantum chemical accuracy[fedik_ml4molprops_2022,fooshee2018deep,musaelian2023allegro]. These ML models, after training on relevant datasets, can be orders of magnitude faster than traditional quantum mechanical simulations while still delivering reasonably accurate predictions[musaelian2023allegro,qiao_orbnet_2020,christensen2021orbnetdenali,batzner2022nequip,tavakoli2025chemically,tavakoli2022quantum]. However, the use of these ML models is often limited to the chemistry of their training data.

Foundation models for machine-learned interatomic potentials (MLIPs), such as UMA[wood2025uma], have demonstrated strong accuracy across a broad range of molecular systems, including those with varying charge and spin states, but with crucial limitations. First, the predictive performance of these models is mainly bounded by the chemical domains represented in their training data. Their strong generalization therefore reflects interpolation across sufficiently diverse training distributions, rather than a principled ability to extrapolate to regimes that are sparsely sampled or absent from the data. As a result, such models may fail in the extreme or unusual cases where physical guidance is most needed. Second, current models often do not natively incorporate environmental effects, such as implicit solvation or external electric fields, which can substantially alter molecular electronic structure and reactivity.

Realistic application to all chemical systems with a generalization capability on par with QM calculations requires taking additional degrees of freedom into account: spin, charge, and environmental effects[faglioni2016battery2,das2024battery3]. Integrating spin effects enables applicability to open-shell systems, where unpaired electrons significantly impact the electronic structure and reaction pathways[borden_diradical_2017,tavakoli2023rmechdb,tavakoli2024ai,kang2024geometry,kang2024orbnetspin]. For example, transition metal complexes and radicals require open-shell calculations to capture spin-related effects accurately[maurer2021organometallic,gutowski2013radicals]. Integrating charge effects allows for the accurate modeling of charged species, which are essential in many chemical systems[he2019dftbattery1]. In biochemistry, charged states play a crucial role in enzyme catalysis and protein-ligand interactions, where protonation states and ionic interactions influence binding affinities and reaction mechanisms[dudev2014competition,roberts2015specific]. Further, integrating environmental effects, such as solvation and external electric fields, enhances the modeling of molecular interactions under realistic conditions. These factors can alter electronic properties and reaction energetics, influencing processes such as molecular solvation and field-assisted catalysis[liang2017efield,zhao2011solvation,wei2000solvation].

Incorporating QM information into molecular representations enhances molecular modeling, resulting in highly accurate predictive models.[qiao_orbnet_2020,qiao_orbnet_equi_2022,karandashev2022qml,cheng2019mobml1,lu2022mobml2,cheng2019mobml3,cheng2022mobml4,cheng_MOBML_2022,venturella2025mlgf]. Compared to purely data-driven or even geometry-based models, QM-informed models offer improved generalization and data-efficiency, mitigating the high computational cost of data generation for training predictive models[qiao_orbnet_equi_2022,rama_delta_learning_2015,karandashev2022qml,fooshee2018deep,irwin2022chemformer]. For example, OrbNet achieved chemical accuracy on relative conformer energies for molecules much larger than those of the training set, demonstrating strong size extrapolation[qiao_orbnet_2020]. Furthermore, leveraging QM-informed features—a strategy we callorbital learning—enables models to more effectively capture underlying physics by implicitly encoding spin, charge, and environmental effects[qiao_orbnet_equi_2022,cheng2019mobml1,fabrizio2022spa,briling2024spahm,venturella2025mlgf].

Figure 1:Schematic of the OrbitAll framework.(A)Illustration of inputs and outputs. OrbitAll can process all possible molecular systems for quantum mechanical calculations with different total spin, total charge, and environments. With this molecular system information, OrbitAll predicts molecular properties.(B)Overall architecture of OrbitAll. First, the molecular system’s information is processed by a semi-empirical QM method, spGFN1-xTB[neugebauer_spgfnxtb_2023], to create the spin-polarized orbital features set,T=(Fα,Fβ,Pα,Pβ,S,Hcore)\textbf{T}=(\textbf{F}^{\alpha},\textbf{F}^{\beta},\textbf{P}^{\alpha},\textbf{P}^{\beta},\textbf{S},\textbf{H}_{\text{core}}). This set of orbital features can represent all molecular systems for quantum mechanical simulations. Second, the orbital features, which are SE(3)-equivariant, are processed by an E(3)-equivariant GNN. The target we predict with the GNN is the delta-label[rama_delta_learning_2015],Δ​y\Delta y, the difference between the low-level approximation obtained from spGFN1-xTB,yspGFN1-xTBy_{\text{spGFN1-xTB}}, and the high-level label,yy.(C)Illustrations of the diagonal reduction and the message passing operations using the QMMs. The diagonal reduction encodes the block-diagonals of QMMs into atom-wise encodings (representations). The messages are created by convolutions of the atom-wise representations with the off-block diagonals, then are aggregated with attention to update the atom-wise representations. Further details can be found in Section4.Our approach:We introduce “OrbitAll”, a geometry- and physics-informed deep learning framework for representing all molecular systems with arbitrary charge, spin, and environmental effects (Figure1). OrbitAll begins by representing electronic structure using a low-cost semi-empirical QM method to produce SE(3)-equivariant orbital features that encapsulate converged mean-field electronic interactions. These quantum mechanical features are then mapped onto atomic graphs to retain both geometric and electronic information. Then, an E(3)-equivariant graph neural network (GNN) with orbital-based message-passing processes these features to predict target properties. Therefore, both the representations and the neural network are grounded in orbital features, closely aligning with the underlying physics and consequently enhancing generalization and extrapolation capability.

OrbitAll is the first orbital-based deep learning framework that uses a single unified representation to jointly model all molecular systems across varying charge, spin, and environmental conditions. In contrast, most state-of-the-art geometric GNNs support only a subset of these features, often requiring separate extensions or input modifications to address them individually.

We conduct a large-scale training of OrbitAll for a general-purpose model, trained on the OMol25-4M dataset[levine2025openmolecules2025omol25]. This dataset is a 4-million-molecule subset of the larger OMol25 dataset, which consists of 140 million molecules and was used to train the UMA foundation model, the state-of-the-art MLIP[wood2025uma]. Although OrbitAll is trained on a substantially smaller dataset (35 times smaller for molecular data and 130 times smaller for the total) and is smaller in size (50 times smaller) compared to UMA, it remains robust for complex and underrepresented species, including highly charged ions, and outperforms UMA on challenging subsets despite.

We fine-tune OrbitAll on an implicit-solvation dataset to construct a solvent-effect-informed interatomic potential. The unified representation, which enables the incorporation of field effects, allows OrbitAll to process molecules under implicit solvation in different solvents, an approach not directly applicable to conventional MLIPs such as UMA, and those models require expensive explicit solvation to account for solvent effects. Using the nudged elastic band method, OrbitAll accurately captures solvent-dependent transition-state structures and reaction energetics at about 100 times lower cost than explicit-solvation umbrella sampling with UMA.

We further do ablations and evaluations of the performance of energy prediction tasks using carefully selected datasets, including QM9star (for varying spins and charges)[tang2024qm9star]and Hessian QM9 (for varying implicit solvents)[williams2025hessianqm9]. We first compare model accuracy across various molecular property predictions, showing that OrbitAll consistently outperforms all augmented models in total energy and frontier orbital energy estimation. Next, the cost-accuracy tradeoff and data-efficiency of the models are compared. The result reveals that OrbitAll achieves a mean absolute error below chemical accuracy (1 kcal/mol) using only 10% of the data required to train the next-best model, while being∼103−104\sim 10^{3}-10^{4}faster than DFT/B3LYP during inference.

We also examine the generalization capability of the models to molecular systems of much larger sizes using a dataset of polypeptides. We observe that OrbitAll robustly outperforms competing models for these much larger molecules. Additionally, we demonstrate OrbitAll’s ability to predict the molecular energies for various implicit solvents using Hessian QM9, which again exhibits remarkable generalization across different solvent environments.

Summarizing, OrbitAll overcomes many crucial limitations of current state-of-the-art foundational MLIPs, which remain constrained by the chemical domains represented in their training data. Furthermore, OrbitAll can natively incorporate environmental effects, such as implicit solvation or external electric fields, that are missed by current MLIPs, which have substantial influence on molecular electronic structure and reactivity.

2Results

2.1The OrbitAll Framework

Incorporating molecular geometry into predictive models significantly improves the accuracy of molecular property prediction[schutt2018schnet,Gasteiger_dimenet_2020]. Moreover, designing models that respect molecular symmetries—through group equivariance—further enhances both accuracy and data efficiency[thomas2018tensor,batzner2022nequip]. Notable improvement is also achieved by orbital learning, where the use of physics-informed features (orbital features) enhances the accuracy and fidelity of predictions[qiao_orbnet_2020,qiao_orbnet_equi_2022,cheng_MOBML_2022,karandashev2022qml]. The OrbitAll framework extends these methods to spin-polarizable systems, which enables processing all molecular systems. As shown in Figure1, the OrbitAll framework accepts as inputs atomic numbers (Z), atomic coordinates (R), total charge (QQ), total spin (SS), and parameters regarding environmental effects, such as implicit solvation. The spin-polarized orbital features,T=(Fα,Fβ,Pα,Pβ,S,Hcore)\textbf{T}=(\textbf{F}^{\alpha},\textbf{F}^{\beta},\textbf{P}^{\alpha},\textbf{P}^{\beta},\textbf{S},\textbf{H}_{\text{core}}), are generated using spGFN1-xTB[neugebauer_spgfnxtb_2023]or g-xTB[froitzheim2025gxtb]. Here,Fα\textbf{F}^{\alpha}andFβ\textbf{F}^{\beta}denote the Fock matrices for up-spin (α\alpha) and down-spin (β\beta), respectively;Pα\textbf{P}^{\alpha}andPβ\textbf{P}^{\beta}are the corresponding density matrices;Sis the overlap matrix; andHcore\textbf{H}_{\text{core}}is the core Hamiltonian matrix. We refer to these as quantum mechanical matrices (QMMs), which collectively represent the converged mean-field electronic structure of the molecular system.

The orbital features that represent the electronic structure of the input system are perturbed correspondingly to the various conditions, since the low-level QM method (e.g., spGFN1-xTB) can physically capture the major effects. Thus, this representation can distinguish molecules with different spin states (singlets, doublets, triplets, etc.), different charges (neutral, anion, cation), and different environmental effects (uniform external electric fields, various implicit solvents, etc.) in a common representation space.

In this study, the set of orbital features used areT=(Fα,Fβ,Pα,Pβ,S,Hcore)\textbf{T}=(\textbf{F}^{\alpha},\textbf{F}^{\beta},\textbf{P}^{\alpha},\textbf{P}^{\beta},\textbf{S},\textbf{H}_{\text{core}}). Each element of a QMMOgenerated from an operator𝒪^\hat{\mathcal{O}}is obtained by,

(O)A​Bn,l,m;n′,l′,m′=⟨ΦAn,l,m|𝒪^|ΦBn′,l′,m′⟩,(\textbf{O})^{n,l,m;n^{\prime},l^{\prime},m^{\prime}}_{AB}=\langle\Phi^{n,l,m}_{A}|\hat{\mathcal{O}}|\Phi^{n^{\prime},l^{\prime},m^{\prime}}_{B}\rangle,(1) wherenn,ll, andmmare the principal, angular, and magnetic quantum numbers, respectively[szabo_quantum_1989].

OrbitAll employs an E(3)-equivariant neural network framework, UNiTE[qiao_orbnet_equi_2022], which serves as the backbone of the geometric GNN module. OrbitAll adheres to SE(3)-equivariance overall based on the predictable change of QMM elements, which are SE(3)-equivariant. For a roto-translational transformationℛ\mathcal{R},

(ℛ⋅O)A​Bl;l′=Dl​(ℛ)​(O)A​Bl;l′​(Dl′​(ℛ))†,(\mathcal{R}\cdot\textbf{O})^{l;l^{\prime}}_{AB}=\textbf{D}^{l}(\mathcal{R})(\textbf{O})^{l;l^{\prime}}_{AB}\left(\textbf{D}^{l^{\prime}}(\mathcal{R})\right)^{\dagger},(2) whereOA​Bl;l′\textbf{O}^{l;l^{\prime}}_{AB}is a block in the QMMOthat represents the interaction between the atomic orbital of atomAAwith an angular momentumlland the atomic orbital of atomBBwith an angular momentuml′l^{\prime}. The dagger symbol denotes the Hermitian conjugate, andDl​(ℛ)\textbf{D}^{l}(\mathcal{R})is the Wigner-D matrix of degreell, an irreducible representation of SO(3).

Within the OrbitAll framework, the self-consistent field (SCF) method used to create orbital features yields the quantum mechanical properties at that level of theory without any additional cost. This low-level approximation enables delta-learning (orΔ\Delta-learning,Δ\Delta-ML) strategy[rama_delta_learning_2015], where the task changes into predicting the difference of the original label (higher-level theory or experimental) from the lower-level approximation. With a delta-learning objective, the ML model predicts the delta-labelΔ​y\Delta y, defined by:

Δ​y=ytarget−ylow-level,\Delta y=y_{\text{target}}-y_{\text{low-level}},(3) whereytargety_{\text{target}}is the target label for prediction andylow-levely_{\text{low-level}}is the low-level approximation (e.g., spGFN1-xTB) of the quantity. The delta-learning strategy has been empirically shown to improve the accuracy and data-efficiency of several different labels, since the low-level approximation can capture essential electronic interactions and the resulting delta energy surface is easier to learn[qiao_orbnet_equi_2022,rama_delta_learning_2015,ruth2022deltalearning1,chen2023deltalearning2,zhu2019deltalearning4]. Here, for comparison, we refer to direct predictions of target labels as direct-learning to highlight the differences with delta-learning.

2.2Diverse Chemistry and Real World Applications

Large-Scale OrbitAll for Diverse Chemistry

We train a large-scale OrbitAll as a QM-informed molecular foundational model on the OMol25 dataset[levine2025openmolecules2025omol25], which contains approximately 4 million single-point hybrid DFT calculations at theω\omegaB97M-V/def2-TZVPD level of theory. The dataset spans chemically diverse systems, various elements, highly non-equilibrium structures, and complex chemistries such as transition metal complexes. The resulting model, OrbitAll-OMol25-4M, is trained using orbital features generated by g-xTB v2.0.0 as the internal QM method[froitzheim2025gxtb], with g-xTB energies and forces used as the delta-learning baseline.

OrbitAll-OMol25-4M achieves an MAE of 62.72 meV (∼\sim1.45 kcal/mol) and a per-atom MAE of 1.10 meV/atom on the evaluation set, comparable to those reported for state-of-the-art benchmark models[levine2025openmolecules2025omol25]. The model also shows strong generalization across different chemical domains, including biomolecules, electrolytes, and metal complexes (AppendixS3.1).

Most importantly, OrbitAll-OMol25-4M performs robustly in extreme and underrepresented charged regimes. It achieves substantially lower MAEs than UMA-s-1p2 for ions with charges of +5 and above or -4 and below, except at -10 (FigureS9), despite a substantially smaller training set: 35-fold smaller for finite molecules and 130-fold smaller overall. Such robustness is particularly valuable for reactive and chemically diverse systems, where unusual charge states can be important.

Solvent-Aware Reaction Pathway Prediction

A key application of MLIP is the construction of a minimum-energy pathway (MEP). In realistic settings, however, MEP calculations must often account for solvent effects, which remain challenging for existing MLIPs to capture directly. A common workaround is to include explicit solvent molecules using the existing MLIPs, but this is often highly sensitive to solvent configuration, hence requiring sampling methods for reliable predictions such as umbrella sampling. OrbitAll, on the other hand, can encode environmental effects, which provide a more practical solution for solvent-aware MEP construction. We discuss this further in Section2.4.

To build a practical solvent-aware reactive MLIP, we create the T1x-Solv dataset with reaction geometries evaluated under four solvent conditions (vacuum, water, methanol, and toluene), and fine-tune the pretrained OrbitAll-OMol25-4M model on T1x-Solv to leverage learned representations. Further details on the dataset are provided in Section4.2.

Figure 2:Transition-state (TS) search under different solvents.(A)Comparison of implicit and explicit solvation: implicit solvation represents averaged solvent effects as a field, whereas explicit solvation treats individual solvent molecules.(B)Claisen rearrangement of allyl-p-tolyl ether (APTE).(C)Schematic and 3D structure of the transition state, with the forming and breaking bond lengths (rCOr_{\text{CO}}andrCCr_{\text{CC}}) labeled.(D–F)TS search results for APTE in different solvents. “Ref.” denotes reference calculations at theω\omegaB97M-V/def2-TZVPD/SMD level of theory. Reaction barriers (TS) and reaction energies (P) from each method are labeled below the corresponding states. Ref. and OrbitAll results are obtained using NEB, whereas UMA results are obtained using TS-informed umbrella sampling. Overlaid TS structures are shown with Ref. in gray and OrbitAll in blue;rCOr_{\text{CO}},rCCr_{\text{CC}}, and RMSD relative to Ref. are reported in Å in the inset tables.(G)Wall-time comparison across methods. UMA wall time is measured for TS-informed umbrella sampling as an explicit solvation baseline, and Ref. with GPU4PySCF and OrbitAll wall times are measured for NEB with implicit solvation.We use the Claisen rearrangement of allyl-p-tolyl ether as an example to predict the MEP with NEB using OrbitAll and the reference method (ω\omegaB97M-V/def2-TZVPD/SMD). The reaction schematic and a representative transition-state (TS) structure are shown in Figures2(B) and2(C), respectively. For both OrbitAll and the reference, we perform NEB simulations under different implicit solvents; experimental details are provided in Section4.3. We also compare to UMA (UMA-s-1p2)[wood2025uma]to provide an explicit solvation baseline with umbrella sampling. Details of the experiments are provided in Section4.3.

The fine-tuned OrbitAll shows strong performance in predicting the TS energetics and structures of the Claisen rearrangement across all solvent conditions (Figures2(D–F)). The predicted reaction barriers (Δ​E‡\Delta E^{\ddagger}) and reaction energies (Δr​E\Delta_{r}E) closely match the reference values in all solvents. For the TS structures, OrbitAll accurately predicts the bond-breaking/forming distances,rCOr_{\text{CO}}andrCCr_{\text{CC}}, with overall RMSDs of only 0.02–0.03 Å. Notably, OrbitAll captures subtle solvent-dependent structural trends, including the longer TS bond lengths in water and methanol than in toluene. It also reproduces the corresponding energetic trend, predicting a higher reaction barrier in toluene than in water and methanol. Compared with the explicit-solvent umbrella-sampling baseline using UMA, OrbitAll achieves reference-level accuracy at approximately two orders of magnitude lower computational cost, providing a scalable approach for TS search under implicit solvation.

Additionally, OrbitAll accurately predicts the Diels-Alder reaction between cyclopentadiene and methyl vinyl ketone (FigureS11). Remarkably, it captures the higher reaction barrier in toluene than in vacuum, as well as the lower barriers in water and methanol than in vacuum. Furthermore, the generalization experiment shows that OrbitAll can generalize to solvents similar to those in the training dataset (FigureS10), highlighting the need for datasets with broader solvent coverage.

2.3Molecular Systems with Different Spins and Charges

To accurately predict properties of charged and/or open-shell species, spin and charge information must be incorporated, quantities often missing from the input data structure of geometry-based GNN models[batzner2022nequip,schutt2018schnet,schutt2021painn,batatia2022mace,liao2023equiformerv2,gasteiger2022dimenetpp]. These quantities are inherently included in the orbital features of OrbitAll, providing physical information such as electronic distribution and intensity of spin polarization.

We compare OrbitAll’s performance to seven geometric GNNs: SpookyNet[unke_spookynet_2021], SchNet[schutt2018schnet], PaiNN[schutt2021painn], DimeNet++[gasteiger2022dimenetpp], SphereNet[liu2022spherenet], NequIP[batzner2022nequip], and EquiformerV2[liao2023equiformerv2]. To incorporate charge and spin into different geometric GNN models, we augment the node representations of these models with additional channels for spin and charge embeddings, as introduced in[unke_spookynet_2021]. Modified models are labeled “-SC” to indicate inclusion of spin and charge. Models trained on delta-labels are marked “(Δ\Delta)”.

Evaluation on Varying Spins and Charges

We use QM9star[tang2024qm9star], a QM9-based dataset[ramakri_qm9_2014]with neutral open-shell doublet molecules (radicals), closed-shell charged molecules (anions and cations), and neutral closed-shell molecules (neutrals). Each molecular property is calculated at the B3LYP-D3(BJ)/6-311+G(d,p) level of theory. We evaluate OrbitAll’s performance on two tasks: (1) the total energy (EtotE_{\text{tot}}) of molecules with varying spins and charges; and (2) the HOMO and LUMO energy levels of radical species for both spins,α\alphaandβ\beta(EHOMOαE^{\alpha}_{\text{HOMO}},ELUMOαE^{\alpha}_{\text{LUMO}},EHOMOβE^{\beta}_{\text{HOMO}},ELUMOβE^{\beta}_{\text{LUMO}}).

Figure 3:The QM9star (B3LYP-D3(BJ)/6-311+G(d,p)) benchmarks with delta-learning.The numbers of learnable parameters are indicated inside brackets next to model names in the legend.(A)Learning curves of different models for all species and for each category (neutral, radical, cation, anion). Each category’s learning curve of different models for the total energyEtotE_{\text{tot}}is plotted. For training sizeNN, the training set containsN/4N/4neutrals, radicals, cations, and anions each. The dotted black lines indicate chemical accuracy (1 kcal/mol≈\approx43.4 meV).(B)HOMO and LUMO energy levels prediction benchmark of radical species for up-spin (α\alpha) and down-spin (β\beta). The test mean absolute errors (MAEs) of different models trained using 100K radicals are presented for each property (EHOMOαE^{\alpha}_{\text{HOMO}},ELUMOαE^{\alpha}_{\text{LUMO}},EHOMOβE^{\beta}_{\text{HOMO}},ELUMOβE^{\beta}_{\text{LUMO}}).For the total energyEtotE_{\text{tot}}prediction, we train all models with the delta-learning strategy. This enables assessing the advantage of orbital learning separately, by comparing against purely geometry-based models. We evaluate each model using training data of varying sizes. All training and validation datasets are composed of equally distributed species (neutral, radical, cation, anion). The largest training set consists of 400K data (100K for each species), 20K validation data (5K for each species), and∼\sim75K test data (∼\sim15K neutrals, 20K for other species). We keep the species ratio the same in all smaller training and validation sets.

As shown in Figure3(A), OrbitAll outperforms the other models consistently on the combined test set for all training sets of different sizes. Additionally, OrbitAll outperforms all other models for each separate species across all training sets of different sizes except for the neutral subset in training sizes above 200K. OrbitAll is especially valuable in predicting non-neutral species (radical, anion, cation), with smaller differences between theEtotE_{\text{tot}}mean absolute errors (MAEs) of neutral and non-neutral species (TableS4).

Furthermore, delta-learning enhances generalization as well. For QM9star, we observe that training with the direct-learning strategy leads to consistently higher errors for cations than anions across all models (TableS4)–a result that is not observed when training with the delta-learning strategy. One explanation for this is the distribution shifts of the labels, where the delta-labels are easier to predict, especially for charged species (FigureS18).

We also use OrbitAll to predict frontier molecular orbital (FMO) energy levels for the QM9star dataset. For radical species, four FMO energy levels exist:EHOMOαE^{\alpha}_{\text{HOMO}},ELUMOαE^{\alpha}_{\text{LUMO}},EHOMOβE^{\beta}_{\text{HOMO}}, andELUMOβE^{\beta}_{\text{LUMO}}. We evaluate and compare the performance of OrbitAll in predicting all these levels with delta-learning. Since unpaired electrons cause spin polarization and split spatial orbital energy levels by spin, predicting FMO levels for radicals is expected to be more challenging than for closed-shell species. As shown in3(B), similar toEtotE_{\text{tot}}, OrbitAll outperforms other competing models across the different FMO levels, demonstrating the effectiveness of spin-polarized orbital features for predicting such quantities.

We present additional benchmarks and comparisons for predicting the singlet and triplet energies and their gaps of carbene molecules (the QMSpin dataset) in AppendixS4.2. OrbitAll achieves the lowest MAEs in nearly all categories, with singlet/triplet errors and vertical/adiabatic spin-gap MAEs below chemical accuracy. OrbitAll uses a unified joint model, which likely improves robustness across spin multiplicities.

Cost-Accuracy Analysis

To employ delta-learning, QC calculations are needed, which are involved for OrbitAll during feature generation. To develop a description of the cost-accuracy tradeoff, we compare OrbitAll with direct-learned geometric GNNs.

Figure 4:Cost, accuracy, and data-efficiency comparisons.“(D)” indicates direct-learning and “(Δ\Delta)” indicates delta-learning[rama_delta_learning_2015]. Black dotted lines represent chemical accuracy (1 kcal/mol≈\approx43.4 meV).(A)Data-efficiency comparison using learning curves of different models. OrbitAll offers greater data-efficiency, requiring∼\sim10 times less data than the next best model (DimeNet++-SC(D)) to achieve chemical accuracy on the QM9star dataset.(B)Cost-accuracy comparisons between different methods on QM9star. OrbitAll achieves∼\sim1,000 times speedup compared to the QM, B3LYP-D3(BJ)/6-311+G(d,p). All models are trained using the 400K training set of QM9star.Delta-learning typically leads to a reduction of errors with a certain offset to the learning curves. This brings a significant advantage in data-efficiency compared to directly learning the targets, as shown in Figure4(A). By interpolation, we find that OrbitAll requires∼\sim7K training data to achieve chemical accuracy, whereas the next-best model (DimeNet++-SC(D)) requires∼\sim70K training data to achieve chemical accuracy. Furthermore, shown in Figure4(B), OrbitAll achieves about 1,000-fold speedup compared to the original DFT method while still providing a near-chemical accuracy MAE in the polypeptide dataset. This means that OrbitAll can accelerate simulations by several orders of magnitude while retaining accuracy. Although orbital learning introduces additional computational costs, OrbitAll achieves comparable acceleration while retaining its key advantages, such as data-efficiency.

Transferability to Much Larger Molecules than in the Training Set

A key advantage of predictive systems is their applicability to large molecular systems where traditional quantum mechanical simulations become impractical. However, most AI models exhibit degraded performance when applied to molecules outside their training set. To assess the ability of models to extrapolate across molecular systems of such sizes, we created thepolypeptide dataset, consisting of 151 data points covering neutral, radical, cation, and anion species. The models trained on the 400K QM9star training set (Section2.3) are used to predict total energies (EtotE_{\text{tot}}) of the polypeptides, computed at the B3LYP-D3(BJ)/6-311+G(d,p) level of theory (the same level of theory as QM9star).

Figure 5:Size-wise extrapolation performance of different models evaluated using the polypeptide dataset for the total energyEtotE_{\text{tot}}prediction on 151 polypeptides.All models are trained using the 400K training set of QM9star. “(D)” indicates direct-learning and “(Δ\Delta)” indicates delta-learning.(A)Size-wise extrapolation capabilities. The MAEs normalized by the number of heavy atoms (C, N, O) is shown as a function of heavy atom count. The five points on the left side are the normalized MAEs for the QM9star test set. The remaining points on the right side are the normalized MAEs for the polypeptide dataset inference, which are grouped into bins of width 3 (e.g., 16–18, 19–21, …) based on the number of heavy atoms. The MAE and standard error of the mean (SEM) are reported across systems within the bin/point.(B)Proportion of predictions ofEtotE_{\text{tot}}surpassed chemical accuracy (1 kcal/mol≈\approx43.4 meV). OrbitAll records∼\sim1.5 times more predictions within chemical accuracy than the next best model.(C)Cost-accuracy comparisons between different models on the polypeptide dataset. The black dotted line indicates chemical accuracy. OrbitAll achieves near-chemical accuracy MAE (48.6 meV), with∼\sim10,000 times speedup compared to the original method, B3LYP-D3(BJ)/6-311+G(d,p).Figure5(A) shows MAEs normalized by the number of heavy atoms (NheavyN_{\text{heavy}}) acrossNheavyN_{\text{heavy}}ranges for both QM9star and the polypeptide dataset. OrbitAll shows the smallest increase in normalized error, maintaining consistently low MAEs and demonstrating superior extrapolation to larger molecules. As a result, OrbitAll yields the highest fraction ofEtotE_{\text{tot}}predictions within chemical accuracy (Figure5(B)). Figure5(C) compares the cost and accuracy of different models on the polypeptides; as in Figure4(B), only models trained with the direct-learning strategy (excluding OrbitAll) are included. Owing to the larger size and electron count of polypeptides compared to QM9star molecules, the speedup is even greater than that in Figure4(B). Solely recording MAE with near-chemical accuracy (48.6 meV) on average, OrbitAll serves as a method with balanced speed and accuracy without sacrificing either.

2.4Molecular Systems Under Various Environmental Effects

Predicting the response of molecular systems to external perturbations from the environment is necessary for calculating important molecular properties. Without explicit embedding layers (e.g., Ref.[gastegger2021fieldschnet,zhang2026consolv]), OrbitAll naturally incorporates environmental effects into a common representation space, since the orbital features obtained by the SCF procedure are directly perturbed by the environment. Environmental effects, such as external electric fields and solvation via implicit treatment, can be applied within the semi-empirical QM simulation (e.g., GFNnn-xTB, g-xTB).

Molecular Systems in Solvents

We trained a single OrbitAll model to predict the total energiesEtotE_{\text{tot}}of molecules in the four solvent environments: vacuum, water, tetrahydrofuran (THF), and toluene. Orbital features were generated using the ALPB method[sigalov2006alpb,ehlert2021alpb2], chosen after ablation studies with other implicit solvation methods (TableS7)

Figure 6:Evaluations of OrbitAll on Hessian QM9 (ω\omegaB97x/6-31G*).OrbitAll is trained using 128K training data (32K of each solvent).(A)Violin plots of error distributions for different solvents. The dotted black lines indicate chemical accuracy (±\pm1 kcal/mol≈\approx43.4 meV).(B)Mean absolute errors (MAEs) for molecules in different solvents and the average (“All”). The error bars represent the standard errors of mean (SEMs).As shown in Figure6(A), the error distributions for different solvents are almost identical, demonstrating the robust performance of OrbitAll for different solvents. Compared to baselines provided in[williams2025hessianqm9], OrbitAll records about 3-4 times smaller MAE across different solvents. Additionally, the MAEs for different solvents are very close while also being significantly smaller than chemical accuracy (Figure6(B)).

An example of this perturbation is illustrated in FigureS8. In addition, we conducted another set of experiments using uniform external electric fields as an environmental effect. The results are presented in AppendixS4.5.

3Discussion

We present OrbitAll, a physics-informed end-to-end SE(3)-equivariant deep learning framework designed to process all molecular systems by simultaneously accounting for varying charges, spin states, and environmental effects. We demonstrate the accuracy and robustness of OrbitAll by evaluating its performance using different quantum chemistry datasets of all molecular systems. OrbitAll consistently outperforms other machine learning approaches, mainly due to its ability to integrate physically-grounded information across these different conditions–conditions that other models often treat separately, if at all. OrbitAll incorporates physical effects through SCF-derived orbital features, providing a practical route toward more unified quantum-chemical learning frameworks.

Incorporating environmental effects remains particularly challenging, as their diversity and complexity increase input dimensionality and often lead to optimization difficulties. In contrast, OrbitAll naturally captures environmental effects by embedding physics-informed information into molecular representations, as shown by the accurate predictions of reaction energetics and TS structures for reactions under implicit solvents. These capabilities could accelerate realistic-condition reaction prediction without requiring expensive DFT calculations or MLIP models in which solvent effects are added post hoc.

OrbitAll enables property prediction at DFT-level accuracy while achieving a10310^{3}–10410^{4}-fold speedup over DFT. The accuracy on par with high-level QM simulations is achieved with strong data efficiency, which is essential when working with costly quantum chemical data. This data efficiency suggests that orbital-feature-based models such as OrbitAll may be useful for future data-scarce applications, such as experimentally labeled molecular properties.

One limitation of OrbitAll is its reliance on the convergence of the underlying SCF method in the low-level semi-empirical QM calculations. Since the neural network depends on QMMs generated from these simulations, failure of the SCF procedure prevents feature generation and, consequently, prediction. A possible future study to mitigate this limitation is to use more robust alternatives, such as non-self-consistent field methods like GFN0-xTB[pracht2019gfn0xtb].

4Methods

4.1The OrbitAll Framework

Unrestricted Open-shell Features

In QM calculations, electron spin can be treated using either restricted or unrestricted methods. Restricted calculations (e.g., restricted Hartree-Fock, RHF) constrain the spin orbitals within a spatial orbital to have the same energy level whereas unrestricted methods (e.g., unrestricted Hartree-Fock, UHF) account for spin polarization arising from unequal numbers of spin-up and spin-down electrons. Depending on the system, one may choose between restricted and unrestricted open-shell treatments. In this work, we adopt an unrestricted formulation to generate orbital features, enabling the model to represent a wider range of molecular systems, including open-shell and spin-polarized species[szabo_quantum_1989].

For unrestricted open-shell systems, the Pople-Nesbet equations, equivalent to the Roothaan-Hall equation for restricted Hartree-Fock, can be expressed as:

Fα​Cα=SCα​εα,Fβ​Cβ=SCβ​εβ,\textbf{F}^{\alpha}\textbf{C}^{\alpha}=\textbf{S}\textbf{C}^{\alpha}\mathbf{\varepsilon}^{\alpha},\qquad\textbf{F}^{\beta}\textbf{C}^{\beta}=\textbf{S}\textbf{C}^{\beta}\mathbf{\varepsilon}^{\beta},(4) whereα\alphadenotes up-spin,β\betadenotes down-spin,Fα\textbf{F}^{\alpha}andFβ\textbf{F}^{\beta}are the Fock matrices,Cα\textbf{C}^{\alpha}andCβ\textbf{C}^{\beta}are orbital coefficients matrices,Sis the overlap matrix, andεα\mathbf{\varepsilon}^{\alpha}andεβ\mathbf{\varepsilon}^{\beta}are diagonal orbital energy matrices. Also, the density matrices for different spins are given by:

Pμ​να=∑aNαCμ​aα​(Cν​aα)∗,Pμ​νβ=∑aNβCμ​aβ​(Cν​aβ)∗,P^{\alpha}_{\mu\nu}=\sum^{N_{\alpha}}_{a}C^{\alpha}_{\mu a}(C^{\alpha}_{\nu a})^{*},\qquad P^{\beta}_{\mu\nu}=\sum^{N_{\beta}}_{a}C^{\beta}_{\mu a}(C^{\beta}_{\nu a})^{*},(5) wherePμ​ναP^{\alpha}_{\mu\nu}andPμ​νβP^{\beta}_{\mu\nu}are elements of the density matricesPα\textbf{P}^{\alpha}andPβ\textbf{P}^{\beta}, respectively, for orbital basesμ\muandν\nu.NαN_{\alpha}andNβN_{\beta}are the numbers of up-spin electrons and down-spin electrons, respectively.Cμ​ναC^{\alpha}_{\mu\nu}andCμ​νβC^{\beta}_{\mu\nu}are elements of the orbital coefficient matricesCα\textbf{C}^{\alpha}andCβ\textbf{C}^{\beta}.

The unrestricted Hartree-Fock (UHF) formulation generalizes the standard Hartree-Fock approach to handle open-shell systems. Specifically, for closed-shell molecules, the resulting UHF matrices reduce to those of the restricted Hartree-Fock (RHF) method in a way that:

Fα=Fβ=F,\textbf{F}^{\alpha}=\textbf{F}^{\beta}=\textbf{F},(6) Pα=Pβ=12​P,\textbf{P}^{\alpha}=\textbf{P}^{\beta}=\frac{1}{2}\textbf{P},(7) wherePis the density matrix for restricted closed-shell Hartree-Fock[szabo_quantum_1989]. The density matrixPfor restricted closed-shell Hartree-Fock is then equivalent to the total density matrixPT=Pα+Pβ\textbf{P}^{T}=\textbf{P}^{\alpha}+\textbf{P}^{\beta}for unrestricted Hartree-Fock. This equivalence enables a unified input representation that can accommodate both open- and closed-shell systems. It is important to note that this assumes the molecule is a closed-shell singlet, not an open-shell singlet.

The QMMs are SE(3)-equivariant as shown in (2). The set of QMMs for OrbitAll,T=(Fα,Fβ,Pα,Pβ,S,Hcore)\textbf{T}=(\textbf{F}^{\alpha},\textbf{F}^{\beta},\textbf{P}^{\alpha},\textbf{P}^{\beta},\textbf{S},\textbf{H}_{\text{core}}), is used as inputs to the equivariant graph neural network backbone, UNiTE[qiao_orbnet_equi_2022]. These features vary with molecular conditions such as spin, charge, and environmental effects. For example, the electron densities for up-spin and down-spin at positionris given by:

ρα​(r)=∑μ∑ν(Pα)μ,ν​Φμ​(r)​(Φν​(r))∗,\rho^{\alpha}(\textbf{r})=\sum_{\mu}\sum_{\nu}(\textbf{P}^{\alpha})^{\mu,\nu}\Phi^{\mu}(\textbf{r})(\Phi^{\nu}(\textbf{r}))^{*},(8) ρβ​(r)=∑μ∑ν(Pβ)μ,ν​Φμ​(r)​(Φν​(r))∗,\rho^{\beta}(\textbf{r})=\sum_{\mu}\sum_{\nu}(\textbf{P}^{\beta})^{\mu,\nu}\Phi^{\mu}(\textbf{r})(\Phi^{\nu}(\textbf{r}))^{*},(9) whereΦμ​(r)\Phi^{\mu}(\textbf{r})andΦν​(r)\Phi^{\nu}(\textbf{r})are the atomic orbital bases. By definition, they satisfy

∫𝑑r​(ρα​(r)+ρβ​(r))=Nelec,\int d\textbf{r}\,(\rho^{\alpha}(\textbf{r})+\rho^{\beta}(\textbf{r}))=N_{\text{elec}},(10) ∫𝑑r​(ρα​(r)−ρβ​(r))=2​S,\int d\textbf{r}\,(\rho^{\alpha}(\textbf{r})-\rho^{\beta}(\textbf{r}))=2S,(11) whereNelecN_{\text{elec}}denotes the number of electrons, which depends on the system’s charge. As a result, the density matrix inherently encodes both spin and charge information. Environmental effects influence the density matrix more subtly by altering the converged mean-field compared to that of an isolated system. An example of such a perturbation is shown in FigureS8.

Any SCF method solving the Roothaan-Hall or Pople-Nesbet equations can be used in this framework. We use spGFN1-xTB[neugebauer_spgfnxtb_2023]for the experiments in Sections2.3and2.4, and g-xTB for the experiments in Section2.2[froitzheim2025gxtb], to efficiently generate orbital features with unrestricted Hartree-Fock.

Equivariant Graph Neural Network

To build a data-efficient model that robustly predicts SE(3) (rotations and translations) equivariant and/or invariant molecular properties such as forces and dipole moments, we construct a GNN, based on the E(3) equivariant UNiTE framework[qiao_orbnet_equi_2022]. The inherent SE(3) equivariance of the QMMs is maintained through our GNN, which results in an end-to-end SE(3) equivariant framework. The GNN satisfies the necessary rotational and translational symmetries:

ℛ⋅ℱ​(T)=ℱ​(ℛ⋅T),\mathcal{R}\cdot\mathcal{F}(\textbf{T})=\mathcal{F}(\mathcal{R}\cdot\textbf{T}),(12) whereℱ\mathcal{F}characterizes the parameters of the GNN,ℛ\mathcal{R}is an arbitrary roto-translational operation, andℛ⋅\mathcal{R}\,\cdotrepresents applying the roto-translational transformationℛ\mathcal{R}. The geometric GNN learns to map a set of QMMs of any molecular system,T, to an atomic or molecular property, by minimizing the objective:

minℱ⁡ℒ​(y,y^),\min_{\mathcal{F}}\mathcal{L}(y,\hat{y}),(13) whereℒ\mathcal{L}is the loss function specific to the predicted property,yyis the target property, either generated by simulation or estimated by experiments, andy^=ℱ​(T)\hat{y}=\mathcal{F}(\textbf{T}).

Diagonal Reduction and Embedding

All QMMs have the same dimensions of(NAO,NAO)(N_{\text{AO}},N_{\text{AO}}), whereNAON_{\text{AO}}is the number of atomic orbitals in the system. Each row and column corresponds to an atomic orbital. The firstNAOAN_{\text{AO}}^{A}entries of rows and columns correspond to the atomic orbitals of atomAA, the nextNAOBN_{\text{AO}}^{B}to those of atomBB, and so on. As a result, the block-diagonal regions of the QMMs capture intra-atomic interactions, while the off-diagonal blocks encode inter-atomic interactions. These inter-atomic blocks are later used to construct messages in the model, as illustrated in Figure1(C).

From the atomic orbital basis QMMs, we construct an atom-based representation via the diagonal reduction module. This module embeds the block-diagonals of QMMOto the reduced embeddingshAO\textbf{h}^{O}_{A}, as follows:

hA,n​l​p​mO={∑μ,ν(O)A​Aμ,ν​(Q~)A,n​l​mμ,ν,p=+1,0,p=−1,\textbf{h}^{O}_{A,nlpm}=\begin{cases}\sum_{\mu,\nu}(\textbf{O})^{\mu,\nu}_{AA}(\tilde{\textbf{Q}})^{\mu,\nu}_{A,nlm},&p=+1,\\ 0,&p=-1,\end{cases}(14) whereppis the parity andQ~\tilde{\textbf{Q}}is an on-site three-index overlap integrals, defined by:

(Q~)A,n​l​mμ,ν=∫ℝ3𝑑r​(ΦAμ​(r))∗​ΦAν​(r)​Φ~An,l,m​(r),(\tilde{\textbf{Q}})^{\mu,\nu}_{A,nlm}=\int_{\mathbb{R}^{3}}d\textbf{r}\,(\Phi^{\mu}_{A}(\textbf{r}))^{*}\Phi^{\nu}_{A}(\textbf{r})\tilde{\Phi}^{n,l,m}_{A}(\textbf{r}),(15) whereΦ~An,l,m​(r)\tilde{\Phi}^{n,l,m}_{A}(\textbf{r})is an auxiliary Gaussian-type basis as defined in[qiao_orbnet_equi_2022]. Note thatQ~\tilde{\textbf{Q}}is proportional to the Clebsch-Gordan coefficients, ensuring the overall process is equivariant to SO(3).

Since g-xTB and spGFN1-xTB employ different basis sets, we compute and tabulate the three-index overlap integrals separately for each method. In both approaches, atomic orbitals are constructed from primitive Gaussian-type orbitals (GTOs). However, while spGFN1-xTB uses fixed contraction coefficients, g-xTB employs environment-dependent (i.e., charge-dependent) contraction coefficients[froitzheim2025gxtb]. Specifically, for an atomAiA_{i}of element typeAA, an AO basis functionΦAiμ\Phi_{A_{i}}^{\mu}is expressed as a contraction over primitive GTOs:

ΦAiμ​(𝐫)=∑λcAiμ​λ​(qAieff)​χAλ​(𝐫;ζAλ),\Phi^{\mu}_{A_{i}}(\mathbf{r})=\sum_{\lambda}c^{\mu\lambda}_{A_{i}}\left(q^{\mathrm{eff}}_{A_{i}}\right)\,\chi^{\lambda}_{A}(\mathbf{r};\zeta^{\lambda}_{A}),(16) whereλ\lambdaindexes the primitive GTOs associated with elementAA,χAλ​(𝐫;ζAλ)\chi^{\lambda}_{A}(\mathbf{r};\zeta^{\lambda}_{A})denotes a primitive GTO with exponentζAλ\zeta^{\lambda}_{A}, andcAiμ​λ​(qAieff)c^{\mu\lambda}_{A_{i}}(q^{\mathrm{eff}}_{A_{i}})are contraction coefficients that depend explicitly on the effective atomic chargeqAieffq^{\mathrm{eff}}_{A_{i}}.

Because the contraction coefficients vary with the chemical environment throughqAieffq^{\mathrm{eff}}_{A_{i}}, the on-site three-index overlap integrals cannot be tabulated directly in the contracted AO basis. Instead, we tabulate the corresponding quantities in the primitive GTO basis,

(𝐐~)A,n​l​mλ,σ=∫ℝ3𝑑𝐫​(χAλ​(𝐫))∗​χAσ​(𝐫)​Φ~An​l​m​(𝐫),(\tilde{\mathbf{Q}})^{\lambda,\sigma}_{A,nlm}=\int_{\mathbb{R}^{3}}d\mathbf{r}\;(\chi^{\lambda}_{A}(\mathbf{r}))^{*}\,\chi^{\sigma}_{A}(\mathbf{r})\,\tilde{\Phi}^{nlm}_{A}(\mathbf{r}),(17) whereλ,σ\lambda,\sigmalabel primitive GTOs. These primitive-basis integrals depend only on the elementAA(and the auxiliary indices) and can therefore be precomputed.

Given the contraction coefficients for a specific atomAiA_{i}, the corresponding contracted on-site three-index overlap integral is obtained by a straightforward expansion,

(𝐐~)Ai,n​l​mμ,ν=∫ℝ3𝑑𝐫​ΦAiμ​(𝐫)∗​ΦAiν​(𝐫)​Φ~An​l​m​(𝐫)=∑λ,σ(cAiλ)∗​cAiσ​∫ℝ3𝑑𝐫​(χAλ​(𝐫))∗​χAσ​(𝐫)​Φ~An​l​m​(𝐫)=∑λ,σ(cAiλ)∗​cAiσ​(𝐐~)A,n​l​mλ,σ.\begin{split}(\tilde{\mathbf{Q}})^{\mu,\nu}_{A_{i},nlm}&=\int_{\mathbb{R}^{3}}d\mathbf{r}\;\Phi^{\mu}_{A_{i}}(\mathbf{r})^{*}\,\Phi^{\nu}_{A_{i}}(\mathbf{r})\,\tilde{\Phi}^{nlm}_{A}(\mathbf{r})\\ &=\sum_{\lambda,\sigma}\left(c^{\lambda}_{A_{i}}\right)^{*}c^{\sigma}_{A_{i}}\int_{\mathbb{R}^{3}}d\mathbf{r}\;(\chi^{\lambda}_{A}(\mathbf{r}))^{*}\,\chi^{\sigma}_{A}(\mathbf{r})\,\tilde{\Phi}^{nlm}_{A}(\mathbf{r})\\ &=\sum_{\lambda,\sigma}\left(c^{\lambda}_{A_{i}}\right)^{*}c^{\sigma}_{A_{i}}\,(\tilde{\mathbf{Q}})^{\lambda,\sigma}_{A,nlm}.\end{split}(18) Thus, while the contracted on-site three-index overlap integrals(𝐐~)Ai,n​l​mμ,ν(\tilde{\mathbf{Q}})^{\mu,\nu}_{A_{i},nlm}become atom-specific through the environment-dependent contraction coefficients, they can be swiftly evaluated from the element-tabulated primitive integrals(𝐐~)A,n​l​mλ,σ(\tilde{\mathbf{Q}})^{\lambda,\sigma}_{A,nlm}once the contraction coefficients are known (e.g., from the basis specification provided by the electronic-structure driver). The contraction coefficients are readjusted with respect to their Frobenius norm for consistent scales of features.

Reduced embeddings from each QMM are passed through linear layers with learnable weights, producing the initial hidden featureshAt=0\textbf{h}^{t=0}_{A}. The number of channels is listed in TableS2.

Message Passing

At layertt, the equivariant message from atomBBto atomAA,mB​At\textbf{m}^{t}_{BA}, is computed via block convolutions over the off-diagonal QMM blocks, following[qiao_orbnet_equi_2022]. These blocks are projected onto convolution channels, with their number specified in TableS1.

Messages from neighboring atoms are aggregated using multi-head attention. The resulting message for atomAAis:

m~At=∑B⨁i,jmB​At,i⋅αA​Bt,j,\tilde{\textbf{m}}^{t}_{A}=\sum_{B}\bigoplus_{i,j}\textbf{m}^{t,i}_{BA}\cdot\alpha^{t,j}_{AB},(19) whereiiis the convolution channel index, andαA​Bt,j\alpha^{t,j}_{AB}is the invariant attention of the atomAAto the atomBBfor thejj-th attention head. The aggregated messagem~At\tilde{\textbf{m}}^{t}_{A}is then coupled with the node representation of the atomAAat layertt,hAt\textbf{h}^{t}_{A}with the point-wise interaction module defined in[qiao_orbnet_equi_2022], which updates the node representation tohAt+1\textbf{h}^{t+1}_{A}. Further details on the message constructions and the point-wise interaction can be found in AppendixS1.

Atom-wise Decoding and Pooling

The atom-wise decoding layer updates each node using its own features. A point-wise interaction is applied tohAt\textbf{h}^{t}_{A}, producing the more abstract updated representationhAt+1\textbf{h}^{t+1}_{A}.

After passing thetmt_{m}message passing layers andtdt_{d}atom-wise decoding layers, the final representationhAtm+td\textbf{h}^{t_{m}+t_{d}}_{A}is used for the task-specific pooling operation. In this work, two distinct pooling operations are employed as in[qiao_orbnet_equi_2022]: one for predicting the total energy and another for the FMO property.

The predicted total energy,E^tot\hat{E}_{\text{tot}}, is given by

E^tot=∑AWo⋅∥hAtm+td∥+bZA,\hat{E}_{\text{tot}}=\sum_{A}\textbf{W}_{o}\cdot\lVert\textbf{h}^{t_{m}+t_{d}}_{A}\rVert+b_{Z_{A}},(20) whereAAis the atom index,Wo\textbf{W}_{o}is a learnable matrix,bZAb_{Z_{A}}is the element-wise energy bias, i.e., the element-wise shift, for atomAA’s atomic numberZAZ_{A}. The element-wise shifts are initialized from a linear regression of the total energy to atomic numbers.

For predicting FMO energies, pooling with global attention is applied as follows:

aA=Wa⋅∥hAtm+td∥∑BWa⋅∥hBtm+td∥,a_{A}=\frac{\textbf{W}_{a}\cdot\lVert\textbf{h}^{t_{m}+t_{d}}_{A}\rVert}{\sum_{B}\textbf{W}_{a}\cdot\lVert\textbf{h}^{t_{m}+t_{d}}_{B}\rVert},(21) E^FMO=∑AaA​(Wo⋅∥hAtm+td∥+bZA),\hat{E}_{\text{FMO}}=\sum_{A}a_{A}(\textbf{W}_{o}\cdot\lVert\textbf{h}^{t_{m}+t_{d}}_{A}\rVert+b_{Z_{A}}),(22) whereaAa_{A}is a scalar attention to the atomAA, andWa\textbf{W}_{a}is a learnable matrix.

4.2Datasets

The OMol25 Dataset

The OMol25 dataset is a large-scale collection comprising 140 million single-point calculations performed at theω\omegaB97M-V/def2-TZVPD level of theory[levine2025openmolecules2025omol25]. It includes molecular systems spanning a wide range of charge states (0 to±\pm10) and spin configurations (spin quantum numbers from 0 to 5, corresponding to up to 10 unpaired electrons and a maximum spin multiplicity of 11). Covering 83 elements and systems containing up to 350 atoms, the dataset captures substantial chemical diversity and complexity.

In this work, we use the 4M (4 million) subset of OMol25, which is uniformly sampled from the full 140M dataset. When employing g-xTB as the underlying semi-empirical quantum mechanical method, molecules containing lanthanides are excluded from both training and evaluation. This is due to the f-in-core strategy adopted in g-xTB[froitzheim2025gxtb], wherein f-electrons are treated as core electrons rather than valence electrons, leading to ambiguities in spin-state representation. Additionally, all molecules that failed to converge or produced errors during preprocessing are removed from the dataset.

The official test set of the OMol25 dataset is not publicly available with labels. Consequently, we designate the provided validation split, whose labels are accessible, as our in-house test set for model evaluation. To retain a validation set for hyperparameter tuning, we randomly sample 65,536 molecules from the original training split and use them for validation. After this re-partitioning, the dataset comprises 3.85M molecules for training, 65.5K for validation, and 2.70M for testing. The model was trained with the delta-learning strategy, without the element-wise or charge biases.

The T1x-Solv Dataset

Transition1x is a dataset containing 10,073 reactions and their reaction pathways, including reactants, products, and transition states. Each reaction includes both converged and unconverged pathways from NEB calculations; therefore, the number of data points varies across reactions. In total, the dataset comprises 9.6M single-point calculations at theω\omegaB97x/6-31G(d) level of theory.

To build a reactive potential under different solvent conditions, we carefully selected geometries from the Transition1x dataset and used them for single-point calculations at theω\omegaB97M-V/def2-TZVPD level of theory. Specifically, we extracted eight geometries per reaction: (1) the relaxed reactant, (2) the relaxed product, (3) the converged transition state and its two neighboring images, and (4) three randomly selected images. This procedure yielded 80,584 geometries, which were then used for single-point calculations under four solvent conditions: vacuum, water, methanol, and toluene. Consequently, the dataset contains 322,336 data points with DFT-evaluated energies and forces. The single-point calculations were performed using ORCA 6.1.0[neese2025orca], with solvation modeled using SMD[marenich2009smd]. The ‘RIJ-COSX’ and ‘TIGHTSCF’ flags were enabled.

Given the current implementation of g-xTB 2.0.0[froitzheim2025gxtb], we used the generalized Born model with finite dielectric constant (GBE) to generate the QMMs. The data points that SCF failed are removed from the data set for training the OrbitAll model. During fine-tuning of the OrbitAll-OMol25-4M model, we additionally applied solvent-specific element-wise shifts.

Evaluation Datasets and Details

QM9star is a quantum chemistry dataset comprised of about 2 million data points in total:∼\sim120K of neutral singlets (neutrals,Q=0Q=0,S=0S=0),∼\sim435K of+1+1charged singlets (cations,Q=+1Q=+1,S=0S=0),∼\sim721K of−1-1charged singlets (anions,Q=−1Q=-1,S=0S=0), and∼\sim731K of neutral doublets (radicals,Q=0Q=0,S=1/2S=1/2)[tang2024qm9star]. Each datapoint is a small drug-like molecule optimized at the B3LYP-D3(BJ)/6-311+G(d,p) level of theory with subsequent calculations of quantum mechanical properties.

In the QM9star dataset, anions, cations, and radicals are generated by removing a hydrogen atom from neutral “parent molecules.” At the site of hydrogen removal, subtracting an electron yields a singlet cation, adding an electron produces a singlet anion, and leaving the molecule with an unpaired electron results in a radical[tang2024qm9star].

To construct training subsets, we ensured that no parent molecule appears in more than one of the training, validation, or test sets. Validation sets are sized at 10% of their corresponding training subsets, except for the full 400K training set, which uses a fixed 20K validation set. Additionally, each training and validation set of sizeNNis sampled to include equal numbers (N/4N/4each) of neutrals, radicals, cations, and anions.

While training all models (including OrbitAll and the baselines) on the QM9star dataset, we applied charge shifts to account for the overall effect of molecular charge. The charge shift, denoted bybQb_{Q}, is a learnable bias term added during the final pooling stage to represent the average contribution of each charge state.

Hence, the predicted energy of a molecule,E^tot\hat{E}_{\text{tot}}, from the models is given by

E^tot=(∑AE^A+bZA)+bQ,\hat{E}_{\text{tot}}=\left(\sum_{A}\hat{E}_{A}+b_{Z_{A}}\right)+b_{Q},(23) where,AAdenotes the atom (node) index,E^A\hat{E}_{A}is the atomic energy predicted by the GNN,bZAb_{Z_{A}}is the element-wise energy bias (i.e., the shift associated with atomAA’s atomic numberZAZ_{A}), andbQb_{Q}is the charge-dependent energy bias (i.e., thecharge shift). We initialize the charge shifts using the average total energyEtotE_{\text{tot}}, minus the sum of element-wise biases for each charge state in the training set. Notably, applying the charge shift stabilized training and reduced both the frequency and magnitude of outliers for OrbitAll.

QMSpin is a dataset for the closed-shell and open-shell species prediction tasks, which consists of 4.9K singlet-optimized and 7.8K triplet-optimized carbene geometries. The singlet (S=0S=0) and triplet (S=1S=1) energies of each geometry is calculated, resulting in a total of 25.6K energy points. The geometries are optimized via restricted open-shell B3LYP/def2-TZVP, and the energies are obtained at MRCISD+Q-F12/cc-pVDZ-F12[schwilk2020qmspin].

The original geometries of polypeptides were obtained from the PEPCONF dataset[prasad2019pepconf], which consists of neutral (Q=0Q=0,S=0S=0), anion (Q=1Q=1,S=0S=0), and cation (Q=−1Q=-1,S=0S=0) polypeptide molecules. The polypeptide molecules are composed only of elements H, C, N, and O. Since OrbitAll predicts open-shell species properties, we created radical species by either removing an electron from an anion or adding an electron to a cation, to neutralize the charge.

After collecting the ground-state geometry of each molecule from the PEPCONF dataset, which then underwent three consecutive geometry optimizations: first with GFN2-xTB[bannwarth2019gfn2], followed by B3LYP-D3(BJ)/def2-SVP and then B3LYP-D3(BJ)/6-311+G(d,p).

The single point total energiesEtotE_{\text{tot}}of polypeptides are the prediction targets, specifically with 58 neutrals, 46 radicals, 25 cations, and 22 anions. The number of heavy atoms of the QM9star dataset ranges from 1 to 9, whereas it ranges from 16 to 36 in the polypeptide dataset, making it a suitable dataset for evaluating the size extrapolation.

Hessian QM9 is an implicit solvation dataset comprising 41.6K molecules in four solvent environments–vacuum, water, tetrahydrofuran, and toluene–distinguished by their dielectric constants (ϵr\epsilon_{r}). Each datapoint is computed at theω\omegaB97x/6-31G* level of theory using the solvation model based on density (SMD)[marenich2009smd]. For generating QMMs, we use the CPCM method[barone1998cpcm1]implemented in tblite[tblite]with the same dielectric constants used for Hessian QM9.

4.3Experiments

Nudged Elastic Band (NEB)

All NEB simulations are performed using the Atomic Simulation Environment (ASE)[hjorth2017ase]. The FIRE optimizer is used for NEB optimization[bitzek2006fire]. Each simulation uses nine images in total, including the reactant, product, and seven intermediate images. The initial path is first relaxed with NEB until the maximum force (fmaxf_{\text{max}}) is below 0.20 eV/Å, followed by climbing-image NEB (CI-NEB) optimization with anfmaxf_{\text{max}}threshold of 0.05 eV/Å. The reactant and product geometries are optimized primarily using the Broyden–Fletcher–Goldfarb–Shanno (BFGS) method, and structures that do not converge are optimized using the FIRE optimizer. All geometry optimizations are performed with the ASE implementation using a convergence threshold offmax=0.02f_{\text{max}}=0.02eV/Å. The relative energies of the reactant, product, and transition state are then used to calculate the reaction barrier,Δ​E‡=ETS−ER\Delta E^{\ddagger}=E_{\text{TS}}-E_{\text{R}}, and reaction energy,Δr​E=EP−ER\Delta_{r}E=E_{\text{P}}-E_{\text{R}}, whereETSE_{\text{TS}},ERE_{\text{R}}, andEPE_{\text{P}}denote the transition-state, reactant, and product energies, respectively.

Umbrella Sampling

We performed explicit-solvent umbrella sampling with UMA-s-1p2 to establish a strong baseline for explicit-solvation reaction modeling. UMA was selected as a competitive MLIP, providing a favorable reference point for explicit-solvent simulations. Direct TS discovery from unbiased explicit-solvent molecular dynamics (MD) is impractical for the Claisen rearrangement because barrier-crossing events are expected to be rare on accessible MD timescales, considering the high barrier height. We therefore used a TS-informed, or oracle, umbrella-sampling in which transition-state-like windows were initialized from reference structures calculated atω\omegaB79M-V/def2-TZVPD/SMD. This setup gives UMA substantial prior information and should be viewed as a favorable lower-bound estimate of the computational burden required for explicit-solvent reaction profiling.

The initial explicit-solvent configurations for the reactant-, TS-, and product-like umbrella windows are generated using PACKMOL[martinez2009packmol]by embedding the corresponding solute geometries in cubic solvent boxes of size25​Å×25​Å×25​Å25~\text{\AA }\times 25~\text{\AA }\times 25~\text{\AA }, reducing solute–image interactions under periodic boundary conditions. Umbrella sampling is performed with UMA-s-1p2 using the collective variableξ=rCC−rCO\xi=r_{\mathrm{CC}}-r_{\mathrm{CO}}, whererCCr_{\mathrm{CC}}andrCOr_{\mathrm{CO}}denote the forming C–C and breaking C–O bond distances, respectively. We sample 31 umbrella windows with centers ranging fromξ=−1.5​Å\xi=-1.5~\text{\AA }toξ=+1.5​Å\xi=+1.5~\text{\AA }at intervals of0.1​Å0.1~\text{\AA }. A harmonic restraint force constant of10.0​eV/Å210.0~\mathrm{eV}/\text{\AA }^{2}is used for all windows. Each window is simulated in the NVT ensemble at300​K300~\mathrm{K}for15​ps15~\mathrm{ps}, consisting of5​ps5~\mathrm{ps}equilibration and10​ps10~\mathrm{ps}production, with a timestep of0.5​fs0.5~\mathrm{fs}. For windows centered betweenξ=−0.6​Å\xi=-0.6~\text{\AA }andξ=+0.6​Å\xi=+0.6~\text{\AA }, the solvated TS structure is used as the initial configuration, whereas reactant- and product-solvated structures are used for windows withξ>+0.6​Å\xi>+0.6~\text{\AA }andξ<−0.6​Å\xi<-0.6~\text{\AA }, respectively. The biased distributions from the production trajectories are reweighted and combined using the weighted histogram analysis method (WHAM) to reconstruct the one-dimensional potential of mean force alongξ\xi, from which the activation barrier and reaction free energy are estimated.

Wall time measurement

The wall times of the reaction pathway predictions for each method are measured using a single NVIDIA H200 GPU and 32-cores of AMD EPYC 9554 @ 3.1 GHz.

All datasets used in this study will be made available upon publication.

All code used in this study will be made available upon publication.

B.S.K. acknowledges graduate research funding from the California Institute of Technology, support from the Pritzker AI+Science fund and the Eddleman Graduate Fellowship. W.A.G. acknowledges support from NSF(CBET-231117). A.A. acknowledges support from the Bren endowed chair, ONR (MURI grant N00014-23-1-2654), and the Schmidt Sciences AI2050 senior fellow program. This work used the Delta system at the National Center for Supercomputing Applications through allocation DMR160114 from the Advanced Cyberinfrastructure Coordination Ecosystem: Services & Support (ACCESS) program, which is supported by National Science Foundation grants #2138259, #2138286, #2138307, #2137603, and #2138296. B.S.K. acknowledges Robert Kalescky, John Santerre, and Bivin Sadler for their help in organizing the computational resources used in this research through SMU’s O’Donnell Data Science and Research Computing Institute.

References

Appendix S1Architecture

Equivariant Normalization

The equivariant normalization (EvNorm) used in the UNiTE framework is defined by:

(h¯,h^)=EvNorm(h),\displaystyle(\overline{\textbf{h}},\hat{\textbf{h}})=\text{EvNorm({h})},(S24)h¯n​l​p:=∥hn​l​p∥−μn​l​phσn​l​ph,\displaystyle\overline{\textbf{h}}_{nlp}:=\frac{\lVert\textbf{h}_{nlp}\rVert-\mu^{h}_{nlp}}{\sigma^{h}_{nlp}},(S25)∥hn​l​p∥:=∑mhn​l​p​m2+ϵ2−ϵ,\displaystyle\lVert\textbf{h}_{nlp}\rVert:=\sqrt{\sum_{m}h^{2}_{nlpm}+\epsilon^{2}}-\epsilon,(S26)h^n​l​p​m:=hn​l​p​m∥hn​l​p∥+1/βn​l​p+ϵ.\displaystyle\hat{h}_{nlpm}:=\frac{h_{nlpm}}{\lVert\textbf{h}_{nlp}\rVert+1/\beta_{nlp}+\epsilon}.(S27) Here,∥hn​l​p∥\lVert\textbf{h}_{nlp}\rVertis the invariant content ofh,μn​l​ph\mu^{h}_{nlp}andσn​l​ph\sigma^{h}_{nlp}are the mean and variance of∥h∥\lVert\textbf{h}\rVert, respectively,ϵ\epsilonis a stability factor andβn​l​p\beta_{nlp}is a positive learnable scalar that controls the amount of∥hn​l​p∥\lVert\textbf{h}_{nlp}\rVertinformation in the vector part.

Creating Messages by Block-wise Convolutions

The block-wise convolutions characterize the message passing procedure of the UNiTE framework[qiao_orbnet_equi_2022]. Supposing we are using one feature matrixO, a message from the atomAAto the atomBBfor theii-th convolution channel is defined by:

mA​B,νi=∑μ(ρi​(hA))μ​(O)A​Bμ,ν,\textbf{m}_{AB,\nu}^{i}=\sum_{\mu}(\rho_{i}(\textbf{h}_{A}))_{\mu}(\textbf{O})^{\mu,\nu}_{AB},(S28) whereμ=(n1,l1,m1)\mu=(n_{1},l_{1},m_{1})andν=(n2,l2,m2)\nu=(n_{2},l_{2},m_{2})are atomic orbital indices,hA\textbf{h}_{A}is the hidden representation of the atomAA, andρ\rhois a matching layer defined by:

(ρi​(hA))μ=Gather​(Wli⋅(hA)l​(p=+1)​m,n​[μ,ZA]).(\rho_{i}(\textbf{h}_{A}))_{\mu}=\text{Gather}(\textbf{W}_{l}^{i}\cdot(\textbf{h}_{A})_{l(p=+1)m},n[\mu,Z_{A}]).(S29) Here,Wli\textbf{W}^{i}_{l}is a learnable weight matrix for theii-th convolution channel specific to angular momentumll. The “Gather” operation rearranges and maps the input hidden features to the valid orbital indexμ\muusing the principal quantum numbernnof atomic orbitalμ\mufor the atomZAZ_{A}.

The messages are aggregated with multi-head attention as in Equation (17), creating an aggregated messagem~A\tilde{\textbf{m}}_{A}(for the atomAA) that is to be used for updating the hidden representation of the atomAA. The aggregated message is passed to a reverse matching layerρ†\rho^{\dagger}, defined by:

(ρ†​(m~A))l​p​m={Wl†⋅∑μScatter​((m~A)μ,n​[μ,ZA]),p=+1,0,p=−1,(\rho^{\dagger}(\tilde{\textbf{m}}_{A}))_{lpm}=\begin{cases}\textbf{W}^{\dagger}_{l}\cdot\sum_{\mu}\text{Scatter}((\tilde{\textbf{m}}_{A})_{\mu},n[\mu,Z_{A}]),\qquad&p=+1,\\ 0,&p=-1,\end{cases}(S30) whereWl†\textbf{W}^{\dagger}_{l}is a learnable weight matrix specific toll, and the “Scatter” operation flattens the messagem~A\tilde{\textbf{m}}_{A}usingn​[μ,ZA]n[\mu,Z_{A}]. The combined operation maps the message and projects into the same dimension as the atomic hidden representationhA\textbf{h}_{A}.

The Point-wise Interaction

The point-wise interaction module,ϕ\phi, performs normalization and applies non-linearity. It couples two equivariant features,handg, to produce an updated representationh′=ϕ​(h,g)\textbf{h}^{\prime}=\phi(\textbf{h},\textbf{g})through the following operations:

(h¯,h^)=EvNorm(h),\displaystyle(\overline{\textbf{h}},\hat{\textbf{h}})=\text{EvNorm({h})},(S31)fl​p​m=(MLP1​(h¯))l​p⊙(h^l​p​m⋅Wl,pin),\displaystyle\textbf{f}_{lpm}=(\text{MLP}_{1}(\overline{\textbf{h}}))_{lp}\odot(\hat{\textbf{h}}_{lpm}\cdot\textbf{W}^{\text{in}}_{l,p}),(S32)ql​p​m=gl​p​m+∑l1,l2∑m1,m2∑p1,p2fl1​p1​m1⋅gl2​p2​m2⋅Cl1​m1,l2​m2l​m⋅δp1⋅p2⋅p(−1)l1+l2+l,\displaystyle\textbf{q}_{lpm}=\textbf{g}_{lpm}+\sum_{l_{1},l_{2}}\sum_{m_{1},m_{2}}\sum_{p_{1},p_{2}}\textbf{f}_{l_{1}p_{1}m_{1}}\cdot\textbf{g}_{l_{2}p_{2}m_{2}}\cdot C^{lm}_{l_{1}m_{1},l_{2}m_{2}}\cdot\delta^{(-1)^{l_{1}+l_{2}+l}}_{p_{1}\cdot p_{2}\cdot p},(S33)(q¯,q^)=EvNorm(q),\displaystyle(\overline{\textbf{q}},\hat{\textbf{q}})=\text{EvNorm({q})},(S34)hl​p​m′=hl​p​m+(MLP2​(q¯))l​p⊙(q^l​p​m⋅Wl,pout),\displaystyle\textbf{h}^{\prime}_{lpm}=\textbf{h}_{lpm}+(\text{MLP}_{2}(\overline{\textbf{q}}))_{lp}\odot(\hat{\textbf{q}}_{lpm}\cdot\textbf{W}^{\text{out}}_{l,p}),(S35) where⊙\odotsymbol indicates an element-wise product,MLP1\text{MLP}_{1}andMLP2\text{MLP}_{2}are multi-layer perceptron (MLP) layers with depth and activation functions specified in TableS1,Wl,pin\textbf{W}^{\text{in}}_{l,p}andWl,pout\textbf{W}^{\text{out}}_{l,p}are learnable weight matrices acting on each(l,p)(l,p),Cl1​m1,l2​m2l​mC^{lm}_{l_{1}m_{1},l_{2}m_{2}}is the Clebsch-Gordan coefficient, andδji\delta^{i}_{j}is the Kronecker delta function.

Updating Hidden Representations

A message passing layer updates the hidden representationhAt\textbf{h}^{t}_{A}of the atomAAwith the reverse-matched aggregated messageρ†​(m~At)\rho^{\dagger}(\tilde{\textbf{m}}^{t}_{A})by:

hAt+1=ϕ​(hAt,ρ†​(m~At)),\textbf{h}^{t+1}_{A}=\phi(\textbf{h}_{A}^{t},\rho^{\dagger}(\tilde{\textbf{m}}^{t}_{A})),(S36) wherehAt+1\textbf{h}^{t+1}_{A}is the updated hidden representation.

An atom-wise decoding layer decodes atom-wise information with local interactions, which is,

hAt+1=ϕ​(hAt,hAt).\textbf{h}^{t+1}_{A}=\phi(\textbf{h}^{t}_{A},\textbf{h}^{t}_{A}).(S37) When training OrbitAll on the OMol25 dataset, we found that the quadratic scaling induced by the self Clebsch-Gordan coupling in Eq.S33led to instability in training when multiple decoding layers are stacked, particularly given the diversity of the dataset. To mitigate this issue, we skip the Clebsch-Gordan coupling step during training on OMol25. In equation, we set𝐪l​p​m=𝐟l​p​m\mathbf{q}_{lpm}=\mathbf{f}_{lpm}.

Appendix S2Training

S2.1Workflow

During inference time, the model requires running semi-empirical QM calculations using spGFN1-xTB[neugebauer_spgfnxtb_2023]or g-xTB[froitzheim2025gxtb]for each molecule. However, during training, since the dataset consists of a fixed set of molecules that are used repeatedly, we preprocess all molecules by performing the semi-empirical QM calculations in advance. The resulting features are then used directly during the GNN training phase.

Algorithm 1Training Procedure# 1. Creating orbital features

for

xix_{i},

yiy_{i}in

𝒟train\mathcal{D}_{\text{train}}do⊳\trianglerightxix_{i}: molecular information,yiy_{i}: label,𝒟train\mathcal{D}_{\text{train}}: train dataset

Ti←\textbf{T}_{i}\leftarrowxTB(

xix_{i})⊳\trianglerightTi\textbf{T}_{i}: orbital features

Add

Ti\textbf{T}_{i}to

𝒟train\mathcal{D}_{\text{train}} endfor

for

xjx_{j},

yjy_{j}in

𝒟test\mathcal{D}_{\text{test}}do⊳\triangleright𝒟test\mathcal{D}_{\text{test}}: test dataset

Tj←\textbf{T}_{j}\leftarrowxTB(

xjx_{j})

Add

Tj\textbf{T}_{j}to

𝒟test\mathcal{D}_{\text{test}} endfor

# 2. Train and test neural network

ℱ←\mathcal{F}\leftarrowfit(

ℱ,𝒟train\mathcal{F},\mathcal{D}_{\text{train}})⊳\trianglerightℱ\mathcal{F}: UNiTE

y^test←\hat{\textbf{y}}_{\text{test}}\leftarrowpredict(

ℱ,𝒟test\mathcal{F},\mathcal{D}_{\text{test}})⊳\trianglerighty^test\hat{\textbf{y}}_{\text{test}}: test set predictions

S2.2Training Configuration

All orbital features are generated either using the tblite package[tblite]with spGFN1-xTB[neugebauer_spgfnxtb_2023]or the binary executable of g-xTB[froitzheim2025gxtb]. For the experiments in Sections2.3and2.4, we used the default hyperparameters used in the previous OrbNet-Equi work[qiao_orbnet_equi_2022](‘Small’ in TableS1). In contrast, for the large-scale experiments in Section2.2, several hyperparameters are modified to increase model capacity (‘Large’ in TableS1).

We used the Adam optimizer[kingma_adam_2014]with a learning rate schedule consisting of a linear warm-up phase followed by cosine annealing. The corresponding hyperparameters are detailed in TableS1. Additionally, the smoothL1Loss loss function was used[SmoothL1Loss], which is a loss function with a quadratic slope below a certain threshold and a linear slope above the threshold.

Table S1:Hyperparameters for OrbitAll used for experiments.DescriptionSmallLargeNode hidden dimension256256Number of channels for each (l,pl,p)Shown in TableS2Shown in TableS2Number of message-passing update steps48Number of point-wise decoding steps44Number of convolution channels816Number of attention heads816Depth of MLPs22Activation functionSwishSwishNumber of radial basis functions1632Stability factorϵ\epsilonin EvNorm layers0.10.1Max neighbors6464Maximum learning rate0.00050.0008Warm-up ratio0.3330.0835Number of epochs30080Batch size64128Stable decoding layerFalseTrueMax gradient norm-100.0Feature normalization-LayerNormEncoding normalizationBatchNormLayerNormMessage-passing node normalization-LayerNormMessage-passing normalizationLayerNormLayerNormDecoding normalizationBatchNormLayerNormMax gradient norm-100.0Total number of learnable parameters2.1M7.5MTable S2:Number of channels for each (l,pl,p).

S2.3The OMol25 Dataset

For training OrbitAll on the OMol25 dataset, we employ the g-xTB framework[froitzheim2025gxtb]. As the current implementation of g-xTB is distributed as a binary executable, direct access to orbital-level quantities, such as the Fock, density, overlap, and core Hamiltonian matrices, is not readily available. To address this limitation, we reconstruct the orbital features from the wavefunction exported in Molden format, which provides access to the converged electronic structures. Using the parsed molecular orbital coefficients, we then build the required orbital features with the PySCF package[sun2020pyscf].

During the construction of the core Hamiltonian matrices using the PySCF package, we observed that their magnitudes were significantly larger than those obtained from spGFN1-xTB as implemented in the tblite package. This difference resulted in feature matrices with widely varying numerical scales, which introduced instability during training. Therefore, to ensure consistency in feature magnitudes across different orbital features, we normalize the core Hamiltonian by the total nuclear charge of the system. Specifically, for moleculeii, we define

𝐇~core,i=𝐇core,i∑AiZAi,\tilde{\mathbf{H}}_{\text{core},i}=\frac{\mathbf{H}_{\text{core},i}}{\sum_{A_{i}}Z_{A_{i}}},(S38) where∑AiZAi\sum_{A_{i}}Z_{A_{i}}is the total number of protons in the molecule, and𝐇~core,i\tilde{\mathbf{H}}_{\text{core},i}denotes the normalized core Hamiltonian matrix. The normalized core Hamiltonian matrix is used with other orbital features during training and evaluating OrbitAll on the OMol25 dataset.

S2.4The T1x-Solv Dataset

For fine-tuning the pretrained OrbitAll model on the T1x-Solv dataset, we unfreeze all model parameters and train for 150 epochs using the same learning rate as in pretraining, 0.0008. We use the GBE (generalized Born with finite epsilon) implementation in g-xTB v2.0.0[froitzheim2025gxtb].

During inference, because only relative energies are relevant, we use the vacuum element-wise shifts for unseen solvents. As a result, absolute energy predictions for unseen solvents are expected to be inaccurate.

S2.5QM9star

The QM9star dataset consists of neutral singlets (neutrals), neutral doublets (radicals),+1+1charged singlets (cations), and−1-1charged singlets (anions). The molecular properties are calculated at the B3LYP-D3(BJ)/6-311+G(d,p) level of theory[stephens1994b3lyp,becke1988becke88,lee1988lyp], where 36 or 39 properties are available for each molecule depending on its spin state, including internal energy (U0U_{0}), heat capacity (CvC_{v}), HOMO energy, LUMO energy, zero-point vibrational energy (ZPVE), dipole moment (μ\mu), and isotropic polarizability (α\alpha). Although the dataset also contains force labels, the magnitudes are practically too tiny for any evaluation to be performed. Hence, we skip the force predictions for the open-shell or charged molecules task for this work.

Figure S7:(a) SpookyNet embedding style for embedding spin and charge information along with the atomic numbers. We apply the SpookyNet embeddings to all the models except for OrbitAll to incorporate spin and charge information[unke_spookynet_2021]. (b) OrbitAll embedding style for embedding electronic structure information, which contains the spin and charge information inherently.Except for SpookyNet[unke_spookynet_2021]and EquiformerV2[liao2023equiformerv2], all models were implemented and trained using the Geom3D benchmark repository[liu2023geom3d], including SchNet[schutt2018schnet], PaiNN[schutt2021painn], DimeNet++[gasteiger2022dimenetpp], and NequIP[batzner2022nequip]. For EquiformerV2, we used the same hyperparameters as in their QM9 experiments[ramakri_qm9_2014,liao2023equiformerv2].

In all models, element-wise biases (or shifts) were initialized via linear regression using only the corresponding training set in each run. The QM9star dataset consists of molecules with total charges of +1, 0, and -1. Since our strategy for initializing the charge shift does not allow other charges than those in the training set distribution, the extrapolation to out-of-distribution charged molecules (e.g., +2, -2, etc.) is expected to fail. This issue can be fixed either by preparing species with various charges in the training set or by implementing appropriate modeling of the charge biases, i.e., the average charge effects.

During training, validation errors were monitored, and the model checkpoint with the lowest validation error was selected for testing. Validation frequencies followed those specified in the original papers (e.g., SpookyNet was validated every 1,000 mini-batches[unke_spookynet_2021]). Early stopping with a patience of 150 epochs was applied to all models unless otherwise specified. Models that failed to converge were excluded from the reported results.

For training the models other than OrbitAll, we use the SpookyNet embedding style[unke_spookynet_2021]as their embeddings so that the different models can accept spin and charge information. The difference in the embeddings of the graph nodes is shown in FigureS7. The models were trained on one of the A100 or GH200 GPUs.

S2.5.1Ablation Studies

With the QM9star dataset, we conducted ablation studies on OrbitAll. We used the 40K subset training set to assess the effects of the different factors. The results of the ablation studies are presented in TableS3. Details on different methods are described below.

Table S3:Ablation studies on OrbitAll. The numbers are mean absolute errors (MAEs) on the 40K QM9star training subset.##### Layer-wise Aggregation

Layer-wise aggregation in this study refers to the architectural consideration of whether to aggregate contributions from each layer or to use the result from the final layer. If layer-wise aggregation is applied, information from each layer is aggregated, which is pooled by a task-specific operation.

First, we generalize the message passing and decoding procedures. After eachtt-th message passing layer,TdectT_{\text{dec}}^{t}subsequent decoding layers are applied. With an expression ofTdectT^{t}_{\text{dec}}as a lengthttvector,Tdec=[Tdec1,Tdec2,…,Tdect]T_{\text{dec}}=[T^{1}_{\text{dec}},T^{2}_{\text{dec}},...,T^{t}_{\text{dec}}], the OrbNet-Equi architecture with four message passing layers followed by four decoding layers can be expressed asTdec=[0,0,0,4]T_{\text{dec}}=[0,0,0,4].

In this ablation study, we considerTdec=[1,1,1,1]T_{\text{dec}}=[1,1,1,1], i.e., one decoding layer per message passing layer. This setting makes the total number of decoding layers the same as the original OrbNet-Equi architecture[qiao_orbnet_equi_2022].

Empirical Physical Terms

Inspired by SpookyNet[unke_spookynet_2021]and FENNIX[ple2023fennix], we consider adding empirical physical terms, especially the electrostatic correction,EelecE_{\text{elec}}, and the dispersion correction,EdispE_{\text{disp}}. We adopted SpookyNet’s implementation in OrbitAll, where the dispersion correction is the D4 dispersion correction of Grimme et al.[caldeweyher2019d4dispersion]. The calculations of both terms require the prediction of atomic partial charges, which is performed by a neural network.

Atomic partial charges can also be obtained from spGFN1-xTB. Hence, delta-learning can be employed, which means that the neural network predictsΔ​qA=qA+qA,spGFN1-xTB\Delta q_{A}=q_{A}+q_{A,\text{spGFN1-xTB}}for atomAA. The atomic partial charges shall add up to the total charge of molecules, i.e.,

Q=∑AqA,Q=\sum_{A}q_{A},(S39) which applies the same to the atomic partial charges from spGFN1-xTB. Therefore, the sum of delta-partial charges adds up to zero,∑AΔ​qA=0\sum_{A}\Delta q_{A}=0. To enforce this as a hard physical constraint, the predicted delta partial charges,Δ​q^A\Delta\hat{q}_{A}, are obtained by following,

Δ​q^A=q^A−1Natom​(∑Bq^B),\Delta\hat{q}_{A}=\hat{q}_{A}-\frac{1}{N_{\text{atom}}}\left(\sum_{B}\hat{q}_{B}\right),(S40) whereq^A\hat{q}_{A}is the predicted partial charge of the atomAAfrom the neural network, andNatomN_{\text{atom}}is the number of atoms in the molecule. The atomic partial charges are then used for calculations of the electrostatic and dispersion correction terms. The details of implementation can be found in[unke_spookynet_2021].

Attention Renormalization

We apply attention renormalization proposed by Liao et al.[liao2023equiformerv2]. The ablation studies of EquiformerV2 suggested that attention re-normalization can bring some improvements in the performance. Similar to EquiformerV2, an MLP-based attention is used in OrbitAll and OrbNet-Equi. Each attention head at layertt,αA​Bt\alpha^{t}_{AB}, for nodeAAattending to the message from nodeBBis created from a two-layer MLP. When attention renormalization is applied, the features for generating the attention are normalized with a layer norm. Hence, the normalized attention is

αA​Bt=σ​(Linear​(σ​(Linear​(Norm​(x))))),\alpha^{t}_{AB}=\sigma(\text{Linear}(\sigma(\text{Linear}(\text{Norm}(x))))),(S41) whereσ​(⋅)\sigma(\cdot)is an activation function,Linear​(⋅)\text{Linear}(\cdot)is a linear function with learnable weights and biases,Norm​(⋅)\text{Norm}(\cdot)is a normalization layer, andxxis primitive information for attention.

S2.6QMSpin

SpookyNet[unke_spookynet_2021], MOB-ML[cheng_MOBML_2022], and TensorNet[simeon2025tensornet_spincharge]evaluated their models using this dataset with a 20K training set and a 1K validation set.

For calculating the adiabatic spin gaps, molecules with both singlet-optimized and triplet-optimized geometries are required. As in Cheng et al.,[cheng_MOBML_2022], we randomly sampled 1K molecules with singlet- and triplet-optimized geometries (2K geometries in total). Each geometry has singlet and triplet energies, so the test set effectively has 4K data points. Similarly, 250 molecules with singlet- and triplet-optimized geometries (500 geometries, 1K data points in total) were randomly selected to create the validation set. Finally, 20K data points were randomly sampled from all remaining data points to generate the training set.

S2.7Hessian QM9

The Hessian QM9 dataset is a dataset that contains 41.6K molecules in four different solvents: vacuum, water, toluene, and tetrahydrofuran (THF). The molecules are all closed-shell, neutral species. Therefore, this prediction task is possible with the orbital features setT=(F,P,S,Hcore)\textbf{T}=(\textbf{F},\textbf{P},\textbf{S},\textbf{H}_{\text{core}}). However, for consistency throughout this work, we use the orbital features setT=(Fα,Fβ,Pα,Pβ,S,Hcore)\textbf{T}=(\textbf{F}^{\alpha},\textbf{F}^{\beta},\textbf{P}^{\alpha},\textbf{P}^{\beta},\textbf{S},\textbf{H}_{\text{core}}). Note that these are closed-shell species, and thereforeFα=Fβ\textbf{F}^{\alpha}=\textbf{F}^{\beta}andPα=Pβ\textbf{P}^{\alpha}=\textbf{P}^{\beta}.

For feature generation, we use the ALPB method[sigalov2006alpb,ehlert2021alpb2]implemented in the tblite package[tblite], with the same dielectric constants for labels[williams2025hessianqm9]. We decided to use the ALPB method based on the result of the ablation study between different solvation methods as shown in TableS7. For vacuum, no implicit solvation model was used in the calculation. The usage of the implicit solvation model perturbs the electronic structure and the orbital features of molecules. As a visualization, FigureS8displays an example of the perturbation from using the different solvents, where it shows the differences between the density matrices at different solvents to those at vacuum.

We split the dataset into a 32K molecules training set, a 3.2K molecules validation set for training, and the rest (∼\sim6.4K molecules) for testing. Since each molecule has energies for the four solvents, the effective total number of data points is 128K for training, 12.8K for validation, and∼\sim25.8K for testing.

Figure S8:Density matrices difference for a sample molecule in different solvents. For example, the “Toluene - Vacuum” figure illustrates the matrix𝐏Toluene−𝐏Vacuum\mathbf{P}_{\text{Toluene}}-\mathbf{P}_{\text{Vacuum}}. The density matrices are generated using GFN1-xTB with ALPB solvation.

S2.8The Polypeptide Dataset

We sampled the initial atomic coordinates of the polypeptides from the PEPCONF dataset[prasad2019pepconf]. The original PEPCONF dataset consists of∼\sim3.8K relative conformational energy data points calculated at the LC-ω\omegaPBE-XDM/auc-cc-pVTZ level of theory. Since the QM9star dataset labels are calculated at the B3LYP-D3(BJ)/6-311+G(d,p) level of theory, we reoptimized the geometries and calculated the single point energy (EtotE_{\text{tot}}) at the same level of theory.

Geometry optimizations of polypeptides were performed using the Orca 6.0 software[neese2025orca]. The single point calculations for the optimized geometries were performed at the B3LYP-D3(BJ)/6-311+G(d,p)[becke1988becke88,lee1988lyp,stephens1994b3lyp,grimme2010dftd]level of theory using the Psi4 software[smith2020psi4], with the default convergence thresholds, given by10−610^{-6}for both energy and density, and10−1210^{-12}for integrals.

S2.9QM9star in Random Uniform External Electric Fields

From the original QM9star dataset[tang2024qm9star], we collect 10K training, 1K validation, and 2K test data points to study the behavior of OrbitAll in random uniform external electric fields. With the data points, we created 3 distinct datasets. First, the dataset with no electric field applied. The task is hence to predict the unperturbed single-point energies. For this dataset, we collected the single-point energies provided by the original QM9star dataset. Second, the dataset has uniform external electric fields in random directions and magnitudes, with a maximum magnitude of 0.005 atomic units (au). Finally, the dataset has uniform external electric fields in random directions and magnitudes, with a maximum magnitude of 0.0005 au. To obtain the single point energies with uniform external electric fields, we used the psi4 software[smith2020psi4], and created the labels at the B3LYP-D3(BJ)/6-311+G(d,p) level of theory, which is the same level as the QM9star dataset. Note that we excluded radicals (Q=0Q=0,S=1/2S=1/2) and collected only neutrals (Q=0Q=0,S=0S=0), anions (Q=−1Q=-1,S=0S=0) and cations (Q=+1Q=+1,S=0S=0).

S2.10Inference Time

Inference times for all machine learning models were measured on an NVIDIA RTX 4090 GPU, using the mean single-batch inference time as the reported metric. For the QM9star calculation times, we sampled 500 randomly selected molecules using Psi4[smith2020psi4]. For measuring spGFN1-xTB using tblite[tblite]and the QM9star calculation times, an AMD Ryzen 5600H CPU was used. For polypeptide dataset molecules, computations were performed on two 16-core Intel Skylake 2.1GHz CPUs. Notably, for the 500 QM9star molecules, the execution times on the Ryzen 5600H were found to be comparable to those on the Skylake CPUs.

Appendix S3Additional Results

S3.1The OMol25 Dataset

Figure S9:Validation-set MAE grouped by molecular charge, ranging from−10-10to+10+10. The value below each charge label indicates the number of samples in that charge group, e.g., n=828 for charge+10+10.FigureS9shows the MAE across molecular charge states sampled from the validation set. Overall, OrbitAll achieves lower MAEs for highly charged species, particularly for cations with charges of +5 and above, and for anions with charges of -4 or below, except at charge -10. The improvement is especially clear for highly charged cations, where OrbitAll is more robust than UMA.

This trend may reflect the role of the g-xTB baseline and the resulting orbital features. For cations, g-xTB appears to provide a sufficiently accurate electronic-structure baseline, enabling OrbitAll to learn more accurate corrections. In contrast, anions are generally more challenging because their electron densities are more diffuse[lynch2003anions_diffuse], which can be difficult to represent with the minimal basis used in g-xTB. Nevertheless, the orbital features still appear to provide useful information for modeling highly charged species, contributing to the improved performance of OrbitAll across most extreme charge states.

Additionally, OrbitAll achieves MAEs of 34.03 meV for “Biomolecules”, 69.50 meV for “Electrolytes”, and 131.99 meV for “Metal Complexes”. Although direct comparison is not straightforward because our training and evaluation sets differ from those used in prior OMol25-4M benchmarks (e.g., lanthanides are excluded in our setting), these errors are competitive with those reported for leading MLIPs trained on OMol25-4M. For example, GemNet-OC trained on OMol25-4M reports MAEs of 39.58 meV for “Biomolecules”, 56.32 meV for “Electrolytes”, and 148.49 meV for “Metal Complexes”[levine2025openmolecules2025omol25]. Notably, OrbitAll achieves lower errors for biomolecules and metal complexes, while remaining within a similar range for electrolytes. These results suggest that OrbitAll can model diverse and challenging chemical systems, including biomolecules, electrolytes, metal complexes, and non-equilibrium geometries, demonstrating its robustness and scalabil/ity.

S3.2Solvation Reactions

Figure S10:Evaluation of the generalization of OrbitAll to unseen solvent conditions. The top two panels, corresponding to benzene and isopropanol, show cases where OrbitAll generalizes well, whereas the bottom two panels, corresponding to acetone and acetonitrile, show cases where OrbitAll fails to generalize.FigureS10shows the zero-shot predictions of OrbitAll fine-tuned on the T1x-Solv dataset, which contains data in vacuum, water, methanol, and toluene. We observe two cases where the model generalizes successfully and two cases where it does not. The successful cases, benzene and isopropanol, are similar to solvents included in the training dataset. Benzene is similar to toluene, a nonpolar solvent without hydrogen-bonding capability, whereas isopropanol is similar to water and methanol, which are polar hydrogen-bonding solvents. In contrast, the two failure cases, acetone and acetonitrile, are both polar aprotic solvents that do not donate hydrogen bonds. SMD is a more sophisticated implicit solvation model than GBE, the base implicit solvation method used during the g-xTB data processing. In particular, SMD includes solvent-specific parameters that account for hydrogen-bonding effects, whereas GBE does not explicitly include such terms. Therefore, in future work, if the base method supports a more sophisticated implicit solvation model such as SMD, OrbitAll may generalize better to unseen solvents.

Figure S11:Transition-state (TS) search under implicit solvation for the cyclopentadiene (CP) + methyl vinyl ketone (MVK) Diels-Alder reaction.(A)Reaction schematic for the CP-MVK Diels-Alder reaction.(B)Schematic and 3D example of the transition state, with the forming bond lengths (rCC1r_{\text{CC1}}andrCC2r_{\text{CC2}}) labeled.(C)NEB-based transition-state search for the CP-MVK Diels-Alder reaction for different solvents. “Ref.” denotes reference calculations at theω\omegaB97M-V/def2-TZVPD/SMD level of theory. For each solvent, the estimated reaction barriers (TS) and reaction energies (P) obtained from different methods are labeled below each state. Overlaid 3D transition-state structures from each method are shown, with Ref in gray, OrbitAll in blue, and UMA in orange.rCC1r_{\text{CC1}},rCC2r_{\text{CC2}}, and the RMSD relative to the reference structure are reported in Å, in the inset tables.We present additional results for transition state searches under implicit solvation for the Diels-Alder reaction between cyclopentadiene (CP) and methyl vinyl ketone (MVK) in FigureS11. Similar to the results shown in Figure2, OrbitAll accurately predicts both the energetics and TS structures of the reaction under four conditions: vacuum, water, toluene, and methanol. All predicted energetics are near or within chemical accuracy (1 kcal/mol), and the TS structures are highly accurate, with RMSD values of only 0.02 Å relative to the reference structures (ω\omegaB97M-V/def2-TZVPD/SMD).

OrbitAll also accurately captures solvent effects. It correctly predicts thatrCC2r_{\text{CC2}}increases in water and methanol compared with vacuum and toluene. In terms of energetics, it also predicts that toluene increases the reaction barrier relative to vacuum, whereas water and methanol lower the barrier.

Appendix S4Transition-state-informed Umbrella Sampling

For umbrella sampling to be reliable, the sampled distributions of the collective variable from neighboring windows must exhibit sufficient overlap. We empirically find that a harmonic restraint force constant of10​eV/Å210~\mathrm{eV}/\text{\AA }^{2}provides good overlap between adjacent windows and enables stable sampling across the reaction coordinate. Representative window histograms for water are shown in FigureS12, and the resulting potential-of-mean-force profiles alongξ\xifor different solvents are shown in FigureS13.

Figure S12:Distribution histograms of the collective variableξ\xifor the umbrella-sampling windows in water.Figure S13:Potential-of-mean-force profiles along the collective variableξ\xifor different solvents.### S4.1QM9star

We present the numerical values of the experiments performed using the QM9star dataset in TableS4.

Table S4:Mean absolute error (MAE) in meV for the QM9star dataset. For the total energyEtotE_{\text{tot}}prediction task, the models were trained using the 400K training subset of the QM9star dataset, with 100K for each category of species. For HOMO and LUMO prediction tasks of radical species forα\alphaandβ\betaspins, a 100K training set of radical species is used.EtotE_{\text{tot}}EHOMOαE^{\alpha}_{\text{HOMO}}ELUMOαE^{\alpha}_{\text{LUMO}}EHOMOβE^{\beta}_{\text{HOMO}}ELUMOβE^{\beta}_{\text{LUMO}}ParamsAllNeutralRadicalCationAnionRadicalSchNet-SC(D)75.2642.4762.95103.4384.51----617KPaiNN-SC(D)41.6323.4034.6756.6647.51----723KDimeNet++-SC(D)19.417.4614.4230.3422.59----1.4MSphereNet-SC(D)17.937.4913.7227.5620.49----2.0MEquiformerV2-SC(D)15.436.2210.6025.3017.45----9.5MNequIP-SC(D)36.9220.4128.5852.1042.67----1.4MOrbitAll(D)14.757.8911.5720.8017.10----2.1MSchNet-SC(Δ\Delta)33.3118.5729.1542.6839.39151.16129.4379.1957.36617KPaiNN-SC(Δ\Delta)25.3614.0922.4031.6230.67148.37125.0075.1946.23723KDimeNet++-SC(Δ\Delta)13.275.9710.4716.9717.9474.4159.7937.9126.771.4MSphereNet-SC(Δ\Delta)12.945.6510.2216.6917.4990.1580.1849.2230.932.0MEquiformerV2-SC(Δ\Delta)9.863.447.1313.6613.7065.9649.3838.6623.599.5MSpookyNet(Δ\Delta)19.2710.7615.4323.9724.93----3.7MNequIP-SC(Δ\Delta)23.0612.7819.7329.2228.09----1.4MOrbitAll(Δ\Delta)8.504.866.8910.6010.7855.4043.3132.2714.732.1M

S4.2QMSpin

Table S5:Mean absolute error (MAE) in kcal/mol on total molecular energy when trained with 10,000 molecules in total on different training sets of the QMSpin dataset.TableS5shows the results from training two OrbitAll models: (1) one using 10K singlets (Singlets-only) as a training set and 500 singlets as a validation set, and (2) one using 5K singlets and 5K triplets (Singlets + Triplets) as training set and 500 singlets and 500 triplets as a validation set, where all are selected randomly. We ensure geometries are not shared between training, validation, and test sets. Both models are tested against three sets: (1) a combined set of 2,316 remaining singlets and triplets (all), (2) the 2,316 singlets, and (3) the 2,316 triplets. We observe that OrbitAll can learn from the combined electronic structures of singlets and triplets. Specifically, the MAE is significantly improved when using both singlets and triplets as training data.

Table S6:The mean absolute errors (MAEs) for the test set in meV on total molecular energyEtotE_{\text{tot}}when trained with 20K molecules for different models using the QMSpin dataset (MRCISD+Q-F12/cc-pVDZ-F12)[schwilk2020qmspin]. The MAEs are taken from original works[unke_spookynet_2021,simeon2025tensornet_spincharge,cheng_MOBML_2022]. OrbitAll and MOB-ML were trained on delta-labels[cheng_MOBML_2022]. “All” is the combined test set of singlet (S=0S=0) and triplet (S=1S=1). SpookyNet[unke_spookynet_2021]and TensorNet[simeon2025tensornet_spincharge]reported “all” MAEs only. The best MAE for each category is highlighted with boldface.We assess the performance of OrbitAll using the QMSpin dataset[schwilk2020qmspin], which contains open-shell species. The model predicts total energies of singlet (ES=0E^{S=0}) and triplet (ES=1E^{S=1}) carbenes at the MRCISD+Q-F12/cc-pVDZ-F12 level of theory. From these predictions, we compute vertical spin gaps at singlet-optimized (EsingletgapE^{\text{gap}}_{\text{singlet}}) and triplet-optimized (EtripletgapE^{\text{gap}}_{\text{triplet}}) geometries, as well as the adiabatic spin gaps (EadiabaticgapE^{\text{gap}}_{\text{adiabatic}}), as defined in[cheng_MOBML_2022]as follows:

Esingletgap=EsingletS=1−EsingletS=0,Etripletgap=EtripletS=1−EtripletS=0,Eadiabaticgap=EtripletS=1−EsingletS=0,\begin{split}&E^{\text{gap}}_{\text{singlet}}=E^{S=1}_{\text{singlet}}-E^{S=0}_{\text{singlet}},\\ &E^{\text{gap}}_{\text{triplet}}=E^{S=1}_{\text{triplet}}-E^{S=0}_{\text{triplet}},\\ &E^{\text{gap}}_{\text{adiabatic}}=E^{S=1}_{\text{triplet}}-E^{S=0}_{\text{singlet}},\\ \end{split}(S42) where the right-hand terms are total energies, with subscripts “singlet” and “triplet” indicating the multiplicity of the optimized geometry, and superscripts “S=1S=1” and “S=0S=0” indicating the spin state at which the energy is computed. The sampling strategy for the training, validation, and test sets is described in Section4.2.

As shown in TableS6, OrbitAll achieves the lowest MAEs across all categories except the singlet state, with both singlet and triplet errors smaller than chemical accuracy (1 kcal/mol≈\approx43.4 meV). This is significant given that the dataset was computed using the highly accurate but computationally demanding MRCISD+Q-F12/cc-pVDZ-F12 method. In contrast to MOB-ML, which trains separate models for singlet and triplet states using features from RHF and ROHF[cheng_MOBML_2022], OrbitAll employs a single model trained jointly on both, using the significantly cheaper spGFN1-xTB method. This unified training likely contributes to OrbitAll’s robust performance across spin states, whereas MOB-ML struggles with triplet predictions. Additionally, OrbitAll achieves the highest accuracy in predicting both vertical and adiabatic spin gaps, with all MAEs smaller than chemical accuracy. These results underscore the strength of OrbitAll’s unified representation in capturing and generalizing across different spin multiplicities.

Figure S14:Error distribution histogram (left) and absolute error log distribution histogram (right) using OrbitAll trained with 10K singlets and 10K triplets, for predicting singlets and triplets in the test set of the QMSpin dataset.FigureS14shows the error distribution histograms in terms of the errors and the absolute errors. The inference is done on the test set of the QMSpin dataset, with the OrbitAll model trained using 10K singlets and 10K triplets. We observe that OrbitAll is able to describe triplets and singlets almost equally well for the QMSpin dataset.

S4.3Hessian QM9

We additionally tested OrbitAll using different solvation methods in GFN1-xTB for generating orbital features. Specifically, we used the CPCM[barone1998cpcm1,takano2005cpcm2], ALPB[sigalov2006alpb], and GBSA[qiu1997gbsa]implicit solvation methods implemented in the tblite package[tblite]. For this test, we used a subset of the Hessian QM9 dataset split into 10K training, 1K validation, and 2K test molecules, with each molecule associated with four solvent data points: vacuum, toluene, THF, and water. For CPCM, we used the same dielectric constants as those used to generate the SMD-solvated data in Hessian QM9, whereas for ALPB and GBSA, we used the empirical solvent parameters implemented in tblite.

Table S7:Mean absolute error (MAE) for the test set in kcal/mol for total molecular energy when trained on 10,000 molecules, corresponding to 40,000 data points across four solvents, from the Hessian QM9 dataset. The implicit solvation method used to generate the orbital features is indicated in parentheses.Interestingly, the results varied substantially depending on the implicit solvation method used to generate the orbital features. Based on this ablation study, we selected the ALPB method for evaluating OrbitAll on the Hessian QM9 dataset.

S4.4Radical Species

Figure S15:Learning curves of different models for the radicals dataset.

Figure S16:Error distribution histogram (left) and absolute error log distribution histogram (right) using OrbitAll trained with 260K molecules for predicting singlets and doublets in the test set of the radicals dataset from St. John et al.[stjohn2020radicals_dataset].In addition to the evaluations we presented in the main text, we also present the evaluation we performed on the dataset reported by St. John et al.[stjohn2020radicals_dataset]. The dataset consists of∼\sim240K radicals and∼\sim40K closed-shell molecules. Specifically, we train on the SCF energies of each molecule for evaluation. Out of the entire dataset, we randomly sampled 260,000 molecules for training, 2,048 for validation, and the remaining 27,589 for testing. Then, we randomly sample a subset from the training set to create the learning curve.

The result implies that training for doublets is certainly more difficult than training for singlets. Although there were many more doublets in the dataset than singlets, we see that generally, singlets record lower MAEs at all training sizes.

S4.5QM9star in Random Uniform External Electric Fields

Figure S17:Learning curves of OrbitAll for different maximum magnitudes of uniform external electric fields.Predicting the molecular energies perturbed by uniform external electric fields is observed to be significantly more challenging than those without external fields. In FigureS17, we observe the higher offset and flatter slopes for both of the learning curves of different magnitudes of uniform external electric fields. The negative slopes for both cases of maximum magnitudes (0.005 au and 0.0005 au) imply that OrbitAll is capable of learning molecules in arbitrary uniform electric fields. However, the learning curves show that assuming the trend continues in a higher data regime, we require several hundred thousand data points to achieve chemical accuracy.

S4.6Delta-learning

Figure S18:The label energy distribution histograms of different species for the QM9star dataset[tang2024qm9star], for direct-learning (left) and delta-learning (right). The delta-labels are obtained from subtractions by spGFN1-xTB energies[neugebauer_spgfnxtb_2023]. The labels here are the original labels subtracted by the element-wise energy bias initializations.Delta-learning is a widely used strategy in the area of chemical properties prediction tasks[rama_delta_learning_2015], which has been especially successful for energetic properties predictions[ruth2022deltalearning1,chen2023deltalearning2,zhu2019deltalearning4]. Delta-learning is known empirically to reduce test error by a certain factor. Another study has shown that delta-learning can also account for long-range interactions[böselt2021deltalearning3_lr]. Since accuracy also relates to data-efficient training, delta-learning can be particularly useful for highly expensive labels, such as energies at the CCSD(T) level of theory[ruth2022deltalearning1].

FigureS18displays the data distribution pattern of different species in the QM9star dataset[tang2024qm9star]. The distributions of direct-learning labels of different species are more dispersed. However, the distributions of the delta-learning labels appear to be more predictable and aligned, with sharper peaks around their corresponding means. Specifically, the anions and cations distributions become equidistant from the origin (E=0E=0), where the center of distribution of the neutral species is located, with opposite signs. These changes in distributions facilitate generalizations, reducing systematic error patterns.

Appendix S5Limitations and Future Works

Our feature generation via the semi-empirical method in tblite does not support analytical gradients[tblite]. Thus, OrbitAll cannot compute forces via backpropagation through atomic positions, unlike fully differentiable models. Obtaining forces from energy gradients ensures conservative force fields, which are necessary for stable geometry optimization and molecular dynamics[bigi2025darkforcesassessingnonconservative]. Qiao et al. addressed this by computing numerical gradients of the intermediate featuresTwith respect to atomic positions,∂T/∂ri{\partial\textbf{T}}/{\partial r_{i}}[qiao_orbnet_equi_2022]. However, this is computationally expensive and can lead to excessive memory usage if cached. Addressing this with analytical gradients[qiao2020orbnetanalytical]or a differentiable semi-empirical backend[friede_dxtb_2024]is possible in future work.