Learning ab initio phase-field models

Ab initio in both derivation and learning

We develop ab initio phase-field (AIPF) models, grounded in microscopic theory and data. The term ab initio is used in the same spirit as in ab initio molecular dynamics (AIMD). In AIMD,

1. the equations of nuclear motion follow from quantum mechanics under the Born–Oppenheimer approximation;
2. the forces are computed using density functional theory, where the exchange–correlation functional is chosen from a ladder of systematically improvable approximations known as Jacob's ladder.

So what makes a phase-field model ab initio? Hover a column of the figure below.

Derivation and learning

Ab initio in both its derivation and its learning. Both start from quantum mechanics.

Hover a column.

Derivation

1. Quantum mechanics to molecular dynamics under the Born–Oppenheimer approximation
2. Molecular dynamics to a generalized Langevin equation through the Mori–Zwanzig formalism
3. The generalized Langevin equation to the ab initio phase-field model under three assumptions
   • Markovian dynamics
   • Local mobility
   • Gaussian white noise

Learning

1. MLIP-driven molecular dynamics of ab initio accuracy
2. Free energy and mobility parametrized by neural networks, with their forms chosen from the ladders

The ladders

The excess free energy carries all interparticle correlations. From exact statistical mechanics down to phenomenology, each lower rung restricts its form further. Hover or tap a rung to see its functional form.

Exact statistical mechanics

$$\mathcal F=\mathcal F_{\mathrm{exc}}+\mathcal F_{\mathrm{id}}+\mathcal F_{\mathrm{ext}}$$ $$\mathcal F_{\mathrm{id}}=k_{\mathrm B}T\sum_\alpha\int\rho_\alpha\big[\ln(\rho_\alpha\Lambda_\alpha^3)-1\big]\,\mathrm d\mathbf r,\qquad \mathcal F_{\mathrm{ext}}=\int\mathbf V_{\mathrm{ext}}\cdot\boldsymbol\rho\,\mathrm d\mathbf r$$

The ideal-gas part $\mathcal F_{\mathrm{id}}$ and the external-potential part $\mathcal F_{\mathrm{ext}}$ are known exactly, with $\Lambda_\alpha$ the thermal wavelength. The excess part $\mathcal F_{\mathrm{exc}}$ accounts for interparticle interactions, and every rung below is an approximation of it.

Rung 4 · Neural operator

$$\mathbf v_0=\boldsymbol\rho,\qquad \mathbf v_{l+1}=\sigma\big(\mathbf k_l\ast\mathbf v_l+\mathbf W_l\mathbf v_l+\mathbf b_l\big)$$ $$\mathcal F_{\mathrm{exc}}=\int\mathbf k_L\,\mathbf v_L$$

A neural operator, with layer fields $\mathbf v_l$, convolution kernels $\mathbf k_l(\mathbf r-\mathbf r',T)$, channel-mixing weights $\mathbf W_l$, biases $\mathbf b_l(T)$ and a pointwise nonlinearity.

Rung 3 · Nonlocal kernel this work

$$\begin{aligned}\mathcal F_{\mathrm{exc}}([\boldsymbol\rho],T)&=\int f_{\mathrm{loc}}(\boldsymbol\rho,T)\,\mathrm d\mathbf r\\&\quad+\frac12\sum_{\alpha\beta}\iint\rho_\alpha(\mathbf r)\,W_{\alpha\beta}(\lvert\mathbf r-\mathbf r'\rvert,T)\,\rho_\beta(\mathbf r')\,\mathrm d\mathbf r\,\mathrm d\mathbf r'\end{aligned}$$

Here $f_{\mathrm{loc}}$ is a function of the local densities and $W_{\alpha\beta}$ the radial kernel between species $\alpha$ and $\beta$. The kernel extends the square-gradient term of rung 2 to the full wavevector dependence of the pair correlations.

Rung 2 · Cahn–Hilliard

$$\mathcal F_{\mathrm{exc}}=\int\big[f_{\mathrm{loc}}(\boldsymbol\rho,T)+\kappa\,\lvert\nabla\boldsymbol\rho\rvert^2\big]\,\mathrm d\mathbf r$$

The Cahn–Hilliard square-gradient form with coefficient $\kappa$.

Rung 1 · Basis expansion

$$\mathcal F_{\mathrm{exc}}=\int\Big[\sum_n c_n(T)\,\phi_n(x)+\kappa\,\lvert\nabla x\rvert^2\Big]\,\mathrm d\mathbf r$$

A prescribed basis expansion in composition $x$, with basis functions $\phi_n$ such as Landau polynomials or Redlich–Kister terms and coefficients $c_n(T)$.

Two examples

Landau $$\mathcal F_{\mathrm{int}}=\int\Big[\tfrac a2(T-T_{\mathrm c})\big(x_A-\tfrac12\big)^2+b\,\big(x_A-\tfrac12\big)^4+\tfrac\kappa2\lvert\nabla x_A\rvert^2\Big]\,\mathrm d\mathbf r$$

Flory–Huggins $$\mathcal F_{\mathrm{exc}}=\int\Big[w\,x_A(1-x_A)+\tfrac\kappa2\lvert\nabla x_A\rvert^2\Big]\,\mathrm d\mathbf r$$

Exact statistical mechanics

The exact Mori–Zwanzig projection gives a generalized Langevin equation with a memory kernel, which describes how the density at the current time is affected by its past. The mobility is the time integral of this kernel, and every rung below approximates it.

Rung 4 · Nonlocal mobility

The mobility couples every pair of points and depends on the whole density field.

Rung 3 · Local mobility this work

$$\mathbf M_\theta=\mathbf C_\theta^{1/2}\,\mathbf D^{1/2}\,\mathbf L_\theta\mathbf L_\theta^{\top}\,\mathbf D^{1/2}\,\mathbf C_\theta^{1/2}$$

• $\mathbf D=\mathrm{diag}(\boldsymbol\rho)$, the species densities
• $\mathbf L_\theta(\boldsymbol\rho)$, lower triangular, parametrized by a multilayer perceptron (MLP)
• $\mathbf C_\theta(T)=\mathrm{diag}\big(c_{\alpha,\theta}(T)\big)$, Arrhenius factors with one learned activation energy $E_\alpha$ per species: $$c_{\alpha,\theta}(T)\propto e^{-E_\alpha/k_{\mathrm B}T}$$

The mobility depends only on the densities and temperature at the same point. The construction makes $\mathbf M_\theta$ symmetric and positive semidefinite.

Rung 2 · Prescribed density dependence

$$\mathbf M=c(T)\,\mathbf g(\boldsymbol\rho)$$

The density dependence $\mathbf g$ is fixed in advance, and only a few constants are fitted.

Examples

Ideal mobility, as in dynamical density functional theory $$M_{\alpha\beta}=\delta_{\alpha\beta}\,\frac{D_\alpha\,\rho_\alpha}{k_{\mathrm B}T}$$

Degenerate Cahn–Hilliard mobility $$M=M_0\,x(1-x)$$

Rung 1 · Constant mobility

$$M=M_0,\qquad M_0=\text{const}$$

A constant mobility, fixed separately or prescribed, as in standard phase-field models. It sets the relaxation rate but carries no dependence on density or composition.

Results

Binary Lennard-Jones mixture

A preliminary check of the framework on a symmetric binary Lennard-Jones mixture under overdamped Langevin dynamics. Below the critical temperature the mixture demixes spontaneously within the spinodal region. The AIPF model is compared with two rung 1 baselines, Landau and Flory–Huggins, trained on the same data with the same loss and mobility parametrization.

Dynamics

Iron–boron melt

FeB₄ is predicted to be thermodynamically stable, but its synthesis has been reported only under high pressure. One route to solid-state synthesis proceeds through a transient melt at the interface between two solid precursors, which requires the melt at the target composition to be locally stable. A single AIPF model of the Fe–B melt is learned at 0, 5 and 10 GPa, and the curvature of its Gibbs free energy of mixing is quantified by the thermodynamic factor

with $\Gamma>0$ for a locally stable melt and $\Gamma<0$ inside the spinodal region.

By learning from the dynamics, the AIPF model accesses the spinodal region directly, resolving the negative curvature inaccessible to equilibrium fluctuation measurements.

Hydrogen–helium under planetary conditions

Hydrogen–helium mixtures make up most of the mass of gas giants such as Jupiter and Saturn. Inside them, helium is hypothesized to become immiscible in metallic hydrogen and settle as droplets toward the planetary center, a process known as helium rain. Whether and where it occurs depends on the immiscibility boundary, and its evolution on the kinetics of droplet nucleation, growth and transport. A single AIPF model spans 200 to 800 GPa and 2000 to 12 000 K.

Thermodynamics

Helium rain in gas giants

The model captures droplet nucleation and growth in helium-rain simulations of a domain corresponding to 2.2 million atoms, far beyond the scale of atomistic modeling at comparable accuracy.

Dynamics simulation

Citation

If you use this code, please cite the paper.

@article{chen2026aipf,
  title   = {Learning ab initio phase-field models},
  author  = {Chen, Mengyi and Zhong, Peichen and Zhang, Zihan and Li, Qianxiao},
  journal = {arXiv preprint arXiv:2610.01432},
  year    = {2026},
  doi     = {10.48550/arXiv.2610.01432},
  url     = {https://arxiv.org/abs/2610.01432}
}