All-electron machine-learned orbital-free density functional for periodic CHNO systems

1 Introduction

• OF-DFT minimizes an explicit density to obtain the ground-state density and energy, avoiding repeated Kohn–Sham Hamiltonian diagonalization.Transferable kinetic-energy approximations remain difficult, especially for covalent bonding and the non-ground-state densities encountered during optimization.Source: §1 · pp. 1–2
• Small energy interpolation errors alone do not ensure accurate functional derivatives for variational optimization.The relevant test is whether a learned functional can be minimized from an initial density to recover a meaningful density and energy.Source: §1 · p. 2
• Earlier learned periodic OF-DFT demonstrations used local pseudopotentials, while all-electron periodic treatment must resolve core densities and periodic electrostatics.Modern transferable pseudopotentials are generally nonlocal, with energies dependent on the one-particle density matrix rather than density alone.Source: §1 · pp. 2–3
• Generalizing variational all-electron learned functionals from molecules to three-dimensional periodic crystals was an unmet capability.Prior work had established periodic learned OF-DFT with pseudopotentials and all-electron learned functionals for finite molecules.Source: §1 · p. 3
• STRUCTURES-P combines analytic periodic electrostatics with a learned kinetic-plus-exchange-correlation energy, demonstrating variational optimization of explicit all-electron densities in periodic crystals.The functional is trained on non-ground-state energies and density-coefficient gradients, and is minimized at fixed electron number.Source: §1 · pp. 3–4

2.1 A variational all-electron functional for periodic crystals

• One periodic all-electron functional combines analytic electrostatics with an equivariant network learning kinetic and exchange-correlation energy, then optimizes density coefficients at fixed electron number.Training uses energies and density-coefficient gradients from non-ground-state densities; updates begin from a superposition of atomic densities and preserve unit-cell neutrality.Source: §2.1 · pp. 4–6

2.2 Evaluation setup

• The 62,109-cell dataset spans even-electron closed-shell CHNO crystals, and all in-distribution validation and test cells have explicit unsmeared ground-state references.The cells include finite bonded components, chains, sheets, three-dimensional networks, and mixed connectivity classes.Source: §2.2 · p. 6, Figure 2 · p. 6

2.3 In-distribution variational optimization

• After optimization, energy MAE is 1.58 mHa/cell (0.254 mHa/atom), and mean electron-normalized L1 density error falls from 6.73% to 0.413%.The energy MAE is reduced by two orders of magnitude; Table 1 reports accuracy across 6211 in-distribution test structures.Source: Table 1 · p. 7, §2.3 · p. 8
• For all 6211 test structures, median per-structure density-error reduction is 20.3-fold and exceeds unity, while 99.9% reach the regularized projected-gradient criterion.The median density error reaches within 10% of its final value after roughly 560 iterations; the tighter criterion takes a median of 1999 cycles.Source: §2.3 · p. 9

2.4 A unified functional for molecules and periodic crystals

• The same shared-parameter functional optimizes held-out isolated QM9 molecules and periodic CHNO crystals, with median absolute errors of 0.173 and 0.727 mHa/cell, respectively.Joint training preserves similar periodic accuracy but does not improve the CHNO result in the reported comparison.Source: §2.4 · pp. 9–10, Table 2 · p. 10

2.5 Transferability to larger cells

• On larger fixed periodic cells, optimized-energy MAE rises from 0.254 to 0.414 mHa/atom, remaining below 1 mHa/atom; the test covers cells beyond the training size.The larger-cell test has 734 cells with 10–24 atoms and 22–176 electrons; the local-pair regularizer reduces largest-interval wall time from 140 to 74 s per structure, with no speedup on small cells.Source: §2.5 · p. 10

2.6 Structure rankings and thermodynamic stability

• Among 538 structures across 19 compositions, STRUCTURES-P changes 4.7% of PySCF pair rankings, compared with 13.4% for CASTEP, and has 3.4 versus 8.2 meV/atom relative-enthalpy MAE.STRUCTURES-P selects the PySCF minimum in 14 of 19 compositions, versus 5 of 19 for CASTEP.Source: Figure 5 · p. 12, §2.6 · p. 13
• On the tested H–N composition slice, STRUCTURES-P matches PySCF’s minimum structure at 10 of 13 stoichiometries and reproduces hull vertices at H, H3N, and N.The reported hulls describe stability within the evaluated candidate set, not a complete phase diagram.Source: §2.6 · p. 13

3 Discussion

• 0.25 mHa per atom mean energy error and 0.4% of the electron count mean density error across 6211 held-out CHNO crystals, from molecular crystals to three-dimensional networks.The model is minimized under periodic boundary conditions from a superposition of atomic densities.Source: §3 · p. 13
• A single jointly trained functional works for isolated QM9 molecules and periodic CHNO crystals, retains sub-mHa/atom accuracy on larger cells, and preserves useful fixed-geometry rankings.Joint training does not improve periodic accuracy in the matched comparison; relative energies and the H–N convex hull retain some differences from PySCF.Source: §3 · p. 13
• All three methods give H, H3N, and N as H–N hull vertices, while STRUCTURES-P and PySCF select different structures at three stoichiometries.The figure uses method-specific elemental references and a search pressure of about 1 GPa.Source: Figure 6 · p. 14
• The tests do not demonstrate a full crystal-structure search because structures are neither generated nor relaxed with the learned functional.The hull comparison is limited to the evaluated candidate set.Source: §3 · p. 14, Figure 6 · p. 14
• Reported errors measure agreement with available PySCF references, which retain finite-orbital-basis errors for some structures; the study covers closed-shell, even-electron CHNO crystals near the AIRSS dataset pressure.Larger-cell tests extend system size, not chemistry or thermodynamic conditions.Source: §3 · p. 14
• More reliable reference data remain a next step, including systematic k-point convergence checks and more robust orbital-basis representations.Alternative all-electron references would still require density fitting into the Gaussian representation and consistent non-ground-state energy and derivative labels.Source: §3 · p. 15
• Plane-wave pseudopotential energies depend on the one-particle density matrix, so using them in this scalar-density functional requires suitable non-ground-state labels or an additional density-matrix approximation.The passage identifies this as a constraint on incorporating nonlocal pseudopotential terms.Source: §3 · p. 15
• Dense matrix operations for global natural reparametrization and regularization are the implementation’s main scaling bottlenecks.The discussion proposes better-conditioned bases or block-sparse and matrix-free schemes as possible ways to reduce their cost.Source: §3 · p. 15

4.1 Density representation in a Gaussian basis

• The even-tempered Gaussian basis uses PySCF’s aug etb from the 6-31G(2df,p) orbital basis with exponent ratio β = 2.561.The Gaussian radial functions are multiplied by spherical harmonics, with exponents in a geometric progression.Source: §4.1 · p. 15
• The total electron count is constrained during density optimization, which also keeps the unit cell charge neutral for periodic electrostatics.The count is a linear function of the density coefficients; the basis expansion itself does not enforce pointwise non-negativity.Source: §4.1 · p. 16
• For optimized densities evaluated on a level-3 Becke grid, cells with detected negative density were rechecked on level-4 and level-5 Becke grids and positive-weight uniform grids; integrated negative charge was at most 6 × 10^-5 electrons per cell.Source: §4.1 · p. 16
• The atom-centered Gaussian basis resolves near-nuclear density variations with fewer degrees of freedom than a uniform real-space grid.The even-tempered basis includes tight, high-exponent functions centered on nuclei.Source: §4.1 · p. 16

4.2 Periodic electrostatic energies

• The linear density expansion makes Hartree energy quadratic and electron–nucleus energy linear in the density coefficients.These forms use the periodic Coulomb metric and nuclear-attraction vector.Source: §4.2 · p. 16
• The periodic Coulomb metric uses a range-separated kernel, summing short-range Gaussian integrals over finite lattice images and the smooth long-range part over reciprocal-lattice vectors.The method uses analytic Gaussian Fourier transforms for the long-range contribution.Source: §4.2 · p. 17
• The zero-wave-vector component is omitted consistently from electron–electron, electron–nucleus, and nucleus–nucleus terms, and cancels in the total electrostatic energy for neutral cells.The nuclear-attraction term uses a cross metric with PySCF’s all-electron nuclear charge representation; nucleus–nucleus repulsion uses PySCF’s point-charge Ewald energy.Source: §4.2 · p. 17
• Because the density is lattice-periodic, the electrostatic quantities require no Brillouin-zone sampling.Source: §4.2 · p. 17
• The geometry-dependent Coulomb metric and external-potential vector are computed once per structure and reused during density optimization.After precomputation, fixed-geometry electrostatic energy and gradient evaluations reduce to matrix-vector and dot products.Source: §4.2 · p. 17

4.3.1 Structures and preprocessing

• 62,109 eligible structures with retained labels enter the training, validation, and test splits after filtering and preprocessing.The even-electron structures are treated as closed-shell; the in-distribution split contains cells with one to nine atoms, while the independent larger-cell test set contains 734 crystals with 10–24 atoms.Source: §4.3.1 · p. 18

4.3.2 Reference Kohn–Sham calculations

• Reference labels come from all-electron, spin-restricted periodic PBE calculations in PySCF using 6-31G(2df,p) and a 0.10 Å−1 reciprocal-space spacing.A separate ranking-sensitivity check recomputes selected fixed cells at 0.04 and 0.03 Å−1 and tests whether pairwise enthalpy differences change by less than 0.5 meV/atom.Source: §4.3.2 · p. 19

4.3.3 Perturbed densities

• Perturbations are applied to the converged reference Fock matrix, producing 40 replay densities per cell and 2,019,152 density samples overall.Unsmeared references also contribute an explicit ground-state sample; smeared references do not, and the replay densities are not self-consistent zero-temperature ground states.Source: §4.3.3 · pp. 19–20

4.3.4 Density fitting

• The Gaussian-basis density fit enforces the correct electron number and exactly reproduces the Kohn–Sham nuclear-attraction energy in the inference electrostatic convention.For 6,211 test ground states, the mean relative L2 fitting error is 0.219%; this fitting error is distinct from learned-model density error.Source: §4.3.4 · p. 21

4.3.5 Energy and gradient labels

• Energy labels are set so the orbital-free total energy matches the Kohn–Sham energy at each sampled density, while gradient labels use a potential-based target.The finite-basis stored gradient is not an exact derivative of the residual energy label; the unperturbed ground-state gradient vanishes by construction, not as an independent consistency check.Source: §4.3.5 · pp. 21–22

4.4.1 Coefficient transformation

• A global symmetric reparametrization based on the density-basis overlap matrix transforms coefficients and dual quantities while preserving the electron-number constraint.The transformation also maps gradients and Coulomb matrices into the coordinates used for training and optimization.Source: §4.4.1 · p. 22

4.4.2 Functional and training objective

• A periodic eSEN equivariant graph network predicts the learned kinetic and exchange–correlation energy, trained on both residual energies and electron-number-projected functional gradients.The gradient loss supervises the energy-landscape slope around training densities; the linear reference captures dominant average contributions before the network learns residuals.Source: §4.4.2 · pp. 22–23

4.4.3 Training settings and joint sampling

• Joint QM9–CHNO training combines 2,420,000 molecular and 2,019,152 periodic density samples without additional domain weighting.The joint model uses the same energy and projected-gradient loss weights as the CHNO-only model.Source: §4.4.3 · p. 23

4.5 Density optimization

• The total optimization gradient combines the learned-functional derivative with analytic electrostatic gradients, then projects updates onto the fixed-electron manifold.The optimization starts from a normalized superposition of atomic densities and stops at a regularized-gradient threshold or after at most 1.5×10^4 cycles.Source: §4.5 · pp. 23–25
• The local orbital-product metric identifies weakly constrained density directions, and a soft spectral regularizer suppresses small-eigenvalue directions without changing energy at fixed density.For periodic cells, the metric uses orbital pairs within a 2.5 Å minimum-image separation, including same-atom pairs; this local-pair approximation has nearly unchanged energy and density error statistics.Source: §4.5 · pp. 24–25
• The reported convergence criterion does not establish stationarity in every charge-conserving direction because the regularizer suppresses weakly constrained directions.The representative cell in Fig. 3 was selected near test-set medians for energy error, density error, improvement factor, and cycle count.Source: §4.5 · p. 25

4.6 Computational cost and size dependence

• The timing benchmark covers 144 structures across six electron-count intervals, using a single A100 GPU for optimization and one CPU thread per PySCF calculation.The hardware setup makes the PySCF timings practical context rather than a hardware-normalized comparison.Source: §4.6 · pp. 25–26
• 1.88-fold speedup in the largest interval is achieved with local pairs; the three smallest intervals show no speedup.Effective wall time per structure ranges from 8.5 to 74.5 s with local pairs and 8.4 to 140.1 s with full pairs; the largest-interval benefit comes from regularizer construction, not fewer optimization steps.Source: §4.6 · p. 26, Figure 7 · p. 26
• With fixed basis size per atom and bounded local coordination, local-pair counts grow as O(N), and dense Gram-matrix construction costs O(N^3) instead of O(N^4).This does not imply quadratic scaling of the complete calculation, which retains other dense operations; a matrix-free alternative also lacks a guaranteed size-independent iteration count.Source: §4.6 · p. 26
• Across five periodic test evaluations, replacing full pairs with local pairs changes energy MAE by at most 0.0051% and mean density error by at most 0.076%.The trained functionals and reference labels remain the same; only the optimization regularizer changes.Source: §4.6 · pp. 26–27

• AI assistants helped write and refactor project code, including analysis and figure scripts.Source: §4.7 · p. 27
• The assistants also helped draft and revise the manuscript text.Source: §4.7 · p. 27
• The authors specified and verified all methods, results, and conclusions and take full responsibility for the content.Source: §4.7 · p. 27

Supplementary Information

• 20.1%, 30.7%, and 23.0% are the energy MAE reductions for the CHNO-only model, joint-model QM9, and joint-model CHNO, respectively, with increased width.Density error reductions are 20.7%, 12.8%, and 17.9% in those matched comparisons; convergence is unchanged or improved.Source: Supplementary Information · p. 34
• The larger model is best among the two tested widths, but whether further capacity increases improve density optimization remains open.The geometry-only model's energy MAE falls by only 1.3%, from 3.361 to 3.316 mHa/cell, when widened.Source: Supplementary Information · p. 34
• Finer k-point sampling changes the lowest-enthalpy CH2 candidate but preserves the separation of higher- and lower-enthalpy groups; the largest pairwise change is 0.20 meV/atom.The model follows the broad grouping but still selects the original minimum.Source: Supplementary Information · p. 35

Reference energy changes and model errors

• For 516 screened cells, mean and median reference changes are 2.78 and 0.43 meV/atom, versus 1.69 and 0.98 meV/atom for STRUCTURES-P deviations from the original reference.Five cells still change by more than 5 meV/atom between 0.04 and 0.03 Å−1, so the comparison does not establish full k-point convergence.Source: Reference energy changes and model… · pp. 35–36
• The data support proof of concept and broad energy-trend comparisons, but reliable ranking of near-degenerate crystals requires more accurate references with systematic k-point and basis-set convergence checks.Source: Reference energy changes and model… · p. 36

Supplementary Information

• Two of 50 sampled optimized densities had small negative regions; their integrated negative charges were very small relative to their unit-cell electron counts.All 50 optimizations converged; 48 had no detected negative density on the level-3 Becke grid.Source: Reference energy changes and model… · p. 37

Reference energy changes and model errors

• Supplementary Data 1 contains source data for the main-text figures, Table 1, and supplementary figures, plus geometries and assignments for 62,109 CHNO cells and 734 larger cells.Source: Reference energy changes and model… · p. 37