Scaling Density Functional Theory with Gaussian Splatting

Andrés Guzmán-Cordero; Cindy Zhang; Majdi Hassan; Marta Skreta; Kirill Neklyudov; Matija Medvidović

physics.chem-phSep 25, 2026

Scaling Density Functional Theory with Gaussian Splatting

Authors: Andrés Guzmán-Cordero, Cindy Zhang, Majdi Hassan, Marta Skreta, Kirill Neklyudov, Matija Medvidović

Organizations: Mila - Quebec AI Institute · Université de Montréal · Princeton University · Institut Courtois ·

Abstract

Density functional theory (DFT) strikes a practical balance between accuracy and computational cost in many problems of computational chemistry and materials science. However, many DFT calculations are limited by fixed atom-centered basis sets, which dictate how accuracy and cost scale with system size. We propose Gaussian Splatting for Density Functional Theory (GS-DFT), which represents molecular orbitals as a cloud of Gaussians whose positions, shapes, and mixing coefficients are optimized jointly by gradient descent to minimize the energy without training data. Conceptually, GS-DFT is 3D Gaussian splatting with the renderer replaced by quantum mechanics. We introduce two key solver components: adaptive density fitting with screening for efficient evaluation of two-electron integrals, and a regularized differentiable orthogonalization of the molecular orbitals. Empirically, the optimized basis reaches the accuracy of the largest conventional basis sets with a fraction of the parameters, converging systematically in energy, density, and nuclear forces. At equal parameter count, it captures the stretched-bond and anion physics that fixed bases only recover with specialized basis augmentation. The resulting solver exhibits quadratic peak memory scaling in the cloud size, allowing us to simulate systems of up to 2,742 atoms (10,406 electrons) without any modifications at triple-zeta scale using a single four-GPU node.

Figures & tables

Figure 1: Gaussian Splatting for Density Functional Theory (GS-DFT). Molecular orbitals are expanded in a cloud of

M

Gaussians ( splats ) whose means and covariances

\Theta

are optimized together with the coefficients

{\bm{C}}

by variational direct energy minimization. The rightmost panels report both scalings for alanine dipeptide (

N{=}78

electrons). The error against the extrapolated references falls faster per free parameter than for the traditional bases with

B

functions, and peak memory grows as

\mathcal{O}(M^{2})

instead of the

\mathcal{O}(B^{3})

of most traditional codes.

Figure 2: Components of Gaussian Splatting for Density Functional Theory (GS-DFT). (a) Parametrization of Equation 8 , and the [Floating Elliptical Gaussian Orbitals, FSGO] ( Vescelius & Frost, 1974 ) optimization trajectories against converged GS-DFT of the same function count. (b) Adaptive density fitting of Equation 14 with screened, chunked evaluation, plotted as peak memory against free parameters. (c) The sharded scaling , with respect to the best energy, against wall-clock time, on chain of 45 alanine residues (

M=4299

Figure 3: Observables against free parameters. Columns are the energy residual per atom, the total variation of the density per electron, and the maximum deviation of the nuclear forces. Top row: measured against the CBS limit. Bottom row: measured against the ISC limit.

Table 4

Figure 4: (a) Fluoride ion (

\text{F}{\vphantom{\text{X}}}^{\text{\hskip 0.93333pt--\hskip 0.93333pt}}

) basis-set convergence.

\mathbf{\Delta}_{\text{CBS}}=E-E_{\text{CBS}}

is the residual from the CBS limit in milli-Hartree (mHa). (b) LiF dissociation at matched count. The margin over the unaugmented basis grows as the bond stretches, while the margin over the augmented basis stays flat.

Figure 5: Accuracy scaling law for GS-DFT . GTO rungs against the extrapolated complete basis set (CBS) energy, and the splat clouds against the infinite splat cloud (ISC) energy. Limits are obtained by fitting

\Delta=E_{P}-E_{\infty}\propto P^{-\alpha}

against four runs with

P

parameters.

System	Atoms	M	GS-DFT	GTO	$\mathbf{E}_{\text{GS-DFT}}$	$\mathbf{E}_{\text{HF}}$	GS-DFT	GS-DFT
			Peak Mem.				Fits GPU node	Simple forces
Trp-cage	304	6,720	$25.2$	$640^{\ast}$	$-7{,}478.3$	$-7{,}439.0$
Defensin	465	10,390	$39.6$	$1{,}060^{\ast}$	$-13{,}685.4$	$-13{,}624.0$
Insulin	947	21,026	$92.8$	$2{,}400^{\ast}$	$-25{,}523.2$	$-25{,}404.2$
$\beta$ -spectrin CH	1,755	38,794	$212.0$	$4{,}850^{\ast}$	$-44{,}255.4$	$-44{,}045.7$
3IFU (chain A)	2,742	60,308	$479.6$	$8{,}070^{\ast}$	$-69{,}823.2$	$-69{,}631.4$

Table 3: FMO systems. Energy in Ha, memory in GB. One

4\times

H200 GPU node has 564 GB in total.

Appendix figures & tables16 assets

Supplementary material from the paper’s appendix.

Appendix

$N_{\mathrm{occ}}$	5	13	39
	$\text{H}{\vphantom{\text{X}}}_{\smash[t]{\text{2}}}\text{O}$	$\text{C}{\vphantom{\text{X}}}_{\smash[t]{\text{2}}}\text{H}{\vphantom{\text{X}}}_{\smash[t]{\text{5}}}\text{OH}$	Ala dipeptide
Gaussian basis	1.40	0.54	0.18
Splat cloud	1.80	0.69	0.23

Appendix

Table 4: Basis parameters per coefficient parameter,

\eta/N_{\mathrm{occ}}

, with

\eta=7

for the Gaussian basis and

\eta=9

for splats.

Obstacle	Earlier cost	Component of GS-DFT	Cost or bound
(i) Integrals	$\mathcal{O}(M^{4})$ per step	Adaptive density fitting	Algorithm 1 , Proposition A.5
(ii) Gradients	$D+1$ evaluations	Reverse-mode differentiation	$\mathcal{O}(1)$ evaluations
(iii) Scale	$\mathcal{O}(D^{3})$ per step	Adam	$\mathcal{O}(D)$ per step
(iv) Dependence	unbounded gradient	Regularized eigensolve	Proposition A.6

Appendix

Table 5: Obstacles that confined floating-Gaussian models to

M=N_{\mathrm{occ}}

, and the component of GS-DFT that removes each.

	3DGS	GS-DFT
Primitive	anisotropic Gaussian	anisotropic Gaussian
Shape chart	covariance, quaternion and log-scales	precision, quaternion and log-eigenvalues
Weights	opacity and color, nonnegative	orbital coefficients ${\bm{C}}$ , signed
Channels	color channels	occupied orbitals
Represented signal	radiance along rays	orbitals $\psi_{k}$ and density $\rho$
Forward map	projection and alpha compositing	closed-form integrals and quadrature

Appendix

Table 6: Components of 3DGS and their counterparts in GS-DFT.

Operation	Algorithm	Compute	Memory
Overlap and kinetic, ${\bm{S}}{\bm{C}}$ and ${\bm{T}}{\bm{C}}$	Dense row blocks	$\mathcal{O}(M^{2}N_{\mathrm{occ}})$	$\mathcal{O}(MN_{\mathrm{occ}})$
Occupied Gram and retraction, Equation 12	Eigensolve of ${\bm{G}}={\bm{C}}^{\top}{\bm{S}}{\bm{C}}$	$\mathcal{O}(MN_{\mathrm{occ}}^{2}+N_{\mathrm{occ}}^{3})$	$\mathcal{O}(N_{\mathrm{occ}}^{2})$
Nuclear attraction	Retained pairs	$\mathcal{O}(MkN_{\mathrm{at}})$	$\mathcal{O}(Mk)$
Fitted Hartree, Equation 14	Fit vector over retained pairs	$\mathcal{O}(MkN_{\mathrm{aux}}+N_{\mathrm{aux}}^{2}+N_{\mathrm{aux}}^{3}/K)$	$\mathcal{O}(N_{\mathrm{aux}}^{2})$
Fitted exchange, hybrids only	One fit vector per orbital pair	$\mathcal{O}(MkN_{\mathrm{aux}}N_{\mathrm{occ}}^{2}+N_{\mathrm{aux}}^{2}N_{\mathrm{occ}}^{2}+N_{\mathrm{aux}}^{3})$	$\mathcal{O}(N_{\mathrm{aux}}^{2})$
Exchange–correlation	Chunked grid	$\mathcal{O}(N_{g}MN_{\mathrm{occ}})$	$\mathcal{O}(N_{g,\mathrm{c}}M)$

Appendix

Table 7: Leading-order cost of one energy-and-gradient evaluation. The factorization of the auxiliary metric is paid once per refresh and appears amortized over the

K

steps it serves.

Setting	Value	Reference
Optimizer	Adam, clip at global norm $1$ , no weight decay	above
Learning rate	$10^{-2}$ peak, $10^{-4}$ floor, cosine	above
Screening	$\tau=10^{-7}$ on the normalized overlap, see below	Section 4
Pair-list padding	$2.0$	Section C.4
Fit regularizer	$\lambda=10^{-8}$	Proposition A.4
Eigensolve	$\varepsilon=\delta=10^{-4}$	Equation 12

Appendix

Table 8: Settings shared by every GS-DFT run. Settings that vary by experiment are in Table 9 .

Result	Cloud size	Init	Steps	$K$	Grid	Seeds
Observables, Figure 3	ladder of $M$	(i)	$12{,}000$	$10$	$3$	$3$
Accuracy per parameter, Figure 5	ladder of $M$	(i)	$12{,}000$	$50$	$3$	$3$
Functional coverage, Table 2	$M=n_{\mathrm{ao}}$ (cc-pVTZ)	(i)	$48{,}000$	$50$	$3$	$3$
Ionization potentials, Table 2	$M=n_{\mathrm{ao}}$ (def2-QZVPP)	(i)	$12{,}000$	$50$	$3$	$3$
Anion convergence, Figure 4 (a)	ladder of $M$	(i)	$12{,}000$	$50$	$3$	$3$
LiF dissociation, Figure 4 (b)	$M=n_{\mathrm{ao}}$ (cc-pVTZ)	(i)	$4{,}000$	$50$	$3$	$3$

Appendix

Table 9: Settings that vary by experiment.

n_{\mathrm{ao}}

is the function count of the named Gaussian basis.

Code	GPU	Largest system
GS-DFT (ours)		$2{,}742$
D4FT ( Li et al., 2023 )		$9$
MESS ( Helal & Fitzgibbon, 2024 )		–
DQC ( Kasim et al., 2022 )	a	$9$
PySCF ( Sun et al., 2020 )		$153$ b

Appendix

Table 10: Capabilities of differentiable DFT codes, with PySCF as the non-differentiable reference. Relevant bases means agreement with PySCF to within

0.1

mHa for water in cc-pVDZ and cc-pVTZ. The last column is the largest system, in atoms, completed at cc-pVTZ on one node before OOM.

Code	$E$ (Ha)	$\Delta$ (mHa)	Iterations	Time (s)
PySCF	$-76.33344$	–	$7$	$1.9$
D4FT (float64)	$-76.33340$	$+0.038$	$637$	$128$
D4FT (float32, default)	NaN	–	$4{,}000$	$206$
MESS a	$-76.92339$	$-588.7$	$110$	$9.7$
DQC (CPU)	$-76.33345$	$-0.010$	–	$15.4$
GS-DFT ( $M=24$ )	$-76.26686$	$+66.6$	$12{,}000$	$154$

Appendix

Table 11: Water, cc-pVDZ, PBE, level-3 grid.

\Delta

is the energy minus the PySCF energy. Times are wall-clock seconds including compilation. Iterations are SCF cycles for PySCF and optimization steps for the others. a Cartesian basis, compared against Cartesian PySCF (

-76.33470

Ha).

System	Atoms	$n_{\mathrm{ao}}$	Mem	Time/step	Wall (s)	Wall (s)	Time/cycle
			GS-DFT		D4FT	DQC	PySCF
Water	$3$	$58$	$0.18$	$0.10$	$474$	$19.9$	$0.7$
Ethanol	$9$	$174$	$1.37$	$0.11$	$4{,}365$	$29.3$	$1.2$
$(\textrm{Ala})_{2}$	$22$	$468$	$1.49$	$0.15$	–	OOM	$4.6$
$(\textrm{Ala})_{5}$	$53$	$1{,}158$	$1.76$	$0.34$	–	–	$22.9$
$(\textrm{Ala})_{15}$	$153$	$3{,}358$	$6.43$	$2.06$	–	–	$238$

Appendix

Table 12: The cc-pVTZ ladder on one node ( Section D.2 ). Memory is the peak GPU memory in GB. GS-DFT reports time per optimization step and PySCF time per SCF cycle, both in seconds; D4FT and DQC report the wall-clock time of the whole run, since neither exposes a per-iteration cost. GS-DFT runs on four H200 GPUs at every rung. GS-DFT uses

M=n_{\mathrm{ao}}

splats; its time per step is the median over

100

steps after compilation. OOM marks an out-of-memory failure, and a dash marks a system not attempted after an earlier failure or, for PySCF, one that did not finish within the day of wall-clock time it was given.

System	$M$	Regularized	Plain	$\Delta$ (mHa)
$\text{CH}{\vphantom{\text{X}}}_{\smash[t]{\text{4}}}$	$86$	$-40.51841$	$-40.47054$	$+47.9$
$\text{C}{\vphantom{\text{X}}}_{\smash[t]{\text{6}}}\text{H}{\vphantom{\text{X}}}_{\smash[t]{\text{6}}}$	$264$	$-232.09545$	$-232.05306$	$+42.4$
$\text{CO}{\vphantom{\text{X}}}_{\smash[t]{\text{2}}}$	$90$	$-188.48252$	$-188.48038$	$+2.1$
$\text{N}{\vphantom{\text{X}}}_{\smash[t]{\text{2}}}$	$60$	$-109.45544$	$-109.45653$	$-1.1$

Appendix

Table 13: Regularized eigensolve against a plain eigendecomposition, at matched settings.

\Delta

is the energy without the regularization minus the energy with it.

Optimizer	Family	$\Delta$ (mHa)	Optimizer	Family	$\Delta$ (mHa)
Adamax	Adam	$-13.7$	Novograd	Adam	$+228.6$
AdaBelief	Adam	$-9.5$	AdaGrad	second moment	$+302.2$
RAdam	Adam	$-9.4$	signSGD	sign	$+321.3$
LAMB	trust ratio	$-8.3$	SGD	non-adaptive	$+449.3$
AdamW	weight decay	$-4.2$	LARS	trust ratio	$+847.2$
Adam	Adam	$0.0$	Fromage	trust ratio	$+998.8$

Appendix

Table 14: Optimizers on

(\textrm{Ala})_{2}

with

M=468

\Delta

is the final energy minus that of Adam.

		$\text{H}{\vphantom{\text{X}}}_{\smash[t]{\text{2}}}\text{O}$ , $M=48$		$\text{C}{\vphantom{\text{X}}}_{\smash[t]{\text{2}}}\text{H}{\vphantom{\text{X}}}_{\smash[t]{\text{5}}}\text{OH}$ , $M=144$
	Adds	Energy	Spread	Energy	Spread
S0	random, atom-centered	$-76.30164$	$11.7$	$-154.87106$	$4.6$
S1	element exponent ranges	$-76.32034$	$1.3$	$-154.89223$	$1.9$
S2	bond-centered splats	$-76.32044$	$1.3$	$-154.89412$	$6.0$
S3	bond-aligned anisotropy	$-76.32078$	$2.3$	$-154.89656$	$1.7$
S4	charge-proportional allocation	$-76.32853$	$2.1$	$-154.90574$	$0.1$

Appendix

Table 15: Initialization ingredients, added cumulatively. Median energy over seeds after

4{,}000

steps, with the seed spread in mHa.

$K$	Energy (Ha)	Time (s)	ms / step
$1$	$-495.60818$	$9{,}629$	$802$
$5$	$-495.59166$	$6{,}427$	$536$
$10$	$-495.63820$	$5{,}878$	$490$
$50$	$-495.57240$	$5{,}580$	$465$
$200$	$-495.55246$	$5{,}462$	$455$

Appendix

Table 16: Refresh period

K

on alanine dipeptide,

M=468

12{,}000

steps,

\tau=10^{-7}

. Time is the total wall-clock time of the run on one A100.

$\tau$	Energy (Ha)	Pairs	Retained	Time (s)
$0$	$-495.61730$	$203{,}560$	$100\%$	$11{,}025$
$10^{-8}$	$-495.59341$	$102{,}450$	$50\%$	$6{,}310$
$10^{-7}$	$-495.63820$	$94{,}658$	$47\%$	$5{,}878$
$10^{-6}$	$-495.65585$	$86{,}488$	$42\%$	$5{,}503$
$10^{-5}$	$-495.61626$	$76{,}734$	$38\%$	$5{,}019$

Appendix

Table 17: Screening threshold

\tau

on alanine dipeptide,

M=468

12{,}000

steps,

K=10

. Pairs is the number retained at the last refresh, out of

203{,}560

unscreened.

	$\zeta$	$M$	cc-pV $\zeta$ Z	aug-cc-pV $\zeta$ Z	GS-DFT	Spread
$\text{F}{\vphantom{\text{X}}}^{\text{\hskip 0.90417pt--\hskip 0.90417pt}}$	D	$14$	$-99.66562$	$-99.77465$	$-99.08264$	$304.1$
	T	$30$	$-99.74783$	$-99.79967$	$\mathbf{-99.80099}$	$4.7$
	Q	$55$	$-99.77711$	$-99.80743$	$\mathbf{-99.81283}$	$0.3$
	5	$91$	$-99.79765$	$-99.81012$	$\mathbf{-99.81330}$	$0.1$
$\text{OH}{\vphantom{\text{X}}}^{\text{\hskip 0.90417pt--\hskip 0.90417pt}}$	D	$19$	$-75.63289$	$-75.73520$	$-75.54983$	$83.2$
	T	$44$	$-75.70329$	$-75.75426$	$\mathbf{-75.75548}$	$1.8$

Appendix

Table 18: Bare anions at matched function count. Energies in Hartree, restricted PBE.

M

is the function count of the plain cc-pV

\zeta

Z basis. Bold marks a splat energy below the augmented basis of the same cardinal number.

Figure 6: Dissociation energy ladders

E(R)

for LiH (left) and LiF (right), restricted PBE on the ionic

\text{Li}{\vphantom{\text{X}}}^{\text{+}}\text{H}{\vphantom{\text{X}}}^{\text{\hskip 0.90417pt--\hskip 0.90417pt}}

\text{Li}{\vphantom{\text{X}}}^{\text{+}}\text{F}{\vphantom{\text{X}}}^{\text{\hskip 0.90417pt--\hskip 0.90417pt}}

diabat. Plain cc-pVDZ/TZ/QZ (blue), augmented aug-cc-pVXZ (green), and splats sized to the plain cc-pVXZ counts (red).

Density Functional Theory

See all 25 papers

Surrogate Functionals for Machine-Learned Orbital-Free Density Functional Theory
Apr 22, 2026Roman Remme, Fred A. HamprechtSurrogate Modeling Quantum Chemistry
Interdisciplinary Center for Scientific Computing (IWR), Heidelberg University
Learning the Kohn-Sham map with neural operators for quasi-linear scaling density functional theory
Aug 24, 2026Danish Khan, Maurice D. Hanisch, Nikolai Argatoff +3Density Functional Theory Computational Quantum Chemistry
Department of Computing and Mathematical Sciences, California Institute of Technology · Department of Mathematics, ETH Zürich · Division of Chemistry and Chemical Engineering, California Institute of Technology +1
Equivariant Flow Matching for Electron Density Prediction
Oct 2, 2026Chenxing Liang, Chengdong Wang, Yuchao Lin +2Computational Quantum Chemistry Quantum Chemistry
Department of Computer Science and Engineering, Texas A&M University · Department of Materials Science and Engineering, Texas A&M University · Department of Electrical and Computer Engineering, Texas A&M University +2

Apr 22, 2026J/K move · Enter open · S save

Surrogate Functionals for Machine-Learned Orbital-Free Density Functional Theory

Roman Remme, Fred A. Hamprecht

Interdisciplinary Center for Scientific Computing (IWR), Heidelberg University

Surrogate Modeling Quantum Chemistry Adaptive Sampling Computational Materials Science Density Functional Theory Computational Quantum Chemistry

Open paper AskPDF arXiv

We introduce surrogate functionals: machine-learned energy functionals for orbital-free density functional theory (OF-DFT) which are defined not by universal fidelity to a physical reference, but merely by the requirement that density optimization with a fixed procedure yields the true ground-state density. Helpfully, training surrogate functionals requires only ground-state densities, no energies or gradients away from the ground state. We here propose a gradient-descent-improvement loss that guarantees exponential convergence of the density to the ground state, and combine it with an adaptive sampling scheme that concentrates learning around the optimization trajectories actually visited during inference. On the QM9 and QMugs benchmarks, surrogate functionals achieve density errors competitive with or improving upon the state of the art for fully supervised machine-learned OF-DFT, while eliminating the need for the

O(N^3)

orthononormalization step required by prior work, yielding improved runtime scaling for larger systems.

Abstract

Density Functional Theory

Explore similar work

Surrogate Functionals for Machine-Learned Orbital-Free Density Functional Theory

Learning the Kohn-Sham map with neural operators for quasi-linear scaling density functional theory

Equivariant Flow Matching for Electron Density Prediction