Deciphering diffuse scattering with machine learning and the equivariant foundation model: The case of molten FeO

Bridging the gap between diffuse x-ray or neutron scattering measurements and predicted structures derived from atom-atom pair potentials in disordered materials, has been a longstanding challenge in condensed matter physics. This perspective gives a brief overview of the traditional approaches employed over the past several decades. Namely, the use of approximate interatomic pair potentials that relate 3-dimensional structural models to the measured structure factor and its associated pair distribution function. The use of machine learned interatomic potentials has grown in the past few years, and has been particularly successful in the cases of ionic and oxide systems. Recent advances in large scale sampling, along with a direct integration of scattering measurements into the model development, has provided improved agreement between experiments and large-scale models calculated with quantum mechanical accuracy. However, details of local polyhedral bonding and connectivity in meta-stable disordered systems still require improvement. Here we leverage MACE-MP-0; a newly introduced equivariant foundation model and validate the results against high-quality experimental scattering data for the case of molten iron(II) oxide (FeO). These preliminary results suggest that the emerging foundation model has the potential to surpass the traditional limitations of classical interatomic potentials.


Understanding disordered structures.
Finding the link between the measured static structure factor from a disordered material, and the predicted structure derived from individual atom-atom interatomic pair potentials, is a long-standing problem in condensed matter physics [1,2].The suitability of classical interatomic potentials in accurately predicting structure are typically assessed by comparing the Sine Fourier transform of the measured diffraction data, the total pair distribution function (PDF), to the appropriately weighted sum of the partial atom-atom PDF's.This comparison is straight forward in the case of neutron scattering where the coherent scattering lengths are constant with momentum transfer, Q, but requires consideration of the Q-dependent electronic form factors in the case of x-rays.Nevertheless, both neutron and x-ray diffraction data typically serve as rigorous tests of model structures predicted by simulation methods [3].
Ornstein-Zerneike have argued that the total correlation function is a direct effect of molecule 1 on molecule 2, plus an indirect effect of all other molecules [4].The Ornstein-Zerneike equation can be extended using angular dependent pair distribution functions, from equations such as Hypernetted Chain Theory [5] or Percus Yevik theory [6] for hard spheres.Solving the angular dependent pair correlation function can be achieved by a series expansion of spherical harmonics, whereby the potential energy is the sum of pair interactions between atoms within different molecules.Alternatively, the use of approximate pair potentials in molecular dynamics or Monte Carlo simulations can offer approximate model structures that capture the underlying physics and chemical bonding in crystals, liquids and glasses.
A notoriously difficult example of a binary liquid x-ray PDF to measure experimentally as well as model, is shown in figure 1.The experimental data are those previously reported by Shi et al. [7] for molten Fe49O51 at 1673K containing 91%Fe 2+ .Here the qualitative local polyhedral bonding and connectivity in molten FeO is reproduced by several classical Molecular Dynamics (MD) interatomic potentials, including Born-Mayer [8][9][10], Buckingham [11], Interionic [12] and Morse potentials [13].However, none of these interatomic potentials are in quantitative agreement with predicting the dominant interactions associated between Fe 2+ and O 2-.The multivalent nature of Fe bonding depends several factors, including oxygen partial pressure, temperature and composition.Fe-O interactions are of fundamental importance in iron and steelmaking, where oxidation states and coordination numbers can have a direct effect on phase stability, melt viscosity, density and heat capacity.were weighted by the Q-dependent x-ray form factors and the total S(Q) truncated at a Qmax=20Å -1 .The MD models were calculated using the potentials of (a) [9] (b) [11] (c) [12] (d) [13] (e) [10] (f) [8].The R-factor for each MD model is calculated with respect to the experimental data in T(r) between rmin=1Å and Rmax=10Å, where a lower Rfactor indicates better agreement with experiment.Here the D(r) representation is shown rather than T(r) to highlight the structural differences.
Typically, a -squared analysis between the computed and measured structure factor, S(Q), [14,15] or R-factor analysis [3,16] between the model and experimental PDF's, is used to assess the validity of any structural model.Here the differential distribution function D(r) is related to S(Q) through a Sine Fourier transform using the Hannon-Howells-Soper formalism [17], and, where where is the average x-ray scattering form factor squared and the total distribution function T(r) includes the radial bulk density term.Comparisons of the R-factor in real space are made between the minimum and maximum atom-atom interaction distances T(rmin<r<rmax).At distances shorter than rmin, oscillations in the experimental data can arise for number reasons, including systematic errors in the measurement, Fourier transform artifacts due to the finite maximum Q, and in the case of x-ray diffraction, approximations in the electron cloud distribution used to analyze the scattering data.Nonetheless these unphysical oscillations below rmin should oscillate about the density line (see figure 1).For these reasons, Reverse Monte Carlo fits are typically refined against the S(Q) reciprocal space data [18], which has the added benefit of a straight forward comparison with statistical experimental errors.

Standard methods for modelling scattering data.
The development of parametrized classical interatomic potentials routinely used in Molecular Dynamics (MD) and Monte Carlo simulations are often based on crystalline properties, lattices and bond lengths.Popular potential types include Lennard-Jones [19], Stillinger-Weber [20], embedded-atom method [21], CHARMM [22] and AMBER [23] force fields, among many others.Although the large system sizes of these simulations enable the calculation of a systems bulk properties, a main aim of fitting potential functions is to make the potential transferable, such that it can be used to describe a different materials' structure and properties other than the one it was designed for.However, there is an inconsistency in the use of the effective pair potential approximation, since it may not effectively capture quantum mechanical phenomenon such as polarizability or many body effects, etc.While pair-theory or fixed functional form approximations greatly simplify and speed-up the computational calculations, such quantum effects are typically accounted for by using approximations to create effective pair potentials.For more accurate quantum mechanical modeling, and the prediction of both the electronic and atomic structures, the interactions between electrons needs to be considered explicitly.Density functional theory (DFT) with systematic advancements in the exchange-correlation approximations [24] has provided an effective approach for quantum mechanical treatment of materials.DFT based MD simulations have played an important in accurate modeling of condensed phases [25].However, even the DFT based MD simulations are often limited to short simulation times (~ 10's of ps), and only a small number of atoms (100's of atoms) can be treated.
By approaching the problem from a tangential direction, Monte Carlo fits can provide a 3-dimensional structures perfectly consistent with the measured data.Methods such as Reverse Monte Carlo (RMC) [18] and Empirical Potential Structure Refinement (EPSR) [15] iteratively move atoms and/or molecules, with a given probability to avoid local minima, that ultimately improve agreement with the scattering pattern.However, these methods often lack knowledge of the underlying physics, restricting them to structure and density refinements alone.These Monte Carlo approaches are also subject to ambiguities caused by uniqueness, whereby vastly different structural models can be made to adequately fit the diffraction data with the correct density [26,27].The empirical potentials used in EPSR in particular are a means to drive the fitting procedure based on the difference between the model and experiment, and not meaningful in the sense of predicting the thermodynamic properties or atomic dynamics of a material.Nonetheless, EPSR has inspired our current efforts described in this perspective to try and complete the link; using modern computing tools to create and guide realistic interatomic potentials that agree with experiment while maintaining quantum mechanical level accuracy.An overly simplistic overview of the correlation between experimental accuracy and theoretical accuracy of existing methods and the objective of developing machine learned interatomic potentials is illustrated in figure 2.

A new perspective on understanding disorder.
An alternative to addressing the accuracy limitations of effective pair potentials while also improving the spatiotemporal limitation of DFT-MD simulation is to directly machine learn the potential energy surfaces from ab initio datasets.Unlike the empirical interatomic potentials, these machines learning interatomic potentials are inherently high dimensional, which allows them to achieve accuracy as afforded by their training reference method.Many machine learning interatomic potential (MLIP) architectures have been developed since the past decade, with some of the most relevant being the Behler-Parrinello neural network [28], the Gaussian approximation potential (GAP) [29], the Spectral Neighbor Analysis Potential (SNAP) [30], the moment tensor potential (MTP) [31], ANI [32,33], SchNet [34], DeepMD [35], NequIP [36], MACE [37] and Allegro [38].Furthermore, addition of long-range interactions has also been proposed in recent studies using MLIP's, which would allow them to model multiply charged systems or interfaces, etc. [39,40].
The early development of the MLIP's (except for ANI) were predominantly focused on simulations of specific chemical systems or materials under different conditions.A bottleneck to training these MLIP's was related to efficient sampling of the training configurations, as labelling these configurations still required expensive DFT calculations.To solve this challenge, active learning strategies were devised to enable efficient sampling of the most relevant minimum training configurations that would lead to optimal coverage of MLIP for the chemical space of interest.Active learning methods were devised [41] to construct quantum mechanically accurate and transferable machine learningbased models of the potential energy surface for the molecular modelling of materials.The method consists of three main steps: phase exploration, then query a small subset of that explored region for labelling, and training of the interatomic potential.These active learners have been shown to dramatically lower the number of labels needs to train the MLIPs by improving the diversity of training space.Machine learning interatomic potentials (MLIP) can maintain near ab initio accuracy, and enable larger system sizes through linear scaling, together with time scales comparable to classical interatomic potentials.MLIP's based on the Gaussian Approximation Potential (GAP) have been successfully trained to model liquids [42], crystals [43], defects [44], amorphous [45], multi-component materials [46], and molecules [47].
Sivaraman et al. [48] have extended the MLIP methodology by using experimental scattering data measured over a wide range of phase space to drive the procedure.An active learning scheme, initialized by model structures, iteratively improves a MLIP until agreement is obtained with experiment within a specified quantum mechanical accuracy (see next section).Using the experimental data to guide the DFT calculations is particularly important for metastable materials such as supercooled, liquid, amorphous or glassy materials, where the system is not in the lowest energy state.A critical pre-requisite of this approach is that the scattering data are of the highest quality.The appropriate experimental corrections and normalization procedures have been well documented for both neutron [49,50] and x-ray diffraction [51,52].However, disordered materials, unlike crystals, do not have a periodic lattice to verify the underlying structural order.Rather, consistency checks are used verify the accuracy of the data, including that the low-r data in the PDF oscillates around the density line and the (dis-)agreement with the measured S(Q) and the Fourier back-transform gives an assessment of the level of systematic error [49,52].A primary goal of this methodology is to develop interatomic potentials that encompass all condensed matter phases, both solid and liquid, that are transferrable.The workflow for mapping MLIP's is shown in figure 3.

Experimentally driven MLIP's.
Experimentally, carefully analysed x-ray and neutron diffraction data can provide rigorous benchmarks for discerning between different structural models of disordered materials.Progress toward the automated fitting of MLIP's using experimental PDFs and active learning have been made in this regard using two approaches.Firstly, reproducing the high temperature phases of the refractory oxide HfO2 [48], has been made using an automated closed loop via an active learner.Crystal phases based on x-ray and neutron diffraction data are used to initialize and sequentially improve an unsupervised machine-learning model over a predetermined range of phase space.The resulting MD simulations were able to reproduce all the experimental phases with near ab initio precision [48].Particularly, the experiment driven MLIP work provided clarity on the ambiguous high-temperature phase behaviour of hafnia immediately preceding melting, demonstrating agreement between simulations and experiments on the mixed cubic and tetragonal phases near the melting point.In the second approach, the initial structure of ionic liquids has been created using classical MD simulations, and down-sampled by using an active learning algorithm.Subsequently, iterative DFT calculations are performed, and a MLIP potential developed until a specified energy threshold is reached.Atomistic models developed using this method have been applied to molten LiCl [53], NaCl [44] and LiCl-KCl [54].Such approaches provide new physical insights into the temperature-dependent coordination environment of liquids, together with property information including density, selfdiffusion constants, thermal conductivity, and ionic conductivity [55].
In the example of liquid NaCl, a GAP model developed by Tovey et al. [44] was trained with 1000 atomic configurations and obtained a DFT accuracy of within 1.5 meV/atom.The MLIP enabled the analysis of the hightemperature molten salt properties on large systems (~10,000 atoms) and longer time scales (>1 ns), currently inaccessible to ab initio simulations.Here, the GAP model reciprocal and real space functions are compared to subsequently measured partial S(Q) and PDF's obtained from Neutron Diffraction Isotopic Substitution (NDIS) measurements, essentially a double difference method [56], see figure 4.An R-factor analysis yields 3-4% at the partial PDF level, well within the experimental errors, and a considerably more rigorous test of any model than that provided by the total neutron or x-ray pair distribution functions.In addition, the GAP model is in very good agreement with previously published experimental diffusion constants, with a root-mean-square deviation of <0.05×10 −8 m 2 s −1 [44].Consequently, these results demonstrate that GAP models are able to accurately capture the many-body interactions necessary to model both the structure and dynamics ionic systems.However, limitations in DFT reference dataset or the approximation can lead to MLIP's to predict overstructured liquid structures [57], leading to a deviation with respected to experimentally measured melt structures.Matin et al. [58] introduced a novel method to refine machine learning potentials by incorporating experimental observations, specifically focusing on the melt phase of pure aluminium.This method leverages iterative Boltzmann inversion, allowing for the integration of experimental radial distribution function data at specific temperatures to act as a correction to DFT trained MLIP's.The addition of this pair correction fixed over-structuring issues seen in melt structures predicted by MLIP's without correction.Their results [58] also showed an enhancement in the prediction of diffusion constants with the added correction.However, the authors left the fitting of RDF's from different temperatures simultaneously for providing continuous corrections to DFTtrained MLIPs for future work.

High performance computing workflow for MLIP fitting over combinatorial chemical spaces.
The discussion so far has centred on the development and fitting of MLIPs for specific chemical systems.We now shift our focus to recent advancements in modelling or the rapid generation of MLIPs for arbitrary chemical systems, with a particular emphasis on disordered melts.In recent work, Guo et al. [59] introduced a high-performance active learning workflow termed AL4GAP.This workflow is designed to generate compositionally transferable MLIPs over chargeneutral mixtures of arbitrary molten mixtures, spanning 11 cations (Li, Na, K, Rb, Cs, Mg, Ca, Sr, Ba, Nd, and Th) and 4 anions (F, Cl, Br, I).The authors demonstrate the efficacy of this workflow by generating and validating GAP-based MLIP models with DFT-SCAN level accuracy for five systems of increasing complexity: LiCl-KCl, NaCl-CaCl2, KCl-NdCl3, CaCl2-NdCl3, and KCl-ThCl4.These hightemperature melts pose significant characterization challenges due to corrosion and radiation, complex sample environments.The emergence of such a powerful framework will aid challenging x-ray synchrotron and neutron experiments in focusing their resources on the most promising melt compositions and conditions, informed by feedback from MLIP-driven MD simulations.

Emergence of foundational models in materials and chemistry.
The materials project based [60], encompassing ~1.5 million configurations across 89 chemical elements, has set the stage for the next generation of foundational models in materials science, including M3GNet [61], CHGNet [62], and MACE-MP-0 [63].These foundation models hold the promise of development of single large models that can learn from increasingly large chemical spaces and elements.Notably, the MACE-MP-0 model demonstrated its versatility in modelling across various chemistries in solid, liquid, and gaseous states.
In the Figure 1, we showed that FeO melt in a challenging system for many classical interatomic potentials.Here we have performed the original simulation for melting of 800 atoms of FeO using MACE-MP-0 foundation model using [64].The starting structure is created using PACKMOL [65] consisting of 100 atoms of FeO using the starting density consistent with the experimental x-ray PDF [7].The starting structure is melted at 2500K and cooled to a target temperature of 2000K.The final structure is replicated to 800 atoms and simulations are performed in NVT ensemble.The simulations are performed using atomic simulation environment [66].The production simulations are performed for 150ps and the PDF's are computed from the last 125 ps.Additionally, we have also computed the PDF's using the MACE-MP-0 model with dispersion correction (labelled as MACE-MP-0+D3) [67].This simulation of molten FeO was not aimed at comprehensiveness, but to offer insights into emerging opportunities such as foundation model through original simulation results.The results are visualized in Figure .5. The R-factor analysis show that MACE-MP-0 simulation shows an improved agreement with experiment in comparison to every other empirical forcefield reported in Figure 1.We should note that the simulations performed in the figure 1 using the empirical forcefield were carried out using a large simulation cell of ~6000 atoms and we would presume that a larger simulation cell could only further improve the agreement for MACE-MP-0 predicted structure with experiments.In addition, it was observed that inclusion of dispersion interactions did not lead to any improvement in the predicted structure.

MLIP prospects and difficulties.
The development of accurate transferrable interatomic potentials has been a longstanding challenge in condensed matter physics.Using scattering experiments to help construct interatomic potentials is not a new idea [2], but with the advent of supercomputers, machine learning, and active learning methods, it is becoming an achievable goal.The results described in this perspective relate primarily to ionic and oxide systems, but they already indicate that MLIP models can capture the many-body interactions required to develop large scale models with quantum mechanical accuracy.Already, experimentally driven active learning methods in the field of molten salts can significantly lower the barrier to the understanding and design of materials across the periodic table.Similar active learning approaches using MLIP's applied to more complex and multi-component materials are expected to uncover new physics, particularly in disordered phases.In this regard, we argue that using experimental scattering data to drive the MLIP testing and training is essential, as the system under investigation may not necessarily be in the lowest energy state.
By scaling up the number of atoms in the simulation box, not only can material properties be more accurately predicted, but much slower quench rates of the system can also be attained.This is particularly useful in the study of metastable states such as supercooled liquids and glasses.Nonetheless, specific challenges arise when modelling metastable amorphous and glassy forms, that can cause the active learning scheme to fail.The limited sampling of configuration space in high dimensionality systems with significant free energy barriers could inhibit access to the necessary atomic arrangements.Initial attempts to widen range of possible configurations using meta-dynamics approaches [68] in the case of amorphous hafnia have not proved successful [69].However, such schemes to increase the diversity of the local chemical environments on which the model is trained is the likely key to future progress.
For stable systems, based on the promising results for liquid FeO shown here, combining a AL4GAP type approach (to sample a vast range of targeted but diverse chemical space) to fine tune the foundation model at a higher accurate density functional theory approximation, along with direct infusion of input from experimental scattering data is a powerful tool.Moreover, this methodology could pave the way in understanding the underlying physics of structure-property relations in a multiplicity of disordered materials.

Figure 1 .
Figure 1.The measured x-ray differential pair distribution function, D(r), of liquid FeO (circles) compared to classical

Figure 2 .
Figure 2. Simplistic overview of scattering model methodologies.

Figure 3 .
Figure 3. Workflow for mapping machine learned interatomic potentials across phase space.(a) Sample the configuration space.(b) Perform single point DFT for the AL samples and fit the ML model.(c) Enrich the configuration space by using meta-dynamics on the ML based MD.(d) Diffraction experiments (e) Perform rigorous validation of ML driven MD simulation using pair distribution functions.

Figure 4 .
Figure 4. Liquid NaCl GAP model compared to the measured (a) x-ray structure factor and (b) partial pair structure factors determined from neutron diffraction isotopic substitution (NDIS) [56] (c) associated partial differential pair distribution functions and their associated R-factors calculated between rmin=2Å and rmax=10Å.

Figure 5 .
Figure 5.The measured x-ray differential pair distribution function, D(r), of liquid FeO at 1673K (circles) compared to the MACE-MP-0+D3 model with dispersion interactions (dashed blue line) and MACE-MP without (solid red line) at 2000K.