Electron-Density-Informed Molecular Property Prediction Using Surrogate Functionals

1 INTRODUCTION

• Ground-state electron density is a physical descriptor because it uniquely determines the external potential, and thus in principle all ground-state properties.Source: §1 · p. 1
• The framework constructs a learned functional whose minimum is the ground-state density and whose value there equals the ground-state energy.Earlier surrogate-functional work reported improved extrapolation to larger systems.Source: §1 · p. 1
• The paper extends surrogate functionals from electronic energy to arbitrary scalar molecular properties, replacing the energy target with the property label.At inference, a separately trained network predicts the required density from molecular structure, avoiding variational optimization.Source: §1 · pp. 1–2
• Surrogate functionals surpass prior state-of-the-art results on six of seven extensive QM9 properties and outperform geometry-only and density-input models across extensive properties, with a larger advantage on QM9-trained evaluation of QM40.The controlled comparisons use a common backbone.Source: §1 · p. 2
• The authors introduce beSEN, an eSEN modification with Gram-matrix-based bilinear activation that further improves prediction, particularly within the surrogate-functional framework.Source: §1 · p. 2

2 BACKGROUND AND RELATED WORK

• Structure-only models advanced from invariant message passing (Sch¨utt et al., 2018; Gasteiger et al., 2020a; Liu et al., 2022) to equivariant architectures, reducing QM9 energy errors to a few meV relative to reference calculations.Source: §2 · p. 2
• Brockherde et al. (2017) and Bogojeski et al. (2020) use staged electronic-structure mappings, while Shao et al. (2023) extend the approach to the one-electron reduced density matrix.Brockherde et al. predict a density expansion before fitting energy; Bogojeski et al. use a DFT density for a coupled-cluster correction.Source: §2 · p. 2
• DenSNet (Bogojeski et al., 2026) predicts atom-centered Gaussian density coefficients for an energy model, while DensIP (Wing et al., 2026) uses a physics-based intermolecular potential in the second stage.Source: §2 · p. 2
• Orbital-free DFT minimizes an energy functional with respect to density without constructing Kohn–Sham orbitals; surrogate functionals use ground-state labels alone and have shown better generalization than direct energy prediction.Recent machine-learned functionals have reported chemical accuracy for small druglike molecules.Source: §2 · p. 2

3.1 DENSITY REPRESENTATION

• Density is expanded in atom-centered Gaussian basis functions and symmetrically Löwdin-orthonormalized, yielding coefficients that transform as SO(3) irreducible representations.These coefficients serve directly as equivariant message-passing node features.Source: §3.1 · p. 3

3.2 ORBITAL-FREE SURROGATE FUNCTIONALS

• Surrogate-functional training maintains cached density coefficients per molecule and updates them by gradient descent when predicted gradients point toward the ground state.The GDI loss requires each step to contract the distance to the true ground state; an energy-bound loss and ground-state energy and gradient penalties further constrain the functional.Source: §3.2 · pp. 3–4

3.3 SURROGATE FUNCTIONALS FOR MOLECULAR PROPERTIES

• For each scalar property, the method trains a dedicated functional by replacing the energy target with that property’s ground-state value.The construction carries over the training pipeline used for surrogate energy functionals.Source: §3.3 · p. 4

3.4 ARCHITECTURE

• The beSEN architecture modifies eSEN’s message passing, node update, and readout, and combines geometry and density embeddings into initial node features.The same beSEN architecture is used for density prediction and property prediction.Source: §3.4 · p. 4, Figure 2 · p. 4, §3.4.1 · p. 5

3.4.2 MESSAGE PASSING

• beSEN uses post-norm residual message passing and a frozen, first-batch per-channel L2 normalization during training and inference.The normalization omits the explicit metric used in the related BOA normalization.Source: §3.4.2 · p. 5

3.4.3 NODE UPDATE; 3.4.4 READOUT

• For node updates, beSEN omits pre-norm and uses a Gram activation to gate tensorial features; its readout also uses the node-feature Gram matrix.The Gram activation derives gating scalars from a subset of the feature-channel Gram matrix, while the readout feeds the Gram matrix to its MLP.Source: §3.4.3 · p. 5, §3.4.4 · p. 6

3.5 SPEED COMPARISON

• beSEN is slightly faster than eSEN on QM9 despite having 6.9 M versus 2.7 M parameters.The comparison measures forward-pass throughput with batch size 128 on one NVIDIA A100 40GB; beSEN’s additional parameters are concentrated in cheaper per-node stages.Source: Table 1 · p. 6, §3.5 · p. 6, Figure 3 · p. 6

3.6 PREDICTING THE GROUND-STATE ELECTRON DENSITY DIRECTLY

• Directly predicting density coefficients with a second network gives comparable performance while reducing inference to two model forward calls.This avoids the roughly 100 forward calls required for variational minimization and avoids explicit orthonormalization of the density coefficients.Source: §3.6 · p. 6

4.1 QM9 DATASET

• On QM9, the surrogate-functional model beats geometry-only and density-input beSEN on all seven extensive properties, and beats GotenNet-L on six.The seven properties are U0, U, H, G, ZPVE, α, and Cv.Source: §4.1 · p. 7, Table 2 · p. 7
• On U0, beSEN outperforms eSEN in geometry-only, density-input, and surrogate-functional modes, with the largest gap in surrogate mode.The comparison reports MAE in meV on the QM9 test split.Source: §4.1 · p. 7, Table 3 · p. 7
• The direct density predictor reaches 0.032% NMAE on QM9, and using predicted rather than label densities only marginally affects property results.The density error is measured against labels after density fitting.Source: §4.1 · p. 8, Table 4 · p. 8

4.2 GENERALIZING TO LARGER MOLECULES

• On 600 QM40 molecules, the surrogate model outperforms geometry-only and density-input beSEN by a large margin on all extensive properties.All models were trained only on QM9; the QM40 molecules have 10–40 heavy atoms.Source: §4.2 · p. 8, Table 5 · p. 8
• For U0 on QM40, surrogate beSEN maintains accuracy better as molecule size increases, yielding much lower median errors for larger molecules.The comparison plots absolute error per atom against heavy-atom count.Source: §4.2 · p. 9, Figure 4 · p. 9

5 DISCUSSION

• The framework extends surrogate functionals from electronic energy to arbitrary scalar molecular properties, with improved extensive-property performance and generalization to larger molecules.The work also introduces beSEN, which performs especially well with the surrogate-functional framework.Source: §5 · p. 9
• The exact mechanism of the improvement is not yet fully understood and is identified for future investigation.The motivation for surrogate functionals is less obvious for arbitrary properties than for variational electronic energy.Source: §5 · p. 9

A.1 ABLATION EXPERIMENTS

• Every ablated U0 configuration has higher error than the full setup combining all training components.The comparison is on the QM9 test split.Source: §A.1 · p. 13, Table 6 · p. 13
• The ablation varies density input, stationarity, energy-bound and GDI losses, and cache use and initialization.The stationarity constraint is ∇pF(p*) = 0 at the ground state.Source: Table 6 · p. 13
• The cache may be initialized with MINAO or with a ground-state density displaced along an electron-count-conserving direction.For the latter, the L2 displacement radius is drawn log-uniformly over [0.002, 0.3].Source: Table 6 · p. 13

A.2 SIZE EXTRAPOLATION FOR ALL TARGET PROPERTIES

• Figure 5 plots median absolute error per atom against heavy-atom count for all twelve target properties and three beSEN variants.The variants are geometry-only, density-input, and surrogate-functional.Source: Figure 5 · p. 14, §A.2 · p. 14
• The leftmost bin is the QM9 test set, while the remaining bins contain the 600 QM40 molecules.The models were trained on QM9 only.Source: Figure 5 · p. 14
• The plotted QM40 molecules are grouped by heavy-atom count, and the reported errors are divided by total atom count.The figure caption identifies median absolute error per atom as the plotted measure.Source: Figure 5 · p. 14, §A.2 · p. 14

• QM9 and QM40 target properties are taken directly from their respective datasets.Source: §A.3 · p. 14
• Ground-state densities are generated with the authors’ Kohn–Sham DFT calculations at B3LYP/6-31G(2df,p), the level used to generate QM9.Source: §A.3 · p. 14
• Density coefficients use an even-tempered Gaussian basis with exponent ratio 2.5, are obtained by density fitting, and are then orthonormalized.Source: §A.3 · p. 14

A.4.1 OUTPUT PARAMETRIZATION

• Extensive-property outputs use least-squares-fitted element-specific atomic references and divide summed per-atom contributions by a fixed normalization constant.Source: §A.4.1 · p. 15
• Intensive-property outputs mean-pool node features before readout so predictions do not grow with molecular size.Source: §A.4.1 · p. 15
• Dipole moment and electronic spatial extent use dedicated output heads, with atomic charges, dipoles, and scalars predicted from node features.Source: §A.4.1 · p. 15

A.4.2 TARGET SCALING

• Target scaling aims to make surrogate-loss hyperparameters transferable across properties with different units and spreads; extensive-property standard deviations are computed after subtracting the atomic reference.Source: §A.4.2 · p. 15

A.5 HYPERPARAMETERS

• Architecture and training settings are shared by all models, while the surrogate-training settings apply only to surrogate models.Source: §A.5 · p. 16, Table 7 · p. 16, Table 8 · p. 16