Phase field crystal is a pretty cool formalism for simulating phase and microstructure evolution using Cahn-Hilliard or Allen-Cahn equations for conserved or non-conserved fields. One of the biggest issues, as is the case with CALPHAD since they rely on the same types of quantities, is the lack of thermodynamic datasets. So I figured, what if I got coding agents to train a small set ML models using open data, deliberately kept blind to one specific chemistry, could I use that to hand off the required material parameters to a phase-field-crystal (PFC) model that would capture something like a displacive phase transformation. Key thing is to hold out the specific chemical system and also not deterministically input the transformation temperature ... so it has to fall out of the model physics.
I went with the material system hold-out being a Ni-rich Ni-Ti-Hf shape-memory alloy, specifically Ni50.3Ti29.7Hf20 [1], which undergoes a diffusionless B2 (austenite) to B19' (martensite) displacive phase transition. We (the agent and I) landed on it because that dissertation reports the parameters to see how good the PFC simulation is, four transformation temperatures, the parent lattice constant, and the martensite cell.
Setup
There were some strict requirements that everything downstream of composition has to be data-driven but also hold-out the key ternary system we are simulating. The requirements were:
- Pull every public dataset that reports transformation temperatures, lattice geometry, or elastic moduli for NiTi and related shape memory alloys [2–7].
- Delete every record whose exact chemical system is Ni-Ti-Hf.
- Train nine small regressors/classifiers on the curated dataset, validated with grouped cross-validation by chemical system, not random, so a whole family (e.g., keep everything with Zr in it) is hidden at once.
- ML models predict PFC inputs for Ni50.3Ti29.7Hf20 using nothing but its composition and the fact that it starts as B2. Then feed into a 2D amplitude-PFC simulation and run a full cooling/heating cycle.
- Compare simulation results to real reported numbers and compute the error.
The point of step 2 is that the target alloy's own transformation behavior should never leak into the surrogate models, directly or even through a near-identical composition sitting in the same training fold.
Amplitude-PFC formalism
The first quesiton someone familiar with phase-field modeling is why are we using phase-field crystal? And what is this amplitude flavor? The important thing to noe about shape-memory alloys that a displacive transformation is due to shear, not a diffusive shuffle, so that every atom shifts a small and coordinated distance but exist now on a different crystal lattice. With classical phase-field model this doesn't really have anywhere for that to happen since there really isn't a crystal lattice to speak of because we represent a "featureless" order parameter. This is of less of an issue when trying to simulate diffusion-controlled transformations.
So why Amplitude-PFC? Full PFC formalism keeps the lattice, but the density field still oscillates on atomic length and time scales, so this has some of the same stiff-timestep problem as molecular dynamics. The amplitude formulation aims to sidesteps that by not trying to resolve the oscillating density directly but rather tracking only the slowly varying amplitude and phase (periodic phase not material phase) of each reciprocal-lattice mode. This means the crystal lattice periodicity is encoded into the phase factors rather than in the field itself. This makes running a full 2D heating/cooling cycle finish in about a minute instead of hours.
To do this we are using the 2D structural amplitude-PFC reduction on a principal 1-3 plane, reconstructing the B2 parent's atomic density directly from reciprocal modes rather than placing Gaussian atoms, and following the amplitude-expansion route for complex crystal structures in [8]. For the numerical side, discretizing a free energy like eq.\ref{eq:free_energy} into something that actually runs, semi-implicit time stepping, the grid/timestep choices, is the part I(we) leaned on Biner's book for [9].
Before writing that free energy down, here's exactly what material-specific numbers it needs (see eq.\ref{eq:theta_vector}). Nine independent, data-driven surrogate models (described in the next section) hand off a parameter vector
$$ \begin{equation} \boldsymbol{\theta}=\left(T_0,\ H,\ a_0,\ \lambda_1,\lambda_2,\lambda_3,\ \mathcal{C}_M,\ K,\ G\right), \label{eq:theta_vector} \end{equation} $$
with, the thermodynamic center $T_0$, the hysteresis width $H=A_f-M_s$, the parent lattice constant $a_0$, the three principal transformation stretches $\lambda_i$, the low-temperature product class $\mathcal{C}_M$ (B19 vs. B19'), and the bulk/shear moduli $K,G$. $M_s$ itself is never in this parameter vector $\boldsymbol\theta$ on purpose, it only exists later as an output of the PFC run!
Eight of these parameters are the coefficients for the free energy density expressions below, where $T_0$ and $H$ set the driving force (eq.\ref{eq:driving_force}), $a_0$ and $\lambda_i$ set the lattice reconstruction (eq.\ref{eq:phases}), $K$ and $G$ set the elastic stiffness $C_i$. $\mathcal{C}_M$ isn't a material input into the free energy density but just a label carried along purely so the simulated low-temperature structure has something independent for us to check against afterward.
A structural order parameter $m(\mathbf r,t)\in[0,1]$ (0 = parent, 1 = product) and three principal log-strains $\epsilon_i$ evolve under the equation
$$ \begin{equation} \begin{aligned} F=\int_\Omega\Big[\ &Bm^2(1-m)^2+\Delta g(T,\mathbf r)h(m)+\frac{\kappa_m}{2}|\nabla m|^2 \\ &+\sum_{i=1}^3\left(\frac{C_i}{2}\left[\epsilon_i-m\epsilon_i^*\right]^2+\frac{\kappa_i}{2}|\nabla\epsilon_i|^2\right)\Big]\,d\mathbf r, \end{aligned} \label{eq:free_energy} \end{equation} $$
with
$$ \begin{equation} h(m)=m^2(3-2m),\qquad \epsilon_i^*=\ln\lambda_i, \label{eq:h_and_eps} \end{equation} $$
with the non-conserved (Allen-Cahn) dynamics given by
$$ \begin{equation} \begin{aligned} \frac{\partial m}{\partial t} &= -L_m\frac{\delta F}{\delta m}+\xi, \\ \frac{\partial\epsilon_i}{\partial t} &= -L_i\frac{\delta F}{\delta\epsilon_i}. \end{aligned} \label{eq:allen_cahn} \end{equation} $$
The temperature driving force, i.e., the driving force for the transformation, is given by
$$ \begin{equation} \begin{aligned} \Delta g(T,\mathbf r) &= s\left[T-T_0-\delta T_0(\mathbf r)\right], \\ s &= \frac{2B}{3H}. \end{aligned} \label{eq:driving_force} \end{equation} $$
This is just a double-well for free energy density where parent and product each lose metastability at $\Delta g=\pm B/3$, so the surrogate-predicted hysteresis width will set how far apart the decomposition sits. $\delta T_0(\mathbf r)$ corresponds to the small fixed-amplitude quenched-disorder field. This means it seeds the spatial front but doesn't prescribe anything for the transformation temperature.
Limitation
One honest limitation of this model is that it predicts where and how a displacive transition happens, not whether one happens at all, the double well assumes bistability exists and only ensures that it will happen at some temperature. Therefore we cannot discover whether a displacive transformation is possible for some new material.
The parent (i.e., B2) atomic density on the simulated plane comes from local phases defined as
$$ \begin{equation} \begin{aligned} \phi_1(\mathbf r) &= \int q_1(\mathbf r)\,dx, \\ \phi_3(\mathbf r) &= \int q_3(\mathbf r)\,dz, \\ q_i &= \frac{2\pi}{a_0e^{\epsilon_i}}. \end{aligned} \label{eq:phases} \end{equation} $$
There are also two sublattice-occupancy fields built that are represented as finite Fourier-mode expansions,
$$ \begin{equation} \begin{aligned} n_{\mathrm{X}} &= \frac{(1+\cos\phi_1)(1+\cos\phi_3)}{4}, \\ n_{\mathrm{Y}} &= \frac{(1-\cos\phi_1)(1-\cos\phi_3)}{4}. \end{aligned} \label{eq:sublattice_density} \end{equation} $$
$\mathrm{X}$ and $\mathrm{Y}$ are the two interpenetrating simple-cubic sites of the parent structure (B2). For the Ni-rich Ni-Ti-Hf hold-out, $\mathrm{X}$ is occupied by Ni and $\mathrm{Y}$ by mixed Ti/Hf, which shows up when we weight the reconstructed density (i.e. using atomic numbers). The $\epsilon_i$ is the single strain field coupled to both sublattices, i.e., one shared displacive mode with no shear/rotation channel buth rather three orthogonal stretches. This does introduce a modeling limitation in that the in-plane lattice vectors mostly remain at $90^{\circ}$ rather than the reported monoclinic angle for the product phase. Most likely if we added the missing shear order parameter this would improve things. For the rendered atomic-density field, we use atomic number weighting:
$$ \begin{equation} n_Z=28\,n_{\mathrm{X}}+\left(22\,c_{\mathrm{Ti}|Y}+72\,c_{\mathrm{Hf}|Y}\right)n_{\mathrm{Y}}, \label{eq:electron_density} \end{equation} $$
then mean-centered and variance-normalized for color. There is no Gaussian projections.
Surrogate formalism
In the diagram below, the whole shape of this surrogate stage is shown. Hopefully this aids in what the input and outputs are for each surrogate model into the APFC model.
The surrogate for the APFC is actually a set of nine models, mostly independent, each trained on online data and cross-validated by chemical system. None of the models are trained on crystal structure, a phase diagram, or an alloy. The only inputs are elemental fractions $c_e$, with exception being $\lambda_1$'s model. This model also consumes $\lambda_2$'s and $\lambda_3$'s own predictions (more on that below). Everywhere else the nine models are trained and evaluated entirely on their own. The parameter vector $\boldsymbol\theta$ they output is handed off to the amplitude-PFC formalism above (eq.\ref{eq:theta_vector}).
$c_e$ by itself is just a list of fractions. For some of the models we expand that into extra signals by looking up generic elemental properties $p_e$ (atomic number, atomic weight, group, period, electronegativity, metallic radius, melting point, valence-electron count, density) and then taking composition-weighted mean, spread, and range,
$$ \begin{equation} \begin{aligned} \mu_p &= \sum_e c_e p_e, \\ \sigma_p &= \sqrt{\sum_e c_e\left(p_e-\mu_p\right)^2}, \\ R_p &= \max_{c_e>0}p_e-\min_{c_e>0}p_e, \end{aligned} \label{eq:compact_stats} \end{equation} $$
there is also the configurational entropy,
$$ \begin{equation} S_{\mathrm{mix}}/k_B=-\sum_e c_e\ln c_e. \label{eq:entropy} \end{equation} $$
These signals, plus the raw fractions, entropy, and element count, provides a fairly compact, hand-crafted, descriptor from $c_e$ and $p_e$. $K$ and $G$ use it. For the $T_0$ surrogat we run the same fractions through the CBFV/JARVIS featurizer [10–11], which provides much richer descriptor. The surrogates for $a_0$, $\lambda_2$, and $\lambda_3$ are expanded further using seven composition statistics recovered from the structural source table (eq.\ref{eq:elemental_coef_solve}). Finally, $H$ and $\mathcal{C}_M$ stay as raw fractions and $\lambda_1$ is fractions plus the other two stretches. The table below gives the breakdowns.
| predicts | features from $c_e$ | model | rows / systems |
|---|---|---|---|
| $T_0$ | CBFV/JARVIS (~2,401), not the compact $p_e$ recipe | Extra Trees | 429 / 41 |
| $H$ | 39 raw element fractions only | Extra Trees | 156 / 18 |
| $a_0$ | 19 fractions + 7 recovered table stats, eq.\ref{eq:elemental_coef_solve} | Ridge | 100 / 17 |
| $\lambda_1$ | 19 fractions + stacked $\hat\lambda_2,\hat\lambda_3$, eq.\ref{eq:lambda1_stack} | Random Forest | 100 / 17 |
| $\lambda_2$ | 19 fractions + 7 recovered table stats | Random Forest | 157 / 19 |
| $\lambda_3$ | same as $\lambda_2$ | Gradient-boosted trees | 100 / 17 |
| $\mathcal{C}_M$ | 19 raw element fractions only | Logistic regression | 157 / 19 |
| $K$ | compact recipe, 23 fractions + $\mu_p,\sigma_p,R_p$ + entropy (52) | Extra Trees | 766 / 393 |
| $G$ | same as $K$ | Extra Trees | 766 / 393 |
The surrogates themselves were kept pretty vanilla on pupose, tree ensembles, ridge, or logistic, whichever generalized best per target in grouped cross-validation. This is low-data, so a fancier model probably just buys more ways to overfit a few dozen chemical systems. Adding the compact $\mu_p,\sigma_p,R_p$ block on top of raw fractions measurably hurt $H$ in grouped-CV (checked, not assumed). For $T_0$, the regression target was choosen as: $$ \begin{equation} T_0=\frac{M_s+M_f+A_s+A_f}{4}+273.15\ \mathrm{K}, \label{eq:t0_target} \end{equation} $$
Because I wasn't sure that a surrogate would even do better than just some statistical value from the data, a model is only used if group cross-validation clears MAE improvement of at least 20% over the training-fold median and $R^2\geq0.20$ for regressors, balanced accuracy $\geq0.85$ for the product-class classifier, and the target's own chemical-system centroid has to sit inside the 95th-percentile applicability domain of the non-target systems. Fortunately, all nine surrogates pass but not by the same margin ... somewhat interesting.
Three of the nine, $a_0$, $\lambda_2$, $\lambda_3$, got a real boost from re-using information that was already sitting in the source data, but wasn't initially used. The structural source table ships precomputed composition-statistics columns (average ionization energy, a Voronoi coordination number, molar-volume statistics) that are themselves fraction-weighted averages of unpublished per-element constants. Since they're a linear function of composition, the underlying per-element values are recoverable by ordinary least squares against the training table,
$$ \begin{equation} \hat p_e=\arg\min_{p_e}\sum_i\left(P_i-\sum_e c_{i,e}\,p_e\right)^2, \label{eq:elemental_coef_solve} \end{equation} $$
which then lets the same statistic be evaluated on a brand new composition instead of only on the training rows, this is the "7 recovered composition stats" in the table above. I checked this wasn't overfitting to noise by comparing the recovered Ni/Ti/Hf ionization energies and molar volumes against published values and the match was good. The other boost, and the one place in this whole stage where a model's output becomes another model's input, was for the weakest of the three principal stretches. The three $\lambda_i$ aren't independent (same displacive shuffle, roughly volume conserving), so letting $\lambda_1$'s surrogate see the other two stretches' own held-out predictions as extra inputs,
$$ \begin{equation} \hat\lambda_1=f_\theta!\left(c_e,\ \hat\lambda_2^{(-k)},\ \hat\lambda_3^{(-k)}\right), \label{eq:lambda1_stack} \end{equation} $$
where the $(-k)$ superscript means "predicted by a sub-model that never saw this row's fold during training." That out-of-fold constraint is what keeps it leak-free so a $\lambda_1$ test row's stacked features never come from a sub-model that had already seen that row's label. In truth could have just had a surrogate that predicts all three stretches at once!
Results
Surrogate Qualities
The full grouped cross-validation (grouped-CV) parity including the regression targets, R², MAE, and RMSE printed on each panel:
|
| Figure 1. Grouped cross-validation parity for the 8 regression surrogates, R-squared/MAE/RMSE annotated per panel |
The split in Figure 1 is stark. $a_0$, $\lambda_2$, $\lambda_3$, $K$, and $G$ all land at $R^2=0.77$-$0.90$. $\lambda_1$ improves to $R^2=0.65$ once it's allowed to see the other stretches (eq.\ref{eq:lambda1_stack}). $T_0$ and $H$ sit at $R^2=0.40$ and $0.23$ and I could not move either one, more on that below.
Phase Field Crystal Predictions
Once every surrogate passed, $\boldsymbol{\theta}$ (eq.\ref{eq:theta_vector}) is frozen, and then running the cooling/heating cycle (eq.\ref{eq:allen_cahn}) gives the hysteresis loop in Figure 2:
|
| Figure 2. Transformed fraction versus temperature, cooling and heating branches, Ms/T50/Mf and As/T50/Af crossings marked |
The "blind" run puts martensitic start temperature, $M_s$ at 457.4 K against a reported 429.2 K [1] (+28.2 K), martensitic finish temperature $M_f$ at 454.9 K against 401.2 K (+53.8 K). For the austenitic start temperature $A_s$ at 493.0 K against 433.2 K (+59.9 K), and austenitic finish temperature $A_f$ at 494.4 K against 462.2 K (+32.2 K). Every one of the temperatures runs "hot" and in the same rough range as the surrogates' own held-out error. I think this is what I'd want to see, i.e., the hold-out error is not some wild outlier relative to the validation numbers but rather consistent with them. For $a_0$ it comes out at 3.1005 Å against a reported 3.092 Å [1].
Slicing the actual density field through the transformation, cooling first (Figure 3) we see:
|
| Figure 3. Cooling branch density, atom-site, and transformed-fraction slices from hot parent to cold product |
and the for heating (Figure 4):
|
| Figure 4. Heating branch density, atom-site, and transformed-fraction slices from cold product back to reheated parent |
A sharp front nucleates and sweeps the domain within a couple of kelvin either direction, and the lattice deforms anisotropically in a pattern that is driven by $\epsilon_i^*=\ln\lambda_i$ as specified in eq.\ref{eq:h_and_eps}. The key though pint is that the model was never handed the transformation temperature.
I wanted to see how sweeping the Hf concentration fixed a fixed Ni concentration of 50.3 at.% would predict the transformation temperatures. Keep in mind Ni-Ti-Hf composition data is held-out of any surrogate model fitting. In Figure 5, $M_s$ rises from 305.7 K at Hf=0% to 535.8 K at Hf=25%, and $A_f$ from 341.7 K to 590.8 K. The profile is a monotonic rise that aligns with Hf's effect in this alloy family [12]. The Hf=0% point is used as a check since binary Ni-Ti data [12] interpolated to Ni=50.3 at.% gives $M_s\approx314$ K and $A_f\approx353$ K, both within about 10 K of the blind prediction. Only the Hf=20% and Hf=0% points have data to check against, so the rest of the trend hasn't been verified against real data.
|
| Figure 5. Predicted Ms and Af versus Hf content at fixed Ni, holdout Ms marked at the one point with an independent label, open triangles are real binary Ni-Ti literature values at Hf=0 |
Overall I think this is prety cool effort and result and makes me want to think more again about phase field methods.
Is this the best we can do?
Honestly, for $a_0$, $\lambda_1$, $\lambda_2$, $\lambda_3$, $K$, and $G$, not entirely, in an absolute sense. There were some gains we found (eqs.\ref{eq:elemental_coef_solve}, \ref{eq:lambda1_stack}, plus widening the elastic-modulus element coverage from 15 to 23 elements) were all cases of unused information already sitting in the data tables we were using. There's probably more of that kind of thing left if I looked harder.
For $T_0$ and $H$ I'm more confident we hit a real wall given the constraints set. For $T_0$ I tried feature selection and PCA dimensionality reduction, neither beat the full descriptor. For $H$ we tried training on individual measurements instead of per-composition medians, weighting by replicate count, swapping its raw fractions for the same compact descriptors that helped everything else, and a lattice-compatibility-theory-motivated feature (predicted $\lambda_2$'s deviation from 1). None of it moved either number outside grouped-CV noise. My honest read is that both are missing processing-history information, heat treatment, grain size, testing protocol, that composition alone cannot supply. So its bascially data starved.
Not a Guarantee
The prompt below is a summarized prompt of several sessions related to working through this side project on having codex/cursor identify the proper available data, surrogate models, phase-field model, and data augmentation. Then how we went about simulating and analyzing the results. With that said I can't guarantee that this will reproduce the exact same results but it should work to initiate the session and get things moving. If you want all the raw data send me an email.
A condensed recipe, written so an agent (or a patient human) could rebuild
this from scratch using only public data and the equations above, no code
from me required.
AI reproduction prompt (click to expand)
References
