Equivariant learning of a transferable three-dimensional classical density functional

Equivariant learning of a transferable three-dimensional classical density functional

Bingqing Cheng

arXiv:2608.13506v1 [cond-mat.stat-mech] 13 Aug 2026

Summary

This paper introduces Equi-cDFT, an Energy-first, EQUIvariant framework for learning the excess free-energy functional Fexc[ρ, T] of a classical fluid as an extensive scalar on a three-dimensional Cartesian grid. The central challenge addressed is that classical density functional theory (cDFT) offers a reusable variational description of liquids, but its central excess free-energy functional is generally unknown, and learned approximations have largely remained restricted to planar or lower-dimensional settings. The authors show that this functional can be learned directly from fully three-dimensional equilibrium density fields while preserving spatial symmetry and variational consistency, without free-energy or chemical-potential labels.

Theory and architecture

The framework learns Fexc[ρ, T] as a function of temperature and a three-dimensional density field represented on a regular grid with voxel volume ΔV = (ΔL)3. The local environment χg contains the density within a finite neighborhood of grid point g. The excess functional is approximated as:

Fexc[ρ, T] = ΔV Σg ρg aexc(χg, T)

using the same local free-energy map aexc at every grid point, making the functional extensive and transferable between system sizes.

The local three-dimensional environments are represented using symmetry-adapted features based on the Cartesian atomic cluster expansion (CACE) adapted to a density grid. For a target grid point rg, the Cartesian moments are:

Al(rg) = Σq≤qcut ρg+q qx lx qy ly qz lz

where l = (lx, ly, lz) specifies Cartesian powers. These moments are covariant under the cubic point group Oh, which consists of the 3! axis permutations combined with 23 independent axis reflections, giving 48 operations in total. For K = (l1,..., lν), the invariant feature is obtained by summing the transformed products over Oh:

BK(ν)(rg) = ΣR∈Oh Πm=1 ν slm(R) A Rlm (rg)

The invariant features, together with the center density and temperature, are passed to a shared scalar readout to obtain aexc(χg, T) at every grid point.

Derivative learning from local chemical-potential balance

The excess free energy is learned through its local density gradient, computed by automatic differentiation:

cg(1)([ρ], T) = - (β/ΔV) ∂Fexc[ρ, T]/∂ρg

where c(1) is the one-body direct correlation functional. For an equilibrium density field in either the canonical or grand-canonical ensemble, the dimensionless local chemical potential is:

µloc g = ln(Λ3ρg) + βVext,g - cg(1)([ρ], T)

The equilibrium condition requires this quantity to be independent of position. Equi-cDFT uses this balance condition directly as the learning signal. For a canonical training density field, the constant chemical potential is an unknown Lagrange multiplier, which is eliminated analytically as the spatial mean of the predicted local chemical potential. The per-field loss is:

L∆µ = Σg (µloc g - µ̄)2

where both the sum and the average are taken over the reliably sampled region of the density field. Canonical balance leaves one unavoidable gauge: the loss is unchanged under Fexc[ρ, T] → Fexc[ρ, T] + b(T)N. This added term is constant under fixed-N density variations and has zero second functional derivative, leaving fixed-N predictions and c(2) unchanged while shifting the absolute chemical potential by b(T).

Results

Equation of state, phase behavior, and interfaces

The learned model behaves as a thermodynamic functional rather than merely reproducing equilibrium fields under external potentials. The training loss contains no explicit supervision for bulk response functions, pressure, phase coexistence or free interfaces.

The two-body direct correlation function c(2) is obtained via differentiation:

c(2)(rg, rg′; [ρ], T) = (1/ΔV) ∂c(1)(rg; [ρ], T)/∂ρ(rg′)

For a homogeneous fluid, the static structure factor is obtained via the Ornstein–Zernike relation:

S(k) = 1 / (1 - ρ ĉ(2)(k))

The resulting low-wavevector structure factors are compared with direct homogeneous-fluid MD across three representative states, reproducing both the magnitude and wavevector dependence of S(k).

The zero-wavevector response determines the inverse isothermal compressibility through ∂βP/∂ρ = S(0)(-1) = 1 - ρĉ(2)(0; ρ, T). Integrating from the ideal-gas limit gives the complete equation of state. The dashed predictions follow the MD-derived EOS. At subcritical temperatures T = 0.8 and 1.0, the learned homogeneous isotherms develop a van der Waals loop, including a mechanically unstable region with ∂P/∂ρ < 0.

The learned functional recovers the vapor–liquid phase diagram. Coexistence densities are obtained directly from fixed-N slab solutions through T = 1.05. The extrapolated critical point is TccDFT = 1.09 and ρcDFT c = 0.31 σ−3, agreeing well with the previous direct-simulation estimate (1.08, 0.32). The learned functional also captures the progressive broadening of the liquid–vapor interface as the critical region is approached across four temperatures.

Solvent-mediated colloid bridging

The mean force between two colloids mediated by the solvent is obtained directly from the equilibrium solvent density through the classical Feynman–Hellmann identity:

⟨F⟩D = -ΔV Σg ρD,g ∂Vext,g(D)/∂D

where ρD,g is the equilibrium density at grid point g. The learned functional predicts the interaction through its equilibrium density alone; negative force denotes attraction and positive force repulsion.

The density fields expose the physical origin of the interaction. At D/σ = 5.0, the solvent-depleted regions surrounding the two cores merge into a continuous neck, producing attraction. At intermediate separation the bridge ruptures and the opposing solvation layers reorganize, generating a pronounced repulsive maximum. At large separation the two density perturbations become independent and the force decays to zero. Equi-cDFT follows the strong attraction at D/σ = 4.5, the sign reversal and repulsive maximum near D/σ = 6.5, and the decay to zero by D/σ = 9–10, while resolving the interaction on a much denser separation grid than sampled directly by MD.

The mean force is evaluated in three ways: direct MD provides the reference ensemble average; applying the Feynman–Hellmann quadrature to the gridded MD density reproduces that force; replacing it with the independently minimized Equi-cDFT density leaves the force curve essentially unchanged. The effective interaction is therefore an emergent prediction of the learned free-energy landscape rather than a separately fitted observable.

Adsorption in a soft porous gyroid

The gyroid makes the test especially demanding: its continuously curved, interconnected pore has no unique wall-normal coordinate and is topologically unlike the external fields used for training. The soft gyroid wall drives the density from nearly complete depletion near the excluded region to dense, connected domains in the pore interior. At the fixed loading shown, Equi-cDFT reproduces this full three-dimensional morphology, including the curved depletion layer and the spatially varying accumulation throughout the interconnected pore. The agreement extends from the strongly perturbed interfacial region into the pore interior.

Adsorption equilibrium requires the confined fluid and homogeneous reservoir to have the same chemical potential. For each measured mean loading, only the particle number is supplied to the fixed-N Equi-cDFT calculation; the reference density and imposed reservoir chemical potential are withheld. After a single bulk calibration required by the canonical gauge, the same shift is applied to the confined branch. Eliminating chemical potential between the independently calculated bulk and confined branches produces the adsorption isotherm without fitting a separate adsorption model.

At T = 1.2, above the bulk critical point, filling is continuous rather than a capillary phase transition. It is nevertheless strongly nonlinear: the pore remains weakly populated in a dilute reservoir and fills rapidly as the reservoir enters the dense-fluid regime.

Benchmark evaluation

The held-out test set contains 1,592 complete fields at T = 0.7, 1.1 and 1.5, temperatures omitted entirely from training and validation. For the forward test, every held-out density field is supplied to the model, and automatic differentiation gives c(1) and hence the external-field prediction. The inverse test removes the reference density and supplies only Vext, T and N, using an adaptive fixed-N Euler solver starting from a uniform density.

Transfer to larger canonical systems is tested by applying the model trained on L = 8σ fields to L = 12σ boxes discretized on 243 grids at the same three held-out temperatures, spanning N = 32–1568 and whole-box mean densities ρσ3 = 0.02–0.91.

Grand-canonical chemical-potential differences are assessed using a GCMC test set containing 22 complete fields at each of T = 1.0, 1.2 and 1.4. The raw prediction is the masked spatial mean of ln(ρΛ3) + βVext - c(1), and reference and predicted values are centered independently at each temperature, corresponding to one additive gauge per temperature and no statewise fitting.

Discussion

Equi-cDFT combines three advances: (1) it learns an explicit scalar excess free-energy functional rather than a direct potential–density map or an independent approximation to c(1); (2) it represents each local three-dimensional density environment with features invariant under the cubic point group Oh, ensuring that the predicted scalar is unchanged and its functional derivatives transform equivariantly under grid rotations and reflections; (3) it uses local chemical-potential balance as the training signal, so canonical equilibrium density fields constrain one common free-energy landscape without requiring free-energy or chemical-potential labels.

The benchmarks show both accuracy and generalization. At temperatures excluded from training and validation, Equi-cDFT accurately recovers external fields from densities and reconstructs densities by fixed-N minimization. Without retraining, it transfers to larger cells and from canonical training to grand-canonical chemical-potential differences. The recovery of structure factors, the compressibility-route equation of state, phase coexistence and interfacial broadening provides a complementary test of emergent response and thermodynamics, because none of these quantities entered the loss.

The largest remaining sensitivities occur in long-wavelength finite-cell response and near criticality; the reported critical point is a mean-field continuation rather than a resolution of critical fluctuations. The present model also uses a fixed real-space discretization: the larger-cell benchmark tests transfer in volume at unchanged grid spacing, while transfer across spatial resolution and multiscale representations remain open directions.

The present LJTS functional is interaction-specific, but the construction suggests several direct extensions: dynamical cDFT simulations, mixtures with species-resolved input channels, ionic and polar fluids requiring an additional long-ranged branch, and reference density fields generated by machine learning interatomic potentials trained on electronic-structure energies and forces.

Methods

The production model represents the excess free energy as the sum of a pointwise baseline and a finite-range lattice-CACE correction:

Fexc,θ[ρ, T] = ΔV Σg ρg [aloc θ(ρ̃g, T̃) + aCACE θ(ρ̃g, Bg, T̃)]

where ρ̃g and T̃ denote network inputs normalized by their corresponding dataset means. The production representation uses a spherical neighbor stencil with a cutoff of three grid spacings, containing the 123 integer offsets satisfying q2 ≤ 32. The center is excluded from the moment sums and supplied separately to the readout, leaving 122 noncentral neighbors. A single constant radial channel and the 20 raw Cartesian monomials with total degree l ≤ 3 are used. Products of the Cartesian moments are retained through correlation order ν = 2 and averaged over the 48 signed axis permutations of the cubic grid. This gives 15 invariant CACE features: two at first order and thirteen at second order. The CACE readout receives 17 inputs: the center density, the 15 invariants and temperature. Its multilayer perceptron has dimensions 17 → 32 → 16 → 1. The pointwise baseline receives the local density and temperature and uses a 2 → 32 → 16 → 1 network. Both networks use SiLU hidden activations and linear scalar outputs, for a total of 1,762 trainable parameters.

Training data were generated for a one-component fluid interacting through the truncated-and-shifted Lennard–Jones (LJTS) potential with cutoff 2.5σ. All training simulations used the canonical ensemble in cubic periodic cells of side L = 8σ. The external potentials were independently randomized superpositions of periodic one-, two- and three-dimensional Gaussian components. Density fields were represented on a 163 grid with spacing 0.5σ. Particle numbers ranged from N = 8 to 464, corresponding to whole-box mean densities from 0.02 σ−3 to 0.91 σ−3. The final corpus contains 10,957 complete fields and spans the reduced temperatures T = 0.625 to 1.8. All retained configurations are liquid; trajectories exhibiting crystalline ordering were removed.

All 1,592 fields at T = 0.7, 1.1 and 1.5 were assigned to the test set to measure interpolation across temperature. Within every remaining source file, 10% of the fields were assigned to validation, giving 8,427 training fields and 938 validation fields.

Training uses complete-field batches of two, Adam optimization with an initial learning rate of 10−4, a reduce-on-plateau scheduler halving the learning rate after three epochs without validation improvement down to a minimum of 10−6, and runs for 200 epochs. The checkpoint with the lowest validation local-chemical-potential loss is used throughout. On a single NVIDIA L40 GPU, the 200-epoch production fit required approximately 1 h 40 min.

The model, automatic functional derivatives and equilibrium solvers are implemented in Python using PyTorch. The implementation is publicly available at https://github.com/BingqingCheng/equicdft. On a four-thread Apple M2 CPU, evaluation of c(1) required 7–110 ms for the 163–403 grids considered, while representative fixed-N equilibrium solves required 0.3–8 s.

Based on this paper, here are the specific improvements you can make to AI systems and what the improved systems can do:

Improvement: Replace loss functions based on labeled outputs with physics-based balance conditions (e.g., local chemical-potential uniformity) that are computed via automatic differentiation of a learned scalar functional. This eliminates the need for expensive labels (free energies, chemical potentials) while enforcing variational consistency.

Improved AI capability: The system can learn thermodynamic functionals directly from equilibrium density fields without any energy or chemical-potential labels, making it applicable to systems where such labels are experimentally or computationally inaccessible.

Improvement: Use group-averaged invariant features (here, under the cubic point group Oh) built from covariant Cartesian moments of local fields. This ensures the learned scalar is invariant and its functional derivatives transform equivariantly under grid rotations/reflections.

Improvement: Represent the target property as an integral over local contributions (here, Fexc = ΔV Σ ρg aexc(χg, T)) with a shared local map, making the functional extensive and size-transferable.

Improvement: Train only on first-order derivatives (c(1)) but extract higher-order derivatives (c(2)), thermodynamic response functions (compressibility), phase coexistence, and interfacial properties via automatic differentiation of the learned functional.

Improvement: Train in the canonical ensemble (fixed N) but predict grand-canonical properties (chemical-potential differences, adsorption isotherms) by analytically eliminating the unknown Lagrange multiplier and applying a single additive gauge correction per temperature.

Improvement: Include temperature as an explicit input to the functional, enabling interpolation across thermal states without retraining.

Improvement: Replace direct molecular dynamics or Monte Carlo with a differentiable neural functional that evaluates c(1) in 7–110 ms and solves equilibrium densities in 0.3–8 s on a CPU.

Improvement: Use grid-based Cartesian moments with finite cutoffs, making the representation applicable to any 3D scalar field (not just particle densities), including charge densities, magnetization, or concentration fields.