PaperlayerView exact revision on arXiv

Source-linked AI summary

Learning ab initio phase-field models

Mengyi Chen, Peichen Zhong, Zihan Zhang, Qianxiao Li

arXiv:2610.01432v1cond-mat.stat-mechcond-mat.mtrl-scics.LG

TL;DR

Atomistic simulations are limited in accessible length and time scales, while conventional equilibrium approaches struggle to characterize unstable Fe–B melts. This work derives a mesoscopic model from molecular dynamics and learns its free energy and mobility, predicting pressure-dependent Fe–B stability and helium-droplet evolution.

  • Problem

    Atomistic simulations have limited access to the length and time scales of helium-rain dynamics, and equilibrium methods struggle to determine Fe–B free-energy curvature inside an unstable melt.

    Source: Iron–boron melt · p. 4, Hydrogen–helium mixture · p. 5

  • Method

    The framework reduces a projected molecular-dynamics equation under Markovian, local-mobility, and Gaussian-noise assumptions, then jointly learns free energy and mobility from dynamics.

    Source: Ab initio phase-field models · p. 3, Discussion · p. 7

  • Results

    The model predicts FeB4 is spinodally unstable at ambient pressure but outside the spinodal region at 10 GPa, and simulates helium-droplet nucleation, growth, and redistribution.

    Source: Iron–boron melt · p. 5, Hydrogen–helium mixture · p. 6

  • Takeaways & Limitations

    The framework points toward mesoscopic thermodynamic and dynamical predictions in regimes where experiments are difficult and direct molecular-dynamics simulations are prohibitively costly.

    Source: Discussion · p. 7

    The present model resolves density fields only, so helium accumulation under gravity occurs through diffusion rather than hydrodynamic droplet settling.

    Source: Discussion · p. 7

Abstract

from arXiv · show

Simulating microstructure evolution requires quantum-mechanical accuracy and mesoscopic reach in length and time scales, a combination that no current method achieves. Classical phase-field models provide this reach, but their accuracy is limited by phenomenological free energies and mobilities. Here we develop a framework for learning ab initio phase-field models, where the mesoscopic equation is not postulated but derived from a Mori-Zwanzig projection of molecular dynamics onto species-density fields under explicit assumptions. The nonlocal free energy and mobility left unspecified by this equation are parametrized by neural networks and learned from short molecular dynamics trajectories generated with machine-learning interatomic potentials of ab initio accuracy. We demonstrate the framework on an iron-boron melt and on hydrogen-helium mixtures under planetary conditions. For iron-boron, the model shows that the melt at the FeB$_4$ composition is spinodally unstable at ambient pressure but stabilized at 10 GPa, offering a thermodynamic rationale for why FeB$_4$ has been synthesized only under high pressure. For hydrogen-helium, the model predicts the immiscibility boundary and captures droplet nucleation and growth in helium-rain simulations of a column corresponding to 2.2 million atoms, far beyond the scale of atomistic modeling at comparable accuracy. Trained across compositions and conditions, such models could provide a mesoscopic counterpart to ab initio molecular dynamics.

Ab initio phase-field models

The framework derives a practical mesoscopic equation by projecting molecular dynamics onto band-limited species-density fields and making explicit assumptions about memory, mobility, and noise. Neural networks learn the remaining free-energy and mobility terms from short trajectories, enabling thermodynamic and dynamical predictions.

Ab initio phase-field models

  • An exact Mori–Zwanzig projection onto band-limited species-density fields yields a generalized Langevin equation with memory-dependent response and fluctuating forces.The density fields retain long-wavelength variations associated with slow mesoscopic dynamics.Source: Ab initio phase-field models · p. 2
  • Markovian dynamics, local-density-dependent mobility, and conditionally Gaussian fluctuating forces reduce the generalized Langevin equation to a practical mesoscopic model.The noise covariance is fixed by the fluctuation–dissipation relation.Source: Ab initio phase-field models · p. 3
  • The excess free energy uses a local density function and nonlocal pair kernel, while the local mobility is separable and may include an Arrhenius temperature factor.The rung-3 free-energy form balances nonlocal expressiveness with data efficiency and a stronger physical prior than rung 4.Source: Ab initio phase-field models · p. 3
  • Short MD trajectories from homogeneous mixtures and composition slabs train the model on phase separation and interface relaxation, with measured density evolution matched to the equation’s deterministic part.The density-evolution interval is chosen to exceed the microscopic memory time while remaining short enough that the drift changes little; equilibrium fluctuations further constrain mobility and free-energy curvature.Source: Ab initio phase-field models · pp. 3–4

Results

In a binary Lennard-Jones mixture, the AIPF model reproduces MD phase-diagram and domain-growth behavior while outperforming Landau and Flory–Huggins baselines. The results are presented as a preliminary test of the framework, alongside the motivation to understand melt conditions relevant to FeB4 synthesis.

Results

  • The binary Lennard-Jones mixture demixes spontaneously within the spinodal region below its critical temperature, while separation between binodal and spinodal boundaries requires nucleation.The phase diagram is derived from the learned bulk free energy; its common tangent gives binodal compositions and negative curvature identifies the spinodal.Source: Results · p. 4
  • The AIPF model agrees closely with MD on the phase diagram and spinodal-decomposition domain size, outperforming Landau and Flory–Huggins baselines.The model is also compared with MD and both baselines for equilibrium composition profiles under external potentials and interface relaxation in slabs.Source: Results · p. 4
  • FeB4 is predicted thermodynamically stable, but reported synthesis has occurred only under high pressure, motivating attention to the precursor melt’s stability.A proposed route forms the target crystal by nucleation from a transient melt at the interface between solid precursors.Source: Iron–boron melt · p. 4

Iron–boron melt.

The AIPF model learns the Fe–B melt’s free-energy curvature from dynamics, including in the unstable spinodal region that conventional equilibrium approaches struggle to resolve. It predicts pressure- and temperature-dependent stability, including stabilization of FeB4 at 10 GPa.

Iron–boron melt.

  • Conventional approaches either require difficult-to-determine nonideal mixing entropy or yield nonnegative curvature estimates that approach zero only with limited resolution.They provide results at discrete state points, so continuous phase boundaries require fitting a prescribed free-energy form.Source: Iron–boron melt · p. 4
  • Learning from phase-separation dynamics gives the model direct access to negative free-energy curvature in the spinodal region, which equilibrium fluctuation measurements cannot resolve directly.A single model is learned for 0, 5, and 10 GPa using the MACE potential; pressure enters implicitly through species densities and temperature.Source: Iron–boron melt · pp. 4–5
  • At 10 GPa, FeB4 at 1500 K lies outside the predicted spinodal region, unlike at 0 and 5 GPa; increasing temperature also stabilizes the melt.The model agrees with MD spinodal classifications and predicts the sign of Γ for states whose MD stability classification is unresolved.Source: Iron–boron melt · p. 5, Figure 4 · p. 5

Hydrogen–helium mixture.

The learned model predicts hydrogen–helium immiscibility and simulates helium rain at mesoscopic scales. Saturn’s profiles intersect the predicted immiscibility region, while Jupiter’s lie above but close to it.

Hydrogen–helium mixture

  • A single model trained from MD trajectories at 200–800 GPa predicts the immiscibility boundary across helium fraction, pressure, and temperature.The figure caption says predictions extend to 150–900 GPa.Source: Hydrogen–helium mixture · p. 5, Figure 5 · p. 6
  • Helium rain is expected where a planetary interior profile lies below the predicted binodal.The simulation uses xHe = 0.089 and 400 GPa, with temperature rising from 3700 K at the top to 3900 K at the base.Source: Figure 5 · p. 6
  • At xHe = 0.089, Saturn’s interior profiles intersect the predicted immiscibility region, while Jupiter’s lie above but close to the boundary.The authors state that more accurate MLIPs and planetary thermal profiles could help determine whether helium demixes in Jupiter.Source: Hydrogen–helium mixture · p. 6, Figure 5 · p. 6
  • For 1 ns, the model simulates helium droplet nucleation, growth, and redistribution in a column corresponding to 2.2 million atoms.Under enhanced gravity, helium accumulates below and is depleted above; depletion shifts the rain zone downward.Source: Hydrogen–helium mixture · p. 6

Discussion

The discussion frames AIPF as a learned, data-efficient route to mesoscopic thermodynamics and dynamics, while identifying limits in the current model and directions for extending it. Broader training could yield reusable mesoscopic models for conditions where experiments are difficult and direct MD is costly.

Discussion

  • Pressure stabilizes the Fe–B melt but promotes demixing in the H–He mixture.The framework is presented as accurately predicting mesoscopic thermodynamics and dynamics from microscopic data.Source: Discussion · p. 7
  • Trajectories starting from unstable mixtures or slabs sample diverse density profiles, enabling short-trajectory training and access to spinodal states.The discussion contrasts this with equilibrium approaches, which cannot sample the spinodal region and can reach it only by extrapolation.Source: Discussion · p. 7
  • Current helium-rain simulations resolve density fields only, so helium accumulates under enhanced gravity by diffusion rather than hydrodynamic droplet settling.Including momentum density could capture how flow shapes phase separation and transport.Source: Discussion · p. 7
  • The Markovian assumption may be inadequate when memory persists on the timescale of resolved fields, as in some supercooled liquids and glasses.The proposed extension is to learn memory kernels and correlated noise.Source: Discussion · p. 7
  • Training across elements, compositions, and conditions could produce reusable mesoscopic foundation models analogous to CALPHAD databases.The proposed models would use general free-energy functionals rather than parametric bulk free energies.Source: Discussion · p. 7

Methods

The methods derive stochastic phase-field dynamics by projecting molecular dynamics onto species-density fields, then learn free energy and mobility from atomistic data with thermodynamic and fluctuation constraints. They also specify the weighted training objective and system-specific regularization choices.

Theoretical derivation

  • Mori–Zwanzig projection yields a generalized Langevin equation for band-limited species-density modes; short-memory, local-mobility, and white-noise assumptions reduce it to stochastic phase-field dynamics.The resulting dynamics admits the Boltzmann distribution defined by the free-energy functional as its invariant density.Source: Methods · pp. 7–8

Training data

  • Training used Fe–B cubes of 3456 atoms spanning xB = 0.10–0.98 at 0, 5, and 10 GPa, and H–He cubes spanning xHe = 0.05–0.98 at 200, 400, 600, and 800 GPa.H–He spanned 2000–12 000 K; slabs of 2000 to 4700 atoms each paired two compositions between pure hydrogen and pure helium.Source: Methods · pp. 8–9

Model parametrization

  • The learned Helmholtz functional combines an exact ideal-gas term with learned local and nonlocal excess terms, while a density- and temperature-dependent mobility is constrained to be positive semidefinite.The excess free energy uses a learned local density and radial interaction kernel; the mobility construction ensures symmetry and positive semidefiniteness.Source: Methods · p. 9

Thermodynamic predictions

  • At fixed temperature and pressure, the Fe–B and H–He phase diagrams follow from Gibbs free energy of mixing: common tangents give binodals and negative composition curvature defines spinodals.The total density varies with composition at prescribed temperature and pressure, and the binodal and spinodal boundaries meet at the critical point.Source: Methods · p. 10

Dynamical predictions

  • For dynamical predictions, the learned AIPF model is integrated using a pseudospectral method on a periodic grid, with the method applied to the model’s evolution on that grid.Source: Methods · p. 10

Training objective

  • Fluctuation measurements anchor mobility and free-energy curvature because drift fitting alone leaves their scale ambiguous and becomes less informative near critical slowing down.The model matches mobility estimates and structure factors; for NPT Fe–B and H–He data it also matches prescribed pressure through an equation-of-state anchor.Source: Methods · p. 11

Methods

  • The objective combines Ldyn with weighted mobility, structure, bulk, and penalty loss terms.L = Ldyn + λM LM + λS LS + λbulk Lbulk + λP LP.Source: Methods · p. 11
  • For the Lennard-Jones mixture, λS = λP = 0; Fe–B and H–He receive additional regularization terms, with forms and weights listed in Supplemental Material, Section S3.Source: Methods · p. 11

Data availability

The paper makes its trained models, figure data, and data-generation pipeline available online, while large trajectories and datasets are available on reasonable request.

Data availability

  • The trained AIPF models, figure data, and pipeline for generating molecular-dynamics trajectories and training datasets are available on GitHub.The trajectories and datasets themselves are available from the corresponding author on reasonable request because of their size.Source: Data availability · p. 11
  • The trajectories and datasets are available from the corresponding author on reasonable request because of their size.Source: Data availability · p. 11
  • The hydrogen–helium and iron–boron MLIPs were taken from Wang et al. and Zhang et al. [41], respectively.Source: Data availability · p. 11

Appendix A: Additional figures

Additional figures compare model predictions with molecular-dynamics results and show thermodynamic, equation-of-state, and computational-scaling diagnostics for the studied mixtures.

Appendix A: Additional figures

  • Figure A1 compares MD and AIPF, Flory–Huggins, and Landau composition profiles under hexagonal and checkerboard potentials at T = 1.20, conditions excluded from training.Each row also shows the imposed potential and profiles along y = Ly/2.Source: Figure A1 · p. 14
  • At T ≈ 1.17 and below, Landau coexistence compositions fall outside [0, 1], making subsequent ψ growth a clamping artifact rather than physical interface relaxation.Figure A2 also presents MD composition fields and x-averaged profiles at t = 0 and t = 2 × 10^4 τ.Source: Appendix A · p. 15, Figure A2 · p. 15
  • Figure A3 plots the iron–boron mixing free energy against boron fraction along learned isobars at 0, 5, and 10 GPa, with temperatures from 1200 to 2600 K.The corresponding spinodal regions are identified by the thermodynamic factor in Fig. 4a,c.Source: Figure A3 · p. 16
  • Figure A4 compares GPU hours per simulated nanosecond and peak GPU memory against atom count for hydrogen–helium at xHe = 0.50, 400 GPa, and 9000 K.It compares AIPF runs using a 5 fs time step with MD using a 0.2 fs time step; the dotted line indicates linear scaling.Source: Figure A4 · p. 17
  • Figure A5 shows hydrogen–helium mixing free energy versus helium fraction along learned isobars, with common tangents determining binodal compositions.Its thermodynamic-factor panel marks ideal mixing at Γ = 1 and the spinodal boundary at Γ = 0; negative values identify the spinodal region.Source: Figure A5 · p. 18
  • Figure A6 shows the learned Lennard-Jones bulk free energy with binodal and spinodal compositions, and marks negative composition curvature in the spinodal region.Source: Figure A6 · p. 19
  • Figure A7 compares hydrogen–helium number-density predictions from AIPF with molecular-dynamics reference states across four pressures and six isotherms.The AIPF predictions are obtained by solving the Euler relation at each composition and temperature.Source: Figure A7 · p. 20
  • Figure A8 compares learned thermodynamic-factor and Onsager-mobility predictions with MD measurements at homogeneous states spanning 6500–12 000 K and 200–800 GPa.The cross-species mobility MHHe is negative at every state and is shown in magnitude.Source: Figure A8 · p. 21

S1. DERIVATION OF THE MESOSCOPIC MODEL

The Mori–Zwanzig projection gives an exact generalized Langevin equation for resolved observables, specialized here to band-limited species-density fields. A Markovian approximation then makes the model tractable, while the exact relations themselves require neither that approximation nor Gaussian noise.

so that, by eq. (4), P2

Within the cutoff these are the Fourier coefficients of the microscopic number density P

so that ρ = Πρ.

  • Because momentum reversal makes the density-mode time derivative odd while the resolved density and equilibrium distribution are even, its projected drift vanishes.The independent resolved variables are the real and imaginary parts of nonzero modes in a wavevector half-space.Source: so that ρ = Πρ · p. 25

Substituting into eq. (6) yields

  • The memory kernel and fluctuating force account for eliminated momenta and density modes above the cutoff; the force has zero conditional mean and satisfies a conditional fluctuation–dissipation relation.These relations are exact under the stated conditions and do not require Markovian dynamics or Gaussian noise.Source: Substituting into eq. (6) yields · pp. 25–26

C. From the GLE to the mesoscopic model

  • The exact equation is made tractable by assuming the memory kernel decays faster than the resolved fields and replacing it with its time-integrated mobility.The approximation freezes the density over the memory time and neglects the remaining short-lived tail.Source: C. From the GLE to the mesoscopic model · p. 26

The fluctuation–dissipation relation in eq. (17), together with the delta approximation in eq. (18), gives

Under the fluctuation–dissipation relation and a short-memory approximation, the integrated mobility can be related to current correlations. Local mobility and Gaussian white-noise assumptions then yield a conservative real-space mesoscopic model.

The fluctuation–dissipation relation in eq. (17), together with the delta approximation in eq. (18), gives

  • The integrated mobility is represented through time-integrated correlations of orthogonally propagated currents, with wavevector factors enforcing mass conservation.The corresponding real-space mobility operator acts on chemical-potential gradients through a two-point tensor kernel.Source: The fluctuation–dissipation relation… · p. 27
  • The local-mobility assumption restricts mobility to a local, isotropic function of density and temperature, while leaving the free energy unrestricted by locality.The local restriction makes integrated mobility easier to identify from trajectory data.Source: The fluctuation–dissipation relation… · p. 27

Proof. Equation (13) gives ∂ργ/∂ˆρk′,β = V −1δγβeik′·r. Consequently,

where i, j are Cartesian components, and set

Using Lemma 1, the three assumptions give

kBT

The model’s Itō dynamics has a specified drift and diffusion, which determine its Fokker–Planck equation. Its invariant density is the equilibrium marginal.

kBT

Proof. The integral of a periodic divergence is zero, so

  • The Fokker–Planck equation for the real mode coordinates follows from this drift and diffusion, with the mobility’s zero divergence used in the derivation.The proof invokes the vanishing divergence established by Lemma 1.Source: Proof. The integral of a periodic… · p. 30

S2. NUMERICAL METHODS

The numerical methods use pseudospectral density-field evolution, implicit–explicit time integration, and stochastic fluxes. They also specify how external potentials enter equilibrium calculations across canonical and grand-canonical settings.

A. Dynamics

  • The two-species solver evaluates nonlocal free-energy contributions as Fourier-space matrix–vector products and computes flux divergences spectrally.Local chemical potentials are differentiated from the local free energy; mobility multiplies chemical-potential gradients in real space.Source: A. Dynamics · p. 31

Time integration. We use a first-order implicit–explicit (IMEX) scheme based on the reference linear operator

  • The IMEX scheme treats a reference linear operator implicitly and the remainder explicitly, requiring only an independent 2×2 linear solve per Fourier mode.The added and subtracted linear terms cancel, so the reference linearization does not replace the nonlinear dynamics.Source: Time integration. We use a first-order… · p. 31
  • Stochastic fluxes preserve cross-species noise correlations through a Cholesky factor of local mobility, and their zero mode is set to zero to conserve each species’ number.The noise amplitude is evaluated at the start of each step using an Itō-convention Euler–Maruyama discretization.Source: Time integration. We use a first-order… · pp. 31–32
  • For the Lennard-Jones mixture, a scalar scheme evolves the composition field using a Cahn–Hilliard gradient-energy stabilizer as its reference linear operator.The effective gradient-energy coefficient comes from the learned kernel’s small-wavevector expansion or the baseline gradient coefficient.Source: Time integration. We use a first-order… · p. 32

The external potential enters Fθ through the coupling term

S3. EXPERIMENTAL DETAILS AND FURTHER RESULTS

The Lennard-Jones mixture model is trained on slab dynamics and homogeneous-cube anchors, with quench trajectories reserved for validation. Its coarsening analysis estimates domain growth from the composition spectrum.

A. Symmetric binary Lennard-Jones mixture

The local energy u0,θ is an even Taylor polynomial in xA −1

  • The dynamical loss compares a windowed density time derivative with model drift averaged over the same interval.For Lennard-Jones trajectories, the window spans 25 saved frames and drift is evaluated at five frames per window.Source: The local energy u0,θ is an even… · p. 33
  • Bulk anchors constrain only the zero-wavevector free-energy curvature at equimolar composition, using Scc(0) extrapolated from the concentration structure factor.Mobility is estimated from density-mode increments after deconvolving the coarse-graining filter.Source: The local energy u0,θ is an even… · p. 34
  • Coarsening length is derived from the inverse first moment of the composition spectrum, with the same estimator applied to molecular dynamics and AIPF fields.The training loss combines dynamics, mobility, and bulk terms, without a finite-wavevector static or isobar anchor.Source: The local energy u0,θ is an even… · p. 34

where S(k, t) = | ˆψ(k, t)|2. For molecular dynamics, ˆψ = 1

The supplementary section details coarsening analysis and the iron–boron and hydrogen–helium model-training procedures. It also reports a largest Lennard-Jones growth-exponent discrepancy and bounds what can be inferred from finite-box fits.

where S(k, t) = | ˆψ(k, t)|2. For molecular dynamics, ˆψ = 1

  • The fitted growth exponent is an effective, window-dependent value rather than an asymptotic law; distinguishing candidate laws requires a box several times larger.The slope is selected after the initial transient and before finite-size saturation.Source: where S(k, t) = | ˆψ(k, t)|2. For… · p. 35
  • At T = 1.10, the AIPF growth exponent is 0.161 ± 0.020, compared with 0.293 ± 0.024 from molecular dynamics, the largest discrepancy across five quench temperatures.Agreement is closest at T = 1.20.Source: where S(k, t) = | ˆψ(k, t)|2. For… · p. 35

B. Iron–boron melt

  • For the iron–boron melt, the learned free energy combines ideal and excess contributions, while the density-dependent mobility is constrained to be positive semidefinite.The model is trained on 0, 5, and 10 GPa datasets; states are classified using the concentration structure factor and trajectory behavior.Source: B. Iron–boron melt · p. 36, The mobility is parametrized as Mθ =… · p. 36

The mobility is parametrized as Mθ = C1/2

The loss is L = λdyn

  • The iron–boron objective combines dynamical, mobility, static, bulk, pressure, and kernel-shape terms, with no convexity penalty on the local free energy.This leaves the local free energy able to become nonconvex in the boron-rich region where instability is observed.Source: The loss is L = λdyn · p. 37

C. Hydrogen–helium mixture

The mobility is parametrized as Mθ = D1/2LθL⊤

  • Hydrogen–helium anchors use only single-phase cube runs above pressure-specific temperature thresholds, selecting 98 runs in total from 403 cube runs.The remaining cubes contribute only to the dynamical loss; all 484 equation-of-state rows contribute to the isobar anchor.Source: The mobility is parametrized as Mθ =… · p. 39

The second penalty constrains the composition dependence of the thermodynamic factor,

  • A single stochastic simulation produces the helium-rain panels, using a periodic 13.4 × 13.4 × 35.8 nm domain whose displayed half corresponds to approximately 2.2 million atoms.The initial state is metastable throughout, with a temperature gradient from 3700 K to 3900 K and imposed gravity.Source: The second penalty constrains the… · p. 41

Related papers

Loading 2610.01432v1…