Predicting neutron experiments from first principles: a workflow powered by machine learning

Abstract

Machine learning has emerged as a powerful tool in materials discovery, enabling the rapid design of novel materials with tailored properties for countless applications, including in the context of energy and sustainability. To ensure the reliability of these methods, however, rigorous validation against experimental data is essential. Scattering techniques—using neutrons, X-rays, or electrons—offer a direct way to probe atomic-scale structure and dynamics, making them ideal for this purpose. In this work, we describe a computational workflow that bridges machine learning-based simulations with experimental validation. The workflow combines density functional theory, machine-learned interatomic potentials, molecular dynamics, and autocorrelation function analysis to simulate experimental signatures, with a focus on inelastic neutron scattering. We demonstrate the approach on three representative systems: crystalline silicon, crystalline benzene, and hydrogenated scandium-doped BaTiO₃, comparing the simulated spectra to measurements from four different neutron spectrometers. While our primary focus is inelastic neutron scattering, the workflow is readily extendable to other modalities, including diffraction and quasi-elastic scattering of neutrons, X-rays, and electrons. The good agreement between simulated and experimental results highlights the potential of this approach for guiding and interpreting experiments, while also pointing out areas for further improvement.

---

1 Introduction

Advancements in materials science are pivotal for technological progress, driving innovations in energy storage, electronics, and catalysis. Computational methodologies, particularly density functional theory (DFT), have become essential tools in materials discovery by predicting materials properties and guiding experimental efforts.^1–5 The integration of machine learning (ML) with these computational techniques has further accelerated the discovery of novel materials, e.g., by enabling rapid screening of vast chemical spaces.^3,4,6–9 ML has also facilitated the development of machine-learning interatomic potentials (MLIPs), which allow accurate and efficient atomic-scale simulations, bridging the gap between empirical potentials and first-principles methods.^10–12

However, the predictive power of these computational approaches necessitates rigorous experimental validation. Scattering experiments, such as neutron, X-ray, and electron scattering, provide critical insights into the structure and dynamics of materials but require precise simulations to interpret the data accurately.¹³ Bridging the gap between computational predictions and experimental observations remains a significant challenge in the field. Predictive simulations could also significantly enhance experimental planning and execution by ensuring that data acquisition is optimized for maximum information gain while reducing the likelihood of inconclusive or ambiguous results.^14–16 Furthermore, such simulations can support the preparation of beamline proposals, providing quantitative justifications for instrument time requests by demonstrating expected signal strengths and resolving power. As experimental facilities increasingly integrate computational tools into their workflows, predictive capabilities are poised to play a crucial role in streamlining the experimental process, improving the overall efficiency of materials characterization, and ultimately accelerating scientific discoveries.

In response to these challenges, we here describe a comprehensive workflow that integrates DFT calculations, MLIPs in the neuroevolution potential (NEP) format, MD simulations using GPUMD,¹⁷ the computation of autocorrelation functions via dynasor,^18,19 and their convolution with atomic form factors, instrument resolution functions and kinematic constraints. This enables instrument-specific predictions of scattering data from first-principles, allowing direct comparisons between simulations and experimental measurements. MLIPs are instrumental to this workflow, as they enable the large-scale MD simulations required to properly converge the density auto-correlation functions, and thus the experiment predictions.

The computational efficiency of NEPs as implemented in GPUMD allows our workflow to be applied to systems containing tens of thousands of atoms, simulated over several nanoseconds. To the best of our knowledge, this allows the access of both larger system sizes and longer simulation times than other workflows predicting inelastic neutron scattering (INS) spectra with MLIPs.²⁰ Furthermore, the workflow is fully implemented in Python, offering direct and easy integration in existing computational workflows. We demonstrate the efficacy of this workflow by applying it to three example systems, including elemental Si, crystalline benzene, and hydrogenated Sc-doped BaTiO₃, showcasing both its potential for guiding experimental design and accelerating the discovery of new materials as well as its current limitation. We focus specifically on simulating INS experiments, but the general workflow can be easily used to simulate other experimental modalities, including diffraction as well as quasi-elastic and inelastic scattering of neutrons, X-rays, and electrons.

2 Methods

The workflow that we demonstrate for predicting neutron scattering experiments from first principles consists of three steps. The first step is the construction of MLIPs based on the NEP framework (Fig. 1a; Subsection 2.1). These accurate and efficient MLIPs enable the second step of the workflow, which are large-scale MD simulations (Fig. 1b; Subsection 2.2). The MD trajectories that result from the second step are then used in the third step of the workflow, in which we compute the dynamic structure factor that is then weighted by neutron scattering lengths as well as the instrument resolution function and kinematic constraint (Fig. 1c; Subsection 2.3). The weighted dynamic structure factor can then be compared directly to experimental data.

Fig. 1 Workflow for simulating neutron scattering experiments from first principles. (a) The first step of the workflow comprises constructing machine-learning interatomic potentials (MLIPs) using an iterative cycle combining both data generation and model validation. Here, the training is facilitated by the GPUMD and calorine packages. (b) The final MLIP is used in the second step to run large-scale MD simulations using the GPUMD package. (c) The dynamic structure factor is computed from the MD trajectories using the dynasor package. The dynamic structure factor is weighted by species-dependent scattering lengths, and broadened with an instrument-specific resolution function in order to predict the outcome of a particular neutron scattering experiment.

2.1 Construction of the machine-learned potentials

The workflow starts with training of a MLIP, or using an already trained model. Here, we used three different MLIPs based on the NEP framework,^17,21,22 one for each of the systems studied in this work. For Si, the model published in ref. 17 was used, while we constructed new models for crystalline benzene and hydrogenated Sc-doped BaTiO₃ (BaTi_1−xSc_xO₃H_x) using the iterative procedure described in ref. 23 utilizing the GPUMD¹⁷ and calorine packages²⁴ (Fig. 1a).

The training set was initially composed of strained and scaled structures, based on ideal structures using reference data from DFT calculations (Subsection 2.4). In the case of benzene, the initial dataset also included dimer configurations to ensure that intermolecular interactions are captured in the training dataset. An initial model was trained on all available data. The dataset was then augmented with structures from several iterations of active learning. To this end, we trained an ensemble of five models by randomly splitting the data into training and validation sets, which was subsequently used to estimate the model uncertainty. MD simulations were then carried out between 10 and 200 K and at pressures ranging from 0 to 10  GPa for benzene, and from 300 to 2000 K and at pressures ranging from −1 to 10 GPa for Sc-doped BaTiO₃ using the respective current generation of NEP models. The NEP model for Sc-doped BaTiO₃ was trained on structures from the extended temperature range from 300 to 2000 K to increase the robustness of the model, by ensuring that the training data set contains a varied set of configurations. Structures encountered at the target temperature 15 K are well within the interpolative regime of the NEP model (Section S1 in the ESI†). The ensemble of models was used to select structures with a high prediction uncertainty, quantified by a range of predictions over the ensemble, for which we computed reference energies, forces, and stresses via DFT. These configurations were subsequently included when training the next-generation NEP model. The training set for benzene consisted of 798 unique benzene structures, corresponding to a total of 94 470 atoms. For Sc-doped BaTiO₃, the training set contained 2280 unique structures, corresponding to a total of 138 438 atoms. Structures were generated and manipulated using the ASE²⁵ and hiphive packages.²⁶

We obtained the final NEP models after 13 iterations for benzene and 6 iterations for Sc-doped BaTiO₃. The final models were trained on all available data. For the benzene model, the resulting average root mean square errors (RMSEs) over the ensemble are 1.1(8) meV per atom for the energies and 63(30) meV Å⁻¹ for the forces, and 8.4(16) meV per atom for the virials. Note that the relatively large uncertainty in the predicted forces is due to one of the ensemble models being an outlier. The corresponding average coefficients of determination on the same folds are R² = 0.9997(3), R² = 0.9949(51), and R² = 0.9988(6) for energies, forces, and virials, respectively. The RMSEs and R² scores for the final benzene model were 8.510 meV per atom and R² = 0.9998 for the energies, 59 meV Å⁻¹ and R² = 0.9962 for the forces, and 9.3 meV per atom and R² = 0.9987 for the virials (Section S2 in the ESI†). For the Sc-doped BaTiO₃ model, the ensemble RMSEs were 6.6(13) meV per atom for the energies, 186(34) meV Å⁻¹ for the forces, and 31(5) meV per atom for the virials. The respective coefficients of correlation (R²) were R² = 0.999 98(1), R² = 0.9765(58), and R² = 0.9964(9) for energies, forces, and virials, respectively. The RMSEs and R² scores for the final Sc-doped BaTiO₃ model were 6.2 meV per atom and R² = 0.9999 for the energies, 172 meV Å⁻¹ and R² = 0.9792 for the forces, and 30 meV per atom and R² = 0.9963 for the virials (Section S3 in the ESI†).

The resulting NEP models along with the reference data used for training are available via zenodo as specified in the Data Availability statement.

2.2 Molecular dynamics

In the second step of the workflow we perform MD simulations using the MLIPs from the first step for large supercells (Fig. 1b). The resulting MD trajectories are later used to compute the dynamic structure factor as detailed in the next section.

For Si, a supercell comprising 38 × 38 × 38 primitive cells for a total of 438  976 atoms was simulated at 300 K, 900 K, 1200 K, and 1500 K, with equilibration of the system in the NPT ensemble and production for 1 ns in the NVE ensemble, with a timestep of 2 fs. The atomic positions were written to file every 14 fs in order to accurately resolve the fastest vibrations in the system when computing the dynamic structure factor.

Crystalline benzene was simulated in a supercell containing a total of 57 024 atoms. The benzene system was equilibrated in the path-integral molecular dynamics (PIMD) ensemble^27,28 to avoid the significant underestimation of the cell volume in the classical NPT ensemble at low temperatures. The MD simulations were conducted at 127 K to strike a balance between computational cost and the number of PIMD beads (see Sections S4 and S5 in the ESI†). Production runs were then performed for 1 ns in the NVE ensemble. A time step of 0.5 fs was used, and the positions were written every 3 fs. Ten independent MD runs were performed to improve the statistics of the computed dynamic structure factor.

Hydrogenated supercells of Sc-doped BaTiO₃ were constructed for various Sc concentrations in the range 16% to 70 % in both the cubic and hexagonal phase. The supercell contained ≃40 000 atoms. Equilibration was performed in the PIMD ensemble at a temperature of 15 K, and production was carried out for 350 ps in the thermostated ring-polymer MD ensemble²⁹ with a timestep of 0.5 fs. This approach captures nuclear quantum effects on the frequencies,²⁸ but it should be noted that the phonon occupation statistics are still classical. Both equilibration and production runs used 32 PIMD beads, for an effective system size of ≃13 00 000 atoms, limiting the length of the production run compared to Si and benzene because of the increased computational cost.

The specific supercell sizes used in this work were chosen in order to strike a balance between computational cost and convergence of the dynamic structure factor with regards to the number of q-points commensurate with the supercell (see Section S6 in the ESI† for a convergence study of the supercell size for crystalline benzene, as well as an extended discussion on supercell size and commensurate q-points).

2.3 Auto-correlation functions and instrument-specific kinematic constraints

The central quantity analyzed in the third step of the workflow is the dynamic structure factor, S(q,ω). S(q,ω) is directly proportional to the intensity measured in scattering experiments, and can be readily extracted from MD simulations. While the procedure has been described in detail in ref. 18 and 19 we briefly summarize it here for completeness. Let n(r,t) denote the particle density defined asr_i(t) is the position of particle i at time t, and N is the total number of particles. The particle density can now be Fourier transformed in space,with the autocorrelation function of n(q,t) yielding the intermediate scattering function F(q,t),where the brackets denote an ensemble average. The intermediate scattering function can then be Fourier transformed in time to yield the dynamic structure factor,

The dynamic structure factor in eqn (4) can be further generalized for multi-component systems. Different atomic nuclei scatter neutrons, X-rays, and electrons with varying intensity, which can be taken into account by weighting the partial dynamic structure factor for species α and β accordingly. In the case of neutrons, the partial dynamic structure factor should be weighted by the scattering lengths, b_α and b_β,

The dynamic structure factor in eqn (5) was computed from the MD trajectories using the dynasor package^18,19 in the third step of the workflow (Fig. 1c). q-points and time lags were selected to match the accessible range of the simulated neutron scattering instruments. Specifically, for Si a Brillouin zone path was sampled connecting the high-symmetry points Γ, X, K, and L. The path was sampled in 52 different Brillouin zones, randomly selected from the first zone up to |q| = 12 Å⁻¹ for a total of 6136 q-points. Randomly selected q-points up to a magnitude |q| = 14 Å⁻¹ and |q| = 18 Å⁻¹ were sampled for benzene and Sc-doped BaTiO₃, respectively, yielding 2116 and 2601 q-points, respectively. Gaussian broadening with a width of 0.01 Å was then applied to each q-point, followed by averaging over spherical shells in |q| to produce S(q,ω).

Instrument-specific resolution functions and kinematic constraints were applied to the calculated spectra using the euphonic package³¹ with the resolution functions defined in the resins package.³² The resolution functions used here are Gaussians with energy-dependent width; the functions for TOSCA and Lagrange are based on implementations in AbINS, and the functions for MAPS and ARCS are based on PyChop.³³ (The instrument functions for both AbINS and PyChop are distributed in Mantid.^34,35) Note that the true resolution functions are four-dimensional and non-Gaussian, but these 1-D approximations are used routinely in INS simulations. The kinematic constraints have their origin in the instrument geometry and transformation from time-of-flight measurements to (q,ω) space. In the simulations they are applied as a mask to data computed directly in the (q,ω) space.

Finally, a quantum correction factor was applied to all dynamic structure factors, in order to correct for the classical phonon statistics generated by the MD simulations. Specifically, we applied the following correction factor based on first-order Stokes-Raman scattering,^36,37

2.4 Density functional theory calculations

To generate reference data for the construction of the MLIPs (Section 2.1) we performed non-spin polarized DFT calculations using the projector augmented wave method^38,39 as implemented in the Vienna ab initio simulation package^40–42 with a plane wave energy cutoff of 520 eV using the vdW-DF-cx exchange correlation functional⁴³ for benzene and the r²SCAN functional⁴⁴ for Sc-doped BaTiO₃. The Brillouin zone was sampled with automatically generated Γ-centered k-point grids with an approximate spacing of 0.25 Å⁻¹ and the partial occupancies in each orbital were set using Gaussian smearing with a width of 0.1 eV. The DFT data are available via zenodo as specified in the Data Availability statement.

2.5 Inelastic neutron scattering experiments on crystalline benzene

For validation of the predictions for crystalline benzene, inelastic neutron scattering experiments were performed at the TOSCA neutron spectrometer^45,46 at the ISIS Neutron and Muon Source. The liquid sample was placed in a 1 mm thick standard flat TOSCA aluminum cell which was then briefly submerged into liquid nitrogen. As soon as the sample solidified it was quickly transferred into the TOSCA closed cycle refrigerator and allowed to further cool to the cryostat base temperature below 10 K. The short INS measurements (approximately 8 minutes per spectrum, i.e., total exposure of 20 μA h) were performed as part of a cooling run at a rate of 3 K min⁻¹, with the initial spectrum taken at a starting temperature of 127 K and followed by other measurements at a starting temperature of 103 K, 75 K, 46 K, and 24 K. The longer INS measurement (approximately 2 hours, i.e., total exposure of 285 μA h) was performed at the base temperature of 10 K, giving a superior spectral signal-to-noise ratio. The raw data, i.e., time-of-flight events, were reduced using Mantid.^34,35

2.6 Post-processing and plotting

The NEP models and calculated correlation functions were post-processed and analyzed using Python scripts, utilizing the NumPy,⁴⁷ Pandas,^48,49 and SciPy⁵⁰ packages. Plots were generated using matplotlib,⁵¹ with color maps from perfect-cmaps.⁵² Atomic structures were visualized and analyzed using OVITO.⁵³

3 Results

3.1 Anharmonicity in Si

We begin by applying the workflow outlined in the methodology section to simulate an INS experiment on crystalline Si at 300 K reported in ref. 30, that was carried out at the ARCS wide range angular spectrometer (BL-18) at the Spallation Neutron Source (SNS) at Oak Ridge National Laboratory (Fig. 2a). The simulation is made instrument-specific by applying the resolution function and kinematic constraint for ARCS to the simulated dynamic structure factor. The q and energy range measured by ARCS is relatively broad (Fig. 3), and thus we computed the dynamic structure in multiple different Brillouin zones to accurately sample the full range allowed by the instrument. In total, 52 Brillouin zones were sampled. The dynamic structure factors S(q,ω) were then weighted by a factor ∝ 1/|q|² since to first order the scattered intensity grows as |q|². Furthermore, the Debye–Waller factor, exp(−q²U/3), was corrected for in each of the Brillouin zones before the zones were averaged together. This expression for the Debye–Waller factor assumes an isotropic displacement of the atoms in all Cartesian directions with U = 〈u²〉 being the mean squared displacement in the system.⁵⁴U was estimated to be 0.013 Ų from a 100 ps MD simulation of the Si system at 300 K, otherwise following the same protocol as the other simulations of Si in this work.

Fig. 2 (a) Simulated INS dispersion of Si from MD for the ARCS spectrometer at the Spallation Neutron Source at ORNL, with the harmonic phonon dispersion calculated using the NEP model overlaid (turquoise lines). (b) The simulated intensity at the X-point, as well as (c) the first peak in the intensity at the X-point in the third Brillouin zone for 300 K, 900 K, 1200 K, and 1500 K (red lines). Note that the dispersion in (a) and intensity in (b) is aggregated over 52 different Brillouin zones, while in (c) the results for only a single Brillouin zone are shown. The experimental data (black circles) is from ref. 30. The MD simulation captures both the anharmonicity and multi-phonon effects present in the experimental data, as well as the mode softening as the temperature is increased. The multi-phonon effects manifest as non-zero intensity between the acoustic and optical branches at the X-point.

Fig. 3 Kinematic constraints for the four neutron instruments simulated in this work. ARCS (BL-18) is a wide-range spectrometer at the Spallation Neutron Source at the Oak Ridge National Laboratory (USA). TOSCA and MAPS are spectrometers at the ISIS Muon and Neutron Source at the STFC Rutherford Appleton Laboratory (UK). Finally, IN1 Lagrange is a spectrometer at the Institute Laue–Langevin in Grenoble (France). Note that all spectrometers differ in q and energy range and resolution, owing to their respective kinematic constraints and resolution functions.

We validate our results by comparing them to the experimental data measured on the ARCS spectrometer in ref. 30. Specifically, we consider the intensity at the high-symmetry X-point (Fig. 2b), where the simulated intensity has been multiplied by an extra factor of ω². Our results are in quantitative agreement with experiments, with the centroids of the phonon mode peaks agreeing well. The relative intensity between the different phonon mode peaks in the experiment are not entirely reproduced in the simulation, which could be due to a missing correction factor or experimental variability.

We note, in particular, the nonzero intensity measured in the experiment and captured by the simulation in the region 30 meV to 50 meV. This scattered intensity corresponds to multi-phonon effects, which are inherently captured by MD simulations. Furthermore, the effect of thermal expansion as the temperature is varied is also directly included by the MD simulations, where specifically the low-energy mode at the X-point is softened as temperature is increased (Fig. 2c). One can observe that some of the predicted mode energies are slightly shifted compared to experiments, by approximately 1 meV. Given that the MLIP accurately reproduces the harmonic phonon dispersion from DFT around the X-point, it is most likely due to the underlying exchange-correlation functional (Section S7 in the ESI†).

At 1500 K it moreover appears that in the region around 40 meV the simulated and experimental intensities differ. This discrepancy could be due to the Brillouin zone in which the simulations have been conducted. Here, we show the X-point in the third Brillouin zone, i.e., q = [0.5, −1.0, −1.5], while in the experimental ref. 30 the exact q-point is not specified. In fact, the intensity at the X-point varies substantially depending on the Brillouin zone, especially the intensity of the multi-phonon shoulder around 30 meV (Fig. 4). For the comparison shown in (Fig. 2c), we selected the Brillouin zone for which the simulated spectrum best reproduces the experimental data, based on the mean-squared error calculated over the spectrum.

Fig. 4 (a) Simulated intensity at the X-point at 1500 K in crystalline Si for each of the 52 Brillouin zones accessible by the ARCS spectrometer. Note that the intensity of the multi-phonon shoulder at 30 meV varies greatly with the Brillouin zone. (b) The kinematic constraint for the ARCS spectrometer at SNS, which limits the range of non-zero intensities in (a) depending on the |q| for the X-point in each Brillouin zone.

In MD, the dynamics of the system described by the potential model is captured at the classical level, including high-order phonon effects, thermal expansion, and full anharmonicity. Efforts have been made in recent years to include the effects of anharmonicity on top of harmonic models, including but not limited to using higher-order force constants,²⁶ temperature-dependent effective potentials,⁵⁵ anharmonic lattice models,⁵⁶ and the self-consistent harmonic approximation.⁵⁷ However, a harmonic model is inherently limited in describing such complicated dynamic events.

3.2 Corrections in crystalline benzene

We now turn to simulating an INS measurement of crystalline benzene at 127 K at the TOSCA spectrometer at the ISIS Neutron and Muon Source (UK) in order to study the effects of the resolution function and quantum correction in more detail (Fig. 5a).

Fig. 5 (a) Simulated INS spectra for crystalline benzene, at increasing levels of refinement, compared to an experimental spectrum obtained at 127 K at the TOSCA spectrometer at the ISIS Neutron and Muon Source (UK). The first level of accuracy is the raw simulated spectrum from MD with only scattering lengths applied (raw spectrum). Correcting for the resolution function and kinematic constraint of the TOSCA spectrometer yields a marked increase in accuracy, with further improvement when additionally applying a quantum correction factor to compensate for the classical statistics in MD simulations, especially for low energies around 10 meV. The spectra have been individually scaled to match the experimental spectrum as closely as possible above 50 meV. (b) Phonon dispersion and density of states compared to the simulated INS spectrum at 127 K, corrected for the kinematic constraint of the TOSCA spectrometer and with quantum statistics.

The scattered intensity from benzene is dominated by incoherent scattering from hydrogen, owing to the exceptionally large incoherent scattering length of hydrogen. We thus study the INS spectrum directly. The dynamic structure factor can be integrated over |q| in order to obtain Comparing the raw simulated spectrum with the experimental data, we find that the simulated spectrum captures the peaks corresponding to different modes but the relative intensity between them is not reproduced. Furthermore, the low-energy peak at 10 meV is not captured. The reason for these discrepancies is to a large extent due to the kinematic constraint and resolution function of the TOSCA spectrometer (Fig. 3). The two detector banks of TOSCA map out two lines in q–ω space, where high (low) frequencies correspond to large (low) q. By sampling along these q–ω lines and convoluting the resulting spectrum with the resolution function of the instrument, the agreement improves notably.

However, the ratio in intensity between the high and low-energy regions is still not reproduced. The main reason for this discrepancy is the classical statistics of MD simulations, which we correct for with the quantum correction factor according to eqn (6). Applying both kinematic constraint and quantum correction yields a simulated spectrum that is in near-quantitative agreement with experiments. The remaining difference to experiments is a redshift of the simulated spectrum by approximately 25 meV. We attribute this redshift to the DFT functional used to train the NEP model, as well as weak intermolecular interactions not being fully captured by the NEP model. A more detailed discussion comparing experiments and first-principles calculations to the predictions from the NEP model can be found in the ESI† (Section S8).

We can further elucidate the simulated INS spectrum by comparing it to the phonon dispersion according to the underlying MLIP (Fig. 5b). The simulated INS spectrum differs notably from the phonon density of states in terms of intensity, owing to the kinematic constraint of the TOSCA spectrometer, and the quantum correction. Furthermore, the full anharmonicity included in the MD simulation in combination with the TOSCA resolution function yields a broadening of the peaks in the simulated INS spectrum.

In summary, this study of crystalline benzene highlights the importance of considering the resolution and kinematic constraints of the specific instrument, as well as correcting the statistics from classical MD simulations, when aiming for quantitative predictions of neutron scattering experiments.

3.3 Hydrogen dynamics in hydrogenated Sc-doped BaTiO₃

Finally, we turn to a more complicated system, in the form of hydrogenated Sc-doped BaTiO₃ (BaTi_1−xSc_xO₃H_x) where x is the doping fraction of the tetravalent site (Ti, Sc). Perrichon et al. have performed a detailed INS study of the hydrogen dynamics in this system at three different spectrometers: the TOSCA and MAPS spectrometers at the ISIS Neutron and Muon source as well as IN1 Lagrange at the Institut Laue–Langevin.⁵⁸ The experiments were carried out at temperatures below 20 K. INS spectra were then obtained by averaging the dynamic structure factor S(q,ω) up to a magnitude of |q| = 12 Å⁻¹. Our simulations presented in this section were averaged over q up to the limit of the kinematic constraint for MAPS, |q| = 18 Å⁻¹, in order to obtain better statistics.

Sc-doped BaTiO₃ undergoes a phase transition from a hexagonal structure to a cubic perovskite structure as the Sc concentration increases. On MD time scales, both structures are, however, at least metastable over the entire composition range, which (in contrast to experiment) allows us to sample structure and composition independently (Fig. 6).

Fig. 6 Hydrogen dynamics in hydrogenated Sc-doped BaTiO₃ (BaTi_1−xSc_xO₃H_x) for various concentrations of dopants, compared with experimental results measured at three different neutron spectrometers: IN1 Lagrange, TOSCA, and MAPS. Experimental data are from ref. 58. (a) IN1 Lagrange is a wide-range spectrometer that probes the dynamics in the region from 0 to 500 meV. (b) TOSCA and (c) MAPS on the other hand can be used to study the region around 100 meV and 500 meV, respectively, utilizing the higher energy resolution they offer. Simulated spectra using our workflow (labeled MD) corrected for quantum statistics and the kinematic constraint of the specific instrument are shown for all concentrations of dopants considered experimentally, with the simulated structure both in the cubic and hexagonal structure. Experimentally, only one of these phases is stable for a given concentration of Sc, but the energy of the two phases are sufficiently close that both phases are stable on the timescales of the MD simulations. The experimental data has thus been duplicated in the upper and lower rows of plots, where the compositions that are stable for each phase are indicated by the black lines. Additionally, simulated spectra based on harmonic phonons obtained via AbINS are included for comparison.

The simulated spectra using our workflow and the experimental spectra agree well in the full energy range 0 to 500 meV. The peak at 125 meV corresponds to O–H vibrations, and is best described by the hexagonal phase for low Sc concentrations and by the cubic phase for high Sc-concentrations (Fig. 6a). However, the overtone peak at 250 meV is underestimated in both simulated phases. This discrepancy could be due to the quantum correction factor in eqn (6), which is only valid for first-order scattering. This is further supported by the simulated spectra obtained using AbINS, which accurately capture the intensity of the 250 meV overtone peak. The latter method handles multi-phonon effects perturbatively and includes quantum effects but does not account for anharmonicity, which explains the sharper first-order features compared to the MD-based simulations.

The results from TOSCA highlight the 125 meV feature further (Fig. 6b). For low Sc concentrations the simulated spectrum using our workflow for the hexagonal structure agrees well with experiments, although the simulated spectrum is redshifted by approximately 25 meV. The simulated spectrum for the cubic structure agrees better with experiments as the Sc concentration is increased, which is in line with the hexagonal to cubic phase transition with increasing Sc concentration. We can thus clearly distinguish the spectra for the two phases of Sc-doped BaTiO₃, as the Sc-doping is varied.

Finally, the MAPS spectrometer probes the high-energy region between 300 meV to 600 meV. The feature in the experimental spectra at 450 meV corresponds to stretching of the O–H bond according to Perrichon et al., with the peak at 550 meV assigned as a combination mode of the O–H wag mode at 120 meV and the O–H stretch mode at 450 meV. The fundamental vibrational peak at 450 meV is captured by the simulations, although with a slight blueshift of 10 meV. However, the intensity for the combination mode is not reproduced by the simulation, neither using the MD-based workflow nor AbINS. In this case, we can further elucidate the nature of the combination modes at 550 meV using AbINS (Fig. 7). In these harmonic incoherent-approximation simulations, the intensity in that region is mainly composed of fourth-order scattering events and above. Such high-order phonons require a higher-order correction factor in order for the statistics to come out correctly using the MD-based workflow. However, applying a higher-order correction factor is not straightforward, as one would have to know a priori in which region of the spectrum to apply the correction, and the order of the higher-order scattering process.

Fig. 7 Contribution to the total dynamic structure factor from different scattering orders for hydrogenated Sc-doped BaTiO₃ with 16% Sc in the cubic phase, simulated for the MAPS spectrometer using AbINS. (a) Fundamental modes with higher orders corresponding to overtones and combination modes in panels (b–d). The color scale is the same for all subpanels, which leads to clipping of the intensity in panel (a). Specifically note that the feature at 550 meV originates from fourth-order scattering processes and above.

3.4 Discussion

The demonstrated workflows enable predictions of neutron scattering experiments, here in the form of INS spectra, from first principles. Starting from an atomic representation of a material, we develop MLIPs using the NEP framework for performing accurate MD simulations for systems comprising at least tens of thousands of atoms over at least a few nanoseconds, from which the dynamic structure factor can be computed using the dynasor package. We have focused on crystalline materials in the examples above, but the workflow is directly applicable to disordered systems as well, including liquids as well as amorphous or biological materials. The predictions are made instrument-specific by applying resolution functions and kinematic constraints, and are additionally corrected to account for the classical statistics inherent to the MD simulations.

The predictions show remarkable quantitative agreement with experiments. Almost all experimental features in the form of vibrational peaks are faithfully reproduced in the predicted spectra, including their relative intensities after applying correction factors for the quantum statistics. Remaining differences between the spectra, such as systematic blue- or redshifts of certain peaks in the case of crystalline benzene, can be most likely attributed to the difficulties in capturing the weak intermolecular interactions with the MLIP, as well as the DFT functional used for training the MLIP. However, the intensity for some higher-order phonon processes, such as in Sc-doped BaTiO₃, are not reproduced faithfully compared to experiments at the present level. In principle these discrepancies can, however, be corrected for by applying higher-order correction factors to recover the correct statistics for overtones and combination modes.

Such a correction would not represent a general-purpose ab initio approach as it is only possible in regions where these features can clearly be identified and separated, such that a single correction can be applied. For an unknown system, identifying such modes in a MD simulation is problematic, although some information can be gained from harmonic calculations, e.g., using the AbINS algorithm in Mantid. However, these corrections only affect the relative intensity of these peaks, while their positions are directly obtained from the MD simulations. Applying higher-order quantum corrections is thus not strictly necessary in order to give a reasonable prediction of a neutron scattering experiment. In general, we suggest applying the lowest order correction for the whole spectrum to be sufficient for the purpose of guiding neutron scattering experiments.

The workflow as presented in this work relies on MLIPs in order to run accurate and efficient MD simulations. Classical force fields could be used but the results might be of limited accuracy, especially in systems involving both bonded and non-bonded interactions. However, training a MLIP or selecting an appropriate force field for a system of interest constitutes a bottleneck in the workflow, requiring domain knowledge and effort. Foundation models trained on large parts of the periodic table, such as MACE-MP-0 or CHGNet among others,^59–63 offer an appealing alternative to creating bespoke MLIPs or using existing force fields. These foundation models can either be used out-of-the-box, or fine tuned with a small number of structures from DFT to yield an accurate model with comparatively low effort. It should, however, be noted that these models are often computationally much more demanding than either NEP models or classical force fields. Recently, a foundation model based on the NEP framework, NEP89, has been released, which combines the computational efficiency of NEP with a competitive transferability and accuracy compared to previously published foundation models out-of-the-box.⁶⁴ Additionally, NEP89 can efficiently be finetuned for a specific application, using only a limited amount of training data from DFT. Using finetuned models based on NEP89 would thus directly alleviate the main bottleneck of our workflow, developing the MLIP, whilst retaining the high computational efficiency. Foundation models further improves the ease of use of our workflow, and enables researchers to quickly predict scattering signatures for systems without existing MLIPs, for the purpose of guiding experiments or materials discovery.

4 Conclusion

In this study, we have presented a workflow that enables predictions of neutron scattering experiments from first principles, by combining DFT calculations, MLIPs, density autocorrelation functions from MD simulations as well as instrument resolution functions and kinematic constraints. We envision this workflow to be of great use in the context of materials discovery, offering an avenue for generating simulated experimental signatures for novel materials that can be directly compared to neutron, X-ray, and electron scattering experiments. ML in the form of MLIPs plays a central role, as the latter enable the accurate MD simulations and the extensive sampling that are the foundation of the present workflow. By integrating these components into a cohesive pipeline, our approach bridges the gap between theory and experiment, facilitating a more efficient feedback loop in the design and characterization of new materials. Ultimately, this workflow stands to accelerate materials analysis and discovery processes by providing high-fidelity, simulation-based insights that are directly aligned with experimental observables.

Data availability

The MLIP models, training data and databases of DFT calculations for benzene and hydrogenated Sc-doped BaTiO₃, as well as the reduced INS data for crystalline benzene are available on zenodo at https://doi.org/10.5281/zenodo.15283532. The development of the GPUMD package is hosted at https://github.com/brucefan1983/GPUMD and its documentation can be found at https://gpumd.org. The calorine package is hosted at https://gitlab.com/materials-modeling/calorine, its documentation is provided at https://calorine.materialsmodeling.org, and releases are available at https://doi.org/10.5281/zenodo.7919206. The dynasor package is hosted at https://gitlab.com/materials-modeling/dynasor, its documentation is provided at https://dynasor.materialsmodeling.org, and releases are available at https://doi.org/10.5281/zenodo.10012241.

Author contributions

Eric Lindgren: conceptualization, methodology, software, formal analysis, investigation, resources, data curation, writing – original draft, writing – review & editing, visualization, project administration Adam J. Jackson: conceptualization, software, investigation, data curation, writing – original draft, writing – review & editing Erik Fransson: software, writing – original draft, writing – review & editing Esmée Berger: software, writing – original draft, writing – review & editing Goran Škoro: resources, investigation, writing – review & editing Svemir Rudić: resources, investigation, writing – review & editing Rastislav Turanyi: software Sanghamitra Mukhopadhyay: supervision, funding acquisition Paul Erhart: conceptualization, methodology, investigation, resources, data curation, writing – original draft, writing – review & editing, supervision, funding acquisition.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

We are grateful to Adrien Perrichon for helpful discussions and for providing the original data and analysis scripts for the neutron measurements on hydrogenated Sc-doped BaTiO₃ (BaTi_1−xSc_xO₃H_x). We gratefully acknowledge funding from the Swedish Foundation for Strategic Research via the SwedNESS graduate school (GSn15-0008) and the Swedish Research Council (No. 2020-04935 and 2021-05072) as well as computational resources provided by the National Academic Infrastructure for Supercomputing in Sweden at NSC, PDC, and C3SE partially funded by the Swedish Research Council through grant agreement no. 2022-06725, as well as the Berzelius resource provided by the Knut and Alice Wallenberg Foundation at NSC. We are grateful to the Science and Technology Facilities Council (STFC), and to the ISIS Neutron and Muon Source in particular, for the provision of beamtime via TOSCA Xpress Access route (RB1990290). We acknowledge use of the computing resources STFC SCARF and SCD Cloud.

References

1. K. Shahzad, A. I. Mardare and A. W. Hassel, Sci. Technol. Adv. Mater.:Methods, 2024, 4, 2292486.
2. C. Li, L. Bao, Y. Ji, Z. Tian, M. Cui, Y. Shi, Z. Zhao and X. Wang, Coord. Chem. Rev., 2024, 514, 215888.
3. C. Chen, D. T. Nguyen, S. J. Lee, N. A. Baker, A. S. Karakoti, L. Lauw, C. Owen, K. T. Mueller, B. A. Bilodeau, V. Murugesan and M. Troyer, J. Am. Chem. Soc., 2024, 146, 20009–20018.
4. S. Fujii, J. Hyodo, K. Shitara, A. Kuwabara, S. Kasamatsu and Y. Yamazaki, Sci. Technol. Adv. Mater., 2024, 25, 2416383.
5. G. Hautier, C. C. Fischer, A. Jain, T. Mueller and G. Ceder, Chem. Mater., 2010, 22, 3762–3767.
6. G. Huang, F. Huang and W. Dong, Chem. Eng. J., 2024, 492, 152294.
7. K. Guo, Z. Yang, C.-H. Yu and M. J. Buehler, Mater. Horiz., 2021, 8, 1153–1172.
8. T. Lookman, P. V. Balachandran, D. Xue and R. Yuan, npj Comput. Mater., 2019, 5, 1–17.
9. S. Wu, Y. Kondo, M.-a. Kakimoto, B. Yang, H. Yamada, I. Kuwajima, G. Lambard, K. Hongo, Y. Xu, J. Shiomi, C. Schick, J. Morikawa and R. Yoshida, npj Comput. Mater., 2019, 5, 1–11.
10. M. H. Müser, S. V. Sukhomlinov and L. Pastewka, Adv. Phys.:X, 2023, 8, 2093129.
11. O. T. Unke, S. Chmiela, H. E. Sauceda, M. Gastegger, I. Poltavsky, K. T. Schütt, A. Tkatchenko and K.-R. Müller, Chem. Rev., 2021, 121, 10142–10186.
12. J. Behler, Chem. Rev., 2021, 121, 10037–10072.
13. Z. Chen, N. Andrejevic, N. C. Drucker, T. Nguyen, R. P. Xian, T. Smidt, Y. Wang, R. Ernstorfer, D. A. Tennant, M. Chan and M. Li, Chem. Phys. Rev., 2021, 2, 031301.
14. G. Ehlers, M. L. Crow, Y. Diawara, F. X. Gallmeier, X. Geng, G. E. Granroth, R. D. Gregory, F. F. Islam, R. O. Knudson, F. Li, M. S. Loyd and B. Vacaliuc, Instruments, 2022, 6, 22.
15. A. Borgschulte, J. Terreni, E. Billeter, L. Daemen, Y. Cheng, A. Pandey, Z. Łodziana, R. J. Hemley and A. J. Ramirez-Cuesta, Proc. Natl. Acad. Sci. U. S. A., 2020, 117, 4021–4026.
16. Y. Q. Cheng, A. I. Kolesnikov and A. J. Ramirez-Cuesta, J. Chem. Theory Comput., 2020, 16, 7702–7708.
17. Z. Fan, Y. Wang, P. Ying, K. Song, J. Wang, Y. Wang, Z. Zeng, K. Xu, E. Lindgren, J. M. Rahm, A. J. Gabourie, J. Liu, H. Dong, J. Wu, Y. Chen, Z. Zhong, J. Sun, P. Erhart, Y. Su and T. Ala-Nissila, J. Chem. Phys., 2022, 157, 114801.
18. E. Fransson, M. Slabanja, P. Erhart and G. Wahnström, Adv. Theory Simul., 2021, 4, 2000240.
19. E. Berger, E. Fransson, F. Eriksson, E. Lindgren, G. Wahnström, T. H. Rod and P. Erhart, Dynasor 2: from Simulation to Experiment through Correlation Functions, 2025, http://arxiv.org/abs/2503.21957.
20. T. M. Linker, A. Krishnamoorthy, L. L. Daemen, A. J. Ramirez-Cuesta, K. Nomura, A. Nakano, Y. Q. Cheng, W. R. Hicks, A. I. Kolesnikov and P. D. Vashishta, Nat. Commun., 2024, 15, 3911.
21. Z. Fan, Z. Zeng, C. Zhang, Y. Wang, K. Song, H. Dong, Y. Chen and T. Ala-Nissila, Phys. Rev. B: Condens. Matter Mater. Phys., 2021, 104, 104309.
22. Z. Fan, J. Phys.: Condens. Matter, 2022, 34, 125902.
23. E. Fransson, J. Wiktor and P. Erhart, J. Phys. Chem. C, 2023, 127, 13773–13781.
24. E. Lindgren, M. Rahm, E. Fransson, F. Eriksson, N. Österbacka, Z. Fan and P. Erhart, J. Open Source Softw., 2024, 9, 6264.
25. A. H. Larsen, J. J. Mortensen, J. Blomqvist, I. E. Castelli, R. Christensen, M. Dułak, J. Friis, M. N. Groves, B. Hammer, C. Hargus, E. D. Hermes, P. C. Jennings, P. B. Jensen, J. Kermode, J. R. Kitchin, E. L. Kolsbjerg, J. Kubal, K. Kaasbjerg, S. Lysgaard, J. B. Maronsson, T. Maxson, T. Olsen, L. Pastewka, A. Peterson, C. Rostgaard, J. Schiøtz, O. Schütt, M. Strange, K. S. Thygesen, T. Vegge, L. Vilhelmsen, M. Walter, Z. Zeng and K. W. Jacobsen, J. Phys.: Condens. Matter, 2017, 29, 273002.
26. F. Eriksson, E. Fransson and P. Erhart, Adv. Theory Simul., 2019, 2, 1800184.
27. M. Parrinello and A. Rahman, J. Chem. Phys., 1984, 80, 860–867.
28. P. Ying, W. Zhou, L. Svensson, E. Berger, E. Fransson, F. Eriksson, K. Xu, T. Liang, J. Xu, B. Song, S. Chen, P. Erhart and Z. Fan, J. Chem. Phys., 2025, 162, 064109.
29. M. Rossi, M. Ceriotti and D. E. Manolopoulos, J. Chem. Phys., 2014, 140, 234116.
30. D. S. Kim, O. Hellman, J. Herriman, H. L. Smith, J. Y. Y. Lin, N. Shulumba, J. L. Niedziela, C. W. Li, D. L. Abernathy and B. Fultz, Proc. Natl. Acad. Sci. U. S. A., 2018, 115, 1992–1997.
31. R. Fair, A. Jackson, D. Voneshen, D. Jochym, D. Le, K. Refson and T. Perring, J. Appl. Crystallogr., 2022, 55, 1689–1703.
32. R. Turanyi, A. Jackson and J. Wilkins, Pace-Neutrons/Resins: Python Library for Resolution Functions of Inelastic Neutron Scattering Instruments, 2025, https://github.com/pace-neutrons/resins, accessed 2025-04-10.
33. K. Dymkowski, S. F. Parker, F. Fernandez-Alonso and S. Mukhopadhyay, Phys. B, 2018, 551, 443–448.
34. O. Arnold, J. C. Bilheux, J. M. Borreguero, A. Buts, S. I. Campbell, L. Chapon, M. Doucet, N. Draper, R. F. Leal, M. A. Gigg, V. E. Lynch, A. Markvardsen, D. J. Mikkelson, R. L. Mikkelson, R. Miller, K. Palmen, P. Parker, G. Passos, T. G. Perring, P. F. Peterson, S. Ren, M. A. Reuter, A. T. Savici, J. W. Taylor, R. J. Taylor, R. Tolchenov, W. Zhou and J. Zikovsky, Nucl. Instrum. Methods Phys. Res., Sect. A, 2014, 764, 156–166.
35. M. Almakki, R. Applin, R. Backman, R. Baust, J. Borreguero, R. Boston, A. Bridger, J. Clarke, A. Diaz-Alvarez, R. Farooq, C. Finn, S. Foxley, D. Ganyushin, J. Haigh, T. Hampson, D. Ioannide, A. J. Jackson, W. P. Jayasundara Abeykoon Wickramasingha, D. Le, M. Lewis, Z. Morgan, M. Patrou, G. Pereira, P. F. Peterson, K. Qianli Ma, A. Savici, S. Schomann, C. Sears, K. TacTac, K. Travis, R. Waite, M. Walsh, R. Whitfield, J. Yusuf, C. Zhang and Y. Zhang, Mantid 6.12.0: Manipulation and Analysis Toolkit for Instrument Data, 2025,  DOI:10.5286/SOFTWARE/MANTID6.12.
36. P. Rosander, E. Fransson, N. Österbacka, P. Erhart and G. Wahnström, Phys. Rev. B: Condens. Matter Mater. Phys., 2025, 111, 064107.
37. Light Scattering, ed. Solids I. I., M. Cardona and G. Güntherodt, Springer, Berlin, Heidelberg, 1982, vol. 50.
38. P. E. Blöchl, Phys. Rev. B: Condens. Matter Mater. Phys., 1994, 50, 17953–17979.
39. G. Kresse and D. Joubert, Phys. Rev. B: Condens. Matter Mater. Phys., 1999, 59, 1758–1775.
40. G. Kresse and J. Hafner, Phys. Rev. B: Condens. Matter Mater. Phys., 1993, 47, 558–561.
41. G. Kresse and J. Furthmüller, Comput. Mater. Sci., 1996, 6, 15–50.
42. G. Kresse and J. Furthmüller, Phys. Rev. B: Condens. Matter Mater. Phys., 1996, 54, 11169–11186.
43. K. Berland and P. Hyldgaard, Phys. Rev. B: Condens. Matter Mater. Phys., 2014, 89, 035412.
44. J. W. Furness, A. D. Kaplan, J. Ning, J. P. Perdew and J. Sun, J. Phys. Chem. Lett., 2020, 11, 8208–8215.
45. S. F. Parker, F. Fernandez-Alonso, A. J. Ramirez-Cuesta, J. Tomkinson, S. Rudic, R. S. Pinna, G. Gorini and J. Fernández Castañon, J. Phys.: Conf. Ser., 2014, 554, 012003.
46. R. S. Pinna, S. Rudić, S. F. Parker, J. Armstrong, M. Zanetti, G. Škoro, S. P. Waller, D. Zacek, C. A. Smith, M. J. Capstick, D. J. McPhail, D. E. Pooley, G. D. Howells, G. Gorini and F. Fernandez-Alonso, Nucl. Instrum. Methods Phys. Res., Sect. A, 2018, 896, 68–74.
47. C. R. Harris, K. J. Millman, S. J. van der Walt, R. Gommers, P. Virtanen, D. Cournapeau, E. Wieser, J. Taylor, S. Berg, N. J. Smith, R. Kern, M. Picus, S. Hoyer, M. H. van Kerkwijk, M. Brett, A. Haldane, J. F. del Río, M. Wiebe, P. Peterson, P. Gérard-Marchant, K. Sheppard, T. Reddy, W. Weckesser, H. Abbasi, C. Gohlke and T. E. Oliphant, Nature, 2020, 585, 357–362.
48. W. McKinney, Proceedings of the 9th Python in Science Conference, 2010, pp. 56–61.
49. The Pandas Development Team, Pandas-Dev/Pandas: Pandas, Zenodo, 2024, https://zenodo.org/records/13819579, accessed 2025-07-02.
50. P. Virtanen, R. Gommers, T. E. Oliphant, M. Haberland, T. Reddy, D. Cournapeau, E. Burovski, P. Peterson, W. Weckesser, J. Bright, S. J. van der Walt, M. Brett, J. Wilson, K. J. Millman, N. Mayorov, A. R. J. Nelson, E. Jones, R. Kern, E. Larson, C. J. Carey, İ. Polat, Y. Feng, E. W. Moore, J. VanderPlas, D. Laxalde, J. Perktold, R. Cimrman, I. Henriksen, E. A. Quintero, C. R. Harris, A. M. Archibald, A. H. Ribeiro, F. Pedregosa and P. van Mulbregt, Nat. Methods, 2020, 17, 261–272.
51. J. D. Hunter, Comput. Sci. Eng., 2007, 9, 90–95.
52. M. Ulmestrand, Perfect-Cmaps, 2025, https://github.com/m-ulmestrand/perfect-cmaps, accessed 2025-04-10.
53. A. Stukowski, Modell. Simul. Mater. Sci. Eng., 2009, 18, 015012.
54. G. L. Squires, Introduction to the Theory of Thermal Neutron Scattering, Cambridge University Press, Cambridge, 3rd edn, 2012.
55. O. Hellman and I. A. Abrikosov, Phys. Rev. B: Condens. Matter Mater. Phys., 2013, 88, 144301.
56. T. Tadano, Y. Gohda and S. Tsuneyuki, J. Phys.: Condens. Matter, 2014, 26, 225402.
57. L. Monacelli, R. Bianco, M. Cherubini, M. Calandra, I. Errea and F. Mauri, J. Phys.: Condens. Matter, 2021, 33, 363001.
58. A. Perrichon, N. Torino, E. Jedvik Granhed, Y.-C. Lin, S. F. Parker, M. Jiménez-Ruiz, M. Karlsson and P. F. Henry, J. Phys. Chem. C, 2020, 124, 8643–8651.
59. C. Chen and S. P. Ong, Nat. Comput. Sci., 2022, 2, 718–728.
60. B. Deng, P. Zhong, K. Jun, J. Riebesell, K. Han, C. J. Bartel and G. Ceder, Nat. Mach. Intell., 2023, 5, 1031–1041.
61. A. Merchant, S. Batzner, S. S. Schoenholz, M. Aykol, G. Cheon and E. D. Cubuk, Nature, 2023, 624, 80–85.
62. F. Xie, T. Lu, S. Meng and M. Liu, Sci. Bull., 2024, 69, 3525–3532.
63. I. Batatia, P. Benner, Y. Chiang, A. M. Elena, D. P. Kovács, J. Riebesell, X. R. Advincula, M. Asta, M. Avaylon, W. J. Baldwin, F. Berger, N. Bernstein, A. Bhowmik, S. M. Blau, V. Cărare, J. P. Darby, S. De, F. D. Pia, V. L. Deringer, R. Elijošius, Z. El-Machachi, F. Falcioni, E. Fako, A. C. Ferrari, A. Genreith-Schriever, J. George, R. E. A. Goodall, C. P. Grey, P. Grigorev, S. Han, W. Handley, H. H. Heenen, K. Hermansson, C. Holm, J. Jaafar, S. Hofmann, K. S. Jakob, H. Jung, V. Kapil, A. D. Kaplan, N. Karimitari, J. R. Kermode, N. Kroupa, J. Kullgren, M. C. Kuner, D. Kuryla, G. Liepuoniute, J. T. Margraf, I.-B. Magdău, A. Michaelides, J. H. Moore, A. A. Naik, S. P. Niblett, S. W. Norwood, N. O'Neill, C. Ortner, K. A. Persson, K. Reuter, A. S. Rosen, L. L. Schaaf, C. Schran, B. X. Shi, E. Sivonxay, T. K. Stenczel, V. Svahn, C. Sutton, T. D. Swinburne, J. Tilly, C. van der Oord, E. Varga-Umbrich, T. Vegge, M. Vondrák, Y. Wang, W. C. Witt, F. Zills and G. Csányi, A Foundation Model for Atomistic Materials Chemistry, 2024, http://arxiv.org/abs/2401.00096.
64. T. Liang, K. Xu, E. Lindgren, Z. Chen, R. Zhao, J. Liu, E. Berger, B. Tang, B. Zhang, Y. Wang, K. Song, P. Ying, N. Xu, H. Dong, S. Chen, P. Erhart, Z. Fan, T. Ala-Nissila and J. Xu, NEP89: Universal Neuroevolution Potential for Inorganic and Organic Materials across 89 Elements, 2025, https://arxiv.org/abs/2504.21286.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: https://doi.org/10.1039/d5ta03325j

---