Understanding dynamic properties of materials using neutron spectroscopy and atomistic simulation

Recent developments in first-principles lattice dynamics and classical force field based molecular dynamics are revolutionising the field of neutron spectroscopy. Herein we present a short review of these methods, their critical role in the supporting of cutting-edge experiments, and how they are improved by matching experimental data. We begin with a brief overview of how lattice dynamics calculations can be compared to inelastic neutron scattering (INS) and molecular dynamics simulations to both INS and quasi-elastic neutron scattering (QENS). We then provide a series of exemplar applications where lattice dynamics and molecular dynamics have been used in conjunction with neutron spectroscopy to bring significant understanding to topical areas of materials science namely: (i) lattice dynamics and INS for the study of hybrid organic-inorganic perovskites (ii) lattice dynamics and INS for the study of flexible porous solids and (iii) molecular dynamics and QENS for probing molecular behaviour in zeolite catalysis. In all three cases, the understanding gained through the synergy of experiment and computation would have been significantly reduced using either in isolation. Finally, we consider the current state of the art, describing outstanding challenges and suggesting future directions in this exciting and fertile area of physical science.


Introduction
The review is intended to summarise and highlight how modern-day atomistic simulations can be used to extract more meaningful information from Neutron Spectroscopy (NS) experiments and how, vice versa, neutron spectroscopy can be used to inform and improve simulations. To contextualise the concepts discussed in the first section, we progress with a series of scientific case studies, showing how NS has been an integral part of the expansion of each field discussed. NS has enjoyed a vibrant, wide-ranging, and successful history ever since its mainstream inception. NS is particularly advantageous over similar vibrational spectroscopy techniques due to the lack of selection rules, the selectivity for light elements and the ability to probe the full Brillouin zone. It played a significant role in the 1991 Physics Nobel Prize winning theoretical concepts of Flory and de Gennes [1] and had a major contribution to the theoretical developments in the dynamics of liquids, such as the memory function formalism, or mode-coupling theory for the description of glasses [2]. It has even played a role in understanding the nature of magnetism in high-temperature superconductors [3]. It is also worth noting that the 1994 Nobel Prize in Physics was awarded to Shull and Brockhouse for the pioneering development of neutron scattering techniques and its contribution to studies of condensed matter. However, the groundbreaking successes of NS did not, rely on modern-day computational techniques, but instead relied on either simple analytical models or basic chemical intuition. As we have pushed the boundaries of physical chemistry further, the applicability of these traditional models has begun to fall short upon probing more complex systems. Cutting-edge scientific problems demand a sophisticated approach to the study of dynamics at the molecular scale.
Fortunately, computational modelling developments over the past 30 years have advanced tremendously both in the fields of ab-initio and classical force field modelling. In particular, the last few years have seen abinitio molecular dynamics (AIMD) systems increase significantly in size for many of the complex materials of technological interest (allowing for sizes and timescales of ∼200 atoms for over 100 ps) and have resulted in an outstanding agreement between experiment and theory for systems exhibiting complex anharmonic motion [4]. The ability to perform AIMD at relevant system sizes/timescales with the appropriate level of description is essential. However, more often than not, static density functional theory lattice dynamics (DFT-LD) calculations are more than adequate in providing quantitative agreement with inelastic neutron scattering (INS) for molecular vibrations. Where in the past, such predictions were costly, modern-day high-performance supercomputing resources can achieve such calculations on the scale of days. Also, INS has, in recent times, seen a paradigm shift to parametric studies in which a variety of stoichiometries are studied in a single experiment [5].
This paradigm shift has tremendous synergy with the newfound capability to run the corresponding DFT-LD calculations for each system, allowing one to study global chemical trends within the class of systems. While simulations can be used to interpret the experiments, the experimental data is also used to validate and refine computational models, resulting in a virtuous circle of improved models. With computation becoming a necessity, so too is there a necessity to develop robust software that is capable of facilitating thedirect comparison between neutron scattering and simulation data.

Summary of basic principles
Throughout this article, we will deal with a quantity that is central to neutron spectroscopy, the dynamic structure factor, S(Q, E). S(Q, E) is a measure of the spatial and temporal correlations between atoms, where it is most convenient mathematically to work in the reciprocal variables (Q; wave vector and E; energy). Importantly, these two quantities are directly related to the change in momentum and change in energy that occurs during the collision of a neutron with the nucleus of the material under study. The time-of-flight spectrometers, which we will discuss, therefore, require a way in which to measure the change in momentum and energy that occurs during this collision process. These spectrometers can be broadly broken down into two classes of instruments; direct and indirect. Both techniques rely on constraining the kinetic energy of the neutrons at some point along their trajectory. Direct instruments achieve this by selecting a narrow band of kinetic energies for the incoming neutrons, either via the use of a mechanical chopper, or monochromatic reflectance. With knowledge of the kinetic energy of the incoming neutron, they can, therefore, infer the change in energy upon interaction with the sample, from the length of the flight path (distance from the source to sample, plus the distance from the sample to the detector) and the time of arrival of each neutron. Indirect spectrometers, on the other hand, take the approach of constraining the post-scattering kinetic energy of the neutrons. In the example of TOSCA (ISIS Neutron and Muon Source, UK), this constraint is achieved via Bragg scattering from a graphitic analyzer. There is a distribution of incoming neutron energies, yet since only one energy is ultimately detected, each total flight time directly maps to a given change in energy after collision with the sample. Having considered how particle collisions allow us to obtain S(Q, E), we now consider the quantum wave nature of the scattering. This wave nature, in fact, allows us to break S(Q, E) down into a contribution resulting from interference of different scattering centres (coherent scattering), and a contribution lacking interference (incoherent scattering). These two contributions relate directly to the inter-particle correlations, and self-correlation, respectively. Thus the former contains information both on the structure and collective motions, whereas the latter primarily concerns the self-motion (diffusion, rotation, incoherent vibrations). The scattering cross sections for these two contributions differ, and in particular, the incoherent cross-section for the hydrogen atom is significantly larger than all other scattering cross-sections. This is, in fact, a key strength of neutron spectroscopy, as it gives extreme sensitivity to hydrogen motions, unlike many other techniques whose sensitivities are determined by electron density. One other consequence of the large incoherent scattering cross-section of hydrogen is that it becomes the basis for approximating S(Q, E) solely as the incoherent contribution (the incoherent approximation). Dispensing with these inter-particle correlations has allowed the development of many simple models of selfmotion, which may be directly compared to experiment. As we shall discuss, however, systems of contemporary interest are often too complex for such an approach. So we must turn to simulation methods, which allow us to explicitly model this complexity with both incoherent and coherent contributions.

Comparing DFT lattice dynamics to INS
A direct comparison of DFT-LD calculations to INS is not an 'out-of-the-box' feature of most DFT codes. Nuances involving the weighting of each mode with neutron cross-sections, combinations, and overtones of the fundamental modes (and in turn, their relative weighting), along with instrumental factors, must be applied to make a quantitative comparison between theory and experiment. The newly launched software package AbINS [6] provides these capabilities and aims to give comprehensive support to as many DFT codes as possible (currently including Gaussian, CASTEP, CRYSTAL, and VASP). The package is neatly included as an algorithm within the widely used neutron scattering analysis software package 'Mantid' [7], allowing a direct comparison of neutron and in-silico data.
A workflow chart is shown in figure 1, which details how AbINS processes DFT-LD output data. The raw input comes from the relative atomic displacements and frequencies from a given DFT-LD calculation. From these, a semi-empirical powder averaging model is applied, whereby each phonon is treated as an independent quantum harmonic oscillator. This process includes weighting by neutron scattering cross-sections and can be applied to cover any order/combination of the fundamental modes, with increasing orders having a progressively more minor contribution. At this point, a theoretical S(Q, E) is obtained, which to be directly compared with experiment, must still be convolved with the specific instrumental resolution function. The functional form for the resolution function of TOSCA [8] is currently included in AbINS; however, both indirect and direct neutron instruments may easily be included in future versions of the code. Processing with AbINS is an essential part of the data analysis as it rules out experimentally observed peaks which are simply overtones/combinations, and explains the absence of modes that involve lower scattering atomic species. Another recently published software package for INS simulation is OClimax [9], which was developed at Oak Ridge National Laboratory (ORNL). OClimax can interface with a similar variety of DFT codes and can perform INS simulations for various neutron spectrometers, including indirect geometry instruments such as VISION at the Spallation Neutron Source (SNS) and the previously mentioned TOSCA spectrometer.
Beyond this essential step of data processing, various considerations should be taken when comparing lattice dynamics (LD) with NS. LD calculations are computationally expensive, and in particular, going beyond a gamma point phonon calculation increases the computing time significantly. A full dispersion calculation attempts to thoroughly sample the Brillouin zone. Thus it is crucial to understand precisely what is being measured on a typical broadband INS instrument such as TOSCA or VISION. Figure 2 shows a two-dimensional (2D) representation of the momentum transfer in a polycrystalline system, where a circle of a given radius is shown and corresponds to a particular magnitude of momentum transfer. It is seen that as the ring sweeps through each new zone, the path sampled is distinct from the previous zones. Therefore, as the dynamics of each zone are a replica of the first zone, it can be a reasonable approximation to sample the full first zone thoroughly. One may also note that this is only the case provided the circle is of sizeable radius, which is a criterion typically fulfilled for low bandpass spectrometers such as TOSCA and VISION. LD, combined with NS, is a powerful tool not only because it allows one to validate the given DFT resultant structure but is also, in effect, a test on the local potential of the DFT model. The validity of which then allows you to make quantitative predictions of the vibrational entropy, which has been shown to be the stabilizing free energy contribution for many systems [10]. As a final note, many codes now incorporate anharmonic corrections, which may produce a more accurate match to experimental frequencies, but at an increased computational cost.

Comparing MD simulations to INS/QENS
LD calculations, although powerful, only answer questions for a subset of the systems studied via NS. For quasielastic neutron scattering studies of diffusion/rotation or even INS of solids with strong anharmonicity, the harmonic approximation of lattice dynamics is no longer appropriate. For such systems, one option is to perform fully-fledged MD simulations, either via semi-empirical force fields or ab-initio methods. Such methods are particularly appropriate for systems such as fluids, where anharmonicity can be large [11,12]. Here we provide a brief introduction to how the results of MD simulations can be used to analyse NS experiments.
As outlined in section 2, The incoherent dynamic structure factor (S inc (Q, E)), relates to the single-particle temporal and spatial correlations. This quantity is related to another function, the self-part of the intermediate scattering function (F s (Q, t)) by a Fourier transformation in the time domain: In a QENS experiment, the obtained S inc (Q, E) is typically fitted with an elastic component (static atoms), and an inelastic (atoms in motion) component. This elastic component takes the form of a delta function convoluted with a Gaussian function (which represents the instrumental resolution). The inelastic broadening of the spectrum, representing motion, is typically fitted with one or more Lorentzian functions. The half-width at half-maximum (HWHM) of this Lorentzian is related to the time scale of motion at a given scattering angle/ momentum transfer. For rotational motion, the width of the fitted Lorentzian displays no Q-dependence. However, for processes of diffusion, a strong Q dependence is seen, where the qualitative dependence of this width with Q can give insight into the nature of the diffusion. These relationships are discussed in detail in previous sources [13][14][15].
F s (Q, t) in equation (1) retains the time variable and can often be more intuitive for comparisons with simulation, especially for relaxation phenomena. Notably, one must derive the powder average of the function F s (Q, t) (as we are rarely working with single-crystal samples, and much more likely polycrystalline samples) from the MD simulations as shown in equation (2). Here, N is the number of hydrogen atoms, and r i is the position of the hydrogen atom i. If one is probing translational diffusion, the absolute atomic displacement is used, but if one is probing rotation, then we use the atomic displacement with respect to the centre of mass of the molecule, the angular brackets represent the ensemble average.
Once F s (Q, t) has been derived from the simulation, it may then be fitted with an exponential function (or more than one if necessary) representing motion in a specific time domain as shown in equation (3) [16]. The pre-exponential factor C denotes the weighting factor (how much the motion is contributing to the total F(Q, t)). The decay constant Γ represents the rate of the motion and can be considered the same as the Lorentzian HWHM used to fit the S(Q, E). At this point, it is worth noting that once each decay constant is known, one may then consider whether this motion is within the time window sampled by the QENS spectrometer, and discount or prioritise which motions to study further, or even plan future experiments as illustrated further in the case studies.
The constant B(Q) in equation (3) can be considered as an exponential where = ¥ t , and thus it represents the final atomic arrangement in Q space; equivalent to the elastic incoherent structure factor (EISF), another observable directly comparable to experiment.
Regardless of the choice of simulation approach, the resultant output is a set of atomic trajectories which, to be compared to NS results, require significant processing/analysis. Many of the emergent quantities probed (e.g., radial distribution functions, mean squared displacements, and rotational correlations) are not the phenomena that are directly probed by NS. Instead, as discussed previously, NS examines the dynamic structure factor, S(Q, E).
MDANSE [17] is a widely used software package for the analysis of MD and subsequent processing/ comparison to NS data. The S(Q, E) may be calculated on an arbitrary discrete grid, which may vary in geometry (spherical, cuboidal, etc.). A key consideration when running an MD simulation for comparison with an NS is that the length and timescales correspond to that covered by the spectrometer in question. Care must be taken that the dimensions of the simulation cell are large enough to capture the low Q limit of the spectrometer. From a time perspective, one must also take care to provide a simulation that is long enough in duration to reach frequencies lower than the energy resolution of the spectrometer in question. The high energy limit of the spectrometer will also dictate the rate at which one should print trajectory snapshots, as producing them more frequently than this will create temporal correlations outside of the energy window accessible to the instrument. It is also wise to provide an in-silico Q grid, which is similar to that of the spectrometer being used in terms of spacing and Q maximum. Once the discrete form of S(Q, E) has been calculated, it may be convolved with a Q dependent resolution function of the instrument. A choice of analytic functions is available, and for QENS applications, a Gaussian is usually an appropriate approximation. Both the cross and self-correlation components of S(Q, E) will be weighted by the coherent and incoherent scattering cross-sections of the atomic species in question. MDANSE provides a decomposition of both the coherent and incoherent dynamic structure factor in terms of all combinations of atomic species, allowing the user to identify the origin of specific features within the total S(Q, E). Here lies the most powerful aspect of MD in comparison to NS, often allowing the user to differentiate incoherent broadenings from loss of coherent scattering for a particular Q, as exemplified in the cited literature [18].
As discussed, generating an experimental/in-silico S(Q, E) is useful for understanding the atomic origin of individual features. But it also has been under-appreciated as a figure of merit for generating empirical classical force fields and improved DFT approximations. A recent study looking at the diffusion of benzene on a nickel catalyst [19], relied on the extraction of the incoherent contribution of S(Q, E) from the benzene, and the subsequent iterative fitting of the force field parameters governing the strength of absorption on the surface. The approach is particularly timely as the empirical potentials for liquids employed by MD are typically designed for interactions in the bulk liquid state, with specific solid surface-liquid/gas force fields being rare. However, potentials that have been fit to QENS spectra are particularly useful as the experiments are probing some very specific dynamical behaviours over a given timescale, enforcing the accuracy of their underlying force fields in this particular environment. With QENS, perfectly symmetrical Lorentzian functions centrally located in the energy range of the instrument are rare. For systems with a particularly weak incoherent fraction, or slow dynamics where only the tails of the Lorentzian protrude from the resolution function, it can be useful to perform a moment's analysis of the scattering function. In particular, the second moment of the experimental data may be compared directly to that of simulation (after processing in MDANSE), thus providing a validation of the simulation model. The simulation may then be used directly to understand the dynamics of the system [18]. In section 7 we provide several further case studies where QENS/MD have been used to understand the structure and dynamics of hydrocarbons in porous materials.

Lattice dynamics of hybrid organic-inorganic perovskites
Hybrid organic-inorganic perovskites (HOIPs) have shot to prominence in the last decade. The material is based on the perovskite crystal structure ABX 3 where octahedra of Bx 6 form a corner-sharing network, with A cations filling the cavities in the structure. In hybrid perovskites at least one of the A, B or X sites is occupied by a molecule, as opposed to a single atom in fully inorganic perovskites [20]. The significant interest in the HOIP systems has been fuelled by their extraordinary success as absorber layers in solar cells [21]. In a few short years, they have reached the point where they can display photovoltaic conversion efficiencies rivalling those of technologies over 50 years old, while necessitating only mild synthesis conditions [22]. This meteoric rise has occurred despite a lack of fundamental understanding of the material characteristics which drive such efficiencies. However, HOIPs suffer from a number of drawbacks, such as the toxicity of lead and inherent instabilities, thus an understanding of their mechanisms of photovoltaic action is necessary for the design of alternatives without these shortcomings. Combined neutron scattering experiments and atomistic simulations have been at the forefront of developing the required level of mechanistic understanding. Here we discuss examples of how they have been instrumental in revealing the effects of molecular motion and disorder, the unusual thermal transport properties, and the factors affecting phase changes in HOIPs.
As mentioned, the molecular orientation and dynamical information that is provided by combined neutron scattering and DFT studies are critical for the emerging understanding of the properties of HOIPs. In particular, the possibility of ordered electric fields from molecular dipoles and polarons (charge carriers dressed with local lattice distortions) are crucial for the operation of these materials in optoelectronic contexts [23]. In a landmark study, neutron powder diffraction provided a complete characterization of the phases of MAPI across a temperature range of 100-352 K; this work revealed increased disorder of the molecules with temperature [24]. Later, molecular cations were found to be positioned slightly off from the symmetric centre of the perovskite cavity, due to hydrogen bonding [25]. This structural representation is supported by DFT calculations [26]. Neutron diffraction revealed the phase transitions around room temperature, while complementary DFT calculations show that directional hydrogen bonding between the cation and the framework drives structure formation [27]. This study showed that octahedral tilting and lone pair distortions are related to the orientational ordering of the amine cation.
By combining QENS with AIMD and Monte Carlo simulations, a deep understanding of the motion of the organic cation was obtained [28]. QENS experiments revealed that the molecular cation in CH 3 NH 3 PbI 3 (MAPI) can re-orientate between faces corners or edges in the crystal-combined with AIMD, it was shown that orientation pointing towards the faces of the cavity is preferred. Monte Carlo simulations of this scenario could then be used to demonstrate the possibility of ferroelectric or anti-ferroelectric domains in these crystals with implications for charge transport. In the context of photovoltaics, the dynamics and characteristics of molecular reorientations become essential. If the molecular motions are rapid, they can provide dynamic screening for charge carriers in the system. On the other hand, slow movements would mean that molecules could serve to give a fixed field, driving charge separation. Recently total neutron scattering combined with DFT calculations and solid-state NMR has provided a comprehensive mapping of molecular reorientation times across a wide temperature range for both MA and FA [29].
Another important feature of the HOIPs is that these materials tend to have extremely low thermal conductivity. This was initially reported for CH 3 NH 3 PbI 3, and the reasons for the low conductivity were postulated to be related to molecular motions [30]. This observation of low thermal conductivity has lead to an interest in the application of HOIPs for heat to electricity (thermoelectric) technologies. Subsequent theoretical calculations probed the origins of low thermal conductivity further, but the results proved to be contentious. Whalley and co-workers report materials with anharmonic phonons giving short quasiparticle lifetimes and mean-free-paths, which result in a low thermal conductivity, based on lattice dynamics calculations with DFT [31]. Notably, these authors draw a distinction with other classical semiconductors (e.g., GaAs) where low energy acoustic modes are responsible for thermal conductivity; in the HOIPs, these modes have ultra-short lifetimes and therefore do not contribute to the thermal transport. The thermal transport is thereby reduced due to the scattering of the phonons of the inorganic cage by the organic cations. A similar situation was described by MD simulations using an empirical interatomic potential [32]. However, another set of MD simulations suggest that organic cations only serve to mediate the vibrations of the inorganic cage and do not act to scatter phonons [33]. The differences were resolved by a combination of INS, QENS, and MD simulations. Li and co-workers published a comprehensive characterisation of the orientations and dynamics of the cations and their interactions with the inorganic cage. INS provided information on the phonons of the inorganic cage, while complementary QENS experiments delivered the details on the organic molecule behaviour [34]. By coupling these observations with MD studies, they were able to show that molecular dipole order plays a dominant role in determining charge and thermal transport properties, favouring an interpretation where thermal transport is influenced by scattering from the organic molecule. This latter study powerfully demonstrates the power of combined neutron scattering/spectroscopy and atomistic modelling to unravel the intricate details of the structure and dynamics in HOIPs, which have a profound effect on their physical properties.
The interpretation and linking of physical properties to the atomistic structure rely on a solid understanding of the crystal structure of a material, as does the ability to design new systems. In the HOIPs, knowledge of the structure-property relationships and the phase evolution of different materials has been the subject of intense study. Early work from Kielsich and Cheetham applied the elegant and straightforward concept of the Goldschmidt tolerance factor to understand how molecular dimensions could be linked to crystal structure [35,36]. While this provides some key indicators for the structure, it has since become apparent that the situation in HOIPs is significantly more complicated. This stems from halide bonds being softer, and a different range of chemical bonding interactions which lead to a more nuanced relationship between structure and property [37,38]. Moreover, it has become apparent that vibrational entropy arising from soft inter-molecular bonds can play a crucial role in determining the structure [39,40].
Beyond establishing the structure at a given temperature, understanding how these structures evolve with temperature is critical for their application. Typically materials display anomalous dielectric behaviour in the region of a phase transition, and this can have profound effects on optical and electronic properties. The perovskites generally have an orthorhombic structure at low temperature, this gradually becomes more symmetric, passing through a tetragonal phase and finally into a high-temperature cubic phase.
Neutron scattering performed on high-resolution instruments reveals that a band of phonons related to octahedral tilting soften and become imaginary at the phase transition, suggesting a displacive mechanism for the phase transition [41]. However, a series of calorimetric measurements have shown a large increase in entropy associated with phase transitions, this is characteristic of an order/disorder phase transition, resulting from increased numbers of equivalent orientations of the organic cation [42]. A broader picture emerges where there is a coupling between the dynamics of the organic and inorganic components that are critical for determining structure and property. Guo and co-workers demonstrated that there is a coupling between the dynamics of organic and inorganic components; this results in an abrupt phase transition when hydrogen bonding is broken and is contrasted to the gradual phase transition in the all-inorganic CsPbBr 3 [43]. Their findings are in line with neutron scattering studies of the phase transition, which find a band of vibrations at −5 meV that shows soft modes at the phase transition [41]. Raman studies that show hydrogen bonding affects octahedral tilting [44]. The overall picture of the phase transition that emerges from these studies is of a mixed displacive/orderdisorder transition. All of these studies emphasize the importance of hydrogen bonding in these systems.
A recent combination of INS and DFT lattice dynamics has highlighted the crucial importance of hydrogen bonding for phase evolution [10]. In this work, Kieslich and co-workers correlated the vibrational spectra collected from INS with calculated phonons, using the ABINS package [6]. Crucially the study reports that hydrogen-bonding, configurational entropy, and vibrational entropy in a model HOIP system are all of the same order of magnitude. This observation could have profound effects on the future design of HOIPs and other systems, where entropy has emerged as a handle through which we can control phase evolution [45].
Although this has been a brief overview of the impact of neutron scattering and atomistic simulations in the field of HOIPs, we believe it highlights the profound synergistic nature of these approaches and how their combination is a powerful tool for understanding functional materials. The insights gained by combining experiment and theory have resolved outstanding questions, revealed new connections, and provided compelling design principles in the field of HOIPs. The successes reported here could be replicated across the full range of functional materials where structure and dynamics on an atomic level drive macroscopic properties.

Lattice dynamics of porous framework materials
Metal-organic frameworks (MOFs) are a diverse class of materials that have attracted a remarkable volume of scientific research. Initially, the primary motivation was on synthesizing a breadth of new materials due to the relatively simple construction from metals or metal clusters connected by organic linkers and in part due to their high levels of porosity and capability for selective guest binding, quantifying the resultant gas adsorption properties [46]. A large percentage of research on MOF materials continues to focus on such areas. However, there has been a growing interest in understanding the structure and properties at the molecular level to investigate a range of promising next-generation applications such as drug delivery, catalysis, and environmental remediation [47].
Vibrational spectroscopy, including inelastic neutron scattering (INS), has been a popular method of analyzing gas and substrate binding in MOF materials, especially hydrogen, and there was recently an excellent review on the topic by Easun et al in [48]. Therefore, this section will focus on work explicitly concerned with the nature of the vibrational modes and how computational methods such as density functional theory (DFT) can be used to explain the collective lattice dynamics [49]. It will also provide a brief overview of some vibrational spectroscopy studies outwith INS that have encouraged recent work focused on the phonon modes of MOFs located in the low-energy (often referred to as the terahertz) region of the vibrational spectrum.
Civalleri et al reported one of the first studies focused on using DFT to calculate the vibrational modes of a MOF material [50]. They computed the vibrational spectra and the electronic structure properties of MOF-5 using the B3LYP hybrid functional [51][52][53], confirming that DFT could be used to obtain accurate data relating to the vibrational dynamics of MOFs. Shortly after, Zhou et al advanced the study of the lattice dynamics of MOF-5 using a combined experimental and computational approach of INS in conjunction with DFT [54]. The computational results indicated that MOF-5 was close to structural instability with a low shear modulus of ∼1.2 GPa and that the structure could yield other phases under pressure. They noted that their computed spectra agreed well with the experimental data and highlighted several interesting vibrational motions, including the softest twisting modes of the organic linker. Recently, these soft phonon modes were confirmed to be linked to a group-subgroup phase transition in MOF-5, resulting in a reduced bulk modulus and a corresponding change of the pressure derivative of the bulk modulus from positive to negative in the lower-symmetry phase [55].
Zhou et al also reported the quantum methyl rotation present in ZIF-8 using QENS [56]. They showed the ground state tunnelling transitions in both the hydrogenated and deuterated framework and indicated that the barrier to internal rotation was low compared to most methylated compounds in the solid-state. They estimated the methyl rotation barrier to be ∼7 meV with a very low energy rotational oscillation, terming the motion as quasi-free methyl rotation. The effect of gas adsorption on methyl tunnelling in ZIF-8 was also recently investigated by Ryder et al, where the shift in the energy of the transition was seen to depend on the specific gas environment [57].
There have been several vibrational studies concerning ZIF-8, mainly due to its significant porosity and promise for storage applications. Hu et al were the first to use infrared (IR) spectroscopy to investigate ZIF-8 in situ at high pressure (up to ∼39 GPa) [58], showing that the structural modifications observed were reversible only for pressure regions below 1.6 GPa and irreversible beyond this region, resulting in a disordered or amorphous phase. ZIF-8 was also shown to exhibit unusual chemical stability with Fairen-Jimenez et al reporting unexpected adsorption of molecules that were larger than the standard pore aperture [59]. This indicated the existence of significant structural flexibility and was investigated via a combined experimental and theoretical approach. The results showed that the ZIF-8 framework underwent a gate-opening mechanism at 1.47 GPa due to a swing effect of the organic linkers, providing increased access to the internal porosity. The same gate-opening mechanism was reported by Ryder et al to be linked to the collective lattice dynamics in the THz region of the vibrational spectra using INS in conjunction with synchrotron far-IR spectroscopy and DFT [60]. The discovery encouraged further study using INS by Casco et al who investigated the effect of different gas molecules on the swinging motion of the imidazolate linkers [61].
The apparent chemical stability of ZIF-8 encouraged further studies, with Hu et al probing its potential for CO 2 storage [62]. The research involved the dosing of CO 2 into the ZIF-8 framework at pressures of 0.8 GPa. Two unique IR-active signals (from the internal and external CO 2 interactions) were identified, providing insight into the interaction sites. The results showed enhanced storage of CO 2 inside the framework, and the storage behaviour was shown to depend significantly on pressure. Soon after, Kumari et al investigated ZIF-8 using Raman spectroscopy [62], showing temperature-induced structural transformations that allowed the framework to adsorb other gases such as N 2 and CH 4 .
A similar study by Aguado et al reported that ZIF-7 could also undergo a kind of gate-opening mechanism but this time more akin to pore aperture breathing [63]. The results showed that the reversible effect arises from a phase transformation upon guest adsorption-desorption. The motion was again linked to the THz vibrations (figure 3) [60], further highlighting the importance and impact of the lattice dynamics on understanding the structural flexibility and elasticity of framework materials. Research regarding the vibrational properties of MOFs has traditionally concentrated on the mid-IR (MIR) region of the vibrational spectra. However, recently there has been a focus on other spectral regions (far-IR and THz) that can reveal interesting properties. As both experimental resolution (at neutron and synchrotron facilities) and supercomputing power increase, these regions can be studied more accurately, allowing for a better understanding of the lattice dynamics. It is worth noting that similar THz frequency vibrational modes had previously been highlighted concerning the structural destabilization of crystalline zeolites [64]. Greaves et al suggested that the Boson peak in the glassy state of a zeolite could be linked to framework features at the nanoscale, namely, connected ring moieties [64]. It was this suggestion that instigated much of the initial research on the THz vibrations of MOFs [49,60].
The initial work on the lattice dynamics of ZIF-4, ZIF-7, and ZIF-8 proved that the THz region phonon modes could be linked not only to gate-opening and pore-breathing mechanisms but also to shear-induced phase transitions and the onset of structural instability [60]. As mentioned previously, the work on ZIF-8 encouraged further study into the gate-opening mechanism upon gas loading [61]. Additional study of the lattice dynamics of ZIF-4 was performed to gain further insight into the phase transitions and amorphization upon thermal stimulus. The work by Ryder et al analyzed the in situ thermal amorphization of ZIF-4 to a-ZIF and the subsequent phase transition to ZIF-zni [65]. They described the nature of the change in the vibrational modes during the amorphisation process and revealed new insights into the effect that temperature has on the Zn-N tetrahedra, suggesting the thermal lability of the phonon softening could be linked to the likelihood of structural amorphization and instability [65]. Recently, Butler et al used temperature-dependent INS to track the change in the lattice dynamics during the transition from the closed-pore to open-pore phase around 140 K highlighting the change in structural entropy [66].
The discovery that such mechanisms could be explicitly linked to the lattice dynamics and phonon modes in the THz region encouraged additional research using INS and other forms of vibrational spectroscopy, in conjunction with computational methods. DFT calculations have been used to reveal the complex nature of the collective THz vibrational modes for numerous MOF materials and allow for detailed correlations to be made with experiments. The motions and mechanisms revealed by the collective lattice dynamics have also been connected to anisotropic elasticity and anomalous mechanical properties in the Cu-based MOF material known as HKUST-1, such as negative Poissonʼs ratios (auxiticity) and negative thermal expansion (NTE) [67]. Several intriguing low energy modes were identified, including trampoline-like deformations and oscillatory motions of the Cu-based paddle-wheel that have both previously been suggested to be the structural cause of the NTE in HKUST-1 [68].  to distinguish between the pressure and temperature responses of the framework motions [69]. They demonstrated the NTE over a 4-400 K range and reported the main contributor to NTE to be the translational transverse motion of the phenyl rings, which resulted in a unit-cell contraction with increasing temperature. They also noted that the lowest-energy mode was a libration motion of the aromatic ring, which did not contribute to NTE. Recently, Ryder and co-workers demonstrated that the lowest energy A 2g symmetry mode (figure 4) was instead linked to the group-subgroup phase transition from phonon softening upon compression in MOF-5 [55]. The work showed that the quasi-harmonic lattice dynamics could reveal the combined temperature and pressure effects on the structure and mechanical stability of MOFs. Using INS in conjunction with synchrotron IR and DFT similar trampoline-like motions suggesting NTE and aromatic ring rotational dynamics were reported for the low-symmetry Zr-based MOF material known as MIL-140A [70]. Of additional impact was the phonon mode related to coordinated shearing dynamics at 2.47 THz, which was linked to the mechanism suggested to destabilize the framework structure, along the exact crystallographic direction of the minimum shear modulus (G min ) of 3.2 GPa [71]. The discovery reinforced the explicit link between framework elasticity and mechanical properties and stability and encouraged further research in the field of using INS to study the complex lattice dynamics of flexible porous materials such as MOFs.

Quasi-elastic neutron scattering studies on sorbate motions in zeolites
The design and optimization of microporous materials, particularly zeolites, is of great interest due to their impact on industries ranging from petrochemicals and gas separation to water treatment and agro-chemistry. However, significant difficulties arise in their study due to all of the processes of interest taking place within the micropores. Neutron based techniques are thus particularly powerful for probing molecular behaviour deep within the porous structure [72,73], especially in the case of aluminosilicate matrices, where their highly penetrating nature and unique 1 H sensitivity means the neutron interacts almost solely with a hydrogenous adsorbate molecule.
QENS is particularly effective for these systems [13][14][15], allowing the measurement of molecular motions taking place over timescales of 2 ps to 100 ns. Depending on the instrument, one may probe long-range diffusion processes throughout the zeolite framework or local motions at the adsorption/active site. The extra strength of these studies is how complementary they are to classical molecular dynamics simulations (see section 4 for more details on this), particularly in a time where increases in computational power have allowed for increased accuracy in the modelling of the framework, sorbate molecule and the interactions between them [74]. Relatively large systems can be modelled explicitly, e.g., 2 to 100 nm supercells along with simulated times reaching hundreds of nanoseconds, which both improves statistics and allows access to longer timescale behaviours not previously observed [75]. The length and timescales probed match those of QENS experiments, and MD is now an important tool for both qualitative verification of QENS observations, but also now for quantitative analysis and reproduction of QENS data. Indeed, while translational and rotational diffusion coefficients may be calculated, and rotational geometries/diffusion mechanisms observed with MD, more revealing still is a comparison of quantitative experimental observables [76,77], such as the intermediate scattering function and its temporal Fourier transform, the dynamic structure factor. These may be calculated directly from the MD atomic displacements to quantitatively compare simulation results with experiment as discussed in detail in section 4, and for zeolite catalysis in this section.
A comprehensive review of the complementary nature of molecular dynamics simulations and QENS studies, and their significant contribution to zeolite science up until the mid-2000s is provided by Jobic and Theodorou [13]. In this article, we focus more specifically on the direct reproduction of QENS observables such as the intermediate scattering function and the dynamical structure factor. These in turn allow us to reproduce line broadenings and the EISF, which represent the rates and geometries of molecular motion respectively. In section 4 we introduced the relationships between such QENS observables and the atomic displacements obtained from a molecular dynamics simulation (adapted from more detailed accounts in references [13] and [14]), here we will present some particularly illustrative case studies where such operations have provided valuable insight to assist in the QENS data analysis.
One set of studies focused on the motion of propane in zeolite NaY where the interest stemmed from the use of faujasite (FAU) zeolites (the topology group to which NaY belongs) in the petrochemical industry, where the mobility of small hydrocarbons is of particular theoretical interest. Sayeed et al [78] were able to decouple the translational motion from that of rotational motions through the use of different QENS instruments. MD simulations were carried out at a loading of 32 propane molecules per unit cell at the experimental temperatures. They reproduced the intermediate scattering function F s (Q, t) from their MD trajectories using equation (2). Exponential functions were then fit to F s (Q, t) as described in equation (3), where the decay constant of each corresponded in the frequency domain to the HWHM (Γ) of the Lorentzian fit to the experimentally measured dynamical structure factor S(Q, E), with the relationship µ -G F Q t Q t , exp s ( ) ( ( ) ). It was noted that more than one exponential was necessary to fit the MD calculated F s (Q, t), suggesting that more than one motion over the time window probed in the simulations was contributing to the intermediate scattering function. A linear sum of three exponentials (where each one represented a fast, intermediate and slow-motion) gave a good fit over the entire time and Q range. The decay constants of the fast, intermediate, and slow exponentials corresponded to a time range of ∼0.3-7, ∼7-70, and ∼70-700 ps respectively. The intermediate motion was the only one in the timescale of the instrument (with the 'fast' motion considered to be shuttling of molecules at an adsorption site, and the longest timescale motion being that of molecules practically immobile over the instrumental timescale due to very slow diffusive behaviour) and was thus chosen for further study. This was found to be the jump diffusion of propane molecules between adsorption sites in the zeolite, which may be quantified with a range of jump diffusion models. In this case, the Chudley and Elliot jumpdiffusion model [79] which employs a fixed jump distance, the Hall and Ross model [80] which allows for a distribution of jump distances, and the Jobic model [81] which allows for delocalisation of the molecule at the site of its residence between jumps were tested to fit the data.
The jump distances of 3-4.5 Å and residence times of 6-9 ps which fit the model best suggested fast jumpdiffusion between adsorption sites as examined through a low-temperature MD run, and corresponded to diffusion coefficients of 2.0−3.6×10 −9 m 2 s −1 which were similar to those calculated through simply plotting the MSD vs time curve, and applying the Einstein relation to obtain the diffusion coefficient.
The experimental data (measured using the quasi-elastic spectrometer at the Dhruva reactor, Trombay) fit well to the Hall-Ross and Jobic jump-diffusion models to give diffusion coefficients of 0.7−4.0×10 −9 m 2 s −1 showing very good agreement both qualitatively and quantitatively between modelling and experiment. The study is, therefore, a clear illustration of how one may sample molecular dynamics simulations appropriately in terms of the timescale measured by the instruments, while also using atomic displacements over this timescale to reproduce QENS observables in a way that agree well with experiment.
Using a different instrument, the same team were able to sample solely the faster timescale of motion, probing the rotation of propane in NaY [82] through the very wide energy window (allowing for fast rotations to be observed) and ∼3.3 meV energy resolution of the Triple axis spectrometer instrument at the Dhruva reactor, Trombay. Thisensured that translational motions would be contained in the resolution function and not contribute to the quasi-elastic broadening). The atomic displacements from the MD simulations were then extracted in terms of the radius vector of each hydrogen atom (rather than the position vector of each propane molecule as in the previous study into translational diffusion), from which the rotational intermediate scattering was calculated. This function is the same as equation (2).However, the coordinates of each proton have the centre-of-mass vectors of the propane molecule to which they are attached subtractedremoving any translational motion from the calculation.
In contrast to the study of translational diffusion, only one exponential was needed to fit F s (Q, t) along with the constant B(Q) as shown in equation (3), which as described in section 4, represents the elastic incoherent structure factor. When B(Q) is plotted as a function of Q, they were able to fit it to the model of isotropic rotation (with a rotational radius of 1.89 Å, in agreement with the average proton distance from the propane centre of mass).
The experimental data could be fit to the same model with the same radius of gyration suggesting that qualitatively the same modes of motion were being sampled by both the QENS experiments and the rotational sampling of MD data, i.e. freely rotating propane in the NaY pores with no most probable orientation. The in order to directly calculate the rotational diffusion coefficient, D r , which was calculated as 0.82×10 12 s −1 for the MD simulations, and 1.05×10 12 s −1 from the experimental data. The study again illustrates both the qualitative and quantitative agreement that theory and experiment can achieve, and the importance of choosing a specific QENS instrument such that only a certain motion will be observed in a given timescale, which may be readily sampled from an MD simulation. There are, however, times where a number of different motions may be taking place over similar timescales, and the deconvolution and direct qualitative/quantitative matching may be more difficult. A recent example of this arose from studies of ammonia in the Chabazite (CHA) and Levynite (LEV) zeolite structures [83,84]. These structures are of interest for commercial use in NO x abatement catalysis and emission control from diesel engines [85], where ammonia is used as a reductant in a small pore zeolite catalyst to turn poisonous NO x gases into N 2 and water. The CHA framework has a 3-dimensional channel structure allowing diffusion in all directions through 8-membered 3.8 Å diameter windows. LEV houses the same windows; however, they are arranged in such a way that diffusion between cages may only take place in 2-dimensions.
The experimental QENS spectra of ammonia in these frameworks showed a significant elastic component at all Q values suggesting that either a localized motion was present or that there were a significant fraction of immobile molecules on the timescale of the instrument (∼2-100 ps). The EISF was fit with the models appropriate to the system, which were diffusion confined to a sphere (representing each cage), isotropic rotation, and 3-site rotation around the C 3 axis. In terms of the EISF fit, there was some ambiguity at higher temperatures as to whether the correct combination of models was that of translational diffusion confined to a sphere with an immobile fraction, or translational diffusion confined to a sphere with a population of molecules undergoing 3-site rotation around a circle.
Molecular dynamics simulations that were originally performed to validate the diffusion coefficients obtained from the QENS experiments were then able to help decipher which motions were in the timescale sampled by the QENS instrument. The dipole correlation function was extracted from the calculations to assess whether the rotational relaxation time of the NH 3 protons around the centre of mass was within the window observable to the QENS instrument. The dipole correlation showed a very fast decay, suggesting that rotation of ammonia molecules was probably too fast to be observed in the time window of the OSIRIS instrument (ISIS Neutron and Muon Source, UK) and would be more likely to create a flat background for each experimental QENS spectrum. On this basis, it was concluded that the motion observed was that of diffusion confined to a sphere of ∼3.5 Å, with a population of molecules that were immobile on the instrumental timescale.
Further confirmation of the motions present came from the direct reproduction of S(Q, E) from the MD atomic displacements. The program MDANSE [17] was used to define the energy window of the OSIRIS spectrometer and also the resolution of the instrument, thus selecting the frequency at which atomic displacements were sampled from the MD trajectory, and the length of the time window factored into the calculation. The resultant S(Q, E) could then be fit similarly to the experimental data to determine the nature and rate of the motions present. An illustration of the fits is shown in comparison with the experimental QENS spectrum for LEV in figure 5(a). One Lorentzian, along with a flat background, is necessary to fit both the experimental and MD calculated spectra, suggesting one motion is dominant in both the experiment and the MD when sampled on the OSIRIS timescale.
The widths of the Lorentzian function fit to the experimental and MD calculated S(Q, E), (indicative of the rates of motion) are, however, significantly different, exhibiting a far more significant broadening in the simulation. This is further illustrated upon plotting the Q-dependence of the HWHM of the widths, as shown in figure 5(b). The calculated widths from the MD simulations also follow the same trend as experiment (the Chudley-Elliot jump-diffusion model) with a matching jump distance to the experiment. This is the case in both zeolites suggesting the ammonia is behaving the same in both frameworks. While the residence time is significantly shorter in the simulation than experiment, it is clear that the same form of motion is being sampled qualitatively between the two techniques. The jump motion is shown through analysing the trajectory plots of the simulation ( figure 5(c)), which show that jump motions between cages.
The MD simulations were run over a 5 ns period, allowing the probing of much longer timescales of diffusion. It was then shown on the larger timescales that ammonia was actually more mobile in the CHA framework than in the LEV framework by a factor of 2 (giving diffusion coefficients of´- . It was concluded that the presence of double the number of 8-ring windows in the CHA framework compared to the LEV framework was responsible, as this doubled the number of opportunities for the ammonia to perform the jump motion observed both experimentally and in the 'picoscale' sample of the MD simulations. In fact, it was shown that 2 diffusion coefficients could be measured in the system, a 'D s pico ' -dictated by the jump motion between spherical cages observed experimentally and through calculation of S(Q, E), occurring at the same rate in both frameworks (~-´- in CHA) due to fundamental differences in channel structure, and opportunities for performing the jump motion.
The study shows very clearly how phenomena observed by simulation may allow for different modes of motion to be discounted based on their rates happening outside the machine resolution while confirming the modes of motion in the instrumental time window. It also shows how one may use the simulation to expand into timescales too slow to be observed by a specific instrument, and uncover consistencies and inconsistencies between what appear to be relatively similar zeolite frameworks, and their structural influence on the behaviour of a sorbed molecule.
It is also possible to differentiate between motions taking place in a simulation (once one has been determined by experiment), and then predict the influence the zeolite composition (Si/Al ratio) is having on either motion. Recent work has studied the effect of zeolite composition (in this case, Brønsted acid site concentration) on the behaviour of phenol in zeolite Beta [86]. Interest in this system stems from the potential use of zeotype systems in the hydrodeoxygenation (HDO) of lignocellulosic biomass, both as supports for nanoparticles and also as direct depolymerization agents through the action of Brønsted acid sites. Phenol, being the most simple model product of the HDO process, offered a suitable starting point to perform QENS and MD simulation studies to understand the fundamental dynamics of such systems to assist in their further design.
The QENS experiments showed a significant population of the molecules were static on the timescale of the instrument from both the significant elastic component in the spectra and the EISF. The isotropic rotation model was found to fit the data best, suggesting that on the timescale of the instrument the phenol molecule was rotating freely in the zeolite Beta channels (∼6.5 Å in diameter), with progressively more of the population becoming mobile over the 393-443 K temperature range. Rotational diffusion coefficients were then extracted in the range of 2.6−3.33×10 10 s −1 , nearly 2 orders of magnitude lower than that of propane in NaY discussed earlier in this section, reflecting the effect of probing a larger molecule in a more restricted pore system. Molecular dynamics simulations performed were first used to quantify the diffusivity of phenol in the acidic zeolite Beta with a Si/Al ratio of 12.5 (replicating that of the experiment) and also a fully siliceous framework with no Brønsted sites present. The diffusion coefficients reflected the lack of H-bonding adsorption sites in the siliceous framework, with D s over a factor of 4 higher at 393 K than those obtained for the acidic structure. The increase in mobility is reflected in the trajectory plots in figure 6 where a typical phenol molecule is able to explore far more of the siliceous zeolite supercell over the 5 ns simulation.
However, the large activation energies of diffusion calculated (∼30 kJ mol −1 ) in the acidic structure suggest that any diffusive processes would likely take place outside of the time window sampled by the QENS experiment (∼2-50 ps), and thus a higher resolution instrument with a smaller energy window would be necessary to isolate translational diffusion experimentally.
The MD simulations could also be analyzed in terms of the rotational intermediate scattering function in equation (2), where r i is the position of each hydrogen atom with respect to the centre of mass of each phenol molecule. The F Q t , s rot ( ) ( ) plots needed two exponential functions to fit them, one of which had a decay constant which fit comfortably within the spectrometer window, and another of significantly higher energy and thus taking place too quickly to be observable by experiment (potentially due to flipping of the aromatic ring or uniaxial rotation of the hydroxyl group of phenol). Fitting exponentials to the F Q t , s rot ( ) ( ) also allowed for the calculation of the EISF as in equation (3). The EISF for the siliceous system fit the isotropic rotation model very well, while that of the acidic system necessitated the inclusion of an immobile fraction to prevent the model from falling below the data points, as was the case in the experiment. It was therefore clear that when sampling the instrumental time window in the simulations, the motion observed was that of isotropic rotation of phenol in the zeolite channels, with the presence of Brønsted acid sites immobilising a population of phenol molecules. This led to a lower translational diffusion coefficient (and higher activation energy) when the nanosecond timescale was probed in the MD simulations. A further agreement was found between experiment and theory when the rotational diffusion coefficients were extracted from MD calculated the F Q t , s rot ( ) ( ) decay constants, which give D rot values all within a factor of 1.5 of the experimental values.
It is clear from the case studies outlined above that classical MD simulations are now at a standard where they may not only provide qualitative guidance on the motions observed in a QENS experiment, but may deconvolute and reproduce QENS data such that quantitative, direct spectral analysis for systems as complex as molecules diffusing and rotating in a porous material may be carried out. While the emphasis of this section has been on microporous catalysis, the use of porous materials in the fields of gas storage, decontamination, controlled release, and other applications will require such detailed characterisation if they are to be fully optimised, for which the combination of QENS and MD simulations will likely play an important role.

The future: challenges and opportunities
The tools which compare MD/LD to neutron scattering experiments have had their development dictated by the available processing power at the time of their inception. However, there is a much more direct way to draw a comparison, which makes fewer approximations than that of the existing methods. This method is based on the neutron ray-tracing program McSTAS [87]. Recent developments in this software package, now allow the user to define a sample with a given S(Q, E), as calculated from MD or LD. One may then use a so-called instrument definition file for the given instrument on which the experiment was performed, and simulate the direct neutron response of their sample on individual detectors. This method has shown great promise in previous studies [88]. Like any new software tool, for the community to become aware of its benefits, and for a more user-friendly functionality to be developed, there will be an energy barrier of testing/developing on different instruments around the world. The benefits and universality of this method are clear, thus with some time, effort, and promotion, this has immense potential to contribute greatly to the simulation led data analysis for NS experiments. Furthermore, as the systems investigated become increasingly more complex, it will become more necessary to dispense with some of the limiting approximations of current methods.
Both neutron scattering experiments and molecular modelling also produce vast quantities of data. With the increased power of new neutron sources and detectors, this is going to increase. Therefore we are poised to take advantage of recent developments in machine learning and data science [89]. For example, methods such as convolutional neural networks, which were developed for computer vision applications have been shown to interpret experimental spectra [90] and generative adversarial networks show promise for simulations of experimental signals [91]. For neutron spectroscopy to fully embrace the potential of machine learning, a community effort to save, curate, and openly share data will be essential. The neutron scattering community has generally been very progressive and forward-looking when it comes to establishing and applying useful data standards. For example, the implementation of the NeXuS standard [92]. We are confident that neutron spectroscopy will be at the forefront in the 'fourth paradigm' of scientific research.

Disclaimer
This manuscript has been authored by UT-Battelle, LLC, under contract DE-AC05-00OR22725 with the US Department of Energy (DOE). The US Government retains and the publisher, by accepting the article for publication, acknowledges that the US government retains a nonexclusive, paid-up, irrevocable, worldwide license to publish or reproduce the published form of this manuscript, or allow others to do so, for US government purposes. DOE will provide public access to these results of federally sponsored research in accordance with the DOE Public Access Plan (http://energy.gov/downloads/doe-public-access-plan).