Paper deep dive
Neural-Network Solutions to Real-Space Charge Density and Generalization
Yuxuan Zeng, Taoyuze Lv, Zhicheng Zhong
Intelligence
Status: not_run | Model: - | Prompt: - | Confidence: 0%
Entities (0)
Relation Signals (0)
No relation signals yet.
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:The Hohenberg-Kohn theorem establishes that, in principle, the ground state (GS) charge density contains all GS information of a many-electron system, such that all GS observables can be expressed as functionals of the GS charge density. Conventional Kohn-Sham density functional theory requires iterative solution of the self-consistent-field equations at substantial computational cost, motivating the development of deep learning surrogates for electronic structure calculations and, in turn, accelerating computer-aided materials design. Here, we propose \textbf{AIDEN}, an \underline{A}tomic-\underline{I}nteraction \underline{D}ensity \underline{E}quivariant \underline{N}etwork for solving real-space charge density. AIDEN separates the element-dependent one-center density from environment-induced density redistribution and represents the latter through complementary atom- and edge-centered tensor correlations. A continuous low-rank Gaussian decoder then reconstructs the density at arbitrary spatial coordinates while reusing atomic encodings independently of the evaluation grid. AIDEN achieves state-of-the-art accuracy on periodic crystal benchmarks while remaining competitive for molecular systems, and further demonstrates zero-shot transferability across several structurally distinct out-of-distribution case studies. Furthermore, AIDEN provides substantially faster inference than both baseline models and full SCF calculations, enabling efficient charge density reconstruction for large-scale electronic structure calculations.
Tags
Links
- Source: https://arxiv.org/abs/2609.14906v1
- Canonical: https://arxiv.org/abs/2609.14906v1
Trouble viewing inline? Open PDF directly →
Full Text
127,356 characters extracted from source content.
Expand or collapse full text
Neural-Network Solutions to Real-Space Charge Density and Generalization Yuxuan Zeng Email: mailto:flotsg@mail.ustc.edu.cnflotsg@mail.ustc.edu.cn Affiliation: School of Artificial Intelligence & Data Science, University of Science and Technology of China Taoyuze Lv †thanks: Corresponding author. Email: mailto:taoyuze.lv@ustc.edu.cntaoyuze.lv Affiliation: School of Artificial Intelligence & Data Science, University of Science and Technology of China Zhicheng Zhong †thanks: Corresponding author. Email: mailto:zczhong@ustc.edu.cnzczhong@ustc.edu.cn Affiliation: School of Artificial Intelligence & Data Science, University of Science and Technology of China Affiliation: Suzhou Institute for Advanced Research, University of Science and Technology of China Abstract The Hohenberg-Kohn theorem establishes that, in principle, the ground state (GS) charge density contains all GS information of a many-electron system, such that all GS observables can be expressed as functionals of the GS charge density. Conventional Kohn-Sham density functional theory requires iterative solution of the self-consistent-field equations at substantial computational cost, motivating the development of deep learning surrogates for electronic structure calculations and, in turn, accelerating computer-aided materials design. Here, we propose AIDEN, an Atomic-Interaction Density Equivariant Network for solving real-space charge density. AIDEN separates the element-dependent one-center density from environment-induced density redistribution and represents the latter through complementary atom- and edge-centered tensor correlations. A continuous low-rank Gaussian decoder then reconstructs the density at arbitrary spatial coordinates while reusing atomic encodings independently of the evaluation grid. AIDEN achieves state-of-the-art accuracy on periodic crystal benchmarks while remaining competitive for molecular systems, and further demonstrates zero-shot transferability across several structurally distinct out-of-distribution case studies. Furthermore, AIDEN provides substantially faster inference than both baseline models and full SCF calculations, enabling efficient charge density reconstruction for large-scale electronic structure calculations. 1 Introduction The exploration of the vast materials space is rapidly shifting from conventional trial-and-error experiments and direct first-principles calculations toward data-driven discovery, with deep learning playing an increasingly important role. The data underlying this paradigm are often generated from high-fidelity Kohn-Sham density functional theory (KS-DFT) calculations (Jones, 2015) or molecular dynamics (MD) simulations (van Gunsteren & Mark, 1998), which can be performed systematically at scale and generally provide greater internal consistency than experimental measurements. Deep-learning approaches to DFT calculations commonly target the charge density (Qin et al., 2026), Hamiltonian (Tang et al., 2024), and density matrix (Dong et al., 2025), which are connected through the self-consistent field (SCF) cycle (Slater, 1969): ρ()→HKS[ρ]→ϕn,ϵn→ρ(). ρ(r)→ H_ KS[ρ]→ \ _nk, _nk \ →ρ(r). (1) Here, ρ()ρ(r) denotes the real-space charge density, HKS[ρ]H_ KS[ρ] the KS Hamiltonian constructed from the density-dependent effective potential, and D the one-particle density matrix constructed from the KS orbitals and their occupations. Compared with learning the Hamiltonian or density matrix, charge density prediction offers several advantages: 1) Weaker dependence on atomic-orbital basis representations. Changing the atomic-orbital basis alters the matrix representations of both the Hamiltonian and density matrix (Yuan et al., 2026), whereas ρ()ρ(r) is intrinsically a real-space scalar field independent of orbital indexing, facilitating a more unified representation across different basis sets. 2) Direct physical interpretability. The real-space charge density provides direct access to electronic-structure characteristics such as charge transfer (Poli et al., 2020) and chemical bonding (Chopra, 2012), offering a microscopic basis for interpreting material properties. 3) Convenient three-dimensional (3D) representation. When discretized on a real-space grid, ρ()ρ(r) forms a regular 3D tensor and can serve as a complementary electronic-structure descriptor for downstream representation learning (Shuang et al., 2026). To this end, we propose an Atomic-Interaction Density Equivariant Network (AIDEN) to efficiently predict real-space GS charge densities directly from atomic configurations. Inspired by the superposition of atomic densities (SAD) (Van Lenthe et al., 2006), AIDEN decomposes the total charge density into an atom-centered term ρinit() _ init(r) and an environment-dependent term ρenv() _ env(r). The former serves as a learnable element-dependent one-center baseline for the dominant local density distribution, whereas the latter captures environment-induced density variations through higher-order many-body features. In the encoder, AIDEN first employs Cartesian atomic cluster expansion (ACE) (Xu et al., 2026a) to encode higher-order geometric and angular information of local atomic environments, followed by tensor edge cluster expansion (TECE) (Xu et al., 2026b) to equivariantly model and adaptively aggregate direction-dependent neighbor interactions, providing expressive representations of complex local electronic environments. The decoder then expands the atomic features using Gaussian-type orbitals (GTOs) (Huzinaga, 1965) and reconstructs a continuous real-space charge-density field. In addition, the computationally expensive equivariant atomic encoding in AIDEN is performed only once and reused across all spatial query points, yielding faster inference than conventional probe-based methods. Together, these designs provide AIDEN with strong learning capacity and substantially improved efficiency, while enabling its extension to larger-scale out-of-distribution (OOD) systems. 2 Related works Equivariant GNNs. Conventional atomistic GNNs mainly learn rotation- and translation-invariant scalar representations (Xie & Grossman, 2018), whereas equivariant GNNs enforce f(D(g)x)=D(g)f(x)f(D_X(g)x)=D_Y(g)f(x) to preserve geometric transformation laws. Tensor Field Networks (Zaccone, 2026) and e3n (Geiger & Smidt, 2022) established E(3)-equivariant representations based on SO(3) irreducible representations and spherical harmonics. NequIP (Batzner et al., 2022) demonstrated their effectiveness for atomistic modeling, while MACE (Batatia et al., 2022) incorporated atomic cluster expansion to capture higher-order many-body interactions. More recent methods improve efficiency through edge-aligned SO(2) operations, including eSCN (Passaro & Zitnick, 2023) and TECE (Xu et al., 2026b), the latter further introducing edge cluster expansion and radial rotary attention. Charge density learning. Recent methods predict real-space charge densities using spatial queries, regular grids, or atom-centered representations. DeepDFT (Jørgensen & Bhowmik, 2022) predicts densities at arbitrary probe points from local atomic environments, while Deep Charge (Lv et al., 2023) learns symmetry-preserving local representations with strong data efficiency. ChargE3Net (Koker et al., 2024) extends probe-based prediction with higher-order E(3)-equivariant features and scales to more than 10510^5 crystals. Li et al. (2025) instead reconstructs molecular densities on regular 3D grids, whereas EAC-Net (Qin et al., 2026) represents the density as a sum of symmetry-consistent atomic contributions. NeuralSCF (Song & Feng, 2026) further learns the Kohn-Sham density map through neural self-consistent iterations. 3 Problem statement Consider a unit cell Ω=∣∈[0,1)3 =\uA ∈[0,1)^3\ with lattice vectors arranged by rows in ∈ℝ3×3A ^3× 3. A periodic structure is written as =(,Zi,ii=1Na)X=(A,\Z_i,s_i\_i=1^N_ a), where ZiZ_i and i∈[0,1)3s_i∈[0,1)^3 are the atomic number and fractional coordinate of atom i, respectively, and i=iR_i=s_iA is its Cartesian position. Its periodic images lie at i+R_i+nA for ∈ℤ3n ^3. For a cell containing NeN_ e electrons, the physical GS density belongs to Ne=ρ:Ω→ℝ+|∫Ωρ()d3=Ne,ρ()=Ne∫∏a=2Ned3a|Ψ(,2,…,Ne)|2, _N_ e= \ρ: _+\ |\ _ ρ(r)\, d^3r=N_ e \, _X(r)=N_ e _a=2^N_ ed^3r_a | (r,r_2,…,r_N_ e) |^2, (2) and KS-DFT obtains ρ⋆=argminρ∈NeE[ρ] _X = _ρ _N_ eE_X[ρ]. Plane-wave calculations (Dunnington & Schmidt, 2012) represent the real-space density numerically on a regular fast Fourier transform (FFT) grid (Ten Eyck, 1973), ℛ=gg=1NgR_X=\r_g\_g=1^N_ g. Below, ρ denotes the smooth PAW grid density used as the learning target, rather than the formal all-electron density in Eq. 2; NeN_ e denotes the corresponding valence-electron count, with Ne≃|Ω|Ng−1∑gρ(g)N_ e | |N_ g^-1 _g _X(r_g). The separately supplied PAW one-center data and the scope of fixed-density evaluation are detailed in Appx. A.2. AIDEN learns a continuous scalar field ρ^(,) ρ_ θ(r;X). For a Euclidean transformation g=(,)g=(Q,t) with ∈O(3)Q∈ O(3) and g=⊤+gr=rQ +t, the target covariance and lattice periodicity are ρ^(g,g)=ρ^(,) ρ_ θ(gr;gX)= ρ_ θ(r;X) and ρ^(+,)=ρ^(,) ρ_ θ(r+nA;X)= ρ_ θ(r;X). Given training samples nn=1Nd\X_n\_n=1^N_ d with full FFT grids ℛn=ngg=1Ng(n)R_X_n=\r_ng\_g=1^N_ g^(n), the model minimizes the grid-averaged absolute deviation to determine the optimal model parameters θ⋆θ , ⋆=argmin1Nd∑n=1Nd1Ng(n)∑g=1Ng(n)|ρ^(ng,n)−ρn(ng)|. θ = _ θ 1N_ d _n=1^N_ d 1N_ g^(n) _g=1^N_ g^(n) | ρ_ θ(r_ng;X_n)- _X_n(r_ng) |. (3) 4 Model architecture Figure 1: Overview of AIDEN. (a) Using a simple Na-Cl atom pair as an example, the initial local charge density ρinit(,NaCl) _ init(r;NaCl) can be decomposed into separate contributions from the two atomic species, while the atomic system is represented as a periodic graph containing elemental embeddings, interatomic distances, and edge-direction information. (b) The geometry-informed embedding (GIE) module initializes scalar and higher-order equivariant atomic features through coordinated scalar and tensor branches. (c) In the encoder, Cartesian ACE constructs many-body correlations after neighborhood aggregation, whereas TECE builds higher-order source-target correlations in edge-aligned local coordinate frames and further modulates the correlated information through radial rotary attention (RRA). (d) In the decoder, the charge density is composed of a local term ρinit() _ init(r) and an environment term ρenv() _ env(r); the equivariant atomic representations are mapped to local GTO coefficients and can be efficiently evaluated at arbitrary spatial coordinates. 4.1 Periodic geometry and equivariant initialization For every target atom i, the periodic atomic graph contains directed edges e=(j,→i)e=(j,n\!→ i) satisfying e=i−(j+),re=‖e‖2,^e=e/re,0<re<rat. _e=R_i-(R_j+nA), r_e=\|d_e\|_2, d_e=d_e/r_e, 0<r_e<r_ at. (4) If more than NnbrN_ nbr images enter the cutoff sphere, only the nearest NnbrN_ nbr are retained for that target. We write iN_i for the resulting incoming edges. Each edge carries a radial vector (re)∈ℝNBb(r_e) ^N_ B and an irreducible rank-ℓ Cartesian harmonic (ℓ)(^e)Y^( )( d_e). Their exact construction is given in Appx. B.1. The element input is i=Norm[(Zi)⊕(Zi)]∈[0,1]133x_i= Norm[o(Z_i) (Z_i)]∈[0,1]^133, where (Zi)∈0,1118o(Z_i)∈\0,1\^118 is a one-hot element vector, (Zi)∈ℝ15q(Z_i) ^15 contains the elemental descriptors, ⊕ denotes concatenation, and Norm Norm is feature-wise min-max normalization over the 118 elements. GIE initializes the scalar and higher-order tensor fields directly from the one-hop atomic environment. We first summarize the radial environment of atom i by ¯i b_i =1|i|∑e∈i(re),i(0,0)=ℱs[ℱ0(i),i⊕¯i], = 1 |N_i| _e _ib (r_e ), _i (0,0 )=F_ s [F_0 (x_i ),x_i b_i ], (5) i(0,ℓ) _i (0, ) =1|i|ℓirr[∑e∈iℱℓ[(re)]⊙ℱZ(j)]⊗(ℓ)(^e),1≤ℓ≤L. = 1 |N_i|P_ irr \ [ _e _iF_ [b (r_e ) ] _ Z (x_j ) ] ( ) ( d_e ) \, 1≤ ≤ L. (6) For the scalar field, ℱ0(i)F_0 (x_i ) provides the elemental seed of the central atom, while ℱs(⋅)F_ s (· ) applies a FiLM-style environment-conditioned (Perez et al., 2018) scale and shift determined from i⊕¯ix_i b_i. For ℓ≥1 ≥ 1, ℱℓ[(re)]∈ℝCF_ [b (r_e )] ^C assigns order-specific radial amplitudes, controlling how strongly a neighbor at distance rer_e contributes to each channel, whereas ℱZ(j)∈ℝCF_ Z (x_j ) ^C supplies elemental amplitudes that modulate these contributions according to the chemical identity of the source atom. Their channel-wise product therefore determines the magnitude of each edge contribution, while the analytic harmonic (ℓ)(^e)Y ( ) ( d_e ) fixes its angular dependence and transformation law. Here ⊙ acts over the C scalar channels, and ⊗ broadcasts each channel weight over the Cartesian components of the rank-ℓ harmonic. The projection ℓirrP_ irr keeps the accumulated tensor in the symmetric traceless rank-ℓ subspace. Consequently, i(0,ℓ)∈ℝC×3ℓh_i (0, ) ^C× 3 forms an equivariant Cartesian field in which the learned maps determine the channel-wise strength of neighboring contributions without independently parameterizing individual Cartesian components, thereby preserving rotational equivariance. 4.2 Equivariant cluster expansion encoder Before each interaction, every angular order is normalized as ¯i(t,ℓ)=ℓ[i(t,ℓ)] h_i^(t, )=V_ [h_i^(t, )] (Appx. B.2). In our case, one Cartesian ACE interaction forms higher-order correlations after neighborhood aggregation, describing collective atom-centered coordination. One subsequent TECE interaction forms source-target correlations within each edge-aligned frame before aggregation, retaining directional information along individual interatomic axes. These are complementary geometric inductive biases for density redistribution. 4.2.1 Cartesian atomic cluster expansion The Cartesian ACE follows the irreducible tensor construction of TACE (Xu et al., 2026a). Let LT_L be the admissible angular paths defined in Eq. 51, and let ℓ1,ℓ2ℓC_ _1, _2 be the irreducible Cartesian coupling defined in Eq. 50. The edge-resolved one-particle tensor is e(ℓ)=∑(ℓ1,ℓ2,ℓ)∈Lℓ1,ℓ2ℓ(e,ℓ1ℓ2ℓ⊙ℒs(ℓ1)¯j(0,ℓ1),(ℓ2)(^e)),e,ℓ1ℓ2ℓ=ℱACEℓ1ℓ2ℓ[(re)], _e^( )= _( _1, _2, ) _LC_ _1, _2 (w_e, _1 _2 _ s^( _1) h_j^(0, _1),Y^( _2)( d_e) ), _e, _1 _2 =F_ ACE _1 _2 [b(r_e) ], (7) where ℒs(ℓ)L_ s^( ) mixes channels at fixed ℓ and ℱACEℓ1ℓ2ℓF_ ACE _1 _2 produces a C-component radial weight. The invariant neighbor coefficient aea_e and the explicit construction of the order-ν Cartesian correlation polynomial ℬν(ℓ)B_ν^( ) are detailed in Appx. B.2. Define the aggregated atomic field i _i as i _i =ℒself¯i(0)+ℒmsg∑e∈iaee, =L_ self h_i^(0)+L_ msg _e _ia_ep_e, i(1,ℓ) _i^(1, ) =¯i(0,ℓ)+ℒout(ℓ)ℬν(ℓ)[C+σ0ℱg(i(0))]⊙i, = h_i^(0, )+L_ out^( )B_ν^( ) \ [1_C+ _0F_ g ( _i^(0) ) ] _i \, (8) where σ0(x)=(1+e−x)−1 _0(x)=(1+ e^-x)^-1 is the sigmoid function, and ℱgF_ g maps the invariant sector i(0) _i^(0) to channel-wise gating coefficients. The symbol C1_C is the all-ones channel vector and is broadcast over every angular component. Bold symbols without an angular superscript denote the direct sum over ℓ=0,…,L =0,…,L, and every ℒL acts only on channels of equal angular order. 4.2.2 Tensor edge cluster expansion and radial rotary attention We first transform the irreducible Cartesian features into a real-spherical representation and then rotate them into an edge-aligned local frame. Let U be the fixed orthogonal map from irreducible Cartesian tensors to real-spherical components, and let eD_e align the local y axis with ^e d_e. With M denoting the maximum retained magnetic order, restricting to |m|≤M|m|≤ M gives DM=(L+1)+2∑m=1M(L+1−m)D_M=(L+1)+2 _m=1^M(L+1-m) real components. The radial operator e[(re)] _e[b(r_e)] and the second-order local edge correlator ℰE are detailed in Appx. B.3. The complete TECE update is the single composition i(2) _i^(2) =12¯i(1)+†∑e∈ie†ℰ[e[(re)]⊙(e¯j(1)⊕e¯i(1))]e, = 1 2 \ h_i^(1)+U _e _iD_e \,E [ _e[b(r_e)] (D_eU h_j^(1) _eU h_i^(1) ) ]A_e \, (9) i⋆(ℓ) _i ( ) =ℓ[i(2,ℓ)]. =V_ [h_i^(2, )]. (10) Here the inverse transforms e†D_e and †U return the correlated edge features to the global Cartesian representation. The product ⊙ in Eq. 9 acts elementwise over the compact components and channels. The H attention heads, each containing Ch=Ce/HC_h=C_ e/H edge channels, act through e=diag(ae1Ch,…,aeHCh)A_e= diag(a_e1I_C_h,…,a_eHI_C_h). The query e,ℓmh∈ℂChq_e, mh ^C_h and key e,ℓmh∈ℂChk_e, mh ^C_h are channel projections of the unmodulated local target and source tensors, respectively; their m=0m=0 components are real. For ,∈ℂChu,v ^C_h, the head-wise Hermitian contraction is ⟨,⟩ch=∑c=1Chuc∗vc ,v _ ch= _c=1^C_hu_c^*v_c. Radial rotary attention assigns ξeh=τhDMCh(∑ℓ=0L⟨e,ℓ0h,e,ℓ0h⟩ch+∑m=1M∑ℓ=mLRe⟨e,ℓmh,eimφehe,ℓmh⟩ch)+βeh, _eh= _h D_MC_h ( _ =0^L _e, 0h,k_e, 0h _ ch+ _m=1^M _ =m^L Re _e, mh, e im _ehk_e, mh _ ch )+ _eh, (11) where βeh∈ℝ _eh and φeh∈(−π,π) _eh∈(-π,π) are radial bias and phase, while τh>0 _h>0 is a learned inverse temperature controlling the sharpness of the attention distribution for head h. The normalized coefficient aeha_eh is a cutoff-weighted softmax over incoming edges, as detailed in Appx. B.3. The score consists only of invariant inner products of equal local frequencies; the radial phase adjusts their relative alignment without changing the SO(2) transformation law. 4.3 Continuous density decoder The decoder maps i⋆(ℓ)h_i ( ) to rank-ℓ coefficient tensors ikpχ(ℓ)=ℒkpχ(ℓ)i⋆(ℓ)+δℓ0μkpχc_ikp^χ( )=L_kp^χ( )h_i ( )+ _ 0 _kp^χ, where χ∈L,Rχ∈\L,R\, k=1,…,Kk=1,…,K, and p=1,…,NGp=1,…,N_ G. Here, ℒkpχ(ℓ)L_kp^χ( ) is an equivariant channel map, δℓ0 _ 0 is the Kronecker delta, and μkpχ _kp^χ is a scalar bias restricted to ℓ=0 =0. The quantities K and NGN_ G denote the environmental rank and the number of Gaussian radial functions, respectively. This Cartesian formulation is equivalent to expressing both the coefficients and harmonics in the same orthogonal real-spherical basis. For a query position r, we define the supported periodic atom images as ℳ()=η=(i,)|η=−(i+),rη=∥η∥2<rorb,^η=η/rη,rη>0,,rη=0. (r)= \η=(i,n)\ |\ δ_η=r-(R_i+nA),\ r_η=\| δ_η\|_2<r_ orb \, δ_η= cases δ_η/r_η,&r_η>0,\\ 0,&r_η=0. cases (12) At an atomic center, we adopt the continuous convention (0)()=1Y^(0)(0)=1 and (ℓ)()=Y^( )(0)=0 for ℓ>0 >0. The radial functions R~ℓp(r) R_ p(r) follow a corrected even-tempered Gaussian construction, while υZik0pχ _Z_ik_0p^χ denotes an element-indexed learnable coefficient for one-center rank k0=1,…,K0k_0=1,…,K_0. Their definitions and the detailed Gaussian basis construction are given in Appx. B.4. Let ⟨⋅,⋅⟩F ·,· _ F denote contraction over all Cartesian tensor indices. The element-only contribution associated with image η and the environment-dependent field at a query position are then ϕηk0χ()=∑pυZik0pχR~0p(rη),Φkχ()=∑η∈ℳ()∑ℓ,pR~ℓp(rη)3ℓ/2⟨ikpχ(ℓ),(ℓ)(^η)⟩F. _η k_0^χ(r)= _p _Z_ik_0p^χ R_0p(r_η), _k^χ(r)= _η (r) _ ,p R_ p(r_η)3 /2 _ikp^χ( ),Y^( ) ( δ_η ) _ F. (13) The continuous density is reconstructed as ρ^(,)=∑η∈ℳ()∑k0ϕηk0L()ϕηk0R()⏟ρinit()+∑kΦkL()ΦkR()⏟ρenv(). ρ_ θ (r;X )= _η (r ) _k_0 _η k_0^L (r ) _η k_0^R (r )_ _ init (r )+ _k _k^L (r ) _k^R (r )_ _ env (r ). (14) Here ρinit _ init is an element-dependent one-center density independent of the encoded environment, whereas ρenv _ env depends on i⋆(ℓ)h_i ( ). In ρenv _ env, the left and right fields are accumulated over ℳ()M(r) before multiplication, naturally introducing cross terms between distinct atom images without explicit pair enumeration. Both terms are expanded in atom-centered GTOs. In our case, we first perform a short pretraining stage dedicated to initializing the representation network for ρinit _ init. During subsequent full training, ρinit _ init and ρenv _ env are jointly optimized. The complete forward process is provided in Appx. B.5. 5 Experiments 5.1 Dataset settings Although AIDEN is primarily designed for periodic crystals, we further evaluate its applicability to both crystalline and molecular systems. For crystals, we use the ECD dataset (Chen et al., 2025), which contains 140,646 inorganic structures spanning 94 elements. The reference charge densities are mainly computed with VASP (Wang et al., 2021) using PAW-PBE (Perdew et al., 1996), with PBE+U (Li et al., 2020) for selected strongly correlated systems and a plane-wave cutoff of 520 eV; an additional 7,147 samples are recalculated with HSE06. We train AIDEN from scratch on the PBE subset and fine-tune it on the smaller HSE subset. For molecular systems, we use the QM9 charge density dataset (Jørgensen & Bhowmik, 2022), containing 133,885 molecules composed of H, C, N, O, and F, with densities computed using VASP and PAW-PBE at a 400 eV cutoff and Γ -point sampling. These datasets provide complementary periodic and molecular regimes for evaluating generalization across chemical spaces, structural scales, and boundary conditions. For OOD evaluation, we additionally consider liquid water, disordered AlMg, amorphous silicon (a-Si), and twisted bilayer graphene (TBG). The first three systems contain 192, 108, and 64 atoms, respectively, while the largest of three TBG structures contains N=148N=148 atoms. 5.2 Evaluation on benchmark datasets Real-space charge density is inherently a continuous physical quantity. For computational convenience, a common coarse-graining strategy is to sample it on a 3D grid, which typically contains millions to tens of millions of grid points. In practice, however, only a tiny fraction of the available density grid can be sufficient for accurate model fitting, particularly when non-uniform or targeted sampling strategies are employed (Jørgensen & Bhowmik, 2022; Focassio et al., 2023). In our case, we uniformly sample 5,000 grid points from each structure for training and quantify the error using the normalized mean absolute error (NMAE), ερ:=1N∑n=1N1Ne(n)∫Ωnd3|ρ^(,n)−ρ(,n)|. _ρ:= 1N _n=1^N 1N_ e(X_n) _ _n d^3r | ρ(r;X_n)-ρ(r;X_n) |. (15) As shown in Tab. 1, AIDEN achieves SOTA performance under both exchange-correlation functionals on the general periodic-crystal benchmark, outperforming existing methods. On the QM9 molecular dataset, AIDEN ranks second, surpassed only by BOA, which is specifically designed for molecular systems. The lower accuracy on HSE than on PBE is expected. First, the HSE training set is substantially smaller, and the ECD benchmark likewise shows that charge-density prediction error decreases as the proportion of HSE data increases (Chen et al., 2025). Second, PBE is a semilocal GGA functional whose exchange-correlation energy mainly depends on the local charge density and its gradient (Perdew et al., 1996), whereas HSE additionally incorporates screened Hartree-Fock exact exchange (Heyd et al., 2003; Krukau et al., 2006). This makes the density more sensitive to orbital localization and complex electronic environments, increasing the difficulty of learning the structure-to-HSE-density mapping. Fig. 2(a) shows that, at the same L, AIDEN provides substantially faster inference than ChargE3Net and BOA. This advantage further increases with system size as NaN_ a grows (see Fig. 10(d)). Fig. 2(b) compares the dependence of parameter count on L. ChargE3Net actively reduces the number of channels assigned to each irreducible representation as L increases (Koker et al., 2024), causing its total parameter count to decrease; however, whether this allocation is optimal for higher-order representation capacity remains unclear. In contrast, BOA maintains an almost constant parameter count by using a fixed feature-channel dimension while encoding angular dependence in non-parametric GTOs (Klockow et al., 2026). As shown in Fig. 2(c), PBE pretraining yields a modest improvement, suggesting that the two functionals share the dominant density structure, while the remaining discrepancy is more closely related to HSE-specific exchange-correlation corrections. The lower performance of AIDEN than BOA on QM9 may partly result from the stronger task-specific priors used by BOA for molecular systems. BOA employs an element-specific def2-QZVPPD basis (Hellweg & Rappoport, 2015) optimized for H/C/N/O/F, together with element-specific radial corrections, and performs message passing directly in the density basis through basis overlaps. Table 1.: Benchmarking charge density learning across different atomic systems and functionals. All reported metrics are quantified by NMAE ερ _ρ (%) (Eq. 15), with lower values indicating better model performance. The SOTA result is highlighted in bold, while the second-best result is underlined. The reference baselines include eqDeepDFT (Jørgensen & Bhowmik, 2022), ChargE3Net (Koker et al., 2024), NeuralSCF (Song & Feng, 2026), SCDP (Fu et al., 2024), ELECTRA (Elsborg et al., 2025), and BOA (Klockow et al., 2026). For the PBE task, two ερ _ρ values are reported for both ChargE3Net and AIDEN, with the left and right entries corresponding to maximum angular order L=3L=3 and L=4L=4. For the HSE task, the two AIDEN results are obtained using the training-from-scratch and fine-tuning strategies. Dataset eqDeepDFT ChargE3Net NeuralSCF SCDP ELECTRA BOA AIDEN PBE 0.799 0.685 / 0.523 - - - - 0.621 / 0.4542 HSE - 1.534 - - - - 1.4307 / 1.3012 QM9 0.284 0.196 0.197 0.178 0.177 0.1339 0.1628 Figure 2: Model efficiency, parameter count, and HSE transfer learning. (a) Per-sample inference time of AIDEN, ChargE3Net, and BOA at different maximum angular orders L. Bar heights indicate the mean values, and the error bars denote one standard deviation. (b) Comparison of model parameter counts. (c) Comparison of the zero-shot inference error of AIDEN and the loss curves under two transfer learning strategies. Opaque curves correspond to validation losses, while semi-transparent curves show training losses. 5.3 Generalization to OOD systems We evaluate zero-shot cross-system transfer of the PBE-pretrained model on these OOD systems, following the system choices of Ref. (Qin et al., 2026), without any system-specific fine-tuning. We compare charge densities and fixed-density energies and forces, and further evaluate non-self-consistent field (NSCF) (Lim & Whitehead, 1967) bands for TBG. For fixed-density property calculations, the predicted grid densities are combined with matched PAW one-center data from the SCF reference, as detailed in Appx. A.2. Figure 3: Visualization of atom-rich charge density cross sections. Rows correspond to Water, AlMg, and a-Si. The first two columns show the DFT reference and AIDEN-inferred ρ()ρ(r); the remaining columns show the normalized pointwise deviation δρ(g):=|ρ^(g)−ρ(g)|/ρ(g) _ρ(r_g):=| ρ(r_g)-ρ(r_g)|/ρ(r_g) of AIDEN, ChargE3Net, and SAD at grid point g. The density and deviation values are indicated by the upper and lower color bars, respectively. All atom-rich cross sections are perpendicular to the x axis. Atoms within 1.5 grid spacings of each plane are counted, yielding H7_7O2_2, Al6_6Mg7_7, and Si5_5 for Water, AlMg, and a-Si, respectively. Detailed settings are provided in Appx. C.4.3. Table 2.: OOD benchmark. Errors are relative to the matched PBE-SCF reference. The best and second-best errors are bold and underlined. Timing ranks compare the two models; SAD times († ) include the complete fixed-density VASP calculation. Here, we define the per-atom energy error as |ΔE|:=|E^−E|/Na| E|:=| E-E|/N_ a, the mean absolute force component as ⟨|Fiα|⟩:=13Na∑i=1Na∑α=x,y,z|Fiα| |F_iα| := 13N_ a _i=1^N_ a _α=\x,y,z\|F_iα|, and the corresponding mean error as ⟨|ΔFiα|⟩:=⟨|F^iα−Fiα|⟩i,α | F_iα| := | F_iα-F_iα| _i,α. System Method ερ _ρ [%] ↓ Time [s] ↓ |ΔE|| E| [meV/atom] ↓ ⟨|ΔFiα|⟩ | F_iα| [eV/Å] ↓ Water AIDEN 1.8764 23.70 3.444 0.11467 ChargE3Net 1.7880 1286.43 6.443 0.12538 SAD 12.3527 1433†1433 485.445 1.11744 AlMg AIDEN 4.3004 13.40 262.809 0.10398 ChargE3Net 4.3274 936.43 276.070 0.03580 SAD 4.8057 1344†1344 49.643 0.15775 a-Si AIDEN 1.3215 4.34 3.213 0.04087 ChargE3Net 1.6391 328.87 4.774 0.04090 SAD 9.9995 307†307 82.107 0.12778 Water, AlMg and a-Si. We evaluate one configuration of each system: a water cell containing 64 molecules, an Al54Mg54 alloy, and a 64-atom a-Si structure. Reference densities and properties are recomputed using PBE-SCF with a 520 eV cutoff; candidate densities are evaluated on the complete grids and held fixed for property calculations (Appx. C.4.2). As shown in Fig. 3 and Tab. 2, AIDEN gives the lowest full-grid ερ _ρ on AlMg and a-Si, while ChargE3Net performs slightly better on water. Both models substantially improve upon SAD on water and a-Si. AIDEN evaluates the full grids in 4.34-23.70 s, a measured 54-76× speedup over the ChargE3Net implementation used here. Its energy errors on water and a-Si are 3.444 and 3.213 meV/atom, compared with 6.443 and 4.774 meV/atom for ChargE3Net; force errors are smaller on water and nearly tied on a-Si. By contrast, AlMg exhibits a different trend. Although both models slightly improve the density NMAE over SAD, their energy errors are substantially larger, while ChargE3Net yields a smaller force error than AIDEN, 0.03580 versus 0.10398 eV/Å, despite nearly identical density NMAEs. Both models overestimate the Mg-centered density peaks by approximately 20%, with the dominant errors localized in the near-core region containing the 2p2p semicore contribution. Interestingly, similar behavior has been reported for learned densities in Al and Si, where small global density errors can produce much larger electrostatic-energy errors because near-nuclear inaccuracies strongly affect the electron-nuclear interaction energy (Lewis et al., 2021). This near-core sensitivity provides a natural explanation for the mismatch between global NMAE and energetic accuracy in AlMg; detailed analyses are provided in Appx. C.4.5. The result therefore highlights the need to evaluate learned densities jointly through density errors and downstream observables. Figure 4: NSCF band structures of TBG. (a) 21.79∘21.79 , Na=28N_ a=28; (b) 13.17∘13.17 , Na=76N_ a=76; and (c) 9.43∘9.43 , Na=148N_ a=148. Yellow solid, blue dashed, and purple dashed curves denote DFT, AIDEN, and ChargE3Net. (d), (e), and (f) show the regions with pronounced errors produced by ChargE3Net in the three TBG band structures, corresponding to the dashed boxes in (a)-(c). (g) Error comparison of TBG band structures computed from the charge densities predicted by AIDEN and ChargE3Net against the DFT reference. NSCF band structure calculation of TBG. We further evaluate TBG at 21.79∘21.79 , 13.17∘13.17 , and 9.43∘9.43 . Fixed-density PAW-PBE calculations use ICHARG=11, a 520 eV cutoff, and the Γ –M–K–Γ path. Fig. 4 shows bands with the DFT EFermiE_ Fermi set to zero and one least-squares rigid shift applied per model and cell. Within EFermi±5E_ Fermi± 5 eV, the aligned band MAEs are 8.42, 7.41, and 11.47 meV for AIDEN, versus 27.35, 10.30, and 12.93 meV for ChargE3Net. These aligned MAEs quantify residual band-shape and dispersion errors after removal of a global reference-energy offset. Full-grid inference takes 13.27-41.43 s, a measured 9.7-15.3× speedup over the evaluated ChargE3Net implementation. Appx. C.3 reports the shifts, raw and tail-error statistics, sampling differences, and timing protocol. 6 Discussion In this work, we introduced AIDEN for accurate and scalable real-space charge density learning. At its core, AIDEN adopts a physically inspired, learnable decomposition that separates an element-dependent one-center baseline from environment-induced charge redistribution. The former captures the dominant local density structure shared across different environments, while the latter describes bonding- and coordination-dependent corrections, providing a natural inductive bias across diverse chemical and structural spaces. This decomposition is integrated with an equivariant encoder and a continuous decoder to preserve the required geometric symmetries while enabling efficient field evaluation. Experiments demonstrate high accuracy on periodic-crystal and molecular benchmarks, effective PBE-to-HSE fine-tuning, and zero-shot transfer to structurally distinct OOD systems. Fixed-density calculations further assess the predicted densities through energies and forces, as well as band structures in TBG. Full-grid inference shows favorable scalability to millions of grid points. Notably, TBG electron counts deviate by only 0.0068-0.0095% without explicit normalization (Appx. C.3), while the AlMg results highlight the complementarity between density NMAE and downstream energy/force errors. Future perspective. AIDEN currently focuses on GS scalar charge densities. Future work may extend it to spin-resolved densities and broader electronic-structure settings, including unified models across exchange-correlation functionals, pseudopotentials, and chemical spaces. Joint learning of charge densities with energies, forces, and other DFT observables may further improve electronic-structure fidelity and generalization. 7 Code availability The code supporting this work is available in the AIDEN repository. The PBE and HSE datasets used for model training are provided by ECDBench, the QM9 dataset is available from the QM9 charge density dataset. 8 Acknowledgements We would like to acknowledge the National Key R&D Program of China (No. 2021YFA0718900), the National Natural Science Foundation of China (Nos. 12374096 and 92477114), and the Jiangsu Funding Program for Excellent Postdoctoral Talent for financial support. Z. Zhong thanks the Suzhou Innovation and Entrepreneurship Leading Talent Program and the Gusu Leadership Program for their support. References Batatia et al. (2022) Ilyes Batatia, David P Kovacs, Gregor Simm, Christoph Ortner, and Gabor Csanyi. MACE: Higher order equivariant message passing neural networks for fast and accurate force fields. In S. Koyejo, S. Mohamed, A. Agarwal, D. Belgrave, K. Cho, and A. Oh (eds.), Adv. Neural Inf. Process. Syst., volume 35, p. 11423–11436, LA, USA, 2022. Curran Associates, Inc. doi: 10.52202/068431-0830. Batzner et al. (2022) Simon Batzner, Albert Musaelian, Lixin Sun, Mario Geiger, Jonathan P Mailoa, Mordechai Kornbluth, Nicola Molinari, Tess E Smidt, and Boris Kozinsky. E (3)-equivariant graph neural networks for data-efficient and accurate interatomic potentials. Nat. Commun., 13(1):2453, 2022. doi: 10.1038/s41467-022-29939-5. Chen et al. (2025) Pin Chen, Zexin Xu, Qing Mo, Hongjin Zhong, Fengyang Xu, and Yutong Lu. ECD: A machine learning benchmark for predicting enhanced-precision electronic charge density in crystalline inorganic materials. In Int. Conf. Learn. Represent., Singapore, 2025. URL https://openreview.net/forum?id=SBCMNc3Mq3. Chopra (2012) Deepak Chopra. Advances in understandingof chemical bonding: Inputs from experimental and theoretical charge density analysis. J. Phys. Chem. A, 116(40):9791–9801, 08 2012. ISSN 1089-5639. doi: 10.1021/jp306169f. Dong et al. (2025) Luqi Dong, Shuxiang Yang, Su-Huai Wei, and Yunhao Lu. Variational machine learning model for electronic structure optimization via the density matrix. Phys. Rev. Lett., 135:256403, Dec 2025. doi: 10.1103/wl9w-8g8r. Dunnington & Schmidt (2012) Benjamin D. Dunnington and J. R. Schmidt. Generalization of natural bond orbital analysis to periodic systems: Applications to solids and surfaces via plane-wave density functional theory. J. Chem. Theory Comput., 8(6):1902–1911, 05 2012. ISSN 1549-9618. doi: 10.1021/ct300002t. Elsborg et al. (2025) Jonas Elsborg, Luca Thiede, Alán Aspuru-Guzik, Tejs Vegge, and Arghya Bhowmik. ELECTRA: A cartesian network for 3d charge density prediction with floating orbitals. In Adv. Neural Inf. Process. Syst., volume 38, p. 31897–31926, San Diego, CA, USA and Mexico City, Mexico, 2025. doi: 10.52202/085713-0947. Focassio et al. (2023) Bruno Focassio, Michelangelo Domina, Urvesh Patil, Adalberto Fazzio, and Stefano Sanvito. Linear jacobi–legendre expansion of the charge density for machine learning-accelerated electronic structure calculations. npj Comput. Mater., 9(1):87, 2023. doi: 10.1038/s41524-023-01053-0. Fu et al. (2024) Xiang Fu, Andrew Scott Rosen, Kyle Bystrom, Rui Wang, Albert Musaelian, Boris Kozinsky, Tess Smidt, and Tommi Jaakkola. A recipe for charge density prediction. In Adv. Neural Inf. Process. Syst., volume 37, p. 9727–9752, Vancouver, Canada, 2024. doi: 10.52202/079017-0310. Gasteiger et al. (2020) Johannes Gasteiger, Janek Groß, and Stephan Günnemann. Directional message passing for molecular graphs. In Int. Conf. Learn. Represent., Addis Ababa, Ethiopia, 2020. URL https://openreview.net/forum?id=B1eWbxStPH. Geiger & Smidt (2022) Mario Geiger and Tess Smidt. e3n: Euclidean Neural Networks, 2022. URL https://arxiv.org/abs/2207.09453. Hellweg & Rappoport (2015) Arnim Hellweg and Dmitrij Rappoport. Development of new auxiliary basis functions of the Karlsruhe segmented contracted basis sets including diffuse basis functions (def2-SVPD, def2-TZVPPD, and def2-QVPPD) for RI-MP2 and RI-C calculations. Phys. Chem. Chem. Phys., 17(2):1010–1017, 01 2015. ISSN 1463-9076. doi: 10.1039/c4cp04286g. Heyd et al. (2003) Jochen Heyd, Gustavo E. Scuseria, and Matthias Ernzerhof. Hybrid functionals based on a screened coulomb potential. J. Chem. Phys., 118(18):8207–8215, 2003. doi: 10.1063/1.1564060. Huzinaga (1965) Sigeru Huzinaga. Gaussian‐type functions for polyatomic systems. I. J. Chem. Phys., 42(4):1293–1302, 02 1965. ISSN 0021-9606. doi: 10.1063/1.1696113. Jones (2015) R. O. Jones. Density functional theory: Its origins, rise to prominence, and future. Rev. Mod. Phys., 87:897–923, Aug 2015. doi: 10.1103/RevModPhys.87.897. Jørgensen & Bhowmik (2022) Peter Bjørn Jørgensen and Arghya Bhowmik. Equivariant graph neural networks for fast electron density estimation of molecules, liquids, and solids. npj Comput. Mater., 8(1):183, 2022. doi: 10.1038/s41524-022-00863-y. Klockow et al. (2026) Manuel V. Klockow, Marc K. Ickler, Peter Lippmann, and Fred A. Hamprecht. A function-centric graph neural network approach for predicting electron densities. In Int. Conf. Learn. Represent., Lisbon, Portugal, 2026. URL https://openreview.net/forum?id=HDdkFjFEZd. Koker et al. (2024) Teddy Koker, Keegan Quigley, Eric Taw, Kevin Tibbetts, and Lin Li. Higher-order equivariant neural networks for charge density prediction in materials. npj Comput. Mater., 10(1):161, 2024. doi: 10.1038/s41524-024-01343-1. Krukau et al. (2006) Aliaksandr V. Krukau, Oleg A. Vydrov, Artur F. Izmaylov, and Gustavo E. Scuseria. Influence of the exchange screening parameter on the performance of screened hybrid functionals. J. Chem. Phys., 125(22):224106, 2006. doi: 10.1063/1.2404663. Lehtola et al. (2020) Susi Lehtola, Lucas Visscher, and Eberhard Engel. Efficient implementation of the superposition of atomic potentials initial guess for electronic structure calculations in gaussian basis sets. J. Chem. Phys., 152(14):144105, 04 2020. ISSN 0021-9606. doi: 10.1063/5.0004046. Lewis et al. (2021) Alan M. Lewis, Andrea Grisafi, Michele Ceriotti, and Mariana Rossi. Learning electron densities in the condensed phase. J. Chem. Theory. Comput., 17(11):7203–7214, 10 2021. ISSN 1549-9618. doi: 10.1021/acs.jctc.1c00576. Li et al. (2025) Chenghan Li, Or Sharir, Shunyue Yuan, and Garnet Kin-Lic Chan. Image super-resolution inspired electron density prediction. Nat. Commun., 16(1):4811, 2025. doi: 10.1038/s41467-025-60095-8. Li et al. (2020) Shikun Li, Yong Li, Marcus Bäumer, and Lyudmila V. Moskaleva. Assessment of PBE+U and HSE06 methods and determination of optimal parameter U for the structural and energetic properties of rare earth oxides. J. Chem. Phys., 153(16):164710, 10 2020. ISSN 0021-9606. doi: 10.1063/5.0024499. Lim & Whitehead (1967) TK Lim and MA Whitehead. Non self-consistent field theory—a new approach in quantum mechanical calculations. Theor. Chim. Acta, 7(1):48–63, 1967. doi: 10.1007/BF00537368. Lv et al. (2023) Taoyuze Lv, Zhicheng Zhong, Yuhang Liang, Feng Li, Jun Huang, and Rongkun Zheng. Deep Charge: Deep learning model of electron density from a one-shot density functional theory calculation. Phys. Rev. B, 108:235159, Dec 2023. doi: 10.1103/PhysRevB.108.235159. Passaro & Zitnick (2023) Saro Passaro and C. Lawrence Zitnick. Reducing SO(3) convolutions to SO(2) for efficient equivariant GNNs. In Proc. Int. Conf. Mach. Learn., volume 202 of Proceedings of Machine Learning Research, p. 27420–27438, HI, USA, 23–29 Jul 2023. PMLR. URL https://proceedings.mlr.press/v202/passaro23a.html. Perdew et al. (1996) John P. Perdew, Kieron Burke, and Matthias Ernzerhof. Generalized gradient approximation made simple. Phys. Rev. Lett., 77:3865–3868, Oct 1996. doi: 10.1103/PhysRevLett.77.3865. Perez et al. (2018) Ethan Perez, Florian Strub, Harm De Vries, Vincent Dumoulin, and Aaron Courville. Film: Visual reasoning with a general conditioning layer. In Proc. AAAI Conf. Artif. Intell., volume 32, LA, US, 2018. doi: 10.1609/aaai.v32i1.11671. Poli et al. (2020) Emiliano Poli, Kwang H Jong, and Ali Hassanali. Charge transfer as a ubiquitous mechanism in determining the negative charge at hydrophobic interfaces. Nat. Commun., 11(1):901, 2020. doi: 10.1038/s41467-020-14659-5. Qin et al. (2026) Xuejian Qin, Taoyuze Lv, and Zhicheng Zhong. EAC-Net: Predicting real-space charge density via equivariant atomic contributions. J. Chem. Theory Comput., 22(9):4813–4821, 2026. doi: 10.1021/acs.jctc.6c00283. Shuang et al. (2026) Liang Shuang, Haocheng Wang, Jiayi Song, Shuquan Ye, and Ben Fei. Ed-dit: Physics-guided diffusion pretraining for transferable molecular representations from electron density, 2026. URL https://arxiv.org/abs/2608.03260. Slater (1969) J. C. Slater. The self-consistent field for crystals. Int. J. Quantum Chem., 4(S3B):727–746, 1969. doi: https://doi.org/10.1002/qua.560040737. Song & Feng (2026) Feitong Song and Ji Feng. Neural network self-consistent fields for density functional theory. npj Comput. Mater., 12(1):289, 2026. doi: 10.1038/s41524-026-02110-0. Tang et al. (2024) Zechen Tang, He Li, Peize Lin, Xiaoxun Gong, Gan Jin, Lixin He, Hong Jiang, Xinguo Ren, Wenhui Duan, and Yong Xu. A deep equivariant neural network approach for efficient hybrid density functional calculations. Nat. Commun., 15(1):8815, 2024. doi: 10.1038/s41467-024-53028-4. Ten Eyck (1973) L. F. Ten Eyck. Crystallographic fast Fourier transforms. Acta Crystallogr. Sect. A, 29(2):183–191, Mar 1973. doi: 10.1107/S0567739473000458. van Gunsteren & Mark (1998) Wilfred F. van Gunsteren and Alan E. Mark. Validation of molecular dynamics simulation. J. Chem. Phys., 108(15):6109–6116, 04 1998. ISSN 0021-9606. doi: 10.1063/1.476021. Van Lenthe et al. (2006) J. H. Van Lenthe, R. Zwaans, H. J. J. Van Dam, and M. F. Guest. Starting SCF calculations by superposition of atomic densities. J. Comput. Chem., 27(8):926–932, 2006. doi: https://doi.org/10.1002/jcc.20393. VASP Software GmbH (2026a) VASP Software GmbH. CHGCAR, 2026a. URL https://vasp.at/wiki/index.php/CHGCAR. VASP Wiki, accessed September 9, 2026. VASP Software GmbH (2026b) VASP Software GmbH. ICHARG, 2026b. URL https://vasp.at/wiki/index.php/ICHARG. VASP Wiki, accessed September 9, 2026. Wang et al. (2021) Vei Wang, Nan Xu, Jin-Cheng Liu, Gang Tang, and Wen-Tong Geng. VASPKIT: A user-friendly interface facilitating high-throughput computing and analysis using vasp code. Comput. Phys. Commun., 267:108033, 2021. doi: https://doi.org/10.1016/j.cpc.2021.108033. Xie & Grossman (2018) Tian Xie and Jeffrey C. Grossman. Crystal graph convolutional neural networks for an accurate and interpretable prediction of material properties. Phys. Rev. Lett., 120:145301, Apr 2018. doi: 10.1103/PhysRevLett.120.145301. Xu et al. (2026a) Zemin Xu, Wenbo Xie, and P. Hu. Spectral/Spatial Tensor Atomic Cluster Expansion with Universal Embeddings in Cartesian Space, 2026a. URL https://arxiv.org/abs/2509.14961. Xu et al. (2026b) Zemin Xu, Wenbo Xie, and P. Hu. Edge cluster expansion with radial rotary attention for interatomic potentials, 2026b. URL https://arxiv.org/abs/2607.10664. Yuan et al. (2026) Zilong Yuan, Zechen Tang, Honggeng Tao, Xiaoxun Gong, Zezhou Chen, Yuxiang Wang, He Li, Yang Li, Zhiming Xu, Minghui Sun, Boheng Zhao, Chen Si, Chong Wang, Wenhui Duan, and Yong Xu. Deep-learning density functional theory hamiltonian in real space. Phys. Rev. Lett., 137:046401, Jul 2026. doi: 10.1103/mbhs-vlby. Zaccone (2026) Alessio Zaccone. Tensor Product Representation, p. 177–188. Springer Nature Switzerland, Cham, 2026. ISBN 978-3-032-27518-9. doi: 10.1007/978-3-032-27518-9_13. Appendix Appendix A Supplementary theory A.1 Kohn-Sham density functional theory and charge density self-consistency We briefly summarize the KS-DFT and SCF procedure underlying the reference charge densities used throughout this work (Jones, 2015; Slater, 1969). As introduced in Sec. 3, the GS density ρ⋆ _X minimizes the total electronic energy over NeD_N_ e. Within KS-DFT, the interacting many-electron problem is mapped onto an auxiliary system of non-interacting electrons that reproduces the same GS density. Up to the ion-ion contribution, which is independent of ρ, the corresponding energy functional is written as E[ρ]=Ts[ρ]+Eext[ρ]+EH[ρ]+Exc[ρ], E_X[ρ]=T_ s[ρ]+E_ ext^X[ρ]+E_ H[ρ]+E_ xc[ρ], (16) where TsT_ s is the kinetic energy of the non-interacting KS system, EextE_ ext^X is the interaction with the external potential generated by the atomic structure, EHE_ H is the classical electrostatic electron-electron interaction, and ExcE_ xc collects the remaining exchange and correlation effects. Taking the functional derivative with respect to the density gives the density-dependent KS Hamiltonian HKS[ρ]=−12∇2+vext()+vH[ρ]()+vxc[ρ](),vH[ρ]=δEHδρ,vxc[ρ]=δExcδρ. H_ KS[ρ]=- 12∇^2+v_ ext^X(r)+v_ H[ρ](r)+v_ xc[ρ](r), v_ H[ρ]= δ E_ Hδρ, v_ xc[ρ]= δ E_ xcδρ. (17) Here and below we use atomic units. Importantly, both vHv_ H and vxcv_ xc depend on the density itself. The KS problem is therefore nonlinear even though the Hamiltonian is linear once a particular input density is fixed. For a periodic system, solving the KS equations at each sampled k point gives HKS[ρ]ϕn()=ϵnϕn(),ρ()=∑nwfn|ϕn()|2, H_ KS[ρ] _nk(r)= _nk _nk(r), ρ(r)= _k _nw_kf_nk | _nk(r) |^2, (18) where w_k and fnf_nk are the Brillouin-zone weights and orbital occupations, respectively. The occupations satisfy ∫Ωd3ρ()=∑nwfn=Ne, _ d^3rρ(r)= _k _nw_kf_nk=N_ e, (19) consistent with the density space NeD_N_ e introduced in Eq. 2. The same construction can be expressed through the one-particle density matrix appearing in the SCF cycle of Sec. 1. For orbitals expanded in an arbitrary finite basis () χ_k(r) with coefficient matrix ()C(k), we have ()=()()†(),ρ()=∑w†()()(). (k)=C(k)f(k)C (k), ρ(r)= _kw_k χ_k (r)D(k) χ_k(r). (20) It should be noted that the numerical representation of D changes with the orbital basis, whereas the reconstructed real-space scalar field ρ()ρ(r) does not. This distinction is also one of the motivations for directly learning the charge density in AIDEN. Although ρ⋆ _X is defined variationally, practical KS calculations usually obtain it through a fixed-point SCF iteration rather than by directly optimizing every real-space grid value. Starting from an input density ρ(t) _X^(t), we construct HKS[ρ(t)]H_ KS[ _X^(t)], solve Eq. 18, and use the occupied orbitals to obtain a new output density ρ~(t+1) ρ_X^(t+1). A simple density-mixing step can be written as ρ~(t+1)=[ρ(t)],ρ(t+1)=(1−α)ρ(t)+αρ~(t+1), ρ_X^(t+1)=K_X [ _X^(t) ], _X^(t+1)=(1-α) _X^(t)+α ρ_X^(t+1), (21) where K_X denotes one KS density update and α is the mixing coefficient. More elaborate Pulay- or quasi-Newton-type mixing schemes replace the second operation in practical electronic-structure codes, but the underlying fixed-point problem remains unchanged: ρ⋆=[ρ⋆]. _X =K_X [ _X ]. (22) The cycle is repeated until the input and output densities, total energies, or their corresponding residuals satisfy the prescribed convergence threshold. Interestingly, the charge density therefore plays two roles simultaneously: it determines the effective KS Hamiltonian and is reconstructed from the eigenstates of that same Hamiltonian. This mutual dependence is precisely the origin of self-consistency. Following the notation of Eq. 1, the complete iteration can be summarized as ρ(t)()→HKS[ρ(t)]→ϕn(t+1),ϵn(t+1)→(t+1)→ρ~(t+1)()→ρ(t+1)(). _X^(t)(r)→ H_ KS [ _X^(t) ]→ \ _nk^(t+1), _nk^(t+1) \ ^(t+1)→ ρ_X^(t+1)(r)→ _X^(t+1)(r). (23) A.2 PAW density target and fixed-density evaluation The formal density in Eq. 2 motivates learning a real-space field. Numerically, AIDEN predicts the smooth FFT-grid density block of the VASP/PAW representation. This target depends on the chosen PAW datasets and valence partition, including semicore electrons where present. Its integral gives the valence-electron count, and it should not be identified with the complete all-electron density. The CHGCAR file additionally stores PAW one-center occupancies, which carry information needed for fixed-density restarts (VASP Software GmbH, 2026a). AIDEN does not predict these occupancies. Its learnable one-center branch ρinit _ init is an element-dependent contribution to the grid field, distinct from the PAW one-center data. For the model-based OOD property calculations, the predicted grid density is used without charge renormalization and combined with one-center data from the matched PBE-SCF reference. The two models therefore share the same reference conditioning within each structure. SAD instead retains its own atomic one-center data. With ICHARG=11, VASP holds the supplied density fixed during electronic minimization (VASP Software GmbH, 2026b). The reported energies, forces, and spectra therefore correspond to fixed-density VASP calculations conditioned on the matched PAW one-center data. In particular, the forces are the fixed-density VASP estimates, not derivatives through AIDEN with respect to atomic positions. Comparisons with SAD also include its different one-center data, as discussed for AlMg in Appendix C.4. A.3 Theory of superposition of atomic densities We briefly introduce the physical motivation of the SAD, which is closely related to the density decomposition adopted in AIDEN. The central approximation of SAD is to regard a many-atom system, at zeroth order, as a collection of isolated atoms and construct its initial electron density by adding the densities of the constituent atoms (Van Lenthe et al., 2006). For an isolated atom of species Z, let ρZatom()=ρZatom(||) _Z atom(r)= _Z atom(|r|) denote its spherically averaged atomic density. For the periodic structure =(,Zi,ii=1Na)X=(A,\Z_i,s_i\_i=1^N_ a), the corresponding SAD density can be written as ρSAD(,)=∑i=1Na∑∈ℤ3ρZiatom[−(i+)]. _ SAD(r;X)= _i=1^N_ a _n ^3 _Z_i atom [r-(R_i+nA) ]. (24) Each contribution depends only on the atomic species and the displacement from its center, so ρSAD _ SAD contains no explicit information about the surrounding chemical environment. The physical basis of this approximation is that a substantial fraction of the real-space density is already determined by the element-dependent atomic electronic structure. In particular, the strong density variation around each nucleus is largely inherited from the corresponding atomic shells. SAD therefore provides a physically meaningful one-center baseline before bonding and other interatomic effects are taken into account. This idea has long been used to initialize SCF calculations, where atomic densities are first generated separately and then combined to construct the molecular density matrix (Van Lenthe et al., 2006). A closely related independent-atom approximation can also be made at the potential level: in the superposition of atomic potentials (SAP), the effective one-particle potential of the full system is approximated by a sum of atomic effective potentials (Lehtola et al., 2020). SAD and SAP operate on different physical quantities, but both rely on the same zeroth-order picture that the dominant local electronic structure can be inherited from isolated atoms. Figure 5: Illustration of the components of the AIDEN-predicted charge density ρ^() ρ(r). (a) The ρ^init(,Z) ρ_ init(r;Z) branch, whose shape depends only on the elemental species Z, is first initialized through pretraining and subsequently optimized jointly with ρ^env() ρ_ env(r). (b) The ρ^env() ρ_ env(r) branch, which is responsible for expressing complex spatial charge-density patterns and can also be interpreted as a refinement of ρ^init(,Z) ρ_ init(r;Z); accordingly, ρ^env() ρ_ env(r) may take signed values and is shown here with dashed lines. We use a fictitious NaCl diatomic unit cell for illustration. Blue, orange, and green denote Na, Cl, and the charge density contributed by the local-environment branch, respectively. The displayed values are actual inferences from a pretrained AIDEN model with L=4L=4, and the final total charge density satisfies ρ^(,NaCl)>0 ρ(r;NaCl)>0. The SAD density is not, however, the self-consistent GS density of the interacting system. Once the atoms are brought together, the KS potential depends on the complete environment and the density is redistributed through chemical bonding, hybridization, polarization, screening, and charge transfer. We can therefore write ρ⋆()=ρSAD(,)+Δρ(), _X (r)= _ SAD(r;X)+ _X(r), (25) where Δρ _X denotes the environment-induced redistribution relative to the independent-atom reference. When the atomic reference and the self-consistent density contain the same number of electrons, this redistribution only rearranges charge in real space, ∫ΩΔρ()d3=0. _ _X(r)\, d^3r=0. (26) Accordingly, SAD should be interpreted as an initial physical reference rather than an approximation to the fully converged density itself. In the original SAD construction, the summed atomic density matrix is generally non-idempotent and is used to construct a starting Fock or KS operator before the usual SCF iterations recover self-consistency (Van Lenthe et al., 2006). This picture directly motivates the decomposition used in AIDEN. In our case, we do not hard-code a fixed SAD density. Instead, we represent the predicted density as ρ^(,)=ρinit()+ρenv(), ρ_ θ(r;X)= _ init(r)+ _ env(r), (27) where ρinit _ init is a learnable element-dependent one-center contribution and ρenv _ env introduces the environment-dependent correction. The former inherits the independent-atom intuition of SAD, while the latter accounts for the density redistribution that cannot be described by isolated atomic densities. Since ρinit _ init is learned jointly with the complete model after a short pretraining stage, it should not be identified with a fixed SAD reference. Rather, AIDEN generalizes the SAD picture into a learnable decomposition in which both the atomic baseline and the environment-induced contribution are optimized directly from self-consistent reference densities. A.4 Spherical-harmonic representation Spherical harmonics provide a natural angular basis for representing 3D geometric information under rotations. For angular order ℓ , the corresponding irreducible representation of SO(3) contains 2ℓ+12 +1 independent components indexed by m=−ℓ,…,ℓm=- ,…, . In AIDEN, however, we do not explicitly store the conventional complex spherical harmonics YℓmY_ m. Instead, we use their equivalent irreducible Cartesian representation throughout the Cartesian encoder. Given a unit direction d, the rank-ℓ Cartesian harmonic used in our implementation is (ℓ)(^)=(2ℓ−1)!!ℓ!ℓirr[^⊗ℓ],(0)=1, ^( )( d)= (2 -1)!! !P_ irr [ d ], ^(0)=1, (28) where ℓirrP_ irr projects the tensor product onto the fully symmetric traceless rank-ℓ subspace. Thus, although (ℓ)Y^( ) is stored using 3ℓ3 Cartesian entries, permutation symmetry and the traceless constraints leave exactly 2ℓ+12 +1 independent degrees of freedom. Under a rotation ∈SO(3)R∈ SO(3), it transforms as (ℓ)(^)=⊗ℓ(ℓ)(^)Y^( )(R d)=R Y^( )( d), while inversion gives (ℓ)(−^)=(−1)ℓ(ℓ)(^)Y^( )(- d)=(-1) Y^( )( d). Thus the angular dependence and its natural parity are fixed analytically, and no individual Cartesian component needs to be learned. Figure 6: Illustration of the real spherical-harmonic basis for angular orders ℓ=0,1,2 =0,1,2. Each order contains 2ℓ+12 +1 independent components indexed by m=−ℓ,…,ℓm=- ,…, , corresponding to the irreducible SO(3) subspace represented by the symmetric-traceless Cartesian tensors in AIDEN. The two colors indicate the positive and negative signs of the angular basis functions. The same irreducible tensor can alternatively be expressed in a compact real-spherical basis. Let ℓU_ denote the fixed orthogonal change of basis used by the model. We then have (ℓ)=ℓvec[(ℓ)]∈ℝ2ℓ+1,(ℓ)=vec−1[ℓ†(ℓ)],=⨁ℓ=0Lℓ. ^( )=U_ vec [Y^( ) ] ^2 +1, ^( )= vec^-1 [U_ y^( ) ], = _ =0^LU_ . (29) So, the Cartesian and real-spherical forms do not represent different physical information; they are simply two bases of the same angular irreducible space. The Cartesian form is convenient for the tensor products, contractions, and symmetric-traceless projections used in GIE and Cartesian ACE, whereas the compact spherical form is convenient for the edge-aligned SO(2) operations in TECE and for density decoding. The transformation U is fixed by the angular algebra and contains no learnable parameters. The resulting hierarchy is intuitive. For ℓ=0 =0, there is only one isotropic scalar component; ℓ=1 =1 contains three vector-like components; and ℓ=2 =2 contains five quadrupolar components, corresponding to the familiar s-, p-, and d-like angular patterns illustrated in our schematic. More generally, each angular order ℓ contains 2ℓ+12 +1 components, while the channel index c=1,…,Cc=1,…,C simply provides multiple learned copies of the same representation and does not change its rotational transformation law. The positive and negative lobes in the schematic indicate the sign of the real angular basis functions rather than different atomic densities or physical orbitals. In this way, learned channel amplitudes encode the chemical environment, while the analytic harmonic basis supplies the required angular structure, SO(3) transformation law, and natural parity (−1)ℓ(-1) . A.5 Equivariance of AIDEN We characterize the equivariance of AIDEN under Euclidean transformations by first considering proper rotations and translations and then extending the symmetry characterization to inversion. Let g=(,)∈SE(3)g=(Q,t)∈ SE(3) with ∈SO(3)Q∈ SO(3) and g=⊤+gr=rQ +t. Under the row-vector convention of Sec. 3, the transformed periodic structure ggX has lattice ′=⊤A =AQ and Cartesian atomic positions i′=i⊤+R_i =R_iQ +t, with the corresponding fractional coordinates understood modulo lattice translations. Accordingly, each periodic image transforms as i+↦(i+)⊤+R_i+nA (R_i+nA)Q +t. Since AIDEN constructs geometric features exclusively from relative displacements, the global translation t cancels identically. We may therefore focus first on the rotational part of the transformation. It is useful to note that every ∈O(3)Q∈ O(3) can be written as =Q=PR, where ∈SO(3)R∈ SO(3) and ∈,−P∈\I,-I\. Thus, together with translations, the E(3) E(3) transformation law can be characterized by proper rotations and spatial inversion. Irreducible Cartesian representation and encoder equivariance. For a rank-ℓ Cartesian tensor T, let ℓ() D_ (Q) denote the natural action of Q on all Cartesian indices, i.e. [ℓ()]a1…aℓ=∑b1,…,bℓQa1b1⋯QaℓbℓTb1…bℓ[ D_ (Q)T]_a_1… a_ = _b_1,…,b_ Q_a_1b_1·s Q_a_ b_ T_b_1… b_ , with 0()=1 D_0(Q)=1. We further write ()=⨁ℓ=0Lℓ() D(Q)= _ =0^L D_ (Q). Restricted to the symmetric traceless rank-ℓ subspace, ℓ D_ realizes the (2ℓ+1)(2 +1)-dimensional irreducible representation of SO(3) SO(3) used throughout AIDEN. Since ⊤=Q Q=I, this representation preserves both the Frobenius inner product and the corresponding tensor norm. The irreducible projection ℓirrP_ irr commutes with rotations: symmetrization is unaffected by applying the same rotation to every Cartesian index, while traces are preserved by ∑aQabQac=δbc _aQ_abQ_ac= _bc. It follows that the irreducible Cartesian harmonics and the Cartesian coupling satisfy (ℓ)(^⊤)=ℓ()(ℓ)(^),ℓirrℓ()=ℓ()ℓirr, ^( )( dQ )= D_ (Q)Y^( )( d), _ irr D_ (Q)= D_ (Q)P_ irr, (30) ℓ1,ℓ2ℓ(ℓ1()1,ℓ2()2)=ℓ()ℓ1,ℓ2ℓ(1,2). _ _1, _2 ( D_ _1(Q)T_1, D_ _2(Q)T_2 )= D_ (Q)C_ _1, _2 (T_1,T_2). (31) The second identity follows directly from the structure of ℓ1,ℓ2ℓC_ _1, _2 , which contains only tensor products, contractions with the Euclidean metric, and irreducible projection. Since ℬν(ℓ)B_ν^( ) is recursively assembled from the same couplings together with fixed-ℓ channel maps, induction over the correlation order gives ℬν(ℓ)[()]=ℓ()ℬν(ℓ)[]B_ν^( )[ D(Q)T]= D_ (Q)B_ν^( )[T]. For an edge e=(j,→i)e=(j,n→ i), Eq. 4 transforms as e′=e⊤d_e =d_eQ , while re′=rer_e =r_e and ^e′=^e⊤ d_e = d_eQ . Thus the radial representation (re)b(r_e), the distance-defined neighborhood iN_i, and the elemental features ix_i are unchanged. The scalar GIE branch is therefore invariant. For ℓ≥1 ≥ 1, all learned factors multiplying (ℓ)(^e)Y^( )( d_e) in Eq. 6 are invariant channel amplitudes, so Eq. 30 immediately yields i(0,ℓ)(g) _i^(0, )(gX) =ℓ()i(0,ℓ)(),0≤ℓ≤L. = D_ (Q)h_i^(0, )(X), 0≤ ≤ L. (32) Every map ℒ(ℓ)L^( ) used by the encoder mixes channels only within a fixed angular order and hence commutes with ℓ() D_ (Q); correspondingly, biases are confined to the invariant ℓ=0 =0 sector. The normalization ℓV_ behaves in the same way, because its denominator is constructed from sums of squared Cartesian components and is therefore rotationally invariant, while its learnable gain acts only on channels. Hence ℓ[ℓ()]=ℓ()ℓ[]V_ [ D_ (Q)h]= D_ (Q)V_ [h]. We next consider Cartesian ACE. Given equivariant normalized inputs, the radial weights e,ℓ1ℓ2ℓw_e, _1 _2 in Eq. 7 remain invariant, whereas the source feature and Cartesian harmonic transform according to ℓ1 D_ _1 and ℓ2 D_ _2, respectively. Using Eq. 31, we obtain e(ℓ)(g)=ℓ()e(ℓ)()p_e^( )(gX)= D_ (Q)p_e^( )(X). In particular, the ℓ=0 =0 component e(0)p_e^(0) is invariant, and so is the scalar neighbor coefficient aea_e. We therefore have i(g) _i(gX) =()i(). = D(Q) _i(X). (33) The gate σ0ℱg[i(0)] _0F_ g[ _i^(0)] is generated entirely from the invariant sector and consequently acts as a scalar on every angular component. Together with the equivariance of ℬν(ℓ)B_ν^( ) and the residual connection in Eq. 8, this gives i(1,ℓ)(g) _i^(1, )(gX) =ℓ()i(1,ℓ)(). = D_ (Q)h_i^(1, )(X). (34) Figure 7: E(3)-equivariance test of the model-predicted charge density for GaAs. (a) Charge density ρ()ρ(r) predicted by AIDEN for zincblende GaAs (mp-2534) on a 32×32×3232× 32× 32 grid. (b)-(d) Normalized pointwise deviation δρ():=|ρg()−ρ()|/ρ() _ρ(r):=| _g(r)-ρ(r)|/ρ(r), where ρg _g denotes the prediction after applying the transformation g and ρ()ρ(r) denotes the reference charge density at ∈ℝ3r ^3 in (a). The tested transformations are (b) translation by =(3,−5,7)/32t=(3,-5,7)/32, (c) proper rotation about =(−1,−2,1)/6n=(-1,-2,1)/ 6 with det=+1 =+1, and (d) inversion →−r→-r with det=−1 =-1. The transformed predictions remain numerically consistent with the reference. Notably, inversion is not a symmetry operation of non-centrosymmetric zincblende GaAs, providing a direct test beyond the crystal space-group symmetries. SO(3) equivariance of TECE and radial rotary attention. TECE performs its nonlinear edge interaction after converting the global irreducible Cartesian representation into an edge-aligned real-spherical representation. Let () S(Q) denote the corresponding real-spherical representation. Since U is a fixed orthogonal change of basis between two realizations of the same SO(3) irreducible representation, it satisfies the intertwining relation () D(Q) =(). = S(Q)U. (35) Equivalently, on the irreducible subspace we have †()=()†U S(Q)= D(Q)U . To describe the edge-aligned frame, let e∈SO(3)G_e∈ SO(3) denote the spatial rotation represented by e=(e)D_e= S(G_e), chosen such that the local y axis is aligned with ^e d_e. After applying the global rotation Q, let e′G_e denote the frame associated with ^e⊤ d_eQ . The two frames describe the same physical edge in local coordinates and can therefore differ only by a rotation around the aligned y axis. Defining e=e′e⊤∈SO(2)yH_e=G_e QG_e ∈ SO(2)_y and e=(e)S_e= S(H_e), we obtain e′() _e S(Q) =ee. =S_eD_e. (36) For the local source and destination features in Eq. 9, Eqs. 34-36 therefore imply es′=eesz_e s =S_ez_e s and et′=eetz_e t =S_ez_e t. For m>0m>0, it is convenient to identify each real-imaginary pair with the complex variable zℓm=zℓmRe+izℓmImz_ m=z_ m Re+ iz_ m Im, for which the residual local rotation takes the simple form zℓm↦eimαezℓmz_ m e im _ez_ m, while the m=0m=0 sector remains unchanged. The radial operator e[(re)] _e[b(r_e)] depends only on rer_e and applies the same radial coefficient to the real and imaginary components of a given m sector, so it commutes with the local SO(2) action. Likewise, the linear maps inside TECE are block diagonal in magnetic order and use the same channel transformation for the real and imaginary components at each m>0m>0. The second-order correlator ℰE then follows the usual SO(2) frequency-addition rule, zm1zm2 z_m_1z_m_2 ↦ei(m1+m2)αezm1zm2,zm1z¯m2↦ei(m1−m2)αezm1z¯m2. e i(m_1+m_2) _ez_m_1z_m_2, z_m_1 z_m_2 e i(m_1-m_2) _ez_m_1 z_m_2. (37) Accordingly, the corresponding products contribute only to sectors with mout=m1+m2m_ out=m_1+m_2 or mout=|m1−m2|m_ out=|m_1-m_2|. Their edge-dependent coefficients are generated from the invariant m=0m=0 sector and radial quantities and thus behave as SO(2) scalars. In this way, the complete local correlator satisfies ℰ(e)=eℰ()E(S_ez)=S_eE(z). We now turn to RRA. At fixed (ℓ,m,h)( ,m,h), the projected query and key acquire the same residual phase, e,ℓmh′=eimαee,ℓmhq_e, mh = e im _eq_e, mh and e,ℓmh′=eimαee,ℓmhk_e, mh = e im _ek_e, mh. Since φeh _eh depends only on radial features and the channel contraction in Eq. 11 is Hermitian, the common local phase cancels: Re⟨eimαee,ℓmh,eimφeheimαee,ℓmh⟩ch Re e im _eq_e, mh, e im _eh e im _ek_e, mh _ ch =Re⟨e,ℓmh,eimφehe,ℓmh⟩ch. = Re _e, mh, e im _ehk_e, mh _ ch. (38) Thus both ξeh _eh and the normalized coefficient aeha_eh are invariant. The attention operator eA_e acts only on channel indices and therefore commutes with eS_e. Let ev_e denote the resulting local edge output after radial modulation, ℰE, and eA_e. We then have e′=eev_e =S_ev_e. Using Eq. 36, e′†e=()e†D_e S_e= S(Q)D_e , and hence †e′†e′ D_e v_e =()†e†e. = D(Q)U D_e v_e. (39) Neighbor summation and the residual connection in Eq. 9 preserve this transformation law. After the final equivariant normalization ℓV_ , we obtain i⋆(ℓ)(g) _i ( )(gX) =ℓ()i⋆(ℓ)(),0≤ℓ≤L. = D_ (Q)h_i ( )(X), 0≤ ≤ L. (40) Therefore, we obtain an SO(3)-equivariant hierarchy of atomic representations, while global translations leave the encoded atomic features unchanged. To complete the E(3) E(3) characterization, we further consider spatial inversion =−P=-I. The irreducible Cartesian harmonics have the natural parity (ℓ)(−^) ^( )(- d) =(−1)ℓ(ℓ)(^). =(-1) Y^( )( d). (41) Accordingly, the Cartesian GIE and ACE representations inherit the parity factor (−1)ℓ(-1) , while all radial quantities remain unchanged. We additionally evaluate the complete transformation numerically on non-centrosymmetric zinc-blende GaAs. The normalized L1L_1 deviations under translation, a generic proper rotation, and inversion are 7.72×10−67.72× 10^-6, 7.27×10−57.27× 10^-5, and 9.51×10−79.51× 10^-7, respectively. In particular, inversion preserves the predicted density to numerical precision even though the transformed GaAs structure is not symmetry-equivalent to the original structure. Since every improper orthogonal transformation can be decomposed into inversion followed by a proper rotation, these results, together with the analytical SO(3) transformation law above, provide numerical evidence for O(3) consistency of the trained AIDEN mapping and hence its E(3) consistency after including translations. Continuous density covariance and lattice periodicity. The decoder maps the final atomic tensors to coefficient tensors according to ikpχ(ℓ)=ℒkpχ(ℓ)i⋆(ℓ)+δℓ0μkpχc_ikp^χ( )=L_kp^χ( )h_i ( )+ _ 0 _kp^χ. Since ℒkpχ(ℓ)L_kp^χ( ) mixes channels only within the same angular order and the bias is restricted to ℓ=0 =0, Eq. 40 gives ikpχ(ℓ)(g) _ikp^χ( )(gX) =ℓ()ikpχ(ℓ)(). = D_ (Q)c_ikp^χ( )(X). (42) For a periodic image η=(i,)η=(i,n) and a simultaneously transformed query position, Eq. 12 gives η′=η⊤ δ_η = δ_ηQ , rη′=rηr_η =r_η, and ^η′=^η⊤ δ_η = δ_ηQ . Thus the radial function R~ℓp(rη) R_ p(r_η) remains invariant, while the angular factor transforms according to Eq. 30. Because ℓ() D_ (Q) is orthogonal, the Cartesian contraction is invariant: ⟨ikpχ(ℓ)(g),(ℓ)(^η′)⟩F _ikp^χ( )(gX),Y^( )( δ_η ) _ F =⟨ikpχ(ℓ)(),(ℓ)(^η)⟩F. = _ikp^χ( )(X),Y^( )( δ_η) _ F. (43) Consequently, every contribution entering Φkχ _k^χ is invariant under the joint transformation of structure and query, giving Φkχ(g,g)=Φkχ(,) _k^χ(gr;gX)= _k^χ(r;X). The one-center field ϕηk0χ _η k_0^χ depends only on the species ZiZ_i and the radial ℓ=0 =0 functions, and therefore obeys the same transformation law. As a result, ρinit(g,g) _ init(gr;gX) =ρinit(,),ρenv(g,g)=ρenv(,), = _ init(r;X), _ env(gr;gX)= _ env(r;X), (44) and the complete continuous density satisfies ρ^(g,g) ρ_ θ(gr;gX) =ρ^(,),∈SO(3). = ρ_ θ(r;X), ∈ SO(3). (45) For inversion, Eq. 41 gives the same scalar contraction because the corresponding angular factors acquire the parity (−1)ℓ(-1) . Combining inversion with an arbitrary proper rotation then extends the scalar-field transformation to improper orthogonal transformations as well. It should also be noted that any pointwise scalar transformation applied after these invariant contractions preserves the same transformation law. Finally, lattice periodicity follows naturally from the periodic-image construction. For any ∈ℤ3m ^3, associate the image η=(i,)η=(i,n) at r with η′=(i,+)η =(i,n+m) at +r+mA. Their relative displacements satisfy η′(+) δ_η (r+mA) =+−[i+(+)]=η(). =r+mA- [R_i+(n+m)A ]= δ_η(r). (46) Thus translating the query by a lattice vector merely relabels the periodic images contributing to ℳ()M(r), without changing any radial or angular quantity entering the decoder. We therefore have ρ^(+,) ρ_ θ(r+mA;X) =ρ^(,). = ρ_ θ(r;X). (47) Taken together, Eqs. 40, 41, 45, and 47 show that AIDEN carries SO(3)-equivariant atomic representations with the natural parity (−1)ℓ(-1) , is invariant to global translations at the representation level, and realizes an E(3)-consistent, lattice-periodic mapping from periodic atomic structures to continuous scalar density fields. Appendix B Detailed implementations of AIDEN This section collects the implementation details needed to reproduce the AIDEN architecture in Sec. 4. We retain only the constructions that define the periodic geometric basis, the two equivariant cluster-expansion interactions, and the continuous GTO decoder; numerical bookkeeping that does not affect the model definition is omitted. B.1 Periodic graph and geometric basis The elemental descriptor (Z)q(Z) contains atomic number, period, group, total and orbital-resolved valence counts, electronegativity, atomic and covalent radii, first ionization energy, electron affinity, dipole polarizability, and GS magnetic moment. It is concatenated with the one-hot vector (Z)o(Z) and feature-wise min-max normalized over the 118 elements, yielding the 133133-dimensional input ix_i used in Sec. 4.1. No structure-dependent information enters this elemental representation. For periodic crystals, all atomic images satisfying the cutoff in Eq. 4 are enumerated before the nearest NnbrN_ nbr images are retained for each target atom. We use Nnbr=100N_ nbr=100 and rat=4År_ at=4\, . Query-atom images are constructed independently with the orbital cutoff rorbr_ orb in Eq. 12 and are not neighbor-capped. For VASP CHGCAR data, FFT-grid points are interpreted in fractional coordinates and mapped to Cartesian positions by g=gr_g=u_gA; the stored density values are divided by |det|| | to obtain the supervised density in e/Å3e/ ^3. The atomic edge basis consists of a zero-order spherical-Bessel radial expansion multiplied by the standard compact polynomial envelope fcf_ c (Gasteiger et al., 2020), together with irreducible Cartesian harmonics, b(r) b_b(r) =2ratsin(bπr/rat)rfc(r),b=1,…,NB, = 2r_ at (bπ r/r_ at)rf_ c(r), b=1,…,N_ B, (48) (ℓ)(^) ^( )( v) =(2ℓ−1)!!ℓ!ℓirr(^⊗ℓ),(0)=1. = (2 -1)!! !P_ irr ( v ), ^(0)=1. (49) Here ℓirrP_ irr denotes the symmetric-traceless projection used throughout the Cartesian encoder. We use NB=8N_ B=8. The same radial vector (re)b(r_e) is shared by GIE, Cartesian ACE, and TECE. B.2 Cartesian atomic cluster expansion For an angular path (ℓ1,ℓ2,ℓ)( _1, _2, ), define κ=(ℓ1+ℓ2−ℓ)/2κ=( _1+ _2- )/2. The Cartesian coupling contracts κ index pairs and projects the remaining tensor to rank ℓ , ℓ1,ℓ2ℓ(1,2) _ _1, _2 (T_1,T_2) =3−κ/2ℓirr(κ(1,2)), =3^-κ/2P_ irr (K_κ(T_1,T_2) ), (50) L _L =(ℓ1,ℓ2,ℓ)∈0,…,L3||ℓ1−ℓ2|≤ℓ≤min(L,ℓ1+ℓ2),ℓ1+ℓ2−ℓ≡0(mod2). = \( _1, _2, )∈\0,…,L\^3\ |\ | _1- _2|≤ ≤ (L, _1+ _2),\ _1+ _2- ≡ 0 2 \. (51) This is the coupling used in Eq. 7; for L=3L=3, the path set contains 2323 admissible triples. The equivariant RMS normalization applied before the interactions and at the encoder output rescales each angular order by a rotation-invariant channel RMS, ℓ[i(ℓ)]c=γℓcic(ℓ)max(C−1∑c′=1C‖ic′(ℓ)‖F2,ϵ). _ [h_i^( )]_c= _ ch_ic^( ) (C^-1 _c =1^C\|h_ic ^( )\|_ F^2,ε ). (52) The learned gain γℓc _ c acts only on channels, so the Cartesian transformation law is unchanged. For the atomic aggregation in Eq. 8, the scalar sector of the edge field defines an invariant logit se=ℱa[e(0)⊕(re)]s_e=F_ a[p_e^(0) (r_e)]. The corresponding destination-wise coefficient is ae=fc(re)ese∑e′∈ifc(re′)ese′. a_e= f_ c(r_e) e^s_e _e _if_ c(r_e ) e^s_e . (53) The same scalar aea_e multiplies every angular order of edge e. To make the correlation polynomial in Eq. 8 explicit, let i=C+σ0ℱg[i(0)]⊙i _i=\1_C+ _0F_ g[ _i^(0)]\ _i. Starting from i[1](ℓ)=i(ℓ)B_i^[1]( )= _i^( ), higher correlation orders are generated recursively and combined as i[q](ℓ) _i^[q]( ) =∑(ℓ1,ℓ2,ℓ)∈Lℓ1,ℓ2ℓ(i[q−1](ℓ1),i(ℓ2)),ℬν(ℓ)(i)=∑q=1νℒcorr,q(ℓ)i[q](ℓ). = _( _1, _2, ) _LC_ _1, _2 (B_i^[q-1]( _1), _i^( _2) ), _ν^( )( _i)= _q=1^νL_ corr,q^( )B_i^[q]( ). (54) We use ν=2ν=2. This order refers to the explicit Cartesian correlation inside the ACE block; the preceding GIE representation already contains one-hop environmental information. B.3 TECE and RRA The orthogonal map U converts each irreducible Cartesian tensor into its real-spherical components, after which eD_e aligns the local y axis with ^e d_e. For |m|≤M|m|≤ M, the compact representation contains DM=(L+1)+2∑m=1M(L+1−m)D_M=(L+1)+2 _m=1^M(L+1-m) real components. With L=M=3L=M=3, all magnetic components are retained and DM=16D_M=16. The source and destination representations are modulated by radial weights before entering the edge correlator. Suppressing channel indices, this operation can be written consistently with Eq. 9 as e=e[(re)]⊙(e¯j(1)⊕e¯i(1)), _e= _e[b(r_e)] (D_eU h_j^(1) _eU h_i^(1) ), (55) where the same radial coefficient is applied to the real-imaginary pair associated with a fixed nonzero m. This preserves the local SO(2) transformation law. The operator ℰE combines a direct branch, an SO(2)-gated branch, and a second-order product branch. Absorbing the fixed-frequency channel projections into the corresponding operators, its structure is ℰ(e)=ℒoutTECE[e+(e)+Π(2)(e,e)3],e=ℒmve,e=ℒmwe. (z_e)=L_ out TECE [ v_e+G(v_e)+ ^(2)(v_e,w_e) 3 ], _e=L_mvz_e, _e=L_mwz_e. (56) For a complex local mode zmz_m, a residual rotation around the edge axis gives zm↦eimϑzmz_m e im z_m. Accordingly, Π(2) ^(2) retains only products with output frequencies m=m1+m2m=m_1+m_2 or m=|m1−m2|m=|m_1-m_2|. Their coefficients are generated from the invariant m=0m=0 sector, so the product branch remains SO(2)-equivariant. The direct, gated, and product branches are then returned to the global Cartesian representation through e†D_e and †U as in Eq. 9. Radial rotary attention is evaluated from the unmodulated local target and source tensors. Splitting the CeC_ e channels into H heads gives Ch=Ce/HC_h=C_ e/H; for m>0m>0, each real-imaginary pair is identified with a complex vector, and the head-wise contraction is ⟨,⟩ch=∑c=1Chuc∗vc ,v _ ch= _c=1^C_hu_c^*v_c. The score ξeh _eh is given in Eq. 11. Its normalized coefficient is aeh=fc(re)eξeh∑e′∈ifc(re′)eξe′h, a_eh= f_ c(r_e) e _eh _e _if_ c(r_e ) e _e h, (57) which is invariant because it depends only on equal-frequency inner products, radial quantities, and the cutoff envelope. We use H=4H=4 and Ce=C=48C_ e=C=48. B.4 Continuous GTO decoder The decoder uses NGN_ G even-tempered Gaussian radial functions for every angular order. Their exponents and corrected radial factors are αp=αmin(αmaxαmin)p−1NG−1,R~ℓp(r)=ℓprℓe−αpr2[1+ζℓp(r)],ℓp=2(2αp)ℓ+3/2Γ(ℓ+3/2). _p= _ ( _ _ ) p-1N_ G-1, R_ p(r)=Z_ pr e^- _pr^2 [1+ _ p(r) ], _ p= 2(2 _p) +3/2 ( +3/2). (58) We use NG=8N_ G=8, αmin=0.15Å−2 _ =0.15\, ^-2, and αmax=256Å−2 _ =256\, ^-2. The element-independent correction ζℓp(r) _ p(r) is predicted from 5050 Gaussian distance features by a small MLP whose final affine layer is initialized to zero, so the decoder starts from the analytic GTO family. The same corrected radial basis is used by the one-center and environment-dependent branches. The coefficient tensors ikpχ(ℓ)c_ikp^χ( ) and the continuous fields ϕηk0χ _η k_0^χ and Φkχ _k^χ are defined directly in Sec. 4.3. Because U is orthogonal on the irreducible subspace, the Frobenius contraction in Eq. 13 is equivalent to contracting the corresponding real-spherical coefficients and harmonics, without introducing an explicit magnetic index in the decoder. We use environmental rank K=8K=8, one-center rank K0=1K_0=1, and orbital cutoff rorb=3År_ orb=3\, . As in BOA, the implementation applies a smooth absolute value to the left factor of each low-rank product before multiplication; this scalar stabilization is omitted from Eq. 14 because it does not alter the equivariant structure of the decoder. The localized GTO basis provides a continuous representation independent of the evaluation grid. The bilinear form in Eq. 14 generates cross-center terms when expanded; with the scalar stabilization, it still couples contributions from different centers. The learned factors are auxiliary fields, not occupied KS orbitals. Likewise, the one-center/environment split is a modeling choice: joint optimization does not impose a unique atomic partition or require the environment-dependent term to integrate to zero. B.5 Forward propagation Given a structure X, the equivariant encoder is evaluated once to obtain the atomic coefficients ikpχ(ℓ)c_ikp^χ( ), after which the continuous decoder evaluates ρ ρ_ θ at arbitrary query positions. Alg. 1 summarizes the resulting forward propagation. Algorithm 1 Forward propagation of AIDEN 1: Structure X and query positions gg=1Nq\r_g\_g=1^N_ q 2: ρ^(g,)g=1Nq\ ρ_ θ(r_g;X)\_g=1^N_ q 3: Construct ii\x_i\_i and the periodic graph ii\N_i\_i with (re),(ℓ)(^e)e\b(r_e),Y^( )( d_e)\_e from Eq. 4. 4: i(0,ℓ)ℓ=0L←GIE()\h_i^(0, )\_ =0^L← GIE(X) using Eqs. 5–6. 5: i(1,ℓ)←ACE(ℓ′[i(0,ℓ′)]ℓ′=0L)h_i^(1, )← ACE\! (\V_ [h_i^(0, )]\_ =0^L ) using Eq. 8. 6: i(2)←TECE(ℓ[i(1,ℓ)]ℓ=0L)h_i^(2)← TECE\! (\V_ [h_i^(1, )]\_ =0^L ) using Eq. 9. 7: i⋆(ℓ)←ℓ[i(2,ℓ)]h_i ( ) _ [h_i^(2, )] using Eq. 10. 8: ikpχ(ℓ)←ℒkpχ(ℓ)i⋆(ℓ)+δℓ0μkpχc_ikp^χ( ) _kp^χ( )h_i ( )+ _ 0 _kp^χ, χ∈L,Rχ∈\L,R\. 9: for g=1,…,Nqg=1,…,N_ q do 10: ℳg←ℳ(g)M_g (r_g) and η,rη,^η∈ℳg\ δ_η,r_η, δ_η\_η _g from Eq. 12. 11: ϕηk0χ(g),Φkχ(g)← _η k_0^χ(r_g),\, _k^χ(r_g)← Eq. 13. 12: ρinit(g)←∑η∈ℳg∑k0ϕηk0L(g)ϕηk0R(g) _ init(r_g)← _η _g _k_0 _η k_0^L(r_g) _η k_0^R(r_g), ρenv(g)←∑kΦkL(g)ΦkR(g) _ env(r_g)← _k _k^L(r_g) _k^R(r_g). 13: ρ^(g,)←ρinit(g)+ρenv(g) ρ_ θ(r_g;X)← _ init(r_g)+ _ env(r_g). 14: return ρ^(g,)g=1Nq\ ρ_ θ(r_g;X)\_g=1^N_ q During training, the ρinit _ init branch is briefly pretrained with the remaining model parameters frozen, after which the complete model is jointly optimized using Eq. 3. Appendix C Supplementary results C.1 Ablation studies We have already examined the effect of the maximum angular order L on the performance of AIDEN in Figs. 2 and 8. Here, we further evaluate the contributions of individual architectural components and the rationale behind their composition. We first remove GIE by retaining only the element-dependent scalar initialization and setting the higher-order initial features to zero. To evaluate Radial Rotary Attention (RRA), we preserve the complete edge correlation operation in TECE but replace its learned multi-head attention weights with the parameter-free cutoff-normalized weights ae(0)=fc(re)/∑e′→ifc(re′)a_e^(0)=f_ c(r_e)/ _e → if_ c(r_e ). This replacement preserves cutoff weighting and neighbor normalization while removing the learned query-key scores, radial phase, and head-dependent adaptive aggregation. We next assess the two physically motivated components of the density decoder. Specifically, we remove the one-center branch ρinit _ init and disable the learnable radial correction of the GTO basis by setting ζℓp(r)=0 _ p(r)=0 in Eq. 58. Finally, to examine the ordering of the interaction layers while keeping their number and parameter count fixed, we compare the full ACE→ architecture with the reversed TECE→ arrangement. Considering computational efficiency, we use a fixed random subset of 2,000 structures drawn from the original PBE training split. All variants are trained for 100,000 steps with L=4L=4 using the same training subset and validation set. For all variants retaining ρinit _ init, we load the same pretrained one-center checkpoint to ensure consistent initialization. The best validation NMAE is evaluated on a fixed subset of 256 structures from the original PBE validation split. The results are summarized in Tab. 3. Table 3.: Module ablation results. Parameter counts and best validation NMAE ερ _ρ for models trained on the same 2,000-structure PBE subset. Relative changes are calculated with respect to the full ACE→ model. Ablation setting Parameters ↓ Best val ερ _ρ [%] ↓ Relative to Full [%] ↓ Full ACE→ 4.459M 1.13151 – TECE→ 4.459M 1.14572 +1.26 w/o GIE 4.393M 1.14622 +1.30 w/o RRA 4.205M 1.16120 +2.62 fixed GTO 4.454M 1.21505 +7.38 w/o ρinit _ init 4.457M 1.53550 +35.70 Removing GIE increases the best validation NMAE by 1.30% relative to the full model, indicating that constructing higher-order equivariant initial representations from one-hop neighborhood geometry before the main interaction blocks facilitates the extraction of local directional information. Replacing RRA with fixed cutoff-normalized weights increases the NMAE by 2.62%. Compared with assigning neighbor contributions solely according to distance, RRA combines equivariant query-key matching, radial rotary phases, and multi-head aggregation to adaptively select edge messages according to the local chemical and geometric environment. The decoder ablations lead to more pronounced performance degradation. Removing ρinit _ init increases the NMAE by 35.70%, demonstrating the importance of the element-dependent one-center density as an atomic-density baseline, which allows the environment-dependent branch to focus on density corrections induced by bonding and coordination. Fixing the GTO basis increases the NMAE by 7.38% while reducing only 4,590 parameters, indicating that the learnable correction ζℓp(r) _ p(r) effectively adjusts the shapes of different angular and radial channels on top of the analytical even-tempered GTO basis to accommodate diverse elemental compositions and local bonding environments. The comparison of interaction order further supports the ACE→ design. Reversing the full model to TECE→ while retaining the same parameter count increases the NMAE by 1.26%. This result suggests that first constructing atom-centered representations with cross-neighbor correlations through the aggregate-then-correlate operation of ACE, followed by directional refinement in edge-aligned local coordinate frames through the correlate-then-aggregate operation of TECE, provides a more effective ordering for charge-density learning. Overall, these ablation results validate the effectiveness of GIE, RRA, one-center density initialization, and learnable GTO radial corrections, while also supporting the ordered ACE→ interaction architecture. Among these components, ρinit _ init and the GTO correction provide the most substantial accuracy gains, whereas the ACE→ ordering effectively combines cross-neighbor correlation modeling with edge-wise directional refinement. C.2 Training details Fig. 8(a-b) shows the loss curves of AIDEN on the PBE crystal dataset and the QM9 molecular dataset. For each dataset, we train models with maximum angular orders L=0,…,4L=0,…,4 while keeping all other training settings unchanged. The models are optimized using Adam, with learning rates of 2×10−42× 10^-4 and 1×10−31× 10^-3 for the encoder and decoder, respectively. We use a batch size of 12 and a gradient-clipping threshold of 0.5. Exponential moving average parameters with a decay rate of 0.995 are used for validation and inference. All models are trained on a single NVIDIA A100-40GB GPU. Figure 8: AIDEN loss curves and violin plots of error distributions. (a), (b) Loss curves of AIDEN on the PBE and QM9 datasets under different L settings. Opaque dashed lines denote the validation set, and semi-transparent solid lines denote the training set. (c), (d) Violin plots of the error distributions of AIDEN on PBE and QM9. Dark blue and cyan respectively denote the training and validation errors. In the embedded box plots, white markers indicate the mean error, the filled boxes span the 25%-75% quantile range, and the upper and lower whiskers are set to the 5% and 95% quantiles. C.3 Twisted bilayer graphene Structures and calculation settings. We consider three unrelaxed commensurate TBG cells, summarized in Tab. 4. All structures have a 20 Å out-of-plane lattice parameter. We reuse the converged PBE-SCF densities, energies, and forces as references and evaluate the complete model density grids without charge renormalization. Each model density is supplemented with the one-center PAW data from the matched reference. Thus, the property comparison evaluates the predicted smooth density under the same PAW conditioning. We use the carbon PAW-PBE potential (08Apr2002), a 520 eV cutoff, Gaussian smearing of 0.05 eV, and non-spin-polarized calculations with fixed atomic positions and symmetry disabled. The fixed-density energy/force and band calculations use ICHARG=11, with electronic convergence thresholds of 10−710^-7 and 10−810^-8 eV, respectively, and PREC=Accurate, LREAL=.FALSE., LASPH=.TRUE., and ADDGRID=.TRUE.. Band calculations additionally use LMAXMIX=2. Energies are the final TOTEN values, and forces are the frozen-density approximate forces returned by VASP. The reciprocal-space path uses fractional coordinates (0,0,0)→(1/2,0,0)→(1/3,1/3,0)→(0,0,0)(0,0,0)→(1/2,0,0)→(1/3,1/3,0)→(0,0,0). We use coarser line-mode sampling for the two larger cells to control memory use. The resulting comparisons contain 180, 18, and 18 k points, respectively, including repeated segment endpoints. It should be noted that the angular trend in the band errors is therefore evaluated at different sampling densities. Table 4.: TBG structures and sampling. The Γ -centered mesh is used for the reference SCF and fixed-density energy/force calculations. The last two columns give the number of points per segment along Γ -M-K-Γ and the number of returned bands. Sampling and band counts are matched across methods within each cell. Twist angle NaN_ a FFT grid NgN_ g SCF mesh Points/segment Bands 21.79∘21.79 28 96×96×30096× 96× 300 2,764,800 9×9×19× 9× 1 60 80 13.17∘13.17 76 160×160×300160× 160× 300 7,680,000 6×6×16× 6× 1 6 168 9.43∘9.43 148 224×224×300224× 224× 300 15,052,800 4×4×14× 4× 1 6 312 Charge density and charge conservation. Tab. 5 reports ερ _ρ using Eq. 15, together with the integrated smooth-grid electron counts. We denote the deviation of the predicted count from the reference by ΔNe N_ e. AIDEN gives lower density errors at all three angles, with ερ=0.2892 _ρ=0.2892-0.2934%0.2934\%, compared with 0.3134-0.3160% for ChargE3Net. Both models slightly overestimate the electron count, while AIDEN has the smaller deviation in every cell. Its relative charge error decreases from 0.00952% to 0.00679% as the cell grows, despite the increase in the absolute count deviation. Although we do not explicitly include a charge-conservation loss term in the training objective in Eq. 3 (as is also common in most related works), AIDEN naturally exhibits near-conservation of the total charge. This behavior is observed not only for the TBG cases but also across the general benchmark datasets (see Fig. 9), and can be largely attributed to the low overall error in charge density prediction. Figure 9: Parity plots of total charge on the general benchmark datasets. Blue and yellow scatter points denote the training and validation samples in (a) PBE and (b) QM9, respectively. The horizontal axis NeN_ e is the total DFT reference charge of each sample evaluated on the full grid, while the vertical axis N^e N_ e is the total charge predicted by AIDEN, obtained by summing the charge density over all grid points. The marginal histograms along the horizontal and vertical axes show the numerical distributions of the total charge. For fixed learned species profiles, integrating the periodic one-center sum gives an extensive contribution ∫Ωρinit()d3=∑iqZi _ _ init(r)\, d^3r= _iq_Z_i, where qZq_Z is the integral of one localized species profile. This branch can therefore supply an atomic baseline for the electron count without learning a global constraint. The remaining environment-dependent contribution is not constrained to have zero integral, and the present results do not isolate the contributions of the two branches. Fig. 10(a)-(c) complements the integrated errors with a spatial comparison of the reference and predicted densities. The density panels use four logarithmically spaced isovalues within their shared positive density range. The pointwise error in panel (c) is expressed relative to the local reference density, with the denominator bounded below by 10−1210^-12 times its maximum magnitude to stabilize the low-density region. Panel (d) then relates the integrated charge deviation to the cost of evaluating the complete density field. Figure 10: Real-space density and property comparison for TBG. (a), (b) DFT and AIDEN density isosurfaces, ρ()ρ(r) and ρ^() ρ(r), shown with a common logarithmic color scale and identical isovalues. (c) Pointwise relative density error δρ _ρ, with isosurfaces at 0.5%, 5%, and 10%. (d) Absolute electron-count deviation |ΔNe|| N_ e| (left axis) and full-grid density inference time (right axis). (e) Total-energy magnitude (left axis) and per-atom energy error |ΔE|| E| (right axis). (f) Mean absolute force component ⟨|Fiα|⟩ |F_iα| (left axis) and force-component MAE ⟨|ΔFiα|⟩ | F_iα| (right axis). Subplots (d)-(f) compare AIDEN and ChargE3Net across all three twist angles; blue and green shades correspond to the left and right axes, respectively. Table 5.: Density and charge conservation for TBG. Electron counts and their signed deviations are reported in units of e; relative deviations use the matched reference count. Twist angle Method ερ _ρ [%] Reference count Predicted count ΔNe N_ e Relative [%] 21.79∘21.79 AIDEN 0.293371 112 112.010667 +0.010667 0.009524 ChargE3Net 0.315961 112 112.035903 +0.035903 0.032056 13.17∘13.17 AIDEN 0.290678 304 304.024022 +0.024022 0.007902 ChargE3Net 0.314475 304 304.098970 +0.098970 0.032556 9.43∘9.43 AIDEN 0.289184 592 592.040188 +0.040188 0.006788 ChargE3Net 0.313427 592 592.178324 +0.178324 0.030122 Energies and forces. Tabs. 6 and 7 give the numerical results corresponding to Fig. 10(e), (f). AIDEN gives |ΔE|=0.753| E|=0.753, 0.747, and 0.741 meV/atom, approximately 28% of the corresponding ChargE3Net errors. The signed total-energy differences are negative for both models at every angle. The force-component MAEs remain approximately 0.021-0.022 eV/Å for AIDEN and 0.040-0.041 eV/Å for ChargE3Net. AIDEN also reduces the component RMSE and the maximum per-atom force-vector error in all three cells. These results show that its improvement in density accuracy is accompanied by consistent improvements in the fixed-density energy and force estimates. Table 6.: TBG total energies and energy errors. The DFT row gives the reused SCF reference. The signed difference E^−E E-E is in eV, while |ΔE|| E| follows the per-atom definition in Tab. 2. Twist angle Method Total energy [eV] E^−E E-E [eV] |ΔE|| E| [meV/atom] 21.79∘21.79 DFT -258.04561259 0 0 AIDEN -258.06669840 -0.02108581 0.753065 ChargE3Net -258.12164598 -0.07603339 2.715478 13.17∘13.17 DFT -700.39986467 0 0 AIDEN -700.45662887 -0.05676420 0.746897 ChargE3Net -700.60567255 -0.20580788 2.707998 9.43∘9.43 DFT -1363.97698248 0 0 AIDEN -1364.08670932 -0.10972684 0.741398 ChargE3Net -1364.36657021 -0.38958773 2.632350 Table 7.: TBG force magnitudes and errors. All entries are in eV/Å. The component RMSE averages squared errors over all atoms and Cartesian components; the last column is the maximum Euclidean norm of the force error over atoms. Twist angle Method ⟨|Fiα|⟩ |F_iα| ⟨|ΔFiα|⟩ | F_iα| Component RMSE Max. vector error 21.79∘21.79 DFT 0.030418 0 0 0 AIDEN 0.050602 0.021058 0.034872 0.062180 ChargE3Net 0.068266 0.041057 0.062076 0.114601 13.17∘13.17 DFT 0.029764 0 0 0 AIDEN 0.049782 0.022210 0.034972 0.065953 ChargE3Net 0.066146 0.041128 0.061317 0.114324 9.43∘9.43 DFT 0.029303 0 0 0 AIDEN 0.049377 0.021770 0.034843 0.067280 ChargE3Net 0.064642 0.040043 0.060362 0.113618 Band structures. For the band metrics in Tab. 8, we compare equal band indices at matched k points and retain reference states within EFermi±5E_ Fermi± 5 eV. This gives 4,450, 996, and 1,663 states for the three cells. Raw errors compare the eigenvalues directly; aligned errors apply one least-squares rigid shift per model and cell before computing the residuals. For the retained state set S, the shift added to the model eigenvalues is s=||−1∑(n,)∈(ϵn−ϵ^n)s=|S|^-1 _(n,k) ( _nk- ε_nk). Aligned errors use ϵ^n+s−ϵn ε_nk+s- _nk, measuring residual band shape after removal of one global offset. Raw errors retain the offset under the calculation’s energy-reference convention; they do not establish vacuum-referenced absolute levels. The DFT band calculation supplies EFermiE_ Fermi. Fig. 4 presents the spectra around this reference. Table 8.: TBG band errors within EFermi±5E_ Fermi± 5 eV. All numerical entries are in eV. The shift is added to the model eigenvalues. MAE, RMSE, and maximum absolute error are evaluated over the same reference-state mask. Raw Aligned Twist angle Method Shift MAE RMSE MAE RMSE Max. abs. 21.79∘21.79 AIDEN -0.179315 0.179315 0.179863 0.008415 0.014021 0.077684 ChargE3Net -0.318982 0.319104 0.324463 0.027352 0.059382 0.327034 13.17∘13.17 AIDEN -0.172100 0.172100 0.173577 0.007410 0.022592 0.512770 ChargE3Net -0.333283 0.333283 0.334107 0.010303 0.023454 0.294968 9.43∘9.43 AIDEN -0.168169 0.171246 0.184834 0.011468 0.076698 1.450108 ChargE3Net -0.319113 0.323720 0.325494 0.012928 0.064134 1.437400 AIDEN has the lower aligned band MAE at every angle. Its raw band MAEs are 171-179 meV (Tab. 8), with the rigid shift removing the dominant reference-energy offset. Its largest reduction in band MAE occurs at 21.79∘21.79 , from 27.35 to 8.42 meV. Interestingly, the ranking of the tail errors is less uniform: ChargE3Net gives a smaller maximum aligned error at 13.17∘13.17 , and a smaller aligned RMSE and maximum error at 9.43∘9.43 . Both models reach maximum individual-state errors of approximately 1.44-1.45 eV in the largest cell. Density inference efficiency. Tab. 9 lists the full-grid density inference times used in Fig. 10(d). Both model timings are measured on a single NVIDIA A100 GPU after warm-up, with CUDA synchronization and model loading excluded; ChargE3Net uses its official full-grid per-sample inference pipeline. AIDEN performs the atomic encoding once per structure and reuses it across grid queries. Its inference time increases from 13.27 to 41.43 s as NgN_ g increases from 2.76 to 15.05 million, whereas ChargE3Net requires 128.81-634.14 s. The measured speedups are 9.71, 10.85, and 15.31×, respectively. The reported speedups refer to these full-grid implementations on the specified hardware. Table 9.: Full-grid charge density inference times for TBG. Times refer to AIDEN’s full-grid inference and ChargE3Net’s official full-grid per-sample inference. The speedup is the ratio of ChargE3Net to AIDEN time. Twist angle AIDEN [s] ChargE3Net [s] Speedup 21.79∘21.79 13.270 128.813 9.71×9.71× 13.17∘13.17 28.933 314.058 10.85×10.85× 9.43∘9.43 41.430 634.141 15.31×15.31× C.4 DFT calculations for water, AlMg, and a-Si C.4.1 Structures and DFT settings We use the water and AlMg structures and one a-Si configuration from the OOD examples of Ref. (Qin et al., 2026), with Na=192N_ a=192, 108, and 64, respectively. The a-Si configuration is frame 11. The supplied files determine the structures and FFT grids; reference densities and properties are obtained from new PBE-SCF calculations. The density grids are 2163216^3 for water and AlMg and 1683168^3 for a-Si. All calculations use VASP with PAW-PBE potentials for H, O, Si, Al, and Mg, a 520 eV cutoff, a 2×2×22× 2× 2 Monkhorst-Pack mesh, and Gaussian smearing of 0.05 eV. The Mg potential is Mg_pv (13Apr2007), with eight valence electrons including 2p63s22p^63s^2. We use non-spin-polarized calculations, disable symmetry, and hold atomic positions fixed. The electronic convergence threshold is 10−810^-8 eV, with PREC=Accurate, LREAL=.FALSE., LASPH=.TRUE., ADDGRID=.TRUE., and LMAXMIX=2. Structure, potentials, grids, and electronic settings are matched across methods for each system. C.4.2 Density and property evaluation AIDEN uses its PBE-pretrained parameters, and ChargE3Net uses the public Materials Project pretrained checkpoint. SAD is generated with ICHARG=12. Each candidate density is subsequently evaluated with ICHARG=11, keeping the density fixed throughout electronic minimization. Model densities are used without charge renormalization and supplemented with the one-center PAW data from the matched PBE-SCF reference; SAD retains its own one-center data. Energies are the final VASP TOTEN values, and forces are the frozen-density approximate forces reported by VASP. C.4.3 Cross sections For Fig. 3, candidate x-normal planes are the FFT planes nearest to atomic centers. We select the plane with the most atoms within 1.5 grid spacings, using proximity to the centers to break ties and periodic distances throughout. The selected zero-based indices are 191, 36, and 97 for water, AlMg, and a-Si, containing 9, 13, and 5 nearby atoms, respectively. AIDEN/ChargE3Net/SAD section NMAEs are 1.887/1.641/12.103% for water, 6.673/6.683/2.459% for AlMg, and 1.302/1.618/9.772% for a-Si. The density panel and error panels use separate logarithmic color scales, with a common error scale across methods for each system. Display clipping is applied only when plotting. C.4.4 Timing Model inference is measured on a single NVIDIA A100 GPU after warm-up, with CUDA synchronization and model loading excluded. Both AIDEN and ChargE3Net use full-grid per-sample inference, with ChargE3Net evaluated through its official inference pipeline. VASP runs use 16 CPU MPI processes with one thread per process. The SAD time measures the complete ICHARG=12 calculation, including diagonalization. The model-to-model speedup ratios compare the evaluated AIDEN and ChargE3Net full-grid implementations. C.4.5 AlMg near-core discrepancy The reference density reaches 15.32 e/Å3 near Mg, whereas the corresponding maxima are 18.75 and 18.40 e/Å3 for AIDEN and ChargE3Net, respectively. Their ten largest pointwise errors all lie within 0.062 Å of Mg centers. Within 0.4 Å of Mg atoms, the local NMAEs are approximately 7.63/7.62/0.37% for AIDEN/ChargE3Net/SAD, compared with 2.72/2.87/46.47% around Al. The dominant discrepancy of the learned models is therefore strongly localized around Mg rather than distributed over the bonding or interstitial regions. This localization is particularly relevant because the Mg_pv PAW potential treats the 2p63s22p^63s^2 shell as valence, introducing a strongly localized 2p2p semicore contribution to the smooth grid density. The different behavior of the density and energy errors can be understood from their spatial weighting. Let δρ()δρ(r) denote the density error. The global NMAE measures its integrated magnitude, whereas the electron-nuclear energy couples the same error to the external potential, δρ()=ρ^()−ρ() δρ(r)= ρ(r)-ρ(r) ,ϵρ=1Ne∫Ωd3r|δρ()|, , _ρ= 1N_e _ d^3r\,|δρ(r)|, δEen=∫Ωd3rvext()δρ() δ E_ en= _ d^3r\,v_ ext(r)δρ(r) ,vext()=−∑IZI|−I|. , v_ ext(r)=- _I Z_I|r-R_I|. Thus, two density errors with comparable contributions to ϵρ _ρ need not contribute equally to the energy. In an all-electron picture, vextv_ ext varies most strongly close to the nuclei, so a spatially compact near-nuclear error can have a disproportionately large energetic effect. In the PAW calculations considered here, the learned target is the smooth grid density rather than the complete all-electron density, while the fixed-density calculation additionally uses atom-centered PAW one-center information. Nevertheless, the same principle remains relevant: the sensitivity of a derived observable to δρ()δρ(r) is highly nonuniform in real space. This behavior closely parallels the analysis of Lewis et al. (Lewis et al., 2021). For periodic Al and Si, they found that density-representation errors below 0.1% could still produce electrostatic-energy errors much larger than suggested by the density metric alone. Their energy decomposition showed that the corresponding Hartree-energy errors were comparatively small and that the dominant contribution arose from the electron-nuclear interaction, which is especially sensitive to density inaccuracies close to the nuclei. They further observed that the sharp near-nuclear density is largely determined by atomic identity and introduced an atom-centered isotropic baseline so that the learned component could focus on chemically induced variations of the density. This is conceptually related to the decomposition used in AIDEN, where ρ^init ρ_ init captures an element-dependent one-center contribution and ρ^env ρ_ env describes environment-induced redistribution. The AlMg results exhibit the same qualitative mechanism. AIDEN and ChargE3Net reproduce most of the full-grid density with similar NMAEs, but their largest errors occur precisely in the Mg near-core region containing the localized semicore contribution. By contrast, SAD gives a much smaller Mg-local NMAE of approximately 0.37%, although its Al-centered and environment-dependent density is substantially less accurate. Since SAD is constructed from atomic densities, it naturally preserves the strongly atom-centered part of the Mg density more faithfully while lacking an accurate description of bonding-induced redistribution. Together with the different PAW one-center data retained in the SAD calculation, this explains why SAD can yield a poorer global density NMAE but a substantially smaller fixed-density energy error. The AlMg case therefore illustrates a general limitation of global density metrics. The NMAE measures the integrated accuracy of the scalar field but does not distinguish spatial regions according to their influence on a particular electronic observable. For AlMg, the dominant model error is concentrated in a small Mg near-core region that carries comparatively little weight in the global density integral but substantially greater energetic sensitivity. The agreement with the mechanism identified in Ref. (Lewis et al., 2021) supports evaluating learned charge densities jointly through global density errors, spatially resolved diagnostics, and downstream electronic-structure observables. Appendix D The use of large language models In this work, we used large language models to assist with language polishing and manuscript editing.