Authors: Mengyi Chen, Peichen Zhong, Zihan Zhang, Qianxiao Li
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
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 FeB4 composition is spinodally unstable at ambient pressure but stabilized at 10 GPa, offering a thermodynamic rationale for why FeB4 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 F splits into the known ideal-gas and external-potential parts and the unknown excess part Fexc . Each lower rung introduces additional restrictions, reducing flexibility and training-data requirements. Rung 4 is a neural operator, with layer fields vl starting from v0=ρ , convolution kernels kl(r−r′,T) , channel-mixing weights Wl , biases bl(T) and a pointwise nonlinearity σ . Rung 3, adopted in this work, restricts the excess free energy to a local free-energy density floc(ρ,T) plus a nonlocal kernel W(∣r−r′∣,T) . 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 x , such as Landau polynomials or Redlich–Kister expansions, with coefficients cn(T) .
Figure 2: Framework for learning ab initio phase-field models. In the simulation stage, molecular dynamics with an MLIP propagates a configuration xt over a short interval δt . In the learning stage, each configuration is coarse-grained into species density fields ρ , from which neural networks with trainable parameters θ output the local mobility Mθ and the free-energy functional Fθ . 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 L(t) at T=1.20 (3 seeds each), compared between MD and the three models, with L(t) defined in Supplemental Material, Section S3 A. Fits to L∝tβ give effective exponents βMD=0.225±0.022 and βAIPF=0.217±0.008 .
Figure 4: Iron–boron melt. a , The predicted thermodynamic factor Γ(xB,T) at 0, 5 and 10 GPa is shown with the MD classification of the sampled states. For each state, ΓMD is obtained from a 50 – 500 ps trajectory by extrapolating the concentration structure factor Scc(k) to k=0 using the Ornstein–Zernike form. States are classified as stable when ΓMD≥1 , Scc(0) is converged to within 10% and composition fluctuations show no growth, and as spinodal when ΓMD≤0.35 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 1500 K, inside the spinodal region at 0 and 5 GPa and outside it at 10 GPa. c , The thermodynamic factor Γ(xB) is shown at three temperatures for each pressure. Dots mark the FeB 4 composition. d , The predictions ΓAIPF are compared with ΓMD 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 800 GPa, with predictions in d,e extending to 150 – 900 GPa. a , Binodal surface in (xHe,P,T) with the critical line Tc(P) . b , Isobaric sections of a at 200, 400, 600 and 800 GPa, with the predicted critical temperatures. Circles denote binodal compositions measured from MD at 800 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 xHe=0.089 , 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 1 ns in a 13.4×13.4×17.9 nm column corresponding to 2.2 million atoms, at xHe=0.089 and 400 GPa (Supplemental Material, Section S3). Temperature increases from 3700 K at the top to 3900 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 T=1.20 . Each row shows, from left to right, the imposed potential Vext , 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 y=Ly/2 . 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 x -averaged profiles ( b ) at t=0 and t=2×104τ in a 9.52×9.52×38.10σAA3 box, with z horizontal. The temperature, T=1.20 , is not included in training. The composition contrast is defined as ψ(t)=[xR(t)−xL(t)]/2 , where xR and xL are the upper and lower fitted plateaus of xA (right and left in b ). c–e , Comparison of ψ(t) 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, xA is clamped to [10−4,1−10−4] at every step and the mean composition is restored. Below T≈1.17 , the Landau model predicts coexistence compositions outside [0,1] . 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, ΔGmix , is shown as a function of boron fraction xB along the learned isobars at 0 , 5 and 10 GPa. Colors indicate temperatures from 1200 to 2600 K in increments of 200 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 xHe=0.50 , 400 GPa and 9000 K. All runs use a single NVIDIA A100 GPU with 40 GB of memory. Molecular dynamics uses a MACE potential trained specifically for hydrogen–helium on vdW-DF data [ 13 ] , with approximately 2×105 parameters. Filled symbols denote the AIPF model simulated with a 5 fs time step and a grid spacing of approximately 0.7 Å. Open symbols, labeled ML potential, denote MD with a 0.2 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, ΔGmix , is shown as a function of helium fraction along the learned isobars at 200 , 400 , 600 and 800 GPa. Colors indicate temperatures from 2000 to 12000 K in increments of 1000 K. Common tangents determine the binodal compositions. b , The corresponding thermodynamic factor Γ ( Eq. 4 ). The dotted line marks ideal mixing ( Γ=1 ), and the solid line marks the spinodal boundary ( Γ=0 ). Negative values identify the spinodal region.
Figure A6: Learned bulk free energy of the Lennard-Jones mixture. a , The learned bulk free energy f(xA,T) is shown with the binodal and spinodal compositions, where xA is the number fraction of species A . b , The composition curvature ∂2f/∂xA2 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 n is shown as a function of helium fraction xHe at 200 , 400 , 600 and 800 GPa. Colors indicate six isotherms from 2000 to 12000 K in increments of 2000 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 6500 – 12000 K at 200 , 400 , 600 and 800 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 2×2 Onsager mobility matrix, measured from short-time density changes. The cross-species mobility MHHe is negative at every state and is shown in magnitude.
Slab
Homogeneous cube
Quench cube
Particles N
3456
4000
6912
Cells
6×6×24
103
123
Temperatures
1.00 – 1.70
1.50,1.60,1.70,1.80
1.10 – 1.30
Production ( τ )
20000
5000
10000
Appendix
Table S1: Microscopic simulations used for the Lennard-Jones mixture. All three datasets use the force-shifted pair potential described in Methods, with ϵAB=0.5ϵAA , total density ρ=1.0σAA−3 and mean composition xA=0.5 . Trajectories are generated by overdamped Langevin dynamics with γt=2.0 and a time step of 2×10−4τ . Configurations are saved every τ . Slab trajectories provide the dynamical training data, with T=1.20 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 55296 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 T=1.10 , 1.15 , 1.25 and 1.30 are shown from top to bottom, complementing the results at T=1.20 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 L(t) , defined in eq. 32 , compared between three independent molecular dynamics (MD) runs and three independent AIPF runs.
0 GPa
5 GPa
10 GPa
State points
145
145
145
Runs
156
159
149
Production (ps)
50 to 500 , per state point
Spinodal / stable / unresolved
12 / 88 / 45
8 / 81 / 56
4 / 72 / 69
Excluded from training
31
36
41
Equation-of-state rows
208
106
128
Appendix
Table S2: Molecular-dynamics datasets for the iron–boron melt. Trajectories are generated for homogeneous cubes of 3456 atoms in the isotropic isothermal–isobaric (NPT) ensemble using the fine-tuned MACE potential of Zhang et al. [ 12 ] . Each system is melted at 2600 K for 10 ps before simulation at the target temperature. The time step is 1 fs, and configurations are saved every 0.1 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 512 atoms, with 5 ps each for melting, equilibration and averaging, and provide the isobar anchor.
Figure S2: Equation of state of the iron–boron melt at 0 , 5 and 10 GPa. Number density n is shown as a function of boron fraction xB along six isotherms between 1200 and 2600 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 442 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
2000 – 4700
3456
512
Ensemble
NPT, z only
NPT, isotropic
NPT, isotropic
Time step (fs)
0.1
0.2
0.2
Production (ps)
40
30
1
Runs at 200 , 400 , 600 , 800 GPa
48 , 56 , 64 , 88
79 , 90 , 101 , 133
121 rows each
Appendix
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 0.02 ps. Slabs are initialized with a double hyperbolic-tangent composition profile containing two interfaces in a 12×12×40 Å box, with pressure control applied only along the interface normal. Cubes are melted for 1.5 ps at 2000 K above the target temperature before production. At 800 GPa, cube simulations span 2000 – 12000 K in increments of 1000 K and helium fractions of 0.05 – 0.98 . At each lower pressure, they span helium fractions of 0.05 – 0.95 at temperatures on both sides of the critical line. Slabs and cubes contribute to the dynamical loss. The 98 cube runs that meet the temperature and single-phase criteria also provide the static, bulk and mobility anchors. The 484 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 6100 , 7800 , 8800 and 8973 K at 200 , 400 , 600 and 800 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 79 , 90 , 101 and 133 cube runs at the four pressures, of which 22 , 22 , 22 and 32 provide anchor data.
Phase-field models play a central role in the continuum description of phase separation, in which the bulk free-energy density and the interfacial thickness parameter determine pattern formation and microstructural evolution. In practice, these constitutive quantities are rarely known a priori and must be inferred from limited dynamical observations. In this work, an extended pseudo-spectral physics-informed neural network (ESPINN) framework is developed for the inverse identification of phase-field models from transient snapshot data. It enables the simultaneous recovery of both the bulk chemical potential and unknown gradient coefficients. Numerical experiments on the one-dimensional Cahn-Hilliard equation demonstrate accurate and statistically stable reconstruction in the noiseless regime, with substantial constitutive information recoverable from even a single snapshot pair. In the presence of noise, reconstruction accuracy degrades gracefully, and increasing the number of snapshots improves robustness by reducing variance across runs. These results establish ESPINN as a data-efficient and physically consistent approach for learning free-energy structure in continuum models of phase separation.
Callum Marsh, Radek Erban, Andreas Munch
Mathematical Institute, University of Oxford, Radcliffe Observatory Quarter, Woodstock Road, Oxford OX2 6GG, United Kingdom
Phase-field modeling provides a powerful approach for predicting microstructure evolution but becomes computationally prohibitive for multicomponent and multiphase systems over large spatial and temporal scales. This work presents an AE-GCN-LSTM surrogate framework for long-horizon forecasting of microstructure evolution in the multicomponent AlCrFeNi high-entropy alloy system containing coexisting BCC and FCC phases. A multi-head autoencoder compresses the four elemental concentration fields and phase-field order parameter into latent representations, which are formulated as graphs for learning their spatial and temporal evolution. The framework accurately forecasts microstructure evolution over horizons extending to 3,000,000 simulation timesteps. Its robustness is systematically evaluated under previously unseen conditions without retraining, fine-tuning, or parameter adaptation. These evaluations include variations in FCC precipitate size and initial position, microstructures containing one, two, and five FCC precipitates, and complex phase interactions involving precipitate merging and splitting. Although trained only on 100 x 100 computational domains containing a single nominal alloy composition, the framework is successfully transferred to larger 256 x 256 and 512 x 512 systems and to previously unseen AlCrFeNi compositions. Across the evaluated configurations, the model preserves the dominant phase morphology and compositional evolution while providing computational speedups ranging from approximately 7200 to 62300 relative to conventional phase-field simulations. These results demonstrate that latent graph-based AE-GCN-LSTM forecasting provides a scalable and computationally efficient surrogate for long-horizon simulation of multicomponent, multiphase microstructures and offers a promising foundation for high-throughput alloy design.
Hamidreza Razavi, Nele Moelans
Department of Materials Engineering, KU Leuven2026 Kasteelpark Arenberg 44 Bus 2450, 3001 Leuven, Belgium
The multi-scale and non-linear nature of phase-field models of solidification requires fine spatial and temporal discretization, leading to long computation times. This could be overcome with artificial-intelligence approaches. Surrogate models based on neural operators could have a lower computational cost than conventional numerical discretization methods. We propose a new neural operator approach that bridges classical convex-concave splitting schemes with physics-informed learning to accelerate the simulation of phase-field models. It consists of a Deep Ritz method, where a neural operator is trained to approximate a variational formulation of the phase-field model. By training the neural operator with an energy-splitting variational formulation, we enforce the energy dissipation property of the underlying models. We further introduce a custom Reaction-Diffusion Neural Operator (RDNO) architecture, adapted to the operators of the model equations. We successfully apply the deep learning approach to the isotropic Allen-Cahn equation and to anisotropic dendritic growth simulation. We demonstrate that our physically-informed training provides better generalization in out-of-distribution evaluations than data-driven training, while achieving faster inference than traditional Fourier spectral methods.