Stable and Accurate Orbital-Free DFT
Powered by Machine Learning

We concentrate here on equilibrium-geometry organic molecules with a mass of up to a few hundred Dalton. While they span only a tiny corner of the entire chemical space, this family is already estimated to comprise anywhere between $10^{24}$ and $10^{60}$ members [47, 48]. Needless to say, this number grows further when taking larger biopolymers such as peptides or nucleotide chains into account.

The QM9 database [30], while restricted to stable molecules and relaxed geometries, is already quite diverse when considering only small organic molecules, see Fig. 4B. Orbital-free DFT as implemented by the STRUCTURES25 functional now affords density optimization which fully converges for each and every of the ca. 13k QM9 test molecules; with the resulting densities deviating from Kohn-Sham ground truth by $5.86\cdot 10^{-4}$ in the metric of the squared error integrated over space, normalized by the number of electrons. The biggest contribution to this deviation comes not from the machine learning model, but stems from the intrinsic error of density fitting (see Fig. S.2 in the supporting information), the process of expressing a Kohn-Sham density in the specific basis $\{\omega_{\mu}\}$ used to expand the orbital-free density in Eq. 2. Mean absolute energy errors of 0.64 mHa relative to the PBE/6-31G(2df,p) KS-DFT reference are well below the barrier of 1.6 mHa, a widely accepted definition of “chemical accuracy”.

Reaching Kohn-Sham accuracy and fully convergent density optimization on the full chemistry embodied by the QM9 dataset is already most auspicious for orbital-free DFT. Yet, for future applications, reliable extrapolation to larger systems is essential: It is here that standard Kohn-Sham DFT becomes intractable. Following prior work [26], we evaluate extrapolation accuracy on the QMugs database [49] comprising substantially larger druglike molecules. To this end, we train on smaller molecules with a mass up to around 200 Da, and evaluate on a subset of 850 test molecules, with a mass up to around 1400 Da. On these, the STRUCTURES25 functional achieves fully convergent density optimization on all molecules, shown in Fig. 4. Remarkably, even though the STRUCTURES25 functional uses message passing only across local neighborhoods up to six Bohr in radius, the mean absolute error per atom is larger, but does not grow with the size of the molecule, across the QMugs database (Fig. 4A). The three largest energy errors were all associated with the trifluoromethoxy group, a moiety which turned out to be represented only in terms of a single molecule in the training data. Larger molecules in the QMugs database highlight the improved scaling of machine-learned OF-DFT in comparison to the Kohn-Sham reference calculations, resulting in nearly an order of magnitude speedup on identical hardware (section S.11 in the supporting information). While density optimization requires many more iterations than a typical Kohn-Sham SCF calculation, this is more than offset by the relative simplicity of each single iteration which does not require matrix diagonalization.

Building on a compact representation of the electron density [27, 28] (Eq. 2), state-of-the-art machine learning architectures [45, 31], automated differentiation [50], efficient training data generation and reparametrization [26] and a strategy to create more balanced training samples [25], orbital-free DFT is now emerging as a viable framework for accurate calculations grounded in the Hohenberg-Kohn theorems.

Still, important limitations persist. The present functional, trained on organic equilibrium-geometry molecules only, does not yet generalize to high-energy geometries (a prerequisite for successful geometry optimization), or across the periodic table. Also, the poor predictions of the trifluoromethoxy group illustrate that with the present size of training data, sufficient examples of chemical moietys must be seen during training. It will be fascinating to see to what extent future models, trained on more extensive data, will “grok” chemistry and generalize across more of chemical space. Indeed, we expect future machine-learned density functionals to work well in those regimes in which Kohn-Sham DFT does. In particular, the viable area should be significantly larger than the kind of chemistries typified by the QM9 and QMugs databases.

For the following discussion, we single out the question of scaling, one principal reason to study OF-DFT in the first place. To become an alternative to linear-scaling KS-DFT [51, 52, 53, 54], OF-DFT will not just need quasi-linear scaling, but also a small prefactor. Let us recall the computational complexity bottlenecks of OF-DFT, beginning with the initialization. While many classical OF-DFT functionals are sufficiently robust to start from cheap initial density guesses, the previous state of the art, the machine-learned M-OFDFT [26] was dependent on MINAO [55, 56] which has cubic complexity in the number of atoms $N$. The robustness of STRUCTURES25 allows starting from our simple dSAD guess (see supporting text) which scales linearly. The next major step in the computational pipeline is natural reparametrization, which still has $\mathcal{O}(N^{3})$ complexity. Finding cheaper alternatives has top priority. Next, message passing on completely connected graphs comes with $\mathcal{O}(N^{2})$ complexity. Our local model operating on a radius graph reduces this to $\mathcal{O}(N)$. The cost of the Hartree term scales quadratically with the number of basis functions when the representation in Eq. 2 is evoked. Approximations such as the fast multipole method can bring this cost down to linearithmic scaling.

Outside orbital-free DFT, direct property prediction using “foundation models” is progressing with great strides[57, 58, 59, 60]. When striving to make valid predictions across a broad range of chemical space, a key question is which approach will prevail: A transductive approach, as in foundation models for property prediction which directly map from molecular geometry to observable; or an inductive ansatz, as in orbital-free DFT, where a universal functional is learned and variational optimization yields the electron density and its energy, from which observables can be derived. We believe the inductive approach will generalize better, but are eager to learn whether this intuition stands the test of time.

Overall, orbital-free DFT is now on the cusp of becoming practically useful in the molecular realm, due in part to the present contributions building on excellent prior work [61, 45, 25, 26, 31].

The model input consists of a molecular graph with atomic positions $\{\mathbf{R}_{a}\}$ and atom types $\{Z_{a}\}$ as well as the coefficients $\mathbf{p}$, representing the electron density in terms of atom-centered basis functions, cf. Eq. 2 in the main text. As is customary in local canonicalization, we compute an equivariant local coordinate system for each atom based on the relative position of adjacent non-hydrogen atoms. The basis functions $\{\omega_{\mu}\}$ are a product of a radial function and spherical harmonics. As a consequence, the density coefficients transform via Wigner-D matrices under 3D rotations. To achieve invariance w.r.t. the global orientation of the molecule, the coefficients are transformed into the respective local frame at each atom. As in M-OFDFT[26], the coefficients undergo a “natural reparametrization” into an orthonormal basis, that is specific to the molecule geometry (see section S.3 for details); and the coefficients are rescaled dimensionwise to standardize the model input and the desired gradient range of the model w.r.t. the input coefficients.

The preprocessed atom-wise coefficients are embedded as node features with a feature dimension of 768. To this end, they are first passed through a shrink gate module

with learnable parameters $\lambda_{\mathrm{in}},\lambda_{\mathrm{out}}$ and subsequently through an MLP. Pairwise distances between nodes are embedded as edge features via a Gaussian basis function (GBF) module, given by

Here, $\eta_{\mathrm{mul}},\eta_{\mathrm{bias}}$ are learnable scalars, which depend on the atom types of sending and receiving node, $k\in\{0,1,\dots,127\}$ and $\mu_{k},\sigma_{k}$ are learnable mean and standard deviation of the $k$-th Gaussian, respectively. We initialize $\eta_{\mathrm{mul}}$ to 1 and $\eta_{\mathrm{bias}}$ to 0, while $\mu_{k}$ and $\sigma_{k}$ are drawn from a uniform distribution in the interval $[0,3]$.

Further, an embedding of the atom number $Z$ and an aggregation of edge features over neighboring nodes $\operatorname{MLP}\left(\sum_{j\in\mathcal{N}(i)}e_{ij}^{k}\right)$ are added to the node features. These are then passed to a message-passing graph neural network. At the core of the architecture, we apply 4 (or 8 for QMugs) Graphormer blocks [45]. For experiments with local architectures, we propagate messages along a graph with a radial cutoff of 6 Bohr and otherwise use the fully connected graph. Before each attention block, we also apply a node-wise layer norm. Finally, we employ an energy MLP which produces an atom-specific energy contribution from the final node features. Combined with an atom-specific contribution (“atomic reference module”, see below) based on the statistics of the data, the individual energies of each molecule are aggregated by summation into the total energy prediction. For the QMugs model, we used a hierarchical energy readout (similar to prior work [74]), where, instead of a single energy MLP after the final layer, we employ an energy MLP after every second transformer block to predict atom-specific energy contributions. In the end, all energy contributions of the individual readouts are summed into the final atom-wise energies. Having readout MLPs also in intermediate layers, where the field of view is still small, enforces the prediction of more local energy contributions whereas MLPs after later layers are expected to capture increasingly non-local effects.

Proper normalization of model inputs and targets is critical for successful training of neural networks. In our pipeline, energy gradients with respect to density coefficients $\nabla_{\mathbf{p}}E_{\mathrm{target}}$ are predicted via back-propagation through our ML functional. This entails that the scales of input density coefficients and gradient labels are linked, as multiplying the former by some factor amounts to rescaling the latter by the inverse factor. Following M-OFDFT [26], we thus employ a number of enhancement modules to constrain energies, gradients and coefficients to favorable ranges.

Natural reparametrization enables the meaningful measurement of density and gradient differences. A difference $\Delta\mathbf{p}$ in density coefficients in the LCAB ansatz (Eq. 2 in the main text) results in a residual density $\Delta\rho(\mathbf{r})=\sum_{\mu}\Delta p_{\mu}\omega_{\mu}(\mathbf{r})$ with an $L_{2}$ norm of

where $\mathbf{W}$ is the density basis overlap matrix. Transforming coefficients by $\mathbf{p}\mapsto\tilde{\mathbf{p}}\coloneq\mathbf{M}^{\top}\mathbf{p}$, where $\mathbf{M}$ is a matrix square root of the overlap matrix, i.e.  $\mathbf{M}\mathbf{M}^{\top}=\mathbf{W}$, leads to an expansion for the density difference in terms of coefficient differences $\Delta\tilde{\mathbf{p}}$ only,

This implies that, after transformation, individual coefficient dimensions now equally and independently contribute to changes in the density. We can solve for $\mathbf{M}$ by considering the eigendecomposition $\mathbf{W}=\mathbf{Q}\boldsymbol{\Lambda}\mathbf{Q}^{\top}$ of the overlap matrix. The solution to $\mathbf{M}\mathbf{M}^{\top}=\mathbf{W}$ has a rotational degree of freedom as any matrix $\mathbf{M}=\mathbf{Q}\boldsymbol{\Lambda}^{\!1\!/2}\mathbf{O}$, with an arbitrary orthogonal matrix $\mathbf{O}$, results in the overlap matrix upon squaring. We here shortly discuss the implications of choosing between two very distinct choices for $\mathbf{O}$. Natural reparametrization is not only carried out on the coefficients $\mathbf{p}$ but also on all related quantities like gradients and basis functions. For the density $\rho(\mathbf{r})$ to remain invariant under reparametrization, basis functions transform via $\omega_{\mu}\mapsto\tilde{\omega}_{\mu}=\sum_{\nu}M^{-1}_{\mu\nu}\omega_{\nu}.$ As this transformation orthonormalizes the basis functions, $\int\tilde{\omega}_{\mu}(\mathbf{r})\tilde{\omega}_{\nu}(\mathbf{r})=\delta_{\mu\nu}$, we stress the similarity to the non-orthogonality problem put forward by Löwdin [75, 76]. This seminal work suggests choosing $\mathbf{O}=\mathbf{\mathbf{Q}^{\top}}$ and therefore $\mathbf{M}_{\mathrm{sym}}\coloneq\mathbf{Q}\boldsymbol{\Lambda}^{\!1\!/2}\mathbf{Q}^{\top}$, the symmetric orthogonalization – often simply called Löwdin orthogonalization [77, 78]. This flavor of reparametrization results in the orthogonalized basis set that is least distant to the original basis set, as measured by the $L_{2}$-norm between basis functions. Furthermore, any rotation of the original basis set commutes with the orthogonalization; that is, the symmetrically reparametrized basis functions transform equivariantly under orthogonal transformations, which includes spatial rotations and permutation of basis functions. The latter property in particular motivates our choice to reparametrize coefficients by $\mathbf{M}_{\mathrm{sym}}$.

Zhang et al. [26] propose $\mathbf{O}=\mathbf{1}$, i.e., $\mathbf{M}=\mathbf{Q}\boldsymbol{\Lambda}^{\!1\!/2}$ in their supplementary. This transformation, called canonical orthogonalization by Löwdin [76], results in highly delocalized orbitals. That is, each individual basis function aims to summarize the information in all untransformed orbitals subject to orthogonality. In our experiments, this reparametrization performed significantly worse than the symmetric version described above. This is because this type of reparametrization not only nullifies any previous transformation into local frames, but also breaks permutation invariance, rendering the energy prediction dependent on the order of atoms in the molecule. However, the code provided by M-OFDFT [26] employs the symmetric reparametrization as we do.

The dimension-wise rescaling and atomic reference modules are implemented as in prior work [26]. Dimension-wise rescaling linearly transforms each density coefficient independently, trading off resulting coefficient and gradient scale for each component. The atomic reference module adds a simple linear fit to the learned part of our functional. It thereby reduces the dynamic range of the predicted energy and effectively centers the gradient labels.

For both QM9 and QMugs, we train the model for 90 epochs with a batch size of 128. Over the course of all epochs, the learning rate is reduced from $7\cdot 10^{-5}$ to $0$ using a cosine annealing schedule [79]. As optimizer, we use AdamW [80] with $\beta_{1}=0.95$, $\beta_{2}=0.99$ and a weight decay factor of $10^{-10}$. We do not use dropout. The training hyperparameters are summarized in Table S.3.

We apply a loss to both the energy $E$ and its gradient $\nabla_{\mathbf{p}}E$. For the energies, a simple $L_{1}$ loss is used to compare the output to the labels. The gradient needs further attention, since it is only known in the hyperplane of normalized densities. To compare the derivative of the predicted energy with the gradient label $\mathbf{g}_{\mathrm{label}}$, we follow M-OFDFT [26] and apply an $L_{1}$ loss on the projected difference,

where $\mathbf{w}$ is the vector of integrals over the density basis with components $w_{\mu}=\int\mathrm{d}^{3}r\,\omega_{\mu}(\mathbf{r})$ and the expression $\left(\mathbf{I}-\frac{\mathbf{w}\mathbf{w}^{T}}{\mathbf{w}^{T}\mathbf{w}}\right)$ corresponds to a projection onto the hyperplane orthogonal to $w$.

For direct comparison with prior work [26], we group the QMugs dataset by the number of heavy atoms into bins of width 5. The first bin contains molecules with $10-15$ heavy atoms and molecules with fewer heavy atoms are discarded from the set. When training on the QMugs dataset, we first combine the QM9 dataset with the first QMugs bin (10-15 heavy atoms) for the initial training set. Following this initial training on the combined dataset, we perform an additional fine-tuning step. This fine-tuning is conducted exclusively on the QMugs subset containing molecules with 10-15 heavy atoms (the first bin). The fine-tuning uses the same hyperparameters as the initial training, except for the learning rate and epoch count. We start with a learning rate of $1\cdot 10^{-5}$, which is also reduced to $0$ following a cosine annealing schedule over 30 epochs. All other training parameters, including the optimizer, loss function, and batch size, remain unchanged. Training on the QM9 dataset alone does not involve this fine-tuning step.

A multitude of established methods for generating initial guesses for the electron density in KS-DFT exist, such as the MINAO initialization [55, 56]. However, while it is cheap compared to the Kohn-Sham iterations because of the minimal basis that it uses, it still scales cubically with system size, and additionally requires density-fitting to transform the guess to the density basis (Eq. 2 in the main text).

A much simpler and linear scaling option is a superposition of atomic densities (SAD). We have datasets with ground state density coefficients at hand, as these are required for training. This allows us to determine average atomic densities with a data-driven approach: We take all instances of each atom type (i.e. chemical element) in the dataset, and take the average of the corresponding coefficients over all these instances. Averages for coefficients corresponding to basis functions with $l>0$ are set to zero. For a given molecule $\mathcal{M}$, concatenating these averages for all atom types in the molecule then yields $\bar{\mathbf{p}}$, our data-driven SAD (dSAD). However, this superposition of atomic densities is not necessarily normalized to the correct electron number. Since the number of electrons stays invariant during density optimization (see section S.6), we need to normalize to generate a valid guess. One approach is to uniformly scale the coefficients linearly to the correct electron number $N_{\mathrm{e}}$, leading to

with the basis integrals $\mathbf{w}$, see section S.4. While simple, this has a major shortcoming: The largest part of the electron density lies close to the cores, and, for elements other than hydrogen, this core density varies only very little between different instances of the same atom type in neutral molecules. Thus, the SAD guess describes the core very precisely. The coefficients of the inner $l=0$ basis functions largely describe this core density and should hence be varied very little in the normalization. However, scaling all coefficients by the same factor to achieve normalization does not respect this. For example, the core density of atomic species with high electronegativity (whose corresponding coefficients, on average, describe a higher number of electrons than their atomic number indicates), would be scaled down and hence underestimated.

This is why we propose a heteroscedastic normalization procedure, which adapts to the variance of the coefficients over the dataset: Coefficients with high variance are scaled more than those with low variance, as they are more likely to be far from the mean. As the mean $\bar{\mathbf{p}}$ and variance $\sigma_{\mu}$ of each coefficient are known, we can formulate this as an weighed least squares optimization problem with a linear constraint, where the weights of the squared deviations from $\bar{\mathbf{p}}$ are given by the inverse squares of the variances (see also Fig. S.1):

with $\Delta N_{\mathrm{e}}=N_{\mathrm{e}}-\mathbf{w}^{\top}\bar{\mathbf{p}}$, the difference between the desired electron number and the one corresponding to the mean coefficients. Introducing a Lagrange-multiplier $\lambda$, we get:

and can solve for $d$ and $\lambda$

to find

This result matches the intuition established above: The correction to each component of the average coefficients $\bar{\mathbf{p}}$ is proportional both to the variance of it over the dataset, and the weight of its corresponding basis function. An illustration of this method compared to simply scaling the guess is shown in Fig. S.1.

One could either apply this normalization per molecule, or per chemical element (in the latter case, $N_{\mathrm{e}}$ denotes the number of electrons corresponding to the atom type). We choose the latter for simplicity, normalizing the electron number per element. The guess is computed before applying natural reparametrization. Comparing our dSAD guess to the MINAO guess on 1000 QM9 molecules in Table S.4, we find that it produces guesses which are slightly closer to the ground state.

After initialization of the electron density using our dSAD guess (section S.5), the learned energy functional enables its iterative optimization in order to find the ground state density and the corresponding energy. To conserve the number of electrons, we must not diverge from the hyperplane of normalized densities $\{\mathbf{p}\colon\mathbf{w}^{\top}\mathbf{p}=N_{e}\}$. We thus project the step $\mathbf{u}$ of the optimizer onto the hyperplane such that the density coefficients are updated according to

We use gradient descent with momentum as our optimizer, with a learning rate of $0.003$ and a momentum of $0.9$ for QM9. For the QMugs dataset, we reduce the learning rate to $0.0015$. These parameters were tuned such that the density optimization shows a fast and robust convergence. The number of iterations is limited to 5000.

Note that while KS-DFT typically requires significantly fewer steps to converge, it relies on much more complex SCF iterations. While OF-DFT requires more steps, each is computationally inexpensive. This results in a significant net gain in efficiency, particularly for larger systems, see section S.11 for details.

The model is never trained on molecules from the test set, and density optimization is performed on the test set only. As Fig. 2B in the main text illustrates, density optimization on the STRUCTURES25 functional converges to gradient norms of $10^{-13}$ Ha for smaller molecules. This level of convergence is more precise than our labels are, due to imperfect density fitting and a finite convergence tolerance in the Kohn Sham calculations of our ground truth, see section S.1.2. We thus stop density optimization when the gradient norm for an entire molecule falls below $10^{-4}$ Ha.

Regarding the definition of chemical accuracy, an error of 1 kcal mol^-1 is a widely used threshold for acceptable energy errors. As there is no widely used definition for chemical accuracy of an electron density, we propose to compare the densities produced by different XC functionals using the same (def2-TZVP) basis set [82]. At the local density approximation (LDA) rung we use BN05[83]. Generalized gradient approximations (GGA) are represented by PBE and B97[84]. As Meta-GGA, incorporating the kinetic energy density, we use R2SCAN[36]. The highest accuracy levels are achieved by the hybrid-GGAs PBE0[85], B3LYP[86] or a meta-hybrid-GGA with PW86 exchange and B95 correlation [87, 88]. The ground state densities of 1000 randomly sampled molecules from QM9 were computed using our default Kohn Sham settings for each of the above functionals. We then measured the density difference between the XC functional PBE[66], which was used for data generation, and this assortment of functionals, with results shown in Table S.5. The density difference between PBE and methods at least at the GGA level is in the range of $7.2\cdot 10^{-4}$ to $1.3\cdot 10^{-3}$ electrons in the L2 norm, while the density difference to the less accurate LDA BN05 is around $4.1\cdot 10^{-3}$ electrons. Given the above deviations from PBE densities, and given that the PBE functional has been used in more than 200 000 publications to date, we here define “chemically accurate densities” at the level of PBE as all those that come within $7.2\cdot 10^{-4}$ electrons in the L2 norm on QM9-sized molecules.

Fig. 4A in the main text, which shows the extrapolation accuracy from small to large organic molecules, reveals three outliers. It turned out that all three contained a trifluoromethoxy group, and that this chemical group was present only in the three conformers of a single molecule in the training set. We thus hypothesized that this group was responsible for the poor predictions. To investigate this hypothesis, we substituted the fluorine atoms of the trifluoromethoxy group by hydrogen, yielding their methoxy derivatives (see Fig. S.3). We re-optimized the geometry at the same level of theory as the original samples (GfN2-xtb using energy and gradient convergence criteria of $5\times 10^{-6}$ Ha and $10^{-3}$ Ha$\,\alpha^{-1}$) [49]. Subsequent OF-DFT calculations employing the STRUCTURES25 functional (trained on QMugs) showed that the energy error per atom dropped from $11.69$ mHa to $0.4$ mHa, from $10.2$ mHa to $0.2$ mHa, and from $7.42$ mHa to $0.47$ mHa respectively, well within the range of the other QMugs test samples. This indicates it was indeed the trifluoromethoxy group which caused the outliers.

To further analyze the phenomenon, we evaluated the accuracy of STRUCTURES25 on smaller molecules containing a trifluoromethoxy group. These exhibited large errors as well, see Table S.6. Finally, the density error for a larger molecule containing the group is visualized in Figure S.4, demonstrating that the error is spatially localized around the group in question. Taken together, these observations support the hypothesis that the lack of sufficiently many training molecules including the trifluoromethoxy group is the root cause of this failure mode. Analyses of this type can help identify and close gaps in the training data.

The representation of the density given in Eq. 2 in the main text in principle allows for regions of negative densities to occur. When using a trained model for density optimization this could lead to problems. Since negative densities are not physical and therefore no training data with negative densities exists, it is unclear if the model can make meaningful predictions for such non-physical densities.

In practice, we report that negative densities are no failure case for our models trained on the $E_{TXC}$ target. When starting from a reasonable initial guess such as dSAD (section S.5), the negative regions of the converged densities are insignificant. This is illustrated for our QM9 model and that of M-OFDFT [26] in Fig. S.5. If negative densities were to become a problem in the future, e.g. for larger molecules or when considering more elements, one could penalize negative densities directly in the density optimization. A penalization term of the form

can be added to the total energy, where $\gamma>0$ is a hyperparameter. Using this penalization term does not significantly change our results, which is why we do not use it for the reported experiments. Such a penalization term has the significant downside of requiring a grid for evaluation, which we otherwise do not need for the $E_{\mathrm{TXC}}$ target.

Zhang et al. [26] require a precise machine-learned guess, “ProjMINAO”, to produce their best results. Their results deteriorate by an order of magnitude when initializing with Hückel [81] densities and become worse still when initializing with a MINAO guess [26, 55, 56]. To gauge the robustness of the STRUCTURES25 functional, we compare density optimizations initialized with our data-driven superposition of atomic densities (dSAD, section S.5) as well as MINAO and Hückel densities. Unlike prior work [26] we refrain from training a model to improve on the initially guessed density.

For the following comparisons, we randomly sample 1000 molecules from our QM9 test set. As a point of reference for the various initial guesses, we report their density error with respect to the ground state in table S.4. Predicted ground state density errors, $\rho^{\mathrm{gs}}-\rho^{\ast}_{\mathrm{guess}}$, from dSAD and MINAO are nearly identical while Hückel fares approximately an order of magnitude worse.

On the same set of molecules, we compare density optimization results starting from the different guesses in Table S.7 and Fig. S.6. Using the identical hyperparameter configuration, both dSAD and MINAO guesses lead to all molecules converging, i.e. achieving gradient norms below $10^{-4}~{\mathrm{Ha}}$, within 5000 optimization steps. Hückel initialization is worse, with only 79% of molecules converging with default settings. Increasing the maximum number of iterations to 20k increases the convergence ratio to 85%. Tuning of the momentum to a value of 0.77 pushes the Hückel convergence ratio to 96.7% within 20k iterations. The ability to converge to good solutions for all molecules from both dSAD and MINAO guesses using the same model is a testament to the generality of the STRUCTURES25 functional, as these two initializations differ considerably (cf. Table S.4).

Starting from different guesses, density optimizations with STRUCTURES25 converges to very similar solutions. Specifically, the predictions on 953 of the 1000 tested molecules differ at least an order of magnitude less to predictions from other initializations than they differ to the ground state label. The remainder exhibit density differences that are of the same order as the distance from the ground truth density.

The design of both the neural network architecture and the training procedure involves numerous choices. This section details our ablation studies, performed to systematically evaluate the impact of key parameters and justify our final model configuration. To minimize the influence of stochasticity, we report the best performance of three independent training runs (seeds) for each configuration, evaluating models based on the energy error unless otherwise specified.

Figure S.8 compares the runtimes of OF-DFT as embodied by STRUCTURES25 and KS-DFT as implemented in GPU4PySCF [89, 90], which can exploit the parallelism of the identical hardware. All computations were run on one Nvidia A100 GPU and 8 cores of an AMD EPYC 7452 processor. For the largest molecule of our QMUGS test set, STRUCTURES25 reaches a speedup of almost an order of magnitude. Importantly for future calculations on larger (e.g., biomolecular) systems, it shows favorable scaling with the size of the system.