Agentic Discovery of Exchange-Correlation Density Functionals

# Agentic Discovery of Exchange–Correlation Density Functionals
Source: [https://arxiv.org/html/2605.05460](https://arxiv.org/html/2605.05460)
1\]ByteDance Seed 2\]Princeton University\\contribution\[\*\]Work done at ByteDance Seed\\contribution\[†\]Corresponding authors

###### Abstract

The development of accurate exchange–correlation \(XC\) functionals remains a longstanding challenge in density functional theory \(DFT\)\. The vast majority of XC functionals have been hand designed by human researchers combining physical insight, exact constraints, and empirical fitting\. Recent advances in large language models enable a systematic, automated alternative to this human\-driven design loop\. This report presents an agentic search system in which an LLM proposes structured functional\-form changes guided by evolutionary history\. The system attempts to improve functional performance through an iterative plan–execute–summarize loop, where improvements are measurable by optimizing functional parameters against a standard thermochemistry dataset, then evaluating performance on a held\-out subset\. The strongest discovered functional,SAFS26\-a\(Seed Agentic Functional Search 2026\), improves upon the gold\-standardω\\omegaB97M\-V baseline by∼9%\{\\sim\}9\\%\. These results also surface a cautionary lesson for AI\-assisted science: models powerful enough to discover genuine improvements are equally capable of exploiting unphysical shortcuts to game the benchmark; domain expertise translated into explicitly enforced constraints remains essential to keeping results scientifically grounded\.

\\correspondence

Yang Sun at , Yixiao Chen at

## 1Introduction

The accuracy of Kohn\-Sham DFT is determined primarily by the quality of the exchange\-correlation functionalExc​\[ρ\]E\_\{\\mathrm\{xc\}\}\[\\rho\]\. Over time, this challenge has produced a dense hierarchy of approximations, from local\-density approximations, generalized gradient approximations, and meta\-GGAs to hybrids, and double hybrids, each trading off accuracy, cost, and generality\. More than 200 XC functionals have been proposed to date, yet there is still no systematic recipe for generating functionals that are reliably accurate across diverse chemical systems\[mardirossian2017thirty\]\.

This lack of a constructive design procedure makes XC functional development a natural target for automated search\. Recent LLM\-driven evolutionary systems have been most successful on algorithmic and engineering tasks, where progress is easily measured and successive iterations can realize concrete improvements\[romeraparedes2024funsearch,novikov2025alphaevolve,yang2023opro,wan2025loongflowdirectedevolutionarysearch\]\. DFT functional discovery, however, poses a fundamentally different challenge\. In algorithmic domains, the starting solution is typically far from optimal, and each evolutionary iteration can identify an objective actionable improvement \(a better heuristic, a tighter bound, a more efficient subroutine\) so that exploitation\-heavy parent selection reliably drives progress\. In functional discovery, by contrast, many apparent improvements are entangled with caveats that scientists have meticulously determined over decades, including conservation laws, numerical stability, interpretability, and transferability\[kaplan2023predictive\]\. Moreover, the XC functional baseline is already a highly optimized human design, leaving few local parameter modifications that do not violate physical constraints or regress broader performance\. Exploitation\-heavy sampling therefore stalls, and repeated micro\-perturbation of a single elite tends to overfit the benchmarks rather than improve genuine generalization\. Incremental refinement of a near\-optimal form is thus not merely ineffective but actively harmful \(see Appendix[G](https://arxiv.org/html/2605.05460#A7)\)\. Genuine progress instead requires firstexploringstructurally novel regions of functional\-form space \(new descriptor combinations, rational enhancement\-factor topologies, and spin\-dependent architectures\), and only thenrefiningwithin a newly discovered promising region\. This exploration\-then\-refinement dynamic accounts for most of the significant score gains observed in our search \(see Appendix[D](https://arxiv.org/html/2605.05460#A4)and[E](https://arxiv.org/html/2605.05460#A5)\)\.

Prior searches over symbolic functional forms can produce competitive functionals\[ma2022evolving\], but semi\-random mutations limit the search to a predetermined set of modifications and only weakly exploit accumulated knowledge of what has succeeded or failed\. Instead, could a symbolic search be performed with an entirely LLM based agent loop, one that can inspect evolutionary history, diagnose plateaus, and propose qualitatively different forms? More importantly, could these agents actually improve upon the most accurate human\-designed functionals? This report focuses on that scientific question\. Our goal is to design an agentic search capable of improving upon the SOTA exchange\-correlation functional forms, analyze the highest\-performing discovered functionals, and extract the design principles and limitations that emerged from the search trajectories\.

### 1\.1Agentic Discovery Framework

![Refer to caption](https://arxiv.org/html/2605.05460v1/x1.png)Figure 1:Sketch of the agentic search framework: a Plan–Execute–Summarize loop over a multi\-island evolutionary population with structured memory management\.Plan–Execute–Summarize Loop\.The agentic search is built on LoongFlow\[wan2025loongflowdirectedevolutionarysearch\], an evolutionary framework that replaces random mutations with a structured cognitive loop\. Each evolutionary iteration passes through three LLM stages: aPlannerthat produces a natural\-language blueprint for the next functional\-form modification, anExecutorthat translates this blueprint into executable JAX code, and aSummarizerthat compares intent against outcome and writes a structured diagnosis into the evolutionary memory\. The exact task definition and full component descriptions are given in Appendices[C](https://arxiv.org/html/2605.05460#A3)and[B](https://arxiv.org/html/2605.05460#A2)\.

Population Structure\.The search maintains a population of candidate functionals across n independently evolving islands, each with a capped archive of top solutions preserving structurally diverse stepping stones\[mouret2015mapelites\]\. Parent selection uses an exploration\-heavy 80/20 split \(80% uniform\-random, 20% from the elite archive\), and periodic migration transfers top solutions between islands, enabling cross\-lineage fusion of independently matured ideas\. This exploration\-heavy regime prioritizes structural diversity over exploitation of the current best, an explicit choice to match the nature of the problem \(Section[2\.3](https://arxiv.org/html/2605.05460#S2.SS3)\)\. Full details are given in Appendix[B\.4](https://arxiv.org/html/2605.05460#A2.SS4)\.

Evolutionary Memory\.Over 400\+ iterations the accumulated search record grows far beyond any single LLM context window\. The framework addresses this through an evolutionary memoryℳt\\mathcal\{M\}\_\{t\}\. In addition to the full parent state, search context is delivered to the planner via two channels: per\-lineage retrieval tools that surface ancestor plans, solutions, and summarizer diagnoses, and a pre\-synthesized global summary of exhausted dead\-end strategies compiled across all islands\. This ensures agents avoid revisiting failed strategies without requiring the full history to fit in a single prompt\. Full details are given in Appendix[B\.5](https://arxiv.org/html/2605.05460#A2.SS5)\.

### 1\.2Functional Search Space And Evaluation

Functional form and search target\.The evolving object is the semi\-local part of a range\-separated hybrid meta\-GGA density functional\[ma2022evolving,mardirossian2017thirty\]\. In a range\-separated hybrid meta\-GGA, exchange correlation is built from the local spin\-densitiesρσ\\rho\_\{\\sigma\}, their gradients\|∇ρσ\|\|\\nabla\\rho\_\{\\sigma\}\|, and the kinetic\-energy densitiesτσ\\tau\_\{\\sigma\}\. The total exchange energy combines short\-range semi\-local exchange, with a fractioncxc\_\{x\}of short\-range exact \(Hartree\-Fock\) exchange and full long\-range exact exchange\. Correlation is likewise decomposed into same\-spin and opposite\-spin semi\-local channels plus a VV10 nonlocal dispersion term\[vydrov2010vv10\]:

Ex\\displaystyle E\_\{x\}=Ex,srmGGA\+cx​Ex,srexact\+Ex,lrexact,\\displaystyle=E\_\{x,\\mathrm\{sr\}\}^\{\\mathrm\{mGGA\}\}\+c\_\{x\}\\,E\_\{x,\\mathrm\{sr\}\}^\{\\mathrm\{exact\}\}\+E\_\{x,\\mathrm\{lr\}\}^\{\\mathrm\{exact\}\},\(1\)Ec\\displaystyle E\_\{c\}=Ec,ssmGGA\+Ec,osmGGA\+Ec,nlVV10\.\\displaystyle=E\_\{c,\\mathrm\{ss\}\}^\{\\mathrm\{mGGA\}\}\+E\_\{c,\\mathrm\{os\}\}^\{\\mathrm\{mGGA\}\}\+E\_\{c,\\mathrm\{nl\}\}^\{\\mathrm\{VV10\}\}\.\(2\)Each semi\-local channel is defined by multiplying the appropriate local density approximation \(LDA\) reference energy density by a dimensionless enhancement factorg​\(w,u\)g\(w,u\); the form of these enhancement factors is the object of the evolutionary search\. The range\-separation parameterω=0\.3\\omega=0\.3, the short\-range exact\-exchange fractioncx=0\.15c\_\{x\}=0\.15, and the VV10 parameters are inherited fromω\\omegaB97M\-V and held fixed throughout; only the three enhancement factorsgxg\_\{x\},gc,ssg\_\{c,\\mathrm\{ss\}\}, andgc,osg\_\{c,\\mathrm\{os\}\}are evolved\.

The raw densitiesρσ\\rho\_\{\\sigma\}enter through the LDA reference energy densities that the enhancement factors multiply\. Since\|∇ρ\|\|\\nabla\\rho\|andτ\\tauare unbounded and not dimensionless, the gradients and kinetic\-energy densities enter the enhancement factors as dimensionless descriptors,

s=\|∇ρ\|/ρ4/3,t=τ/ρ5/3,s=\|\\nabla\\rho\|/\\rho^\{4/3\},\\qquad t=\\tau/\\rho^\{5/3\},\(3\)whereτ=12​∑i\|∇φi\|2\\tau=\\tfrac\{1\}\{2\}\\sum\_\{i\}\|\\nabla\\varphi\_\{i\}\|^\{2\}is the positive\-definite kinetic energy density\. The bounded variables are

w=kσ−tkσ\+t,u=γ​s21\+γ​s2,w=\\frac\{k\_\{\\sigma\}\-t\}\{k\_\{\\sigma\}\+t\},\\qquad u=\\frac\{\\gamma s^\{2\}\}\{1\+\\gamma s^\{2\}\},\(4\)withw∈\[−1,1\]w\\in\[\-1,1\]andu∈\[0,1\]u\\in\[0,1\], wherekσ=310​\(6​π2\)2/3k\_\{\\sigma\}=\\tfrac\{3\}\{10\}\(6\\pi^\{2\}\)^\{2/3\}is the UEG kinetic\-energy prefactor andγ\\gammais a fixed gradient parameter inherited fromω\\omegaB97M\-V\.

Baseline and evaluator\.The initial baseline isω\\omegaB97M\-V, the gold standard range\-separated hybrid meta\-GGA density functional\[mardirossian2016wB97MV,mardirossian2017thirty\]\. In the baseline form, exchange is represented by a compact polynomial inwwanduu, and correlation is represented by low\-order polynomials in spin\-resolved or averaged descriptors\. Candidate descendants are implemented in a JAX\-based codebase that evaluates functional forms on electron densities held fixed at theω\\omegaB97M\-V self\-consistent values, following a non\-self\-consistent evaluation protocol\[ma2022evolving\]\. The empirical parameters of each candidate are fitted by minimizing the training error using a quasi\-Newton optimization algorithm \(L\-BFGS\), with gradients computed via auto\-differentiation through JAX\. Loss is measured through the weighted root mean squared deviation \(WRMSD\) on the MGCDB84 dataset\. To avoid leaking the final test set in the evolution, evaluation uses the validation subset WRMSD to calculate the score\. Following the procedure established forω\\omegaB97M\-V, final performance is measured by the WRMSD computed over the combined validation and test subsets of MGCDB84, excluding the RG10 dataset\. RG10 is excluded because it does not belong to any of the eight thermochemical data types into which the remaining 82 datasets are classified; in total 3547 points make up this final set\[mardirossian2016wB97MV\]\. MGCDB84 training\-split was minimally changed with the transfer of two molecules from the test set to the validation set \(C20C24 subset\)\. This was done because this subset is known to be significantly out\-of\-distribution; thus, its inclusion added a critical evolution signal\.

Because each candidate is evaluated on a density that is not self\-consistent for it, the reported WRMSD values may shift once the density is allowed to relax\. Ma et al\.\[ma2022evolving\]found that self\-consistent GAS22 errors were close to their non\-self\-consistent counterparts, but analogous SCF validation remains a necessary follow\-up\.

Scoring\.Letssdenote a candidate functional form together with its optimized empirical parametersΘ\\Theta\. Parameters are fitted on the MGCDB84 training split, and every candidate is scored by the validation weighted root\-mean\-square deviation

J​\[s​\(Θ\)\]=WRMSDval​\[s​\(Θ\)\],R​\[s​\(Θ\)\]=JtargetJ​\[s​\(Θ\)\]\.J\[s\(\\Theta\)\]=\\mathrm\{WRMSD\}\_\{\\mathrm\{val\}\}\[s\(\\Theta\)\],\\qquad R\[s\(\\Theta\)\]=\\frac\{J\_\{\\mathrm\{target\}\}\}\{J\[s\(\\Theta\)\]\}\.\(5\)HereJtarget=3\.45J\_\{\\mathrm\{target\}\}=3\.45kcal is an arbitrarily chosen target approximately 15% below theω\\omegaB97M\-V validation WRMSD \(4\.02 kcal\)\. To discourage unphysical solutions, the raw score is penalized for violations of the enforced physical constraints \(Section[1\.3](https://arxiv.org/html/2605.05460#S1.SS3)\):

Revlv​\[s​\(Θ\)\]=R​\[s​\(Θ\)\]×0\.9n,R\_\{\\mathrm\{evlv\}\}\[s\(\\Theta\)\]=R\[s\(\\Theta\)\]\\times 0\.9^\{\\,n\},\(6\)wherennis the number of constraints violated by the candidate\. The search therefore favors functional forms whose parameters can be optimized successfully, that generalize to the held\-out validation subset, and that satisfy the enforced physical constraints\. The test split is reserved for final reporting only\.

### 1\.3Verification and Robustness

Beyond numerical robustness, the discovered functionals are screened against three physical constraints \+ one technical constraint that any well\-behaved meta\-GGA exchange\-correlation functional should satisfy\[kaplan2023predictive\]\. Results collected without penalizing violations of these constraints reveal enforcement is not merely precautionary but forms essential guardrails \(see Appendix[F](https://arxiv.org/html/2605.05460#A6)\)\.

1\. Spin symmetry\.The XC energy density must be invariant under exchange of spin channels:εxc​\(𝐫;ρα,ρβ\)=εxc​\(𝐫;ρβ,ρα\)\\varepsilon\_\{\\mathrm\{xc\}\}\(\\mathbf\{r\};\\rho\_\{\\alpha\},\\rho\_\{\\beta\}\)=\\varepsilon\_\{\\mathrm\{xc\}\}\(\\mathbf\{r\};\\rho\_\{\\beta\},\\rho\_\{\\alpha\}\)at every point in space\. This was tested on an O2triplet molecular grid \(whereρα≠ρβ\\rho\_\{\\alpha\}\\neq\\rho\_\{\\beta\}\) by swapping the spin channels and comparing the resulting energy densities pointwise\. The maximum pointwise deviationmax𝐫⁡\|εxcα​β−εxcβ​α\|\\max\_\{\\mathbf\{r\}\}\|\\varepsilon\_\{\\mathrm\{xc\}\}^\{\\alpha\\beta\}\-\\varepsilon\_\{\\mathrm\{xc\}\}^\{\\beta\\alpha\}\|must fall below10−510^\{\-5\}Ha\.

2\. Uniform electron gas \(UEG\) exchange limit\.In the uniform electron gas limit \(∇ρ=0\\nabla\\rho=0,τ=τUEG\\tau=\\tau\_\{\\mathrm\{UEG\}\}\), the short\-range DFT exchange enhancement factor, extracted by dividing out both the LDA exchange and the RSH attenuation factorfSR​\(ρ,ω\)f\_\{\\mathrm\{SR\}\}\(\\rho,\\omega\), must satisfygx​\(s=0,t=kc\)=cx,0g\_\{x\}\(s\{=\}0,\\,t\{=\}k\_\{c\}\)=c\_\{x,0\}, whereaHF=0\.15a\_\{\\mathrm\{HF\}\}=0\.15is the short\-range Hartree–Fock exchange mixing fraction ofω\\omegaB97M\-V, givingcx,0=1−aHF=0\.85c\_\{x,0\}=1\-a\_\{\\mathrm\{HF\}\}=0\.85, andkc=310​\(6​π2\)2/3k\_\{c\}=\\tfrac\{3\}\{10\}\(6\\pi^\{2\}\)^\{2/3\}is the UEG value of the kinetic descriptortt\. This was tested on synthetic unpolarized UEG grids at four representative densities spanning typical valence\-electron regimes\. The maximum deviation\|gx−cx,0\|\|g\_\{x\}\-c\_\{x,0\}\|must remain below 0\.5% of true UEG exchange\.

3\. Uniform coordinate scaling of exchange\.Under uniform coordinate scaling𝐫→𝐫/λ\\mathbf\{r\}\\to\\mathbf\{r\}/\\lambda, the dimensionless descriptorss=\|∇ρ\|/ρ4/3s=\|\\nabla\\rho\|/\\rho^\{4/3\}andt=τ/ρ5/3t=\\tau/\\rho^\{5/3\}are invariant\. Likewise, the exchange enhancement factorgx​\(s,t\)g\_\{x\}\(s,t\)must therefore be unchanged\. This was verified on the H2O molecular grid at three scaling factors \(λ=0\.5,2\.0,5\.0\\lambda=0\.5,2\.0,5\.0\) and comparing the enhancement factor before and after scaling, with the RSH attenuation factor divided out\. To achieve scaling invariance, the energy\-weighted integrated deviation must fall below10−410^\{\-4\}Ha\.

4\. AE18 grid convergence\.Atomization energies of the H\-Ar are computed on two integration grids: a coarse\(99,590\)\(99,\\,590\)and a fine\(250,974\)\(250,\\,974\)Lebedev grid\. Following the grid\-sensitivity protocol established forω\\omegaB97M\-V\[mardirossian2016wB97MV,liang2026coach\], the maximum absolute energy difference across the 18 atomization energies must remain below 0\.015 kcal/mol\. This threshold ensures that the functional is sufficiently insensitive to integration\-grid resolution for practical use, including reliable second\-order properties\[sitkiewicz2022reliable\]\. The behavior of the search when these constraints are not enforced is examined in Appendix[F](https://arxiv.org/html/2605.05460#A6)\.

## 2Results

### 2\.1Overall Performance

Table[1](https://arxiv.org/html/2605.05460#S2.T1)summarizes the best functionals discovered by each LLM backend under each constraint regime\. The strongest constrained functional,SAFS26\-a, reduces the WRMSDtotto 3\.70 kcal/mol \(∼9%\{\\sim\}9\\%improvement overω\\omegaB97M\-V\) while satisfying all enforced physical constraints except grid convergence\. A parallel grid\-constrained search producedSAFS26\-b, which passes all four constraints at a WRMSDtotof 3\.83 kcal/mol \(∼6%\{\\sim\}6\\%improvement\)\. Both top agent functionals were discovered by the Seed 2\.0 backend\[seed2026modelcard\]\. The GPT 5\.4 functionals improve on validation WRMSD but exhibit a widening validation\-test gap indicative of benchmark overfitting \(Section[3](https://arxiv.org/html/2605.05460#S3)\)\. The best unconstrained functional\(SAFS26\-x\) achieves competitive performance but fails all four physical constraints \(Appendix[F](https://arxiv.org/html/2605.05460#A6)\)\.

FunctionalModelConstraintsWRMSDtot% Improve\.ω\\omegaB97M\-V——4\.07—GAS22——3\.777\.4%SAFS26\-xSeed 2\.0None3\.816\.6%SAFS26\-aSeed 2\.0Physical3\.709\.1%SAFS26\-bSeed 2\.0Physical \+ Grid3\.835\.9%AFS26\-aGPT 5\.4Physical3\.904\.3%AFS26\-bGPT 5\.4Physical \+ Grid3\.904\.2%Table 1:WRMSDtotin kcal/mol on MGCDB84 for the best functionals discovered by each LLM backend, grouped by model and constraint regime\. Theω\\omegaB97M\-V baseline and GAS22\[ma2022evolving\]are shown for reference\. Full per\-split results are given in Table[3](https://arxiv.org/html/2605.05460#A1.T3)\. SAFS26 = Seed Agentic Functional Search 2026; AFS26 = Agentic Functional Search 2026\.
### 2\.2Per\-Category Accuracy

Categoryω\\omegaB97M\-VGAS22SAFS26\-aSAFS26\-b\\multirow@dimawidth 0ptNon\-covalent

\{\\left\.\\vbox\{\\vrule height=0\.0pt\\multirow@dima width 0pt\}\\textnormal\{\\rotatebox\[origin=\{c\}\]\{90\.0\}\{\\scriptsize Non\-covalent\}\}\\right\\\{RG102\.722\.172\.762\.66NCED2\.252\.102\.422\.15NCEC6\.355\.815\.965\.95NCD2\.362\.372\.792\.75\\multirow@dimawidth 0ptCovalent

\{\\left\.\\vbox\{\\vrule height=0\.0pt\\multirow@dima width 0pt\}\\textnormal\{\\rotatebox\[origin=\{c\}\]\{90\.0\}\{\\scriptsize Covalent\}\}\\right\\\{IE2\.982\.652\.682\.81ID13\.2510\.9310\.3012\.23TCE3\.974\.163\.633\.78TCD1\.992\.361\.762\.56BH8\.287\.757\.717\.37Table 2:WRMSD in kcal/mol per MGCDB84 reaction category \(validation and test splits only\) forω\\omegaB97M\-V, GAS22, and the two strongest discovered functionals\. Cells are shaded on a continuous green \(best\) to red \(worst\) scale within each row\.Table[2](https://arxiv.org/html/2605.05460#S2.T2)decomposes the aggregate WRMSD improvement into category\-level contributions\. Functional performance is often a tradeoff between different chemical domains\. This is reflected in the per\-category performance, which suggests non\-uniform improvement, especially between GAS22 andSAFS26\-a\. GAS22 leads in the non\-covalent categories: non\-covalent clusters \[easy\] \(NCEC\), non\-covalent dimers \[easy\] \(NCED\), non\-covalent dimers \(NCD\) \[difficult\], as well as rare\-gas dimers, reflecting well\-tuned exchange\-correlation behavior in the low\-density regime that governs intermolecular binding\.SAFS26\-a, by contrast, achieves its largest gains on covalent and thermochemical categories: it substantially outperforms bothω\\omegaB97M\-V and GAS22 on the thermochemistry \[difficult\] \(TCD\), thermochemistry \[easy\] \(TCE\), and isomerizations \[difficult\] \(ID\), all categories governed by short\-range exchange\-correlation in high\-density covalent bonding regions\.

### 2\.3Search Dynamics

![Refer to caption](https://arxiv.org/html/2605.05460v1/score_plot_test.png)Figure 2:Score trajectories over 400 iterations per backend–constraint regime\. Dots are per\-island candidates, the line is the running best validation score, and stars mark the top scoring test set score \(Section[1\.2](https://arxiv.org/html/2605.05460#S1.SS2)\)\. Seed 2\.0 explores more broadly with test stars above its validation best, while GPT 5\.4 clusters tightly with test oscillating at or below the starting point, signaling validation overfitting\.The score trajectories in Figure[2](https://arxiv.org/html/2605.05460#S2.F2)hint at a marked difference in proposal style between model backends\. Seed 2\.0 produces a wide score spread from early iterations, reflecting structurally diverse seed solutions rather than incremental refinement of a single family\. This diversity seeds the population with structurally distinct solutions, giving the multi\-island architecture richer material for recombination and allowing the cumulative best to saturate sooner\. Population\-wide plan analysis reveals GPT 5\.4’s tighter score clustering is a by\-product of ignoring the parent solution in favor of the top performing seed\. Ablation study in Appendix[G](https://arxiv.org/html/2605.05460#A7)reveals this proposal pattern also causes overfitting issues\. Indeed, although GPT 5\.4 performance is tied with the highest score \(lowestWRMSDv​a​l\\text\{WRMSD\}\_\{val\}\) and displays a significantly higher average score, this is completely offset by their performance on the held out test\-set \(see Table[3](https://arxiv.org/html/2605.05460#A1.T3)\)\.

### 2\.4Structural Analysis:SAFS26\-a

The evolutionary trajectory leading toSAFS26\-ais detailed in Appendix[D](https://arxiv.org/html/2605.05460#A4); below we describe the structural changes in the final functional relative toω\\omegaB97M\-V\.

Exchange\.The baselineω\\omegaB97M\-V exchange enhancement is a three\-term polynomial in the bounded descriptors\(w,u\)\(w,u\):

gx=cx,0\+c1​w\+c2​u\.g\_\{x\}=c\_\{x,0\}\+c\_\{1\}w\+c\_\{2\}u\.\(7\)SAFS26\-aretains this polynomial structure but adds one new descriptor: a bounded cross termv=s​t/\(1\+s​t\+ϵ\)v=st/\(1\+st\+\\epsilon\)coupling gradient and kinetic information

The exchange enhancement becomes

gx=Px​\(w,u\)\+cx,v,0​v\+cx,v,1​v2,g\_\{x\}=P\_\{x\}\(w,u\)\+c\_\{x,v,0\}\\,v\+c\_\{x,v,1\}\\,v^\{2\},\(8\)wherePx​\(w,u\)P\_\{x\}\(w,u\)is the baselineω\\omegaB97M\-V polynomial form\.

Correlation\.Inω\\omegaB97M\-V, same\-spin and opposite\-spin correlation are low\-order polynomials in\(w,u\)\(w,u\):

gα​α=Pcss​\(wα,uα\),gβ​β=Pcss​\(wβ,uβ\),gα​β=Pcos​\(wavg,uavg\)\.g\_\{\\mathrm\{\\alpha\\alpha\}\}=P\_\{\\mathrm\{css\}\}\(w\_\{\\alpha\},u\_\{\\alpha\}\),\\qquad g\_\{\\mathrm\{\\beta\\beta\}\}=P\_\{\\mathrm\{css\}\}\(w\_\{\\beta\},u\_\{\\beta\}\),\\qquad g\_\{\\mathrm\{\\alpha\\beta\}\}=P\_\{\\mathrm\{cos\}\}\(w\_\{\\mathrm\{avg\}\},u\_\{\\mathrm\{avg\}\}\)\.\(9\)SAFS26\-areplaces these with bounded rational enhancement factors\. For each correlation channel, two additional bounded cross descriptors are constructed from the von Barth–Hedin spin\-polarization interpolation functionfζ​\(ζ\)=\[\(1\+ζ\)4/3\+\(1−ζ\)4/3−2\]/\(24/3−2\)f\_\{\\zeta\}\(\\zeta\)=\[\(1\+\\zeta\)^\{4/3\}\+\(1\-\\zeta\)^\{4/3\}\-2\]/\(2^\{4/3\}\-2\)\[vonBarth1972\]:

z=fζ​\(ζ\)⋅s⋅t1\+fζ​\(ζ\)⋅s⋅t\+ϵ,x=z⋅s⋅t1\+z⋅s⋅t\+ϵ,z=\\frac\{f\_\{\\zeta\}\(\\zeta\)\\cdot s\\cdot t\}\{1\+f\_\{\\zeta\}\(\\zeta\)\\cdot s\\cdot t\+\\epsilon\},\\qquad x=\\frac\{z\\cdot s\\cdot t\}\{1\+z\\cdot s\\cdot t\+\\epsilon\},\(10\)where same\-spin channels use per\-spinssandttwhile opposite\-spin uses spin\-averaged values\. The enhancement factor is then a bounded rational:

g=N​\(w,u,v,z,x\)D​\(w,u,v,z,x\),g=\\frac\{N\(w,u,v,z,x\)\}\{D\(w,u,v,z,x\)\},\(11\)with numerator

N=P​\(w,u\)\+cv​\[0\]​v\+cv​\[1\]​v2\+tanh⁡\(cz​\[0\]​z\+cz​\[1\]​z2\)\+cx​\[0\]​x\+cx​\[1\]​x2,N=P\(w,u\)\+c\_\{v\}\[0\]\\,v\+c\_\{v\}\[1\]\\,v^\{2\}\+\\tanh\\\!\\bigl\(c\_\{z\}\[0\]\\,z\+c\_\{z\}\[1\]\\,z^\{2\}\\bigr\)\+c\_\{x\}\[0\]\\,x\+c\_\{x\}\[1\]\\,x^\{2\},\(12\)and denominator

D=1\+∑isoftplus​\(di\)​ϕi,D=1\+\\sum\_\{i\}\\mathrm\{softplus\}\(d\_\{i\}\)\\,\\phi\_\{i\},\(13\)where the basis\{ϕi\}\\\{\\phi\_\{i\}\\\}comprises the monomials ofPc​\(w,u\)P\_\{c\}\(w,u\)together with

\{v,v2,z,z2,x,x2\}\.\\\{v,\\;v^\{2\},\\;z,\\;z^\{2\},\\;x,\\;x^\{2\}\\\}\.\(14\)The two channels share this functional form but differ inPc​\(w,u\)P\_\{c\}\(w,u\)\(each inheriting its monomial powers from theω\\omegaB97M\-V baseline\)\.softplus​\(⋅\)\\mathrm\{softplus\}\(\\cdot\)helps ensure non\-negative denominator coefficients, suppressing the risk of sign changes\.

The total number of trainable parameters is 50 \(compared to 12 in theω\\omegaB97M\-V baseline\): 4 exchange parameters \(2 polynomial, 2 cross\-termvv\) and 46 correlation parameters \(11 same\-spin numerator, 11 same\-spin denominator, 12 opposite\-spin numerator, 12 opposite\-spin denominator\)\. In addition, 6 baselineω\\omegaB97M\-V structural parameters are held fixed: the descriptor mappings \(γx\\gamma\_\{x\},γcss\\gamma\_\{\\mathrm\{css\}\},γcos\\gamma\_\{\\mathrm\{cos\}\}\), and the polynomial power templates\.

### 2\.5Structural Analysis:SAFS26\-b

The evolutionary trajectory leading toSAFS26\-bis detailed in Appendix[E](https://arxiv.org/html/2605.05460#A5)\. The final functional differs fromω\\omegaB97M\-V primarily in its correlation channels, with more modest changes to exchange\.

Exchange\.Relative toω\\omegaB97M\-V, two modifications were introduced\. First, a new descriptorv=1/\(1\+α2\)v=1/\(1\+\\alpha^\{2\}\), whereα=\(τ−τW\)/τUEG\\alpha=\(\\tau\-\\tau\_\{W\}\)/\\tau\_\{\\mathrm\{UEG\}\}is the iso\-orbital indicator, was added to the exchange polynomial with two cross termsw​vwvandu​vuv\[Perdew2015SCAN\]\. Second, a multiplicative sigmoid correction was appended:

gx=Px​\(w,u,v\)⋅\(1\+csig,x⋅σx\),g\_\{x\}=P\_\{x\}\(w,u,v\)\\cdot\\bigl\(1\+c\_\{\\mathrm\{sig\},x\}\\cdot\\sigma\_\{x\}\\bigr\),\(15\)whereσ​\(⋅\)=1/\(1\+e−\(⋅\)\)\\sigma\(\\cdot\)=1/\(1\+e^\{\-\(\\cdot\)\}\)denotes the logistic sigmoid andσx=σ​\(ax​\(s2−bx\)\)\\sigma\_\{x\}=\\sigma\(a\_\{x\}\(s^\{2\}\-b\_\{x\}\)\)\. Although this breaks the UEG constraint slightly, the coefficient values are small, meaning\(1\+csig,x⋅σx\)≈1\(1\+c\_\{\\mathrm\{sig\},x\}\\cdot\\sigma\_\{x\}\)\\approx 1falls under the threshold\.

Same\-spin correlation\.The same\-spin enhancement factor combines a polynomial base in\(w,u\)\(w,u\), a sigmoid correction, additiveu⋅wu\\cdot wterms, and a multiplicative rational correction incorporating the Wigner–Seitz radius\(rsr\_\{s\}\), thevvdescriptor, and spin polarization:

gα​α=\[Pcss​\(wα,uα\)⋅\(1\+csig,ss⋅σssα\)\+c1​uα​wα\+c2​uα2​wα2\]⋅\(1\+βss​rs,α\+cv,ss​rs,α​vα\+cζ,ss​ζα1\+Γss​rs,α2\),g\_\{\\alpha\\alpha\}=\\Bigl\[P\_\{\\mathrm\{css\}\}\(w\_\{\\alpha\},u\_\{\\alpha\}\)\\cdot\(1\+c\_\{\\mathrm\{sig,ss\}\}\\cdot\\sigma\_\{\\mathrm\{ss\}\}^\{\\alpha\}\)\+c\_\{1\}u\_\{\\alpha\}w\_\{\\alpha\}\+c\_\{2\}u\_\{\\alpha\}^\{2\}w\_\{\\alpha\}^\{2\}\\Bigr\]\\cdot\\biggl\(1\+\\frac\{\\beta\_\{\\mathrm\{ss\}\}r\_\{s,\\alpha\}\+c\_\{v,\\mathrm\{ss\}\}r\_\{s,\\alpha\}v\_\{\\alpha\}\+c\_\{\\zeta,\\mathrm\{ss\}\}\\zeta\_\{\\alpha\}\}\{1\+\\Gamma\_\{\\mathrm\{ss\}\}r\_\{s,\\alpha\}^\{2\}\}\\biggr\),\(16\)wheregβ​βg\_\{\\beta\\beta\}is defined analogously fromβ\\beta\-spin quantities,σssα=σ​\(ass​\(tα−kc/2\)\)\\sigma\_\{\\mathrm\{ss\}\}^\{\\alpha\}=\\sigma\(a\_\{\\mathrm\{ss\}\}\(t\_\{\\alpha\}\-k\_\{c\}/2\)\)gates on the kinetic descriptor, andζα=ζ​\(1−ζ\)\\zeta\_\{\\alpha\}=\\zeta\(1\-\\zeta\),ζβ=−ζ​\(1\+ζ\)\\zeta\_\{\\beta\}=\-\\zeta\(1\+\\zeta\)are spin\-dependent polarization terms ensuring spin symmetry\.

Opposite\-spin correlation\.The opposite\-spin channel received the most extensive modifications\. The spin\-averaged descriptors used in this channel are computed assavg=\(\(sα2\+sβ2\)/2\)1/2s\_\{\\mathrm\{avg\}\}=\(\(s\_\{\\alpha\}^\{2\}\+s\_\{\\beta\}^\{2\}\)/2\)^\{1/2\},tavg=2​tα​tβ/\(tα\+tβ\+ϵ\)t\_\{\\mathrm\{avg\}\}=2\\,t\_\{\\alpha\}\\,t\_\{\\beta\}/\(t\_\{\\alpha\}\+t\_\{\\beta\}\+\\epsilon\),vavg=\(vα\+vβ\)/2v\_\{\\mathrm\{avg\}\}=\(v\_\{\\alpha\}\+v\_\{\\beta\}\)/2, andrs,avg=\(3/\(4​π​ρ\)\)1/3r\_\{s,\\mathrm\{avg\}\}=\(3/\(4\\pi\\rho\)\)^\{1/3\}\. The base enhancement factor is a polynomial in three descriptors\(wavg,uavg,z\)\(w\_\{\\mathrm\{avg\}\},u\_\{\\mathrm\{avg\}\},z\), wherez=fζ​\(ζ\)z=f\_\{\\zeta\}\(\\zeta\)is the von Barth–Hedin interpolation function defined above andPcosP\_\{\\mathrm\{cos\}\}includes 11 terms with cross powers inzzup to quadratic order:

gα​β=Pcos​\(wavg,uavg,z\)⋅\(1\+βos\+cv,os​vavg\+cζ,os​\|ζ\|1\+Γos​rs,avg2\),g\_\{\\alpha\\beta\}=P\_\{\\mathrm\{cos\}\}\(w\_\{\\mathrm\{avg\}\},u\_\{\\mathrm\{avg\}\},z\)\\cdot\\biggl\(1\+\\frac\{\\beta\_\{\\mathrm\{os\}\}\+c\_\{v,\\mathrm\{os\}\}v\_\{\\mathrm\{avg\}\}\+c\_\{\\zeta,\\mathrm\{os\}\}\|\\zeta\|\}\{1\+\\Gamma\_\{\\mathrm\{os\}\}r\_\{s,\\mathrm\{avg\}\}^\{2\}\}\\biggr\),\(17\)
The total number of trainable parameters is 19 \(compared to 50 forSAFS26\-aand 12 for theω\\omegaB97M\-V baseline\): 8 scalar correction parameters \(2 same\-spin additive, 3rsr\_\{s\}/vv/ζ\\zetafor each of same\-spin and opposite\-spin\) and 11 opposite\-spin polynomial coefficients \(5 same\-spin polynomial coefficients are inherited and frozen from the parent\)\.

## 3Discussion

### 3\.1The Value of Structured Evolution Protocol

A natural question is whether a general\-purpose coding agent, operating in a single context window, could achieve an equivalent result if given the same functional\-design problem and evaluation pipeline\. We argue that two architectural features of the framework, evolutionary memory and population structure, address fundamental limitations of the single\-agent setting\.

Evolutionary memory enables sustained structured search\.Over 400\+ iterations the accumulated evolutionary record \(plans, summaries, code, constraint\-test outcomes, and scores\) grows far beyond the context window of any current LLM\. A single\-thread agent has limited mechanisms to synthesize lessons across this history: even with tool\-augmented context management, the working memory available at each decision point is bounded, making it difficult to maintain a coherent picture of which strategies have been tried, why they failed, and which structural regions remain unexplored \(running multiple such agents in parallel does not resolve this: without shared state, each independently rediscovers the same dead ends\)\. Our evolutionary memoryℳt\\mathcal\{M\}\_\{t\}\(Appendix[B\.5](https://arxiv.org/html/2605.05460#A2.SS5)\) addresses this through three granularities of context: the planner always receives its direct parent’s plan, solution, and summarizer diagnosis as immediate context; retrieval tools allow it to selectively query deeper into its lineage history, surfacing ancestor plans and diagnoses that track the broader research trajectory; and a global dead\-end list, synthesized across all islands, prevents population\-wide repetition of exhausted strategies\.

Population structure prevents premature convergence\.A single\-thread agent naturally gravitates toward exploiting the current best solution\. This is because that is what occupies its context and attention, while suboptimal but structurally distinct solutions are discarded rather than preserved for later refinement\. The population structure \(Appendix[B\.4](https://arxiv.org/html/2605.05460#A2.SS4)\) addresses this through two complementary mechanisms\. First, the population archive retains structurally diverse solutions that are globally suboptimal but occupy useful regions of the search space; exploration\-heavy parent selection provides repeated chances for these ideas to be refined\. Second, island isolation protects nascent ideas through their initially unpromising stages\. The decisive fusion that producedSAFS26\-arequired two qualitatively different ideas—spin\-enhanced exchange and rational bounded correlation—to mature independently over multiple iterations on separate islands before being combined at iteration 197 \(Appendix[D](https://arxiv.org/html/2605.05460#A4)\)\. The ablation study in Appendix[G](https://arxiv.org/html/2605.05460#A7)demonstrates the consequence of removing these mechanisms\. Naively parallelized agent systems can share limited context through shared filesystems or message passing\[anthropic2026agentictrends\], but these coordination mechanisms do not provide the structure needed for reliable ultra\-long horizon searches: deciding which suboptimal ideas to preserve, when to revisit them, or how to introduce ideas across lineages\.

We should note that significant improvements have been made employing parallel and interacting agents\[anthropic2026agentictrends,lin2026selfdriving,carlini2026compiler\]\. Nevertheless, these systems have so far been applied primarily to software engineering tasks where correctness is binary and work decomposes naturally into independent sub\-tasks\. Scientific discovery search poses distinct coordination challenges, including a continuous and deceptive objective landscape, higher risk of benchmark exploitation, and a requirement for strict adherence to domain\-specific rules\. That said, the rapid progress in multi\-agent orchestration is encouraging and we are optimistic that agent architectures will extend their reach from engineering to open\-ended scientific optimization\.

### 3\.2The Risk of Unproductive Shortcuts

The unconstrained results in Appendix[F](https://arxiv.org/html/2605.05460#A6)illustrate a broader hazard of agentic scientific search: without careful design, optimization pressure will reward unphysical shortcuts over scientifically meaningful improvements\. During one unconstrained search,everytop\-scoring functional violated at least three of four basic constraints, and the majority violated all four\. The violations are not subtle numerical drift; they represent physics breaking failures \(antisymmetric spin terms, improper asymptotic behavior\) that were implicitly selected because they lower the training objective\. Without explicit constraint enforcement, benchmark exploitation consistently dominates over physically meaningful innovation\. This echoes a broader pattern in AI\-assisted science: recent work deploying AI for research\-level mathematics and physics similarly finds that publishable\-quality results required sustained human–AI collaboration, with fully autonomous outputs mostly limited to well\-posed problems with verifiable answers\[feng2026autonomous,woodruff2026accelerating\]\.

Even with all physical constraints enforced, however, a subtler form of benchmark gaming persists at longer search horizons\. Inspection of the evolutionary record reveals a recurring pattern: rather than building on their assigned parent solution, agents will source from the current top\-performing functional, freeze its existing parameters, and graft small additive corrections on top to eke out marginal score improvements\. This behavior becomes especially pronounced when the mechanisms designed to maintain creativity and innovation \(multi\-island isolation and exploration\-heavy parent selection\), are removed from the evolution architecture; the ablation study in Appendix[G](https://arxiv.org/html/2605.05460#A7)documents the most extreme instance, a functional that froze 61 parameters from its parent and trained only 4, achieving a record\-low validation WRMSD while generalizing substantially worse than theω\\omegaB97M\-V baseline on the test set; instead of exploring qualitatively different functional forms, the search collapsed into repeated micro\-perturbations of a single elite solution\. Although each perturbation is trained on the training split, it is evolved on its performance on the validation split\. Over many generations, the cumulative effect of hundreds of such incremental additions is to progressively leak the validation set\. The overfitting signature is unambiguous and gets worse with time: when the grid\-constrained search is extended from 400 to 800 iterations, the validation WRMSD continues to improve slightly, but the test WRMSD degrades substantially, a hallmark of data leakage rather than genuine generalization\. This underscores that physical constraints alone are insufficient to prevent all forms of benchmark exploitation\.

### 3\.3Conclusion

This work demonstrates that an agent\-guided evolutionary search can discover exchange–correlation density functionals that improve over a strong human\-designed baseline while respecting enforced constraints\. Starting fromω\\omegaB97M\-V, the search producedSAFS26\-a, which reduces the WRMSDtot∼9%\{\\sim\}9\\%and satisfies three of four constraints, andSAFS26\-b, which passes all four constraints with a∼6%\{\\sim\}6\\%WRMSDtotimprovement\. These results reveal three properties of agent\-driven search that are relevant beyond DFT\. First, structured evolutionary memory enables agents to sustain coherent searches over hundreds of iterations without collapsing into repetition or blowing up context, a capability that single\-conversation coding agents lack because their context is bounded and eventually discarded or compressed\. Second, population diversity is essential, and multiple steps were needed to maintain it: the multi\-island topology, high exploration, and large population archive maintain a pool of structurally distinct “stepping\-stone” solutions that were individually suboptimal but collectively enabled the cross\-lineage fusions responsible for the largest score improvements\. Third, creative high\-variance proposal generation is a strength rather than a liability in population\-based search; over a long span, diversity of ideas matters more than per\-iteration reliability\. This stands in contrast to recent LLM post\-training paradigms centered on verifiable domains such as mathematics, coding, and formal reasoning, where reward signals derive from objectively correct solutions\[guo2025deepseekr1,wei2025swerl\]\. Our results suggest that in discovery\-oriented settings, preserving high\-variance exploratory proposal generation may be more valuable than maximizing short\-horizon correctness\.

Equally important is the central cautionary finding: the same search machinery, without physical constraints, reliably produces functionals that game the benchmark through unproductive mechanisms \(Appendix[F](https://arxiv.org/html/2605.05460#A6)\)\. Every top\-scoring unconstrained functional violated at least three of four basic constraints\. The physical constraints imposed by domain expertise are not incidental to the result: they are what separates scientifically meaningful improvement from benchmark exploitation\. Unless these guardrails are carefully enforced or the task itself is very well defined, the trustworthy role of agent\-based scientific agents remains constrained search under human oversight rather than fully autonomous research\.

## 4Future Directions

As previously stated, validation\-set leakage remains a problem for these models\. The most direct mitigation of this is to train and score on a larger, more chemically diverse benchmark\[liang2025gold\]; a richer data distribution raises the statistical cost of the incremental micro\-perturbations that drove leakage here, making exploitation of any narrow subset less rewarding\. The increased breadth and complexity of such a corpus—spanning transition\-metal energetics, barrier heights, and electric\-field responses in addition to the thermochemistry and noncovalent interactions of MGCDB84—would simultaneously constitute a more stringent test of the agent’s capacity for genuine chemical innovation, requiring discovered functionals to generalize across a broader swath of chemical space\. This is far from the only mitigation factor, penalization scoring based on problematic behavior \(i\.e\. perturbative additions onto frozen parameters\) or semantic similarity may also afford the intended effect\. A natural progression of this pipeline would be to include a judge or rubric agent to systematically grade trial functionals on these more nebulous metrics\.

Second, the results reported here rely on non\-self\-consistent evaluation using electron densities held fixed at theω\\omegaB97M\-V solution\. Computing SCF\-converged energies for the discovered functionals is a necessary follow\-up to confirm that the observed improvements persist once the density is allowed to relax\. More broadly, incorporating periodic rounds of SCF convergence into the evolutionary search itself \(in the spirit of established protocols\[liang2026coach\]\) would ensure that the frozen\-density approximation remains valid as the functional form drifts further from the baseline, preventing the search from inadvertently optimizing for an increasingly inaccurate surrogate\.

Finally, the functional\-evolution task is naturally suited to LLM post\-training through reinforcement learning: the symbolic functional form defines a state, the LLM\-proposed modification is an action, and the WRMSD score provides a continuous reward\. Recent work on test\-time RL\[yuksekgonul2026learningdiscovertesttime\]has shown that updating model weights on a single problem \(rather than prompting a frozen LLM\) can discover stronger solutions by internalizing lessons from the search trajectory, and applying such an approach to functional evolution could complement or replace the prompt\-based evolutionary memory used here\.

## References

## Appendix AFull Train/Validation/Test WRMSD

FunctionalWRMSDtrainWRMSDvalWRMSDtestω\\omegaB97M\-V3\.324\.023\.64GAS223\.273\.743\.19SAFS26\-x2\.913\.953\.49SAFS26\-a2\.763\.653\.41SAFS26\-b2\.753\.783\.44AFS26\-a2\.813\.743\.78AFS26\-b2\.653\.653\.96Table 3:Per\-split WRMSD in kcal/mol on MGCDB84 for each functional\. Bold entries indicate the best value in each column among the discovered functionals\. Theω\\omegaB97M\-V baseline and GAS22\[ma2022evolving\]are shown for reference\.
## Appendix BEvolutionary Framework

### B\.1Planner: lineage\-aware hypothesis generation

The planner is the central component of the agentic framework\. Given a parent solutionsts\_\{t\}, it generates a natural\-language blueprintbbconditioned on the evolutionary memoryℳt\\mathcal\{M\}\_\{t\}:

b∼πθ​\(b∣st,ℳt,ℐplan\),b\\sim\\pi\_\{\\theta\}\\\!\\bigl\(b\\mid s\_\{t\},\\;\\mathcal\{M\}\_\{t\},\\;\\mathcal\{I\}\_\{\\mathrm\{plan\}\}\\bigr\),\(18\)whereℐplan\\mathcal\{I\}\_\{\\mathrm\{plan\}\}is the planner system prompt\. The construction and delivery ofℳt\\mathcal\{M\}\_\{t\}is described in Section[B\.5](https://arxiv.org/html/2605.05460#A2.SS5)\.

In practice, the planner’s blueprints are detailed enough to specify the mathematical form of the proposed change, the initialization of new parameters, and the rationale for why the change should improve accuracy\.

### B\.2Executor: implementation and verification

The executor translates the planner’s blueprint into executable Python code implementing the new functional form:

s′∼πθ​\(s′∣b,st,ℐexec\)\.s^\{\\prime\}\\sim\\pi\_\{\\theta\}\\\!\\bigl\(s^\{\\prime\}\\mid b,\\;s\_\{t\},\\;\\mathcal\{I\}\_\{\\mathrm\{exec\}\}\\bigr\)\.\(19\)Its value lies not in proposing new ideas but in implementation quality: correctly translating mathematical specifications into executable computational components, handling bounded descriptors, and ensuring numerical stability\. Before the candidate is sent to the full training\-and\-evaluation pipeline, the executor runs a local verification loop that catches syntax errors and import failures, preventing low\-quality candidates from consuming expensive evaluation resources\.

When the plan is ambiguous \(for example, specifying unpolarized formulas for a spin\-resolved functional\), the executor must interpret the physics and select the correct spin\-polarized variants\. The trajectories show instances where the executor also debugged numerical instabilities \(e\.g\., replacing manual polynomial constructions with stable JAX built\-ins\) through multiple iterative attempts before producing a working implementation\.

### B\.3Summarizer: causal analysis for future iterations

After evaluation, the summarizer compares the planner’s intent with the actual outcome and generates a structured insight:

z∼πθ​\(z∣s′,R​\(s′\),b,ℐsum\),ℳt\+1←ℳt∪\{z\}\.z\\sim\\pi\_\{\\theta\}\\\!\\bigl\(z\\mid s^\{\\prime\},\\;R\(s^\{\\prime\}\),\\;b,\\;\\mathcal\{I\}\_\{\\mathrm\{sum\}\}\\bigr\),\\qquad\\mathcal\{M\}\_\{t\+1\}\\leftarrow\\mathcal\{M\}\_\{t\}\\cup\\\{z\\\}\.\(20\)Each summary contains an executive summary, data\-driven findings \(score deltas, sibling rankings\), a root\-cause analysis of success or failure, a key insight for future iterations, and identified risks\. Critically, the summarizer is instructed to report facts and diagnose causes, not to suggest strategies: that responsibility belongs to the planner\.

This separation ensures that the accumulated memory is a reliable factual record rather than a set of potentially contradictory recommendations\.

### B\.4Population structure and selection

The population is managed by a multi\-island evolutionary model with two components:

- •Multi\-island topology\.The population is partitioned intoKKislands \(four in this work\), each evolving independently\. Different islands can develop distinct functional\-form “species” in isolation\. At regular intervals, the top elites from each island are migrated to neighboring islands\.
- •Elite archive and exploitation\.Each island maintains a capped set of thekkhighest\-scoring solutions \(elites\)\. Parent selection draws uniformly at random from this elite set, concentrating compute on the most promising regions of the search space\. The exploration rateϵ\\epsilonoverrides this and selects a uniformly random solution from the full population, providing the exploration counterpart\.

### B\.5Evolutionary memory

The evolutionary memoryℳt\\mathcal\{M\}\_\{t\}reaches the planner’s context through two complementary channels:

- •Per\-lineage retrieval\.The planner is equipped with database tools to query its own lineage history: the plans that created each ancestor and the summarizer outputs that recorded what worked or failed\. This providesintent tracking\(understanding the research trajectory of prior generations\) and localcourse correction\(avoiding strategies already tried within its lineage\)\.
- •Global memory synthesis\.Before each planning step, an LLM synthesizer scans the full iteration history across all islands and injects a list ofexhausted dead endsdirectly into the planner’s context\. A strategy is declared exhausted when it has been attempted many times across the population, never improved the score, and its failures share a common root cause\. This global channel is necessary because reconstructing a cross\-island picture through sequential tool calls would exceed practical context limits\.

## Appendix CTask Definition

The following is the exact task description passed to the LLM planner and executor at every iteration of the evolutionary search\.

Task DefinitionYou are a senior computational chemist implementing a new local \(short\-range exchange \+ correlation\) functional for a range\-separated hybrid meta\-GGA\.Goal\.Create a more physically accurate DFT functional, as measured by the WRMSD on a comprehensive thermochemistry benchmark dataset\. Aim to innovate beyond existing architectures rather than replicating them\.What You’re Building\.Thelocal/semi\-local partof a range\-separated hybrid meta\-GGA functional\. Long\-range HF exchange \(ω=0\.3\\omega\{=\}0\.3,cx=0\.15c\_\{x\}\{=\}0\.15\) and VV10 dispersion \(b=6b\{=\}6,C=0\.01C\{=\}0\.01\) arefixedand cannot be modified or trained\. You implement: attenuated short\-range exchange \+ local correlation\.Execution Notes\.•Youmustviewbase\.pybeforewriting any code\.•Re\-use existing utilities where appropriate, but you are not limited to them\.•Your solution file must contain onlyonefunctional class; full training runs automatically\.Available Utilities\.•SeparatedPxExcFunc: base class; exchange defined for polarized spins, correlation handles spin internally\.•spinsum\_and\_zeta\(rho, drho, tau\)→\\to\(rho\_total, zeta, drho\_total, tau\_total\)•spin\_scaling\_exchange\_polarized\(fn\): decorator for exchange on single spin channel\.•spin\_decompose\_correlation\(fn\): wraps correlation to return\(ec\_aa, ec\_bb, ec\_ab\)\.•spin\_fz\(ζ\\zeta\):\[\(1\+ζ\)4/3\+\(1−ζ\)4/3−2\]/\(24/3−2\)\[\(1\{\+\}\\zeta\)^\{4/3\}\+\(1\{\-\}\\zeta\)^\{4/3\}\-2\]\\,/\\,\(2^\{4/3\}\-2\);spin\_phi,spin\_ds,spin\_dx: analogous spin\-scaling functions\.•gga\_s\_from\_drho:\|∇ρ\|/ρ4/3\|\\nabla\\rho\|/\\rho^\{4/3\};mgga\_t\_from\_tau:τ/ρ5/3\\tau/\\rho^\{5/3\};fermi\_kf: Fermi wavevectorkFk\_\{F\}\.•lda\_x\_slater\_pol: Slater exchange \(single spin\);lda\_c\_pw: PW92 correlation;pw92\_g,pw92\_g\_deriv: PW92g​\(rs\)g\(r\_\{s\}\)and its derivative\.•rsh\_exchange\_factor\(rho, omega, polarized\): short\-range attenuation factor \(ω=0\.3\\omega\{=\}0\.3fixed\)\.

## Appendix DPhysical Constraint Search: Discovery ofSAFS26\-a

The evolutionary record does not show smooth monotonic improvement\. Instead, it reveals a small number of reusable functional\-form ideas that appeared on separate islands, stalled independently, and were ultimately combined across lineages\. The best functional,SAFS26\-a, emerged at the intersection of two independently evolved branches\.

### D\.1Discovery trajectory \(iterations 14–197\)

The direct lineage ofSAFS26\-aruns through two independently evolved branches on Island 2, that migrated across islands and ultimately fused at iteration 197\.

Iterations 14–97: Spin\-exchange lineage\.The first lineage began on Island 2 and proceeded through three generations of incremental descriptor additions to theω\\omegaB97M\-V baseline\. The first descendant,ExtendedWB97MV\(iteration 14\), introduced a bounded cross descriptorv=s​t/\(1\+s​t\+ϵ\)v=st/\(1\+st\+\\epsilon\), to all three enhancement factor channels \(exchange, same\-spin correlation, opposite\-spin correlation\)\. This single addition reduced the performance non\-negligibly \(WRMSDval4\.40→\\to4\.09\) kcal/mol\.

Iteration 14 summary\.“Bounded cross\-terms betweenssandttcan add accuracy without breaking constraints\. The newvvdescriptor captures gradient\-kinetic interactions that the independentwwanduudescriptors cannot represent\.”

The second descendant,SpinEnhancedExtendedWB97MV\(iteration 30\), added spin\-dependentzzdescriptors to the correlation channels only, constructed asz=f​\(ζ\)⋅s⋅t/\(1\+f​\(ζ\)⋅s⋅t\+ϵ\)z=f\(\\zeta\)\\cdot s\\cdot t/\(1\+f\(\\zeta\)\\cdot s\\cdot t\+\\epsilon\)wheref​\(ζ\)=\[\(1\+ζ\)4/3\+\(1−ζ\)4/3−2\]/\(24/3−2\)f\(\\zeta\)=\[\(1\+\\zeta\)^\{4/3\}\+\(1\-\\zeta\)^\{4/3\}\-2\]/\(2^\{4/3\}\-2\)is the von Barth–Hedin spin\-polarization interpolation function, bounded in\[0,1\]\[0,1\]with unpolarized and polarized cases returningf​\(0\)=0f\(0\)=0andf​\(±1\)=1f\(\\pm 1\)=1\. This yielded a modest improvement \(WRMSDval4\.09→\\to4\.05\) kcal/mol\. The third descendant,SpinExEnhancedWB97MV\(iteration 97\), extended the same idea to exchange, adding a spin\-dependentzxz\_\{x\}descriptor\. The gain was negligible:

Iteration 97 summary\.“The negligible improvement is attributable to high quality execution of a plan with overestimated expected gains: the short\-range exchange channel was already sufficiently optimized in prior iterations, so adding small spin\-dependent cross terms provided no meaningful accuracy benefit\.”

By iteration 97, this lineage had reached a form\-level plateau at WRMSDval≈\\approx0\.006454\.

Iterations 18–106: Rational\-correlation lineage\.A second, independent lineage independently made progress using a different approach\. Its first descendant,WB97MV\_Rational\(iteration 18\), attempted to replace both exchange and correlation enhancement factors with rational formsP​\(w,u\)/Q​\(w,u\)P\(w,u\)/Q\(w,u\)\. The raw WRMSDvalimproved \(4\.40→\\to4\.17\) kcal/mol, but the trainable denominator terms shifted the exchange enhancement factor away from the required UEG value, causing a constraint failure:

Iteration 18 summary\.“Raw valid WRMSD improved, but UEG exchange constraint failed\. The trainable denominator broke the UEG guarantee\. The constant term of the denominator must be fixed to 1 to preservegx=0\.85g\_\{x\}=0\.85at the UEG limit\.”

This failure was instructive: the summarizer’s causal diagnosis \(that the denominator’s leading term must be fixed, not trainable\) was preserved in evolutionary memory and informed subsequent designs\. The next descendant,TanhBoundedSpinCrossEnhancedWB97MV\(iteration 82\), abandoned the rational exchange form entirely and instead added tanh\-bounded spin descriptors and a new cross descriptorx=z​v/\(1\+z​v\+ϵ\)x=zv/\(1\+zv\+\\epsilon\)to the correlation channels while keeping exchange as a standard polynomial\. All four constraints passed, and the performance from its parent improved \(WRMSDval4\.17→\\to4\.08\) kcal/mol\.

The third descendant,RationalBoundedCrossEnhancedWB97MV\(iteration 106\), reintroduced rational forms, but only to correlation, and with a critical fix: all denominator coefficients were wrapped insoftplus​\(x\)=log⁡\(1\+ex\)\\mathrm\{softplus\}\(x\)=\\log\(1\+e^\{x\}\), guaranteeing non\-negative contributions and a denominatorD≥1D\\geq 1everywhere\.

Iteration 106 plan\.“Replace correlation enhancement factors with bounded rational form\. The numerator retains all existing terms; the denominator uses softplus\-wrapped coefficients so thatD≥1D\\geq 1always, eliminating division by zero risk\. Exchange is unchanged\.”

This produced a WRMSDvalof 3\.81 kcal/mol, the highest score on any single island before the fusion event\. The key design insight, learned from the iteration 18 failure, was that rational forms should be applied to correlation only, with exchange left as a simple polynomial to preserve the UEG constraint algebraically\.

Iteration 197: fusion as the decisive step\.By iteration 197, the two branches were joined through migration across islands: the spin\-exchange lineage’sSpinExEnhancedWB97MVresided on Island 1, and the rational\-correlation lineage’sRationalBoundedCrossEnhancedWB97MVmigrated to from Island 2 to Island 1\. The planner recognized that the two lineages had complementary strengths \(one in exchange, the other in correlation\) and proposed fusing them:

Iteration 197 plan\.“Fusion of the current parent’s spin\-enhanced exchange channel with the top\-performing rational bounded correlation channel from island 1’s best solution, delivering combined improvements from both independent high\-performing lineages\. \[…\] Initializes exactly to the performance level of the top island 1 solution when exchange spin parameters are set to 0, eliminating all regression risk\.”

The fusion producedFusionRationalSpinEnhancedWB97MV\(SAFS26\-a\): exchange from the spin\-exchange lineage \(with the spin\-dependentzxz\_\{x\}descriptor\) combined with correlation from the rational\-correlation lineage \(rational bounded forms with softplus denominators and tanh\-bounded spin cross terms\)\. The result was a WRMSDvalof 3\.65 kcal/mol, a significant score improvement over the direct parent:

Iteration 197 summary\.“Score improved substantially\. \[…\] Fusion of independently optimized, complementary functional components across separate lineages delivers far larger performance gains than incremental descriptor additions to already optimized single channels\.”

The functional remained the global best through the entire remaining 200 iterations of the search \(iterations 198–400\), with no subsequent modification able to improve upon it\.

Agent Trajectory CaveatIteration 97 describes an additional exchange spin\-dependent cross\-term descriptor\. This term is excluded from Section[2\.4](https://arxiv.org/html/2605.05460#S2.SS4)because an oversight led to the effects of the spin becoming nullified\. Because the exchange enhancement factor is evaluated separately on each spin channel’s \(scalar\) density, the spin\-resolution needed to compute\(ζ\)\(\\zeta\)collapses by the time the functional is called, nullifying explicit\(ζ\)\(\\zeta\)\-dependence and effectively yielding \(zx≈vz\_\{x\}\\approx v\)\.

## Appendix EGrid\-Constrained Search: Discovery of a Fully Compliant Functional

The search described in Section[D](https://arxiv.org/html/2605.05460#A4)enforced spin symmetry, UEG exchange limit, uniform coordinate scaling but did not penalize AE18 grid sensitivity during evolution\. A parallel search was conducted with AE18 grid convergence added as a fourth enforced constraint\.

### E\.1Discovery trajectory \(iterations 8–384\)

The direct parent chain ofSAFS26\-blies on Island 0, but the two largest score jumps each drew on solutions from other islands:SigmoidRegimeTunableSwitch\(Island 1, iteration 141, WRMSDval4\.12 kcal/mol\) andEvolveSRMGGA\_ZetaAugmented\_Optimized\(Island 3, iteration 267, WRMSDval3\.94 kcal/mol\)\. One marginal intermediate generation \(iteration 288,rs×vr\_\{s\}\\times vcross terms\) is omitted from the narrative below\.

Iterations 8–80: alpha descriptor, UEG failure, and descriptor partitioning\.The first descendant,WB97MV\_AlphaDescriptor\(iteration 8\), introduced an iso\-orbital indicator descriptorv=1/\(1\+α2\)v=1/\(1\+\\alpha^\{2\}\), \(whereα=\(τ−τW\)/τUEG\\alpha=\(\\tau\-\\tau\_\{W\}\)/\\tau\_\{\\mathrm\{UEG\}\}is the iso\-orbital indicator\) to all three enhancement\-factor channels\. The iso\-orbital indicatorα\\alphais the same quantity whose grid sensitivity problems are well documented in the SCAN functional\[wheeler2010grid\]; however, the functional does not useα\\alphadirectly in interpolation functions with sharp features as SCAN does\. Instead, it mapsα\\alphathrough a smooth inverse decayv=1/\(1\+α2\)v=1/\(1\+\\alpha^\{2\}\), which saturates to zero for largeα\\alphaand is clipped to\[0\.01,1\.0\]\[0\.01,\\,1\.0\]\. This bounded mapping suppresses the effect of large grid\-dependent fluctuations inα\\alphathat occur in low\-density tail regions\.

The validation WRMSDvalimproved relative to the parent seed \(4\.40→\\to4\.34\) kcal/mol, but thevvdescriptor in exchange shifted the enhancement factor away from the UEG limit:

Iteration 8 summary\.“UEG Exchange constraint failed with max\|gx−0\.85\|=1\.06×10−2\|g\_\{x\}\-0\.85\|=1\.06\\times 10^\{\-2\}\. \[…\] Expanding the functional form with new descriptors without explicitly fixing their limiting values for standard physical constraint cases leads to avoidable constraint violations even when raw performance improves\.”

The next descendant,CorrelationOnlyVDescriptor\(iteration 80\), proposed removingvvfrom exchange while retaining it in correlation, directly citing the parent’s failure:

Iteration 80 plan\.“Fixes the parent’s UEG Exchange constraint violation by eliminatingvvterms from the exchange channel, which were the root cause ofgxg\_\{x\}deviating from the required 0\.85 value at the UEG limit\. Retains allvvdescriptor terms in correlation channels, which were proven in the parent iteration to deliver measurable improvement in raw thermochemistry prediction accuracy\.”

All four constraints now passed\. The raw WRMSDvalwas slightly worse \(4\.54 kcal/mol\), but the score improved because the UEG penalty was eliminated\. No subsequent generation reintroduced descriptors to the exchange channel\.

Iteration 192: absorbing the sigmoid regime\-switch architecture from Island 1\.The largest single\-step gain in the lineage came at iteration 192, where the planner absorbed the sigmoid regime\-switch architecture fromSigmoidRegimeTunableSwitch, a solution on Island 1 with WRMSDvalof 4\.12 kcal/mol\. That solution had independently evolved tunable sigmoid switches for exchange and additiveu⋅wu\\cdot wcorrection terms for same\-spin correlation\. The iteration 192 planner combined this architecture withrsr\_\{s\}\-dependent rational correction terms for correlation:

Iteration 192 plan\.“Combining the highest\-performing sigmoid regime switch architecture \(score 0\.8385\) with validatedrsr\_\{s\}\-dependent rational correction terms for correlation channels \[…\] Merges two independently validated performance improvements: tunable exchange regime switches andrsr\_\{s\}\-dependent correlation corrections, targeting complementary error sources \(high\-gradient bonding regions and low\-density correlation effects\) for synergistic improvement\.”

The exchange enhancement factor was augmented with a multiplicative sigmoid correction\(1\+csig⋅σ​\(ax​\(s2−bx\)\)\)\(1\+c\_\{\\mathrm\{sig\}\}\\cdot\\sigma\(a\_\{x\}\(s^\{2\}\-b\_\{x\}\)\)\)and a trainable switch functionfswitch​\(u\)f\_\{\\mathrm\{switch\}\}\(u\)\. Same\-spin correlation received anrsr\_\{s\}\-dependent multiplicative factor\(1\+βss​rs/\(1\+γss​rs2\)\)\(1\+\\beta\_\{\\mathrm\{ss\}\}r\_\{s\}/\(1\+\\gamma\_\{\\mathrm\{ss\}\}r\_\{s\}^\{2\}\)\), whereβss\\beta\_\{\\mathrm\{ss\}\}is a trainable coefficient controlling the linear\-rsr\_\{s\}correction strength andγss=0\.1\\gamma\_\{\\mathrm\{ss\}\}=0\.1is a fixed regularization constant whose sole role is to damp the correction at largersr\_\{s\}, preventing divergence in the low\-density limit \(it is unrelated to the descriptor\-mapping parametersγx\\gamma\_\{x\},γcss\\gamma\_\{\\mathrm\{css\}\},γcos\\gamma\_\{\\mathrm\{cos\}\}\); opposite\-spin correlation received an analogous term\. The score improved substantially \(WRMSDval4\.54→\\to4\.13\) kcal/mol:

Iteration 192 summary\.“Score improved significantly\. \[…\] Combining pre\-validated, complementary architectural modifications targeting separate electronic regimes delivers measurable performance gains without increasing constraint violation risk\.”

Iteration 332: spin\-polarization corrections\.After a marginal intermediate generation that addedrs×vr\_\{s\}\\times vcross terms \(iteration 288\), the iteration 332 planner identified spin polarization as an unexploited structural dimension:

Iteration 332 plan\.“Targets the key unexplored structural gap in the parent architecture: existing correlation terms do not account for spin polarization dependence, which is a known source of error for spin\-polarized systems in thermochemistry benchmarks\. \[…\] Leaves the exchange channel completely unchanged to preserve the parent’s already good UEG constraint performance\.”

Spin\-polarization \(ζ\\zeta\)\-dependent terms were added to the correlation corrections: aζ​\(1−ζ\)\\zeta\(1\-\\zeta\)term for same\-spin channels and an\|ζ\|\|\\zeta\|term for opposite\-spin \(WRMSDval4\.12→\\to4\.04\) kcal/mol:

Iteration 332 summary\.“Score improved\. \[…\] The trainedcss,ζ=−0\.770c\_\{\\mathrm\{ss\},\\zeta\}=\-0\.770indicates that spin\-polarization corrections to same\-spin correlation carry substantial weight\.”

Iteration 384: absorbing zeta polynomial cross terms from Island 3\.The final and second\-largest gain came from absorbing the opposite\-spin zeta polynomial architecture fromEvolveSRMGGA\_ZetaAugmented\_Optimized, a solution on Island 3\. That solution had independently developed an 11\-term polynomial in\(wavg,uavg,z\)\(w\_\{\\mathrm\{avg\}\},u\_\{\\mathrm\{avg\}\},z\)for opposite\-spin correlation, wherez=fζ​\(ζ\)z=f\_\{\\zeta\}\(\\zeta\)is the von Barth–Hedin interpolation function defined above:

Iteration 384 plan\.“Extending the validated parent by merging the zeta polynomial cross terms from the global top\-performingEvolveSRMGGA\_ZetaAugmented\_Optimizedarchitecture exclusively into the opposite\-spin correlation channel, with no changes to exchange or same\-spin correlation channels, plus explicitsig\_coeff\_osclamping\. \[…\] Only modifies the opposite\-spin correlation channel, which is the highest\-impact underoptimized channel per past iteration results\.”

The opposite\-spin polynomial basis was extended to includezzas a third descriptor\. Thesig\_coeff\_osparameter was explicitly clamped to\[0,2\]\[0,2\]\. The resulting functional,WB97MV\_RegimeSwitch\_RSVZetaPolynomial\(SAFS26\-b\), produced the best score in the lineage \(WRMSDval4\.04→\\to3\.78\) kcal/mol:

Iteration 384 summary\.“Score improved\. \[…\] Merging validated features from top\-performing architectures exclusively into specific underoptimized channels delivers incremental performance gains while maintaining full constraint compliance with zero regression risk when new parameters are initialized to zero\.”

Agent Trajectory CaveatIteration 384 describes additional clipped sigmoid term added to the same spin correlation:\(1\+clamp​\(csig,os,0,2\)⋅σ​\(aos⋅savg2\)\)\(1\+\\mathrm\{clamp\}\(c\_\{\\mathrm\{sig,os\}\},0,2\)\\cdot\\sigma\(a\_\{\\mathrm\{os\}\}\\cdot s\_\{\\mathrm\{avg\}\}^\{2\}\)\)\. This term is excluded from Section[2\.5](https://arxiv.org/html/2605.05460#S2.SS5)becausecsig,osc\_\{\\mathrm\{sig,os\}\}optimized below the clipping range, nullifying its effects\.

## Appendix FAblation Study: Dangers of an Unconstrained Search

As noted in Table[1](https://arxiv.org/html/2605.05460#S2.T1), the top unconstrained functional \(SAFS26\-x\) achieves a WRMSDvalcomparable to the constrained search while failing all four physical constraints\. This subsection examineshowthe optimizer exploits unphysical degrees of freedom\.

Inspection of the implementation reveals the specific mechanisms by which the optimizer exploited unphysical degrees of freedom\.

1\. Spin symmetry violation \(3\.5 mHa; tolerance 0\.01 mHa\)\.The opposite\-spin correlation enhancement factor includes an additive polarization term:

`The optimizer converged to alpha\_cos=−0\.174\\texttt\{alpha\\\_cos\}=\-0\.174, which is antisymmetric under spin exchange: swapping ρα↔ρβ\\rho\_\{\\alpha\}\\leftrightarrow\\rho\_\{\\beta\} flips the sign of ζ\\zeta, directly breaking the invariance Exc​\[ρα,ρβ\]=Exc​\[ρβ,ρα\]E\_\{\\mathrm\{xc\}\}\[\\rho\_\{\\alpha\},\\rho\_\{\\beta\}\]=E\_\{\\mathrm\{xc\}\}\[\\rho\_\{\\beta\},\\rho\_\{\\alpha\}\]\. The optimizer retained this term because the antisymmetric correction lowers training loss on open\-shell systems like O2\. 2\. UEG exchange limit violation \(1\.67%; tolerance 0\.5%\)\. The exchange enhancement factor passes the descriptors \(w,u\)\(w,u\) along with a normalized Laplacian through a learned two\-unit sigmoid layer before constructing gxg\_\{x\}: In the uniform electron gas \(s=0s=0, τ=τUEG\\tau=\\tau\_\{\\mathrm\{UEG\}\}\), the descriptors evaluate to w=0w=0, u=0u=0, lapl​\_​norm=0\\mathrm\{lapl\\\_norm\}=0, and the sigmoid units produce σ​\(bx\)≠0\\sigma\(b\_\{x\}\)\\neq 0\. The resulting hidden\-layer residual shifts gxg\_\{x\} away from the required value cx,0=0\.85c\_\{x,0\}=0\.85 by 1\.67%\. In the baseline ω\\omegaB97M\-V, exchange is a simple polynomial that evaluates to cx,0c\_\{x,0\} exactly when w=u=0w=u=0; the neural network breaks this algebraic guarantee\. 3\. Uniform coordinate scaling violation \(12\.2 mHa; tolerance 0\.1 mHa\)\. The exchange network uses the normalized Laplacian as an input feature: Under uniform coordinate scaling 𝐫→𝐫/λ\\mathbf\{r\}\\to\\mathbf\{r\}/\\lambda, the dimensionless descriptors ss and tt are invariant, but lapl​\_​norm\\mathrm\{lapl\\\_norm\} is not: the Laplacian of a molecular density does not scale identically to kF2​ρk\_\{F\}^\{2\}\\rho away from the uniform\-gas limit\. Because the neural network mixes lapl​\_​norm\\mathrm\{lapl\\\_norm\} with \(w,u\)\(w,u\), the entire exchange enhancement factor acquires spurious λ\\lambda\-dependence\. 4\. Grid stability violation \(max \|Δ​Exc\|=0\.100\|\\Delta E\_\{\\mathrm\{xc\}\}\|=0\.100 kcal/mol; toleran1ce 0\.015 kcal/mol\)\. Two architectural choices make the functional acutely sensitive to quadrature resolution, as revealed by comparing the \(99,590\) and \(250,974\) grids on the AE18 atomization set\. First, the exchange enhancement factor depends on a numerical Laplacian: \_laplacian\_from\_drho applies jnp\.gradient \(a finite\-difference stencil\) to each Cartesian component of ∇ρ\\nabla\\rho on the integration grid\. Because the stencil width and point spacing change when the radial or angular grid is refined, lapl​\_​norm\\mathrm\{lapl\\\_norm\} is not a true pointwise density functional; it is an inter\-grid\-point numerical derivative whose value shifts with quadrature resolution\. This shift propagates through the sigmoid network into gxg\_\{x\} and ultimately into the integrated exchange energy\. Second, the same\-spin correlation adds a correction controlled by a steep sigmoid switch in rsr\_\{s\}: The slope of 10 produces a peak derivative of ∼\\sim2\.5 near rs=0\.5r\_\{s\}=0\.5, which falls squarely in the valence\-density regime of AE18 atoms\. Small shifts in ρ\\rho between grids translate into rsr\_\{s\} shifts that are amplified by the switch, toggling the correction on or off in slightly different spatial regions\. Together, the two mechanisms drive the maximum \|Δ​Exc\|\|\\Delta E\_\{\\mathrm\{xc\}\}\| to 0\.100 kcal/mol, nearly seven times the 0\.015 kcal/mol threshold recommended by Refs\. \[mardirossian2016wB97MV, liang2026coach\]\.`

`Appendix G Ablation Study: Creativity as a Critical Ingredient This report argues that genuine improvement in this search required structural exploration rather than incremental exploitation, and credits multi\-island topology and exploration\-heavy parent selection as the important mechanisms for maintaining that diversity \(Section 3\)\. Here, we present the outcome of removing these mechanisms, using only a single island and reverting to a 20/80 exploration split\. The result was a functional that achieved the lowest recorded WRMSDval of any candidate across all runs, yet generalizes catastrophically to the test set\. Functional WRMSDval WRMSDtest Gap ω\\omegaB97M\-V 4\.02 3\.64 −\-0\.38 SAFS26\-a 3\.65 3\.41 −\-0\.24 SAFS26\-b 3\.78 3\.44 −\-0\.34 WB97MVCrossTerm…ActivatedCSS\_v1 3\.64 4\.09 \+0\.45 Table 4: Validation, test, and val–test gap WRMSD in kcal/mol for selected functionals\. WB97MV\_CrossTerm\_Gated\_TauW\_MultiGate\_RatOSC\_ActivatedCSS\_v1 achieves a record\-low validation score while simultaneously producing the widest positive gap of any candidate\. WB97MV\_CrossTerm\_Gated\_TauW\_MultiGate\_RatOSC\_ActivatedCSS\_v1 exhibits the most severe validation–test gap of any candidate in the study\. The reason is clear, supplied with no outside ideas and repeatedly viewing the same small set of high performing solutions, any large structural change is almost certain to degrade performance\. Thus, agents chose ‘safer’ smaller changes\. Among other micro\-exploitations, the most striking was an emergent strategy of systematically freezing high performing parent’s parmeters and grafting a minimal correction on top\. Of the 65 total numeric parameters in the functional, 61 were inherited verbatim from the parent at their fully trained values, and only the 4 elements of c\_css\_act were free to move during L\-BFGS training\. The structural modification multiplied the existing same\-spin enhancement factor by a ratio of two sigmoid activations, a correction which amounts to very little grafted onto an otherwise fixed 61\-parameter form\. Any improvement over the parent’s validation score is thus driven by overfitting to the validation distribution rather than by genuine representational change\. Significance for the creativity claim\. WB97MV\_CrossTerm\_Gated\_TauW\_MultiGate\_RatOSC\_ActivatedCSS\_v1 came from Seed 2\.0, the backend that otherwise produced structurally diverse proposals and the two best\-generalizing functionals in the study\. Its failure therefore cannot be attributed to the intrinsic conservatism of the LLM backend\. Rather, it demonstrates that naively iterating on top\-performing solutions collapses the entire population diversity in a negative feedback loop\. The emergent frozen\-parameter strategy is emblematic of this collapse: it offers a near\-zero\-risk way to guarantee no regression on validation while potentially squeezing out marginal score gains\. Without structural mechanisms that force the search away from the current elite even an otherwise high performing SOTA models still converge to this degenerate regime at sufficient search depth\. class WB97MV\_CrossTerm\_Gated\_TauW\_MultiGate\_RatOSC\_ActivatedCSS\_v1\( SeparatedPxExcFunc\): def \_\_init\_\_\(self, trainable\_params=\("c\_css\_act",\)\) \-\> None: \.\.\. params\.update\(dict\( c\_x\_0 = 0\.85, \# frozen \(UEG constraint\) c\_x\_train = jnp\.array\(\[0\.46467883, 3\.27078968\]\), \# frozen c\_css = jnp\.array\(\[1\.10881579, \-1\.72331564, \-7\.26001241, 16\.62646155, \-1\.15694472\]\), \# frozen c\_cos = jnp\.array\(\[0\.24988053, \-9\.39610915, 6\.10997423, \-9\.37973002, \-15\.28420953, 16\.03908365\]\), \# frozen c\_x\_cross = jnp\.array\(\[\-0\.12446728, \-2\.87756573\]\), \# frozen c\_css\_cross = jnp\.array\(\[ 0\.98494613\]\), \# frozen c\_cos\_cross = jnp\.array\(\[ 1\.40793417,\-12\.8656106 \]\), \# frozen c\_x\_gate = jnp\.array\(\[\-0\.03675481, 0\.50439803\]\), \# frozen c\_css\_gate = jnp\.array\(\[\-2\.29262665\]\), \# frozen c\_cos\_gate = jnp\.array\(\[14\.47816245, 26\.0946234, 0\.45369871\]\), \# frozen c\_x\_tw = jnp\.array\(\[\-0\.26107215, \-0\.12455014\]\), \# frozen c\_css\_tw = jnp\.array\(\[14\.58461232\]\), \# frozen c\_cos\_tw = jnp\.array\(\[\-2\.56354756, 1\.05986868, \-0\.939609 \]\), \# frozen c\_x\_gate\_tw = jnp\.array\(\[ 0\.35455838, 2\.1550685 \]\), \# frozen c\_css\_gate\_tw= jnp\.array\(\[\-0\.25252967\]\), \# frozen c\_cos\_gate\_tw= jnp\.array\(\[ 0\.37069902\]\), \# frozen c\_cos\_tw\_coupled = jnp\.array\(\[1\.97558549, 0\.24267317\]\), \# frozen c\_rat\_num = jnp\.array\(\[\-0\.24245632, 0\.07568424, \# frozen \-0\.58029248, 0\.0038596 , \-0\.13676558, 0\.05632591, 0\.02322058, \-0\.02319411\]\), c\_rat\_den = jnp\.array\(\[\-0\.24487811, \-0\.15191397, \# frozen \-0\.54878029, \-0\.00818291, \-0\.01495359, 0\.01088918\]\), c\_rat\_num\_ext2 = jnp\.array\(\[ 0\.00163673, \-0\.00321903\]\), \# frozen c\_rat\_den\_ext2 = jnp\.array\(\[\-6\.44e\-04, 1\.74e\-05 \]\), \# frozen c\_osc\_elf = jnp\.array\(\[ 0\.82254053, 1\.40289355, \# frozen \-0\.20961353\]\), c\_osc\_kingrad= jnp\.array\(\[ 0\.0, 0\.0, 0\.0\]\), \# frozen \# \-\-\-\- only trainable parameters \-\-\-\- c\_css\_act = jnp\.array\(\[ 0\.0, 0\.0, 0\.0, 0\.0\]\), \# TRAINABLE, zero\-initialized \)\) \.\.\.`