Direct prediction of inelastic neutron scattering spectra from the crystal structure

Inelastic neutron scattering (INS) is a powerful technique to study vibrational dynamics of materials with several unique advantages. However, analysis and interpretation of INS spectra often require advanced modeling that needs specialized computing resources and relevant expertise. This difficulty is compounded by the limited experimental resources available to perform INS measurements. In this work, we develop a machine-learning based predictive framework which is capable of directly predicting both one-dimensional INS spectra and two-dimensional INS spectra with additional momentum resolution. By integrating symmetry-aware neural networks with autoencoders, and using a large scale synthetic INS database, high-dimensional spectral data are compressed into a latent-space representation, and a high-quality spectra prediction is achieved by using only atomic coordinates as input. Our work offers an efficient approach to predict complex multi-dimensional neutron spectra directly from simple input; it allows for improved efficiency in using the limited INS measurement resources, and sheds light on building structure-property relationships in a variety of on-the-fly experimental data analysis scenarios.


Introduction
The collective vibration of atoms in a material is a fundamental 'vital sign' that dictates many of its properties. Experimentally, it is often probed by spectroscopic methods in which a beam of probe particles (photons, neutrons, or electrons) is directed into the sample. The atomic dynamics is then measured by the energy and momentum transfer when the particles are scattered. Inelastic neutron scattering (INS) is a technique with a few unique advantages. For example, it is sensitive to light elements such as H, O, Li, etc, which are often challenging to see with other methods such as x-ray scattering. The high penetration of neutron beams allows the use of complex sample environments and allows neutrons to probe deep into bulk material. INS is also ideal for studying phonons because the energy and momentum of thermal neutrons are well matched to those of phonons, allowing high-resolution measurement of the whole dynamic structure factor S(Q,E), where Q is the momentum transfer and E is the energy transfer. The dynamic structure factor contains rich microscopic information on correlations and interactions. For instance, neutrons interact directly with nuclei, making it possible to calculate INS spectra from lattice or molecular dynamics simulations rigorously. In cases where researchers are probing for magnetic excitations using neutron scattering, phonons and molecular vibrations may be considered spurious signals overlapping with the desired magnetic signal. In this case, having simple planning tools to calculate the phonon spectrum would allow scientists to measure the magnetic excitations more efficiently and effectively. Unlike Raman scattering or Infrared spectroscopy, INS has no selection rules for lattice excitations or molecular vibrations, meaning that all such modes can, in principle, contribute to INS intensities. This makes INS powerful, but also the spectra more challenging to decipher and dependent on comprehensive analysis and modeling. In fact, modeling INS from first-principles calculations (e.g. using density functional theory, DFT) requires significant computing resources and expertise. It is one of the most significant barriers limiting the scientific productivity at INS experimental facilities.
With the rapid advance in high-performance computing and machine learning, a data-driven solution to overcome these barriers becomes increasingly possible and intriguing [1][2][3][4]. Recently, Chen et al [5] demonstrated that one could directly predict the phonon density of states (PDOS) from a material's crystal structure without explicitly calculating phonons via DFT. Fundamentally, the vibrational behavior of a material is indeed determined by its structure (i.e. the atomic coordinates and the elemental species), but the connection is rather complicated and implicit. Thus, a simple prediction seemed almost impossible to establish. The key to success by Chen et al is using an Euclidean neural network (as implemented in the e3nn library) [6,7], a highly efficient way to represent the crystal structure with intrinsic equivariance to three-dimensional (3D) rotation, translation, and inversion. They showed that when the crystal structure represented by e3nn is trained with the corresponding pre-calculated PDOS (represented by a scaler of 51 numbers in the range of [0,1]), the resulting neural network can make reasonably accurate predictions of PDOS. This achievement opens the door for directly predicting other related properties, such as INS spectra, from the crystal structure. There are, however, two additional barriers before this can be applied to INS. First, we need an INS spectra database for training. For a rigorous calculation, one needs the PDOS and complete phonon information, including eigenfrequencies and polarization vectors, directly coupled with performing high-throughput calculations of INS spectra for different experimental scenarios. Second, although the prediction of the one-dimensional (1D) S(E) spectra would be helpful, the two-dimensional (2D) S(Q,E) map contains much more information and is what most INS users are interested in when measuring powder samples (the 4D S(Q,E) for single crystals will be deferred for examination in future studies). A 2D map, however, contains significantly more variables (e.g. for a 300 × 300 map in wave-vector, Q, and energy transfer, E, space, that is 90 000 variables to predict), and thus direct prediction is unlikely to work. In this paper, we demonstrate how the above issues can be solved by implementing a workflow for high-throughput INS spectra calculation, as well as an autoencoder to compress and represent the 2D S(Q,E) in latent space. The hybrid network proposed in this work, combining an autoencoder for data compression and a symmetry-aware neural network for feature prediction in the latent space, can be used for not only INS prediction, but more generally connecting the crystal structure directly with other multi-dimensional properties of the material.

Results
To develop the machine learning model for prediction, a training dataset of INS spectra covering all common elements and typical crystal structures needs to be produced. We have adopted the phonon database by Togo [8][9][10] as it provides the calculated force constants needed to obtain the eigenfrequencies and polarization vectors (via Phonopy) [11]. The OCLIMAX [12] program, which allows streamlined, rapid, and accurate calculations of INS spectra, was used to obtain the 1D S(E) and 2D S(Q,E) from the Phonopy output. Following this workflow, we created an INS database containing 10 023 inorganic crystals. With the DFT-derived force constants available (i.e. the most time-consuming DFT calculations already performed by Togo), the high-throughput workflow allows the construction of such an INS database in hours (in this case, on a computer cluster equipped with 50 compute nodes or 1600 CPU cores). This makes it possible to customize the database training for a specific application (e.g. a specific instrument, sample environment, resolution, Q/E range, or material class). In this study, the 1D spectra were simulated for the VISION spectrometer [13], and the 2D S(Q,E) were simulated in the energy transfer (E) range of 0-150 meV and momentum transfer (Q) range of 0-15 Å −1 . All simulations were done at 0 K (phonon population terms at finite temperatures can be applied later following the Bose-Einstein distribution).
A straightforward extension of the original work by Chen et al [5] is to train a Euclidean neural network to predict 1D S(E) spectra, which is, in principle, the same as the prediction of the PDOS. Indeed, we have achieved similar performance, as shown in figure 1, using essentially the same training parameters. The accuracy of the prediction can be seen in the distribution of the mean squared errors (MSE), as well as the example spectra with MSE ranging from the smallest 25% (top row) to the largest 25% (bottom row). Overall, even in the worst-case scenarios, the general features of the spectra, such as the energy range of fundamental excitations, distribution of phonon bands and band gaps, position and relative intensity of the major peaks, are well predicted. Some discrepancies such as the sharpness of the peaks and some fine structures are observed. These observations are important for the users to decide whether the prediction accuracy is sufficient for their specific applications. In addition to what has been presented in figure 1, different INS features can also be trained as the prediction output. For example, one can train one neural network to predict the fundamental excitations and another neural network for all the higher order combinations and overtones. One can also obtain partial contributions to the total spectra for a specific element or atom. This will be very useful in guiding the decomposition of the experimentally measured total spectra.
While the prediction of 1D S(E) is helpful, to fully take advantage of neutron scattering, which can measure the 2D S(Q,E) for powder samples with high resolution and accuracy, we would like to be able to perform the same kind of prediction, from structure directly to 2D S(Q,E). However, while a 1D S(E) in figure 1 only has 50 variables to predict, a 2D S(Q,E), assuming in a typical case of 300 × 300 bins, needs 90 000 variables to reconstruct. This is beyond what can be achieved with the current neural network. Fortunately, many 2D S(Q,E) maps show certain similarities in the overall pattern. For example, bands and stripes indicating energy levels and phonon branches, common intensity profiles with increasing Q, etc. This means the 90 000 variables are far from uncorrelated, and one should be able to significantly compress the 2D S(Q,E) to a much smaller dimension without losing the essential features. Recently, Samarakoon et al [4] demonstrated that an autoencoder could efficiently compress single crystal diffuse neutron scattering data, a 3D dataset, into a small latent space. Enlightened by this work, here we train an autoencoder to compress the synthetic 300 × 300 S(Q,E) into a 50-dimensional latent space vector, which will later be associated with the crystal structure and another neural network trained for predictions. In fact, in two follow-up publications [14,15], Samarakoon et al also demonstrated the usefulness of this approach in exploring the parameter phase space in magnetic materials. The architecture of the autoencoder is illustrated in figure 2 (details in section 4), together with an example of the reconstruction. The successful reconstruction (including some fine details and coherent scattering features) indicates that we can 'decompose' any 2D S(Q,E) onto 50 'principal components' , which are the unit vectors in the 50-dimensional latent space. With the compressed representations of the 2D S(Q,E), the next step is to connect the crystal structures with the latent space vectors. The trained encoder and decoder will be used separately in this step. Specifically, we run all the 10 023 S(Q,E) through the encoder and convert them into latent space representations. Then we perform another training of a Euclidean neural network to predict the latent space representations from the crystal structure (i.e. a latent space predictor). Here the 50-dimensional latent space vectors are treated the same as the 1D S(E). An important detail in this training is that a relatively large cutoff radius of 6 Å was used to include more neighbors into consideration. This is intended to capture the coherent effects present in some of the 2D S(Q,E), especially for crystals with small unit cells (simple and high symmetry structure) and coherent scattering elements. In the application cycle, we take a crystal structure as input with the autoencoder and the latent space predictor. The latent space predictor will predict the latent space vector. The predicted latent space vector can then be passed through the decoder to reconstruct the predicted 2D S(Q,E). The workflow is illustrated in figure 4, and the predicted latent space vectors for some crystals that are not part of the training dataset are shown in figure 5, in the same style as figure 1, with the latent space vectors presented as 1D spectra.   . Workflow for 2D S(Q,E) prediction. In training cycle 1, an autoencoder is trained to compress the 2D S(Q,E) into a 50-dimensiontal latent space. In training cycle 2, a latent space predictor is obtained by training the latent space vectors with the corresponding crystal structure. In the application cycle, the latent space vector is predicted from the crystal structure, and then the S(Q,E) is reconstructed from the decoder. The right column shows the predicted S(Q,E) from the crystal structure (i.e. reconstructed from the 'predicted' latent space vectors). Therefore, a comparison between the left and middle columns indicates how faithfully the autoencoder compresses and reconstructs the S(Q,E); and a comparison between the middle  autoencoder and the latent space predictor. For most systems, the predictions are fairly accurate, reproducing the most distinctive features and sometimes even the details such as the coherent scattering patterns. It is also observed that the worst predictions are often on systems containing hydrogen (thus having multiple strong and sharp peaks in the spectra) or very strong and unique coherent effects (with unique phonon dispersion patterns that are difficult to learn if such patterns are rare/non-existent in the training data). One possible way to alleviate this issue is to train multiple neural networks for different classes of materials based upon composition and structure, i.e. we sort the data into different categories according to the distinctive features they may have and train a dedicated neural network for each category. For example, hydrogen containing materials can be treated separately because they usually have their distinctive features in S(Q,E).
Thus far, we have been working with synthetic INS data to train and test the model. In the following example we will show a direct comparison with experimental data that demonstrates how the machine learning model is used in a real application. A research team is studying the magnetic properties of Ba 3 ZnRu 2 O 9 , which is not found in our training database. The powder INS data was collected at SEQUOIA [16] at 5 K with an incident energy of 160 meV. Figure 7 compares the S(Q,E) from experiment, machine learning model prediction, and full DFT simulation. While the full DFT simulation reveals more details associated with coherent scattering, the rapid prediction from the machine learning model is similar to the experimentally measured spectrum. We emphasize that the full DFT calculation took days on a computer cluster using hundreds of CPU cores while the machine learning prediction was done on a single graphics processing unit (GPU) in a fraction of a second. Moreover, for materials with more complicated crystal structures, full DFT calculations may not be feasible. The comparison in figure 7 immediately shows anomalous scattering at low Q (below 2 Å −1 ) in the experimental spectra, which is assigned to magnetic correlations. The level of accuracy from the rapid prediction is sufficient for many application scenarios such as experiment planning, ad-hoc decision making, and quick analysis where full DFT simulation is either impractical or unnecessary. One thing can be learned from the rapid prediction in this case is that because of overlapping contributions of magnetic scattering and phonon scattering in the region of 40 meV energy transfer, one may need to adjust the measurement to measure lower ranges of wave-vector transfer to fully quantify the magnetic scattering in this portion of the spectrum.
To further demonstrate that our approach can handle much more complex systems beyond the crystals in the training dataset (or in general beyond crystals that can be calculated by DFT), in figure 8 we show the predicted INS spectrum for a TiO 2 (8 mol%)-SiO 2 glass. Since the structure is amorphous and the distribution of Ti is also random, a large structure model must be used (in this case it contains 5184 atoms with 3456 O, 1590 Si, and 138 Ti, as depicted in figure 8). Phonon calculation with DFT would be impractical for such a complex system. However, our machine learning model can still make rapid predictions in the same way as what it does for other simpler crystals, and the accuracy of the prediction is decent when it is compared with experimental data reported in literature [17]. This example clearly demonstrates how our method can enable and greatly accelerate INS data analysis, especially as additional complex materials are being measured at the neutron scattering facilities.

Discussion
The availability of the synthetic INS database and the capability to perform on-demand, effortless prediction of S(Q,E) are expected to have significant value to the neutron scattering and other spectroscopic research communities. It will advance scientific research using neutron scattering by enabling better designed and more efficient experiments given the extremely limited experimental resources. The ability to rapidly understand entire regions of the phonon spectrum will significantly enhance the cycle of structural design and optimization of materials' properties. The developed framework in this work will be deployed to the data analysis servers at neutron scattering facilities (e.g. analysis.sns.gov at the Spallation Neutron Source) to support the facility users in the stages of experimental planning, execution, and data analysis. The users will be able to upload their crystal structure (e.g. in the standardized format of crystallographic information file (CIF) files) to the server and they will receive the predicted INS spectra immediately. A graphic user interface will be developed to further facilitate visualization for a better user experience. The machine learning model on the server will be retrained and updated periodically as new training data becomes available. Users may also choose to download the trained model to perform predictions on their own computers. To illustrate how tools developed based on this work can help, we assume a user would like to perform an experiment at the SEQUOIA direct geometry chopper spectrometer instrument. Before the measurements, important instrument settings such as the incident neutron energy and chopper frequency need to be decided. Knowing approximately what the spectra would look like will allow users to obtain the optimal setup efficiently well before arriving at the instrument. This software could even be used in beamtime proposals in order to illustrate ranges of energy transfer that will be explored if beamtime is granted to the applicant. During the experiment, the collected data will be examined to determine whether the experiment is succesfully measuring the expected signals from the sample, how strong is the contribution from the instrument background and sample container, and if there are any unusual observations that require more focused measurements. All these can be greatly facilitated with the predicted spectra at hand, and the capability to perform on-demand predictions with additional/modified structure models. Rapid prediction of S(Q,E) will also allow users to disentangle phonons and magnons in the measured total spectra in real-time. This has been a long-standing challenge of signal separation, given the dual neutron-matter interactions (nuclear interactions and spin interactions). It is worth mentioning that after the training stage (which the users do not have to do as the trained models will be provided on the data analysis servers), all these predictions can be achieved on a personal computer almost instantaneously by a generic user even with no prior modeling expertise. This is in stark contrast to the conventional approach of running DFT simulations on a computer cluster or supercomputer, which could take significant computational cost and expertise to perform. The prediction accuracy is close to ab initio results within acceptable error, and will be further improved as we include more data in the training database.
The workflow and methodology are not limited to neutron spectra prediction but have a few immediate and broader implications. First, calculating the phonon contribution of the specific heat is a natural next step in using the predicted PDOS. This enables the selection of materials with low specific-heat for thermal switch and high specific-heat for thermal storage materials. Additionally, several other straightforward extensions such as inelastic x-ray scattering, Raman scattering, infrared spectroscopy, and electron energy loss spectroscopy can be performed in a similar fashion, aiming to connect crystal structure with multi-dimensional experimental data or properties. Third, since vibrational spectroscopy like INS is very sensitive to local vibrational modes associated with defects, the work, which allows fast prediction and analysis of the spectra for a large variety of defective structure models, also offers an alternative avenue to perform non-contact, non-destructive defect characterization and quantification, particularly for point defects that can directly create new local vibrational nodes.

Generation of INS database
The phonon database at Kyoto university, created by Atsushi Togo, contains DFT calculated force constants for over 10 000 inorganic crystals. The 2018-04-17 version [8] database was used to calculate the INS spectra in this work. Phonopy [11] was used to calculate the eigenfrequencies and polarization vectors on a Γ-centered mesh with spacing smaller than 0.1 Å −1 . The T = 0 K INS spectra were then calculated with OCLIMAX [12]. The 1D S(E) spectra were calculated for the VISION spectrometer using its instrument parameters [13]. The 2D S(Q,E) spectra were calculated in the energy range of 0-150 meV (bin size 0.5 meV), Q range of 0-15 Å −1 (bin size 0.05 Å −1 ), using an energy resolution of 1.5 meV. The calculations were performed on a computer cluster with 1600 CPU cores. It took about a day of computing time in total (this implies that regenerating this database for a different instrument, sample environment, E/Q range, or temperature will not be an issue).

Euclidean neural network
A Euclidean neural network very similar to that used by Chen et al [5] was trained to perform the prediction of 1D S(E) spectra and the latent space vector, respectively. In the former case, the 1D spectra calculated from OCLIMAX were processed to have an energy range of 10-1200 cm −1 (1 cm −1 = 0.124 meV), with a bin size of 10 cm −1 (i.e. 120 bins to represent the spectra). In the latter case, the 50-dimensional latent space vectors are the features to predict. A relatively large cutoff radius (6 Å) was used here to include more neighbors around each center atom. This is to capture potential coherent effects in the S(Q,E) which only becomes apparent when the unit cell size is small, the symmetry is high, and all elements in the crystal are strong coherent scatterers. The dataset was categorized as 80% for training, 10% for validation, and 10% for testing.

Autoencoder
A fully connected deep feed-forward neural network was used as the autoencoder. The seven hidden layers contain 90 000, 2500, 500, 50, 500, 2500, 90 000 nodes, respectively. Since the original S(Q,E) has strong intensities in a small fraction of the image and a large dark area in the rest, we applied this cubic function g(x) = 1 − (1 − x) 3 for brightening. This function was chosen because it is invertible, maps [0, 1] → [0, 1], and is monotonous. The latent space vectors are normalized so that the maximum number in the vector is unity. The dataset was categorized as 5% for validation, and 1% for testing, and 94% for training.

DFT calculation of Ba 3 ZnRu 2 O 9
Spin-polarized DFT calculations of Ba 3 ZnRu 2 O 9 (interlayer antiferromagnetic configuration) were performed using the Vienna Ab initio Simulation Package (VASP) [18]. The calculation used projector augmented wave [19,20] method to describe the effects of core electrons, and Perdew-Burke-Ernzerhof [21] implementation of the generalized gradient approximation for the exchange-correlation functional. Energy cutoff was 800 eV for the plane-wave basis of the valence electrons. The lattice parameters and atomic coordinates from literature [22] were used as the initial structure. The electronic structure was calculated on a 6 × 6 × 2 Γ-centered mesh. The total energy tolerance for electronic energy minimization was 10 −8 eV, and for structure optimization it was 10 −7 eV. The maximum interatomic force after relaxation was below 0.001 eV Å −1 . A (U-J) of 2.3 eV was applied to account for the localized 3d orbitals of Ru. A 3 × 3 × 1 supercell was created for phonon calculations. The interatomic force constants were calculated by density functional perturbation theory as implemented in VASP, and the vibrational eigen-frequencies and modes were then calculated using Phonopy [11]. The OCLIMAX software [12] was used to convert the DFT-calculated phonon results to the simulated INS spectra.

Data availability statement
The data that support the findings of this study are openly available at the following URL/DOI [23]: www.nature.com/articles/s41597-022-01926-x.