cond-mat.stat-mechOct 1, 2026
SaveLearning ab initio phase-field models
Organizations: Department of Mathematics, National University of Singapore, Singapore · Department of Materials Science and Engineering, National University of Singapore, Singapore · Institute for Functional Intelligent Materials, National University of Singapore, Singapore
Abstract
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 composition is spinodally unstable at ambient pressure but stabilized at 10 GPa, offering a thermodynamic rationale for why FeB 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.
Figures & tables
Figure 1: Ab initio phase-field models. a , The proposed ab initio phase-field (AIPF) model is ab initio in both its derivation and its learning, which start from quantum mechanics. The derivation proceeds from quantum mechanics to molecular dynamics under the Born–Oppenheimer approximation, then to a generalized Langevin equation through the Mori–Zwanzig formalism, and finally to the ab initio phase-field model under the stated assumptions. For the learning process, the AIPF model is trained on molecular dynamics trajectories generated with a machine-learning interatomic potential (MLIP) of ab initio accuracy. The free energy and mobility of the density field are parametrized by neural networks, with the free-energy functional chosen from the ladder in b . b , Ladder of free-energy functionals from exact statistical mechanics to phenomenology, where splits into the known ideal-gas and external-potential parts and the unknown excess part . Each lower rung introduces additional restrictions, reducing flexibility and training-data requirements. Rung 4 is a neural operator, with layer fields starting from , convolution kernels , channel-mixing weights , biases and a pointwise nonlinearity . Rung 3, adopted in this work, restricts the excess free energy to a local free-energy density plus a nonlocal kernel . Rung 2 approximates the nonlocal term by the Cahn–Hilliard square-gradient form with coefficient . Rung 1 further restricts the local free energy to a prescribed basis expansion in composition , such as Landau polynomials or Redlich–Kister expansions, with coefficients .
Figure 2: Framework for learning ab initio phase-field models. In the simulation stage, molecular dynamics with an MLIP propagates a configuration over a short interval . In the learning stage, each configuration is coarse-grained into species density fields , from which neural networks with trainable parameters output the local mobility and the free-energy functional . At inference, the thermodynamics follows from minimizing the learned functional and the dynamics from integrating the mesoscopic equation.
Figure 3: Binary Lennard-Jones mixture. a , Phase diagrams from the AIPF model and the Landau (upper small panel) and Flory–Huggins (lower small panel) baselines, compared with MD [ 32 ] . All models are trained on the same data using the same loss function and mobility parametrization, differing only in the form of the free-energy functional. b , Domain size at (3 seeds each), compared between MD and the three models, with defined in Supplemental Material, Section S3 A. Fits to give effective exponents and .
Figure 4: Iron–boron melt. a , The predicted thermodynamic factor at 0, 5 and GPa is shown with the MD classification of the sampled states. For each state, is obtained from a – ps trajectory by extrapolating the concentration structure factor to using the Ornstein–Zernike form. States are classified as stable when , is converged to within and composition fluctuations show no growth, and as spinodal when and the melt coarsens. All remaining states are classified as unresolved. b , The spinodal regions predicted by the AIPF model at the three pressures. The star marks the FeB 4 composition at K, inside the spinodal region at 0 and GPa and outside it at GPa. c , The thermodynamic factor is shown at three temperatures for each pressure. Dots mark the FeB 4 composition. d , The predictions are compared with for states classified by MD as spinodal or unresolved. The AIPF model agrees with the spinodal classification from MD and predicts the sign of for unresolved states.
Figure 5: Hydrogen–helium mixture. The AIPF model is trained at 200, 400, 600 and GPa, with predictions in d,e extending to – GPa. a , Binodal surface in with the critical line . b , Isobaric sections of a at 200, 400, 600 and GPa, with the predicted critical temperatures. Circles denote binodal compositions measured from MD at GPa. c , Schematic of helium demixing from metallic hydrogen and settling toward the center of a gas giant. d , Predicted binodal (solid) and spinodal (dashed) at the protosolar composition , compared with published binodals from Schöttler and Redmer [ 33 ] , Morales et al. [ 34 ] and Lorenzen et al. [ 35 ] , and the spinodal from Wang et al. [ 13 ] , computed with the vdW-DF functional. e , The same boundaries are compared with the interior profiles of Jupiter and Saturn. Bands span published isentropes and present-day thermal profiles compiled by Wang et al. [ 13 ] . Helium rain is expected wherever a profile lies below the binodal. f,g , Helium rain simulated for ns in a nm column corresponding to 2.2 million atoms, at and GPa (Supplemental Material, Section S3). Temperature increases from K at the top to K at the base, with enhanced gravity directed toward the planet’s center. f , Opaque surfaces show helium droplets colored by radius. White surfaces indicate helium-rich fluctuations below the droplet threshold. g , The helium fraction is shown on a wedge section of the same domain at the same times.
Appendix figures & tables14 assets
Supplementary material from the paper’s appendix.
Appendix
Figure A1: Equilibrium distribution of the Lennard-Jones mixture under external potentials. a,b , Equilibrium profiles are compared for hexagonal ( a ) and checkerboard ( b ) potentials at . Each row shows, from left to right, the imposed potential , with its magnitude indicated by the vertical axis, the coarse-grained composition field from molecular dynamics, the fields of the AIPF, Flory–Huggins (FH) and Landau models, and profiles from all four along . The fields of all three models are time averages of stochastic simulations started from a uniform composition. The external potentials and this temperature are not included in training.
Figure A2: Interface relaxation in the Lennard-Jones mixture. a,b , MD composition fields ( a ) and corresponding -averaged profiles ( b ) at and in a box, with horizontal. The temperature, , is not included in training. The composition contrast is defined as , where and are the upper and lower fitted plateaus of (right and left in b ). c–e , Comparison of with MD at seven temperatures for the Landau ( c ), Flory–Huggins ( d ) and AIPF ( e ) models. Curves show means over four independent seeds for MD and each model, and bands show one standard deviation. In the stochastic model simulations, is clamped to at every step and the mean composition is restored. Below , the Landau model predicts coexistence compositions outside . In these gray panels of c , the subsequent increase in is a clamping artifact rather than physical interface relaxation.
Figure A3: Mixing free energy of the iron–boron melt. The Gibbs free energy of mixing per atom, , is shown as a function of boron fraction along the learned isobars at , and GPa. Colors indicate temperatures from to K in increments of K. The corresponding spinodal regions are identified by the thermodynamic factor in Fig. 4 a,c.
Figure A4: Comparison of computational scaling. a , GPU hours per simulated nanosecond and b , peak GPU memory as functions of the number of atoms for cubic boxes of the hydrogen–helium mixture at , GPa and K. All runs use a single NVIDIA A100 GPU with GB of memory. Molecular dynamics uses a MACE potential trained specifically for hydrogen–helium on vdW-DF data [ 13 ] , with approximately parameters. Filled symbols denote the AIPF model simulated with a fs time step and a grid spacing of approximately Å. Open symbols, labeled ML potential, denote MD with a fs time step. The dotted line indicates linear scaling with system size.
Figure A5: Mixing free energy and thermodynamic factor of hydrogen–helium. a , The Gibbs free energy of mixing per atom, , is shown as a function of helium fraction along the learned isobars at , , and GPa. Colors indicate temperatures from to K in increments of K. Common tangents determine the binodal compositions. b , The corresponding thermodynamic factor ( Eq. 4 ). The dotted line marks ideal mixing ( ), and the solid line marks the spinodal boundary ( ). Negative values identify the spinodal region.
Figure A6: Learned bulk free energy of the Lennard-Jones mixture. a , The learned bulk free energy is shown with the binodal and spinodal compositions, where is the number fraction of species . b , The composition curvature is negative in the spinodal region.
Figure A7: Comparison of the hydrogen–helium equation of state from the AIPF model and molecular dynamics. Number density is shown as a function of helium fraction at , , and GPa. Colors indicate six isotherms from to K in increments of K. Points denote molecular-dynamics reference states obtained with the machine-learning interatomic potential. Lines show the predictions of the AIPF model, obtained by solving the Euler relation ( Eq. 12 ) at each composition and temperature.
Figure A8: Agreement with the hydrogen–helium training anchors. The predictions of the AIPF model are compared with molecular-dynamics measurements at the homogeneous states used for the static and mobility anchors. These states span – K at , , and GPa. Marker shapes indicate pressure, colors indicate temperature, and dashed lines denote equality. The first panel compares calculated from the learned free-energy curvature with estimates from MD concentration fluctuations. The remaining panels compare the independent components of the symmetric Onsager mobility matrix, measured from short-time density changes. The cross-species mobility is negative at every state and is shown in magnitude.
| Slab | Homogeneous cube | Quench cube | |
| Particles | |||
| Cells | |||
| Temperatures | – | – | |
| Production ( ) |
Table S1: Microscopic simulations used for the Lennard-Jones mixture. All three datasets use the force-shifted pair potential described in Methods, with , total density and mean composition . Trajectories are generated by overdamped Langevin dynamics with and a time step of . Configurations are saved every . Slab trajectories provide the dynamical training data, with held out from both training and model selection. Homogeneous cubes provide the bulk and mobility anchors. Quench cubes are used exclusively to validate spinodal decomposition. Two additional systems of particles are simulated with an underdamped integrator under checkerboard and hexagonal external potentials. These runs sample equilibrium density profiles and are excluded from training.
Figure S1: Coarsening of the Lennard-Jones mixture. Results at , , and are shown from top to bottom, complementing the results at in the main text. a,c,e,g , Composition fields from the ab initio phase-field (AIPF) model at four times. b,d,f,h , Domain size , defined in eq. 32 , compared between three independent molecular dynamics (MD) runs and three independent AIPF runs.
| GPa | GPa | GPa | |
| State points | |||
| Runs | |||
| Production (ps) | to , per state point | ||
| Spinodal / stable / unresolved | / / | / / | / / |
| Excluded from training | |||
| Equation-of-state rows | |||
Table S2: Molecular-dynamics datasets for the iron–boron melt. Trajectories are generated for homogeneous cubes of atoms in the isotropic isothermal–isobaric (NPT) ensemble using the fine-tuned MACE potential of Zhang et al. [ 12 ] . Each system is melted at K for ps before simulation at the target temperature. The time step is fs, and configurations are saved every ps. The state-point grid and selection criteria are described in the text. Retained trajectories provide the dynamical training data, and stable states that satisfy the anchor-selection threshold also provide the static and mobility anchors. Separate equation-of-state calculations use atoms, with ps each for melting, equilibration and averaging, and provide the isobar anchor.
Figure S2: Equation of state of the iron–boron melt at , and GPa. Number density is shown as a function of boron fraction along six isotherms between and K, distinguished by color. Points denote molecular-dynamics reference states obtained with the fine-tuned potential. Lines show the learned isobars, calculated by solving the model pressure equation at fixed composition and temperature. The reference states are included in the isobar anchor during training, so this comparison assesses the fit to that constraint.
| Slabs | Homogeneous cubes | Equation of state | |
| Atoms | – | ||
| Ensemble | NPT, only | NPT, isotropic | NPT, isotropic |
| Time step (fs) | |||
| Production (ps) | |||
| Runs at , , , GPa | , , , | , , , | rows each |
Table S3: Molecular-dynamics datasets for hydrogen–helium. All runs use the MACE potential of Wang, Hamel and Cheng [ 13 ] , a Nosé–Hoover thermostat and a Martyna–Tobias–Klein barostat. Configurations are saved every ps. Slabs are initialized with a double hyperbolic-tangent composition profile containing two interfaces in a Å box, with pressure control applied only along the interface normal. Cubes are melted for ps at K above the target temperature before production. At GPa, cube simulations span – K in increments of K and helium fractions of – . At each lower pressure, they span helium fractions of – at temperatures on both sides of the critical line. Slabs and cubes contribute to the dynamical loss. The cube runs that meet the temperature and single-phase criteria also provide the static, bulk and mobility anchors. The equation-of-state rows provide the isobar anchor and the density paths for the thermodynamic-factor prior.
Figure S3: Coverage of the hydrogen–helium cube datasets at the four training pressures. Each run is shown at its helium fraction and temperature. Rust denotes runs used only for the dynamical loss. Navy denotes runs that also provide the static, bulk and mobility anchors. Anchor selection requires temperatures above , , and K at , , and GPa, respectively, and single-phase behavior over the sampled interval. Runs below these thresholds are excluded from the anchors because the uncertainty in the concentration structure factor is underestimated. The datasets contain , , and cube runs at the four pressures, of which , , and provide anchor data.