Nuclear quantum dynamics in Hexamethylenetetramine and its deuterated counterpart: a DFT-augmented neutron study

Despite being one of the most thoroughly characterised molecular crystals, hexamethylenetetramine (HMT) and its deuterated counterpart (DHMT), are still not fully understood, especially regarding anharmonic and nuclear quantum effects. In this work, an unprecedented combination of experimental techniques, including neutron and x-ray diffraction, inelastic neutron scattering, neutron transmission, and Compton scattering, all augmented ab initio by harmonic lattice dynamics calculations, was applied. The main question that motivated the presented work was the interplay between the phonon anharmonicity and isotope and nuclear quantum effects related to the zero-point energies of proton and deuteron. Signatures of the combined effects of isotopic substitution, temperature, anharmonicity and nuclear quantum effects were found in data from all experimental methods. In the case of neutron and x-ray diffraction, these signatures manifested as systematic discrepancies between the structural and atomic displacement parameters and thermal diffuse scattering obtained from harmonic lattice calculations and their experimental counterparts. To this end, an important effect was found that could not have been explained by the harmonic lattice modelling; the reverse Ubbelohde effect, i.e. the observation that deuteration decreases hydrogen bond length in HMT. In the case of neutron transmission, further discrepancies between theoretical predictions and experimental data were found at cryogenic temperatures. Finally, applying the diabatic theory of the local potential of the intermolecular hydrogen bond in HMT, it was possible to elucidate the degree of anharmonicity of the C–H···N bonds by relating it to the magnitude of the vibrational isotope effect for the C–H bond stretching observed in inelastic and neutron Compton scattering experiments. It was found that the combined nuclear quantum and anharmonic effects of the protons (deuterons) in hydrogen bonds in HMT (DHMT) manifest as systematic discrepancies between the ab initio predictions for the widths of nuclear momentum distributions and the experimental values.


Introduction
HMT was one of the first molecular crystals whose structure was determined by x-ray diffraction (XRD), and one of very few organic compounds known to crystallise in a cubic space group [1,2]. In their pioneering XRD experiment, Dickinson and Raymond were the first to assume that the crystal structure of HMT is composed of molecules of (CH 2 ) 6 N 4 [2]. Since then, owing to its application in the synthesis of many organic compounds [3,4], including pharmaceuticals [5], its usefulness in chemotherapy [6], and its role as a potential constituent of several objects in space and the interstellar medium [7][8][9][10], the crystal and molecular structures of HMT and DHMT have been frequently determined using both XRD and neutron diffraction (ND) [2,[11][12][13][14][15][16][17][18][19][20][21][22].
From the perspective of modelling vibrational properties, HMT has proven to be both a challenging and rewarding system. The HMT molecule is suitable for testing different theoretical predictions and benchmarking force fields using vibrational spectroscopy methods for a number of reasons. First, all the hydrogen atoms are symmetry equivalent in HMT, and second, in contrast to many crystals containing spherical molecules, the HMT crystal has no phase transition below 400 K [23]. Moreover, the HMT molecule can be represented as a very stable and rigid framework formed by an octahedron of carbon atoms and a tetrahedron of nitrogen atoms (see [22] and figure 1).
The direct consequence of the rigidity of the HMT molecule is that the separation of the normal modes into distinct intramolecular and intermolecular components is a good approximation for the HMT crystal, as confirmed by XRD [12] and coherent inelastic neutron scattering (INS) experiments [22]. Quantitatively, the degree of this separation has been established by Raman and infrared measurements, indicating that the lowest intramolecular frequencies are approximately five times higher than the rigid-body intermolecular mode frequencies [24]. The dissection of the total vibrational density of states (total VDOS) in the HMT crystal into intra-and intermolecular modes enabled the analysis of the rigid-body translational and librational amplitudes of the HMT molecules in the crystal in terms of a lattice model in which the translational modes were assumed to obey a Debye distribution and the librations an Einstein distribution [12]. Coherent INS experiments on DHMT crystals enabled the analysis of the dispersion curves of intermolecular modes propagating along various highsymmetry directions using phenomenological models involving interactions between supposedly rigid molecules [22]. A model postulating interactions with both the first and second nearest-neighbour centrosymmetric molecules was found to be successful in describing the VDOS and heat capacity of HMT under the assumption that the deuteration of HMT does not change the intermolecular forces [22].
Despite being one of the most thoroughly characterised molecular crystals, HMT is still not fully understood, especially in terms of nuclear quantum effects (NQEs). NQEs, such as nuclear zero-point energy (ZPE), are responsible for a plethora of effects that manifest themselves through the unusual physical and chemical properties of condensed matter systems and molecules [25][26][27]. In the context of HMT, Kampermann et al first noticed that the anisotropic mean-square displacement tensor, when analysed assuming no coupling between the external (rigid body) and internal vibrations, yields mean square displacements for which rigidbody vibrational ZPE effects are apparent [19]. One of the signatures of the vibrational ZPE effects, in this case, is the fact that methylene protons have internal vibrations whose atom-projected VDOS (apVDOS) is approximately independent of temperature below 200 K [19].
In the case of molecular crystals, NQEs are inextricably intertwined not only with anharmonic effects but also with the characteristics of van der Waals interactions and hydrogen bonds [28][29][30][31][32][33]. Kampermann et al observed that the thermal vibrations in the HMT crystal have an anharmonic component in the temperature range of 120-50 K [19]. A more detailed analysis of the anharmonic effects in the anisotropic atomic mean square displacement tensors of the HMT was conducted by Burgi et al [34][35][36]. Moreover, molecular reorientation in HMT crystals has also been interpreted in terms of intermolecular hydrogen bonds and van der Waals contacts [1,[37][38][39][40][41].
Apart from vibrational ZPE effects, molecular systems such as HMT exhibit non-trivial NQE-related isotope effects. One of them is the normal (reverse) Ubbelohde effect, which is the elongation (contraction) of the hydrogen bond upon deuteration due to the quantum nature of protons and deuterons [42][43][44]. The Ubbelohde effect is yet another manifestation of a more general NQE: the selective strengthening or weakening of hydrogen bonds due to the ZPE effects of specific vibrational modes [26]. As most of the benchmarking of the force fields for HMT has been performed under the assumption that deuteration of HMT does not change the intermolecular forces, the Ubbelohde effect may be responsible for systematic discrepancies in the modelling of vibrational and thermodynamic data of this system.
Another non-trivial NQE-related isotope effect is the equilibrium isotope fractionation of hydrogen and deuterium in molecular systems between two phases (e.g. liquid and vapour water). This quantity can be directly related to the quantum kinetic energy difference between isotopes [26]. Deuterium fractionation upon the formation of HMT is not well understood, especially when it is produced from partly deuterated methanol isotopologues, such as CH 2 DOH and CH 3 OD [8]. Thus, knowledge of the vibrational ZPE values of H(D) in HMT (DHMT) may help understand deuterium fractionation in HMT and its temperature dependence.
The main aim of this work is a critical appraisal of the interplay between the phonon anharmonicity and isotope and nuclear quantum effects related to the zero-point energies of hydrogen and deuteron in HMT and DHMT. To perform this appraisal, we provide extensive benchmarking of experimental data against the predictions based on the harmonic lattice dynamics (HLD) simulations.
Critical for delivering the main aim of this work is the right choice of experimental and computational protocol. The first criterion governing the choice of the experimental techniques for this work has been their 'clear-probe character'. To this end, neutron scattering techniques provide an obvious choice. The deeply penetrating nature of neutrons and the lack of spectroscopic selection rules have largely contributed to the selection of diffraction (ND), neutron transmission (NT), inelastic neutron scattering (INS), and neutron Compton scattering (NCS) as the experimental methods.
VESUVIO is the only neutron beamline in the world that offers five different experimental techniques simultaneously (including the NT) and the only one that can deliver information about nuclear quantum effects via NCS in an isotope-sensitive manner [45,46]. The NT technique, when applied with epithermal neutrons, uniquely available on VESUVIO, can deliver information about the total neutron cross-section in a manner which is independent of the sample density changes due to isotopic substitution [47][48][49]. The NT and NCS data obtained on VESUVIO can be employed to correct the INS spectra for both sample and container self-shielding and to normalize the experimental spectral intensity to an absolute physical scale (barn/energy unit), thereby facilitating the comparison with computer simulations and interpretation [50]. Moreover and most importantly, to the best of the authors' knowledge, there is no experimental and theoretical work so far that would aim at assessing the sensitivity of the NT technique to the combined effects of anharmonicity and NQEs.
Finally, the question of isotopic substitution as a tool for highlighting the combined effects of anharmonicity and NQEs needs to be briefly introduced. The H (D)-sites in HMT (DHMT) are not only expected to play a role in the nuclear dynamics but also to affect the experimental interaction with the measurement probe through effects of multiple scattering, nuclear/atomic absorption and coherent/incoherent scattering ratio. Disentangling the combined effects of anharmonicity and NQEs from the total neutronic responses in the two isotopic systems may thus be hindered by experimental effects, rendering the isotopic substitution an ambivalent tool. In the NCS, unlike in other techniques used in this work, the effects of isotope-dependant nuclear absorption and the effects of coherent/incoherent ratio are immaterial due to the very high incident energy of neutrons and the completely impulsive and incoherent character of the scattering [45,46]. Moreover, the ratio between the single and multiple scattering intensity, shaping the signal correction in the NCS, is fully controlled by the neutron sample transmission in the epithermal incident neutron energy region. On VESUVIO, both NT and NCS are routinely concurrently measured for any given sample, and the ratio of single to multiple scattering is directly inferred from the epithermal plateau region of the neutron transmission curve.
As far as the assessment of the combined nuclear quantum effects and anharmonicity is concerned, the HLD method has been chosen as it provides a robust, easily accessible and rapid benchmarking vehicle. Other methods, such as ab initio molecular dynamics (AIMD) and path integral molecular dynamics (PIMD), provide increasingly better descriptions of the nuclear quantum effects but are much more complex and computationally expensive than HLD.
In this work, we use the fact that the same apVDOSs underlay ND, INS, NT, and NCS data to globally benchmark neutron data collected on HMT and DHMT over seven orders of magnitude in energy and two orders of magnitude in momentum against the results of the HLD simulations. Moreover, we capitalise on the isotopic mass-selective nature of the NCS technique [45,[51][52][53][54][55][56] and compare the experimental values of the ZPE and widths of nuclear momentum distributions of each isotopic species present in HMT and DHMT with their theoretical counterparts calculated based on the underlying apVDOS, thereby providing a more extensive characterisation of the combined effect of anharmonicity and NQEs beyond the usual discussion of protons.

Materials and methods
Sample preparation HMT was purchased from Sigma Aldrich in a powder form of 99% purity (product number H11300) and used directly in this form in the NCS and INS experiments, while single crystals of HMT were grown by evaporation from an aqueous solution. DHMT was prepared in the ISIS Deuteration Laboratory by the reaction of ammonia-d 3 with formaldehyde-d 2 . The powder product was dried under helium and then used to prepare single crystals by growth from a D 2 O solution kept at room temperature under a helium atmosphere. Full deuteration was determined spectroscopically, as described later.

Neutron diffraction
Single-crystal neutron diffraction was performed at the SXD diffractometer installed at the ISIS pulsed neutron and muon source [57], using the Time-of-Flight Laue method. Diffuse scattering data were collected on largesize (~3 × 3 × 4 mm 3 ) single crystals of HMT and DHMT at 300 K in a series of respectively 8 and 10 orientations with exposure times of 17 (HMT) and 8 (DHMT) hours/orientation. The thermal diffuse scattering (TDS) in HMT and DHMT was modelled following the methodology outlined elsewhere [58][59][60][61]. DFT calculation (see subsequent modelling section) and Fourier interpolation were used to compute the phonons for TDS calculations. The diffraction patterns corresponding to first-order thermal diffuse scattering were calculated using phonon polarisation vectors and eigenfrequencies, according to the formalism of Xu and Chiang [61].
Smaller-sized crystals (∼ 1.2 × 1.2 × 1.5 mm 3 ) were used to measure diffraction data of both HMT and DHMT at multiple temperatures in the range of 10 to 200 K. Bragg integrated intensities were obtained for all data using the 3D profile-fitting method implemented in the SXD2001 software [62] and corrected for the Lorentz factor. No absorption corrections for the shape of the crystals were applied. All structures were refined by full-matrix least squares on F 2 using the JANA2020 software [63]. All atoms were refined anisotropically, and the atomic displacement parameters for the 300 K structures were described in the anharmonic approximation using Graham-Charlier expansion coefficients.
X-ray diffraction X-ray diffuse scattering patterns were recorded at room temperature on beamline I12 of the Diamond Light Source [62] equipped with a Pilatus 2M CdTe detector. The crystal was mounted on a spindle and rotated around a vertical axis in steps of 0.25°for a total of 360°, with an exposure time of 2-3 s per frame. Reciprocal space reconstructions were performed using I12 in-house-developed software. The details of image processing are given in the Supplementary Information (SI).

Inelastic neutron scattering
INS experiments were performed using the recently upgraded TOSCA spectrometer at the ISIS pulsed neutron and muon source [63][64][65][66][67]. The powder specimens of HMT and DHMT were loaded into rectangular aluminium sample cells sealed with indium wire and cooled to a temperature of T = 10 K, at which the experiments were conducted. INS data for the samples and empty containers were corrected for incident neutron-dependent attenuation, and the corrected empty container spectra were then subtracted from the corrected sample and container spectra, following Romanelli et al [50]. Only the spectra recorded by the backscattering TOSCA detectors were used for further analysis. The corrected backscattering INS spectrum of the HMT (DHMT) sample was normalised by scaling to the intensity of the INS peak present at approximately 375 cm −1 (312 cm −1 in the deuterated sample). Complete deuteration was ascertained from the intensity of the C-H stretching band. In the protiated sample, this band appeared at approximately 2950 cm −1 ., whereas, in the case of the DHMT sample the band was red-shifted to around 2220 cm −1 (C-D stretching). High-resolution TOSCA data were analysed up to a neutron energy transfer of 4000 cm −1 .

Neutron transmission
Neutron transmission was measured at the VESUVIO spectrometer at the ISIS pulsed neutron and muon source using an established protocol described elsewhere [47,[68][69][70][71][72][73][74][75]. The samples were loaded in rectangular flat aluminium containers oriented perpendicular to the direction of the neutron beam. Aluminium spacers were inserted between the container faces to adjust the thickness of each sample. Spacers of thickness 0.5 mm (1.5 mm) for HMT (DHMT) were chosen. For each sample, such a setting guaranteed that the product of the sample thickness, sample number density, and total free scattering cross-section of all nuclear species present in the sample, weighted by the sample stoichiometry, yielded an asymptotic value of the sample scattering power (transmission) as close as possible to the canonical value of ten per cent (ninety per cent). HMT (DHMT) sample shapes corresponded to cuboids with a thickness of 0.5 mm (1.5 mm) and a cross-sectional area of their faces of 40 cm 2 , a value sufficient to cover the entire circular beam profile with a diameter of 5 cm [76].
The HMT sample was measured at 10, 50, 120, 160, 200, and 300 K, and the DHMT sample was measured at 300 K. The temperature values were stabilised using a closed-circuit refrigerator (CCR) mounted inside the VESUVIO beamline. For the HMT (DHMT) sample at 300 K, data were collected for approximately 34 h (40 h), corresponding to an integrated proton current within the ISIS synchrotron of 3067 μAh (3608 μAh). For the HMT sample, the integrated proton current values were 2435, 2976, 2977, 3517, and 3067 μAh for the experiments conducted at VESUVIO at 10, 50, 120, 160, and 200 K, respectively. In all cases, each subsequent hour of beam time corresponded to an additional value of 180 μAh for the integrated proton current.
Sample transmission was recorded using the time-of-flight (TOF) technique. The TOF values of neutrons transmitted through the sample were obtained using standard neutron beam monitors of type GS20 equipped with 6 Li-doped scintillators. With respect to the position of the VESUVIO moderator, the incident monitor was placed at 8.57 m, the transmission monitor at 13.45 m and the sample at 11.00 m, respectively [76].
For the calculation of total cross-section data from neutron transmission measured in the TOF domain, a bespoke VESUVIO Transmission algorithm [73], implemented in the MantidPlot computational environment, was employed. The algorithm is based on the Beer-Lambert law and works under the assumption that the effects of multiple scattering in sample transmission on VESUVIO are completely negligible, given the instrument geometry [47,73]. Moreover, the algorithm allows obtaining the sample's transmission and total-cross-section data from their counterparts recorded for the sample and container system by normalising the sample-pluscontainer spectra to the empty-container spectra [47].
The negative logarithm of the measured transmission curve, taken at the asymptotic region of the incident neutron energy above 3 eV, was normalised to the value of the respective constant free scattering cross section per formula unit (for details of the procedure, see [47]).

Neutron compton scattering
Neutron Compton scattering experiments were conducted using the VESUVIO spectrometer at the ISIS pulsed neutron and muon source [45,56,[76][77][78][79][80]. The measurements were performed concurrently with the neutron transmission experiments described above; thus, only NCS-specific details are described here.
Isotopic mass-resolved NCS spectra were recorded in the neutron TOF domain using the double-difference (resonance) configuration for the backscattering (forward scattering) VESUVIO detectors. The emptyinstrument background at VESUVIO is flat when viewed by the forward scattering detectors but has a nonvanishing smooth profile when viewed by backscattering detectors. In both cases, the background was temperature independent within the experimental accuracy. For this reason, two sets of NCS experiments were conducted for backscattering detector data analysis of the HMT and DHMT samples. First, the samples were measured in the CCR as a function of temperature. Subsequently, a single experiment was conducted for the empty CCR at room temperature, and the data from both sets of experiments were subtracted.
The background-corrected NCS raw data were subjected to traditional data treatment with multiple scattering (MS) and gamma background (GB) corrections applied in a self-consistent manner [45,46,53,56,78,80,81]. The data were fitted to the TOF domain using linear combinations of the recoil peaks. Each recoil peak is described by the underlying nuclear momentum distribution (NMD) of a given atomic species present in the HMT and DHMT samples. An additional recoil peak represents the contribution of aluminium in the sample container. All NMDs were assumed to be described by Gaussian functions with standard deviations s (hereinafter referred to as the NMD widths) with additional final state effect (FSE) corrections taken into account [45,46,53,56]. In the fitting, stoichiometric fixing [82] was employed for the integral intensities of the recoil peaks of all atomic species, with the exception of aluminium. The aluminium NMD width was fixed at 14 ± 1 Å −1 , which was obtained from a series of calibration experiments performed at VESUVIO for an empty aluminium sample container [83]. Bespoke NCS data treatment routines implemented in the MantidPlot computational environment have been employed [84,85].

Modelling
The phonon-related properties of crystalline HMT and DHMT were modelled using ab initio harmonic lattice calculations (HLD). In the first step of the modelling, geometry optimisation was performed, followed by the calculation of the apVDOSs.
Calculations were performed using CASTEP version 19.1 [86][87][88], which implements the plane-wave/ pseudo-potential formulation of the density functional theory (PW-DFT). The exchange-correlation and van der Waals interactions were approximated by the Perdew-Burke-Ernzerhof (PBE) functional [89,90] and Tkatchenko-Scheffler's (TS) semi-empirical dispersion correction scheme [91], respectively. A plane-wave energy cut-off of 1012 eV was chosen for precise convergence with on-the-fly generated norm-conserving pseudopotentials from the NCP19 library. The total energy convergence tolerance in electronic minimisation was 1 × 10 −8 eV atom −1 , and a Monkhorst-Pack grid with a maximum k-point spacing of 0.05 Å −1 was used.
Following geometry optimisation, the HLD calculations were performed by diagonalisation of the dynamical matrices, computed using density-functional perturbation theory (DFPT). The apVDOSes were obtained from the CASTEP HLD output using the DensityOfStates algorithm in MantidPlot [85,92].
To account for the temperature effect, the apVDOSs obtained from the HLD simulations were scaled using Boltzmann population factors and then used to compute both the neutron spectroscopic [93][94][95] and diffraction observables [60], as outlined below.

Inelastic neutron scattering
In addition, in the case of the simulation of the INS spectra, the underlying phonon dispersion and apVDOSs of individual atomic species present in HMT and DHMT were obtained from the HLD simulations. The details of these simulations are described above, and only INS-relevant details are mentioned here. Following the protocol established in a previous study [93,95,96], the phonon eigenvalues and eigenvectors obtained from the HLD calculation were used as inputs for the calculation of the INS spectra of HMT and DHMT using Oclimax software [97]. In the Oclimax inputs, polycrystalline powder spectra were simulated from HLD in the 0 K limit, assuming an incoherent approximation for both HMT and DHMT and a maximum order of excitation of ten. The fundamental modes and two-quantum events were calculated analytically using the two-quantum events obtained by considering all possible combinations of the fundamental modes. Starting from the three-quantum events, higher-order excitations (multiphonons) were calculated using numerical convolution. Phonon wings were not simulated. The neutron energy transfer range used for the simulations was 8-5000 cm −1 . INS spectra recorded by both forward and backward scattering detectors (scattering angles of 45°and 135°, respectively) were simulated, but only the backscattering simulations were used for comparison with their counterparts recorded at TOSCA.

Neutron transmission
Simulations of the total neutron cross-section as a function of the incident neutron energy were compared with their counterparts obtained from neutron transmission measurements on VESUVIO. For the simulations, apVDOSs obtained from the HLD simulations described above were used as the underlying source of information. The apVDOSs were used as components of the input files for the simulation conducted using NCRYSTAL software (version 2.7.3) [98,99]. The total cross-section data were generated for the incident neutron values spanning a grid of points spaced in logarithmic progression between 0.1 meV and 10 eV. Total (incoherent and coherent) and isolated incoherent approximations were employed to calculate the total neutron scattering cross-section data for both HMT and DHMT. In the case of coherent approximation, apart from the apVDOSs, the input files contained geometrical information (unit cell specification including the fractional coordinates of atoms). The neutron absorption cross-section data were also simulated in NCRYSTAL and added to the simulated total neutron cross-section data for comparison with the experiments. Following the methodology adopted in [47], nuclear resonances were not taken into account in the calculation of the absorption cross-sections as they are not observed from H, D, C, and N atomic species for incident neutron energy values between 0.1 meV and 10 eV. The neutron absorption cross-section data for a given atomic species E , a ( ) s tabulated as a function of the incident neutron energy E, were then assumed to adopt a simple form is the value of the neutron absorption cross-section at the incident neutron energy E 0 = 25.3 meV, taken from values tabulated in the literature [100].

Neutron compton scattering
The second moment of purely Gaussian NMD T 2 ( ) s can be expressed using the following formula [45,56,101,102] where m is the atomic mass of a given nucleus, g ( ) w is its apVDOS tabulated as a function of the frequency of the atomic vibration , w k B the Boltzmann constant, and T is the absolute temperature. Using the second moment of NMD of a given nucleus, the value of the nuclear kinetic energy can be obtained using the following formula: To calculate the apVDOSs of different isotopic species, we assumed that the interatomic forces do not depend on the isotopic mass (and, thus, depend only on the electronic structure). Following this assumption, the apVDOSs of individual atomic species present in DHMT were obtained from the CASTEP phonon calculations performed for the HMT counterpart and postprocessed using the Phonons routine distributed with CASTEP [87,88].
HMT sublimes at 553 K, and because of this exceptional thermal stability, the interactions between HMT molecules in the crystal have been the focus of many studies [1]. The analysis of the HMT crystal packing was performed by Becka and Cruickshank [11]. It was found that each HMT molecule is surrounded by eight others, along 〈111〉, at 6.08 Å between centres (298 K) and by six others, along 〈100〉, at 7.02 Å [11]. The separations between the centres of the nitrogen and hydrogen atoms are right for van der Waals contacts if the radius of hydrogen is taken as 1.3 Å, which is typical in benzene, naphthalene and anthracene [11]. The contacts along 〈100〉 allow a CH 2 group at the top of one molecule to nestle at 90°across the CH 2 group at the bottom of the next molecule. The short H-H distances within a single molecule are 1.77 Å between the two hydrogens of a CH 2 group and 2.40 Å between hydrogens in adjacent groups [11]. These are again appropriate to a van der Waals radius of about 1.3 Å for hydrogen [11]. On the whole, the authors argued that the cohesion of the crystal structure must be described in terms of the contacts between one molecule and its 14 (8 plus 6) neighbours [11]. However, an independent proton nuclear magnetic resonance (H-NMR) study confirmed the results of earlier XRD studies [40,41] and concluded that intermolecular C-H···N hydrogen bonds are the source of both the stability of the HMT crystal and the barrier to the rotation of the HMT molecules inside the crystal, coinciding with the latent heat for crystal sublimation [38]. This result was later confirmed by Eletr using 14 N nuclear quadrupole resonance (NQR) [37]. In addition, Hirshfeld surface analysis [1] showed that both intermolecular C-H···N hydrogen bonds and van der Waals contacts might be responsible for HMT crystal cohesion.
The phonon dispersion curves resulting from the CASTEP HLD calculation based on the underlying optimised unit cell structures of HMT and DHMT are shown in figure 2. The primitive unit cell of the HMT contains one C 6 H 12 N 4 molecule; thus, a phonon dispersion relation is expected to consist of 66 branches in the general direction in reciprocal space [18]. Owing to the separation of the inter-and intramolecular modes and the fact that there is only one HMT molecule in the primitive cell, there are only six dispersive, low-frequency phonon branches due to the intermolecular modes expected for the HMT crystal [18]. The branches below 12 meV, visible in the upper (lower) right panels in figure 2 for HMT (DHMT), represent correlated rotational and translational motions of HMT and DHMT molecules within their respective crystals. Because of degeneracies along specific directions within the first Brillouin zone (1BZ), imposed as a result of the space-group symmetry of elemental bcc structures, the actual number of phonon branches for the intermolecular modes is smaller than six. At the 1BZ zone centre, a longitudinal accoustic (LA) mode is clearly visible as a solid green line in the right panels in figure 2. Going away from the 1BZ zone centre, the energy of the LA mode is clearly higher than the energies of the two transverse accoustic (TA) modes, depicted as solid red and black lines. For a highly symmetric body-centered cubic system, phonon frequencies are triply degenerate at the high-symmetry H points, which also can be clearly seen in figure 2. At this point, the two TA modes have the same energies as the optical mode depicted by the blue (brown) line for HMT (DHMT), and the LA mode has the same energy as the two optical modes depicted as brown and yellow (blue and yellow) lines for HMT (DHMT).
Starting from ca. 40 meV, dispersionless, high-frequency internal modes populate the phonon dispersion chart (see figure 2), reflecting the localised nature of the vibrations of individual atomic species present in the HMT and DHMT molecules.
There is a wide gap between approximately 181 meV (145 meV) and 371 meV (269 meV) in the phonon dispersion chart of the HMT (DHMT) crystal. This gap separates high-frequency symmetric and antisymmetric CH 2 (CD 2 ) stretching modes from the deformation, rocking, wagging and stretching modes at lower frequencies [103,104]. The longitudinal (LO) and transverse optic (TO) modes, lying below the gap, are characterised by relatively small amounts of LO-TO splitting [104]; thus, they are not visible in the charts plotted in figure 2.
The total and atom-projected VDOSes resulting from the same simulation are shown in figure 3. The highfrequency symmetric and antisymmetric CH 2 (CD 2 ) stretching modes have eigenvectors corresponding to large displacements of the hydrogen atomic species and relatively smaller displacements of the carbon atoms, as reflected in the highest-frequency modes present in their respective apVDOSs of the HMT (DHMT). The deformation, rocking, wagging, and stretching modes at lower frequencies involving hydrogen, carbon, and nitrogen, as well as the intermolecular modes involving the correlated rotations and translations of the entire HMT (DHMT) molecules, are expected to be separated from the highest frequency CH 2 (CD 2 ) stretching modes by large gaps. In the lower frequency region of the apVDOSs of hydrogen, a dominant feature is associated with the peaks centred in the region of 150-180 meV, which correspond to CH 2 wagging and CH 2 deformation modes [103,104].
In the case of the apVDOSs of carbon, lower-intensity peaks are visible in the same frequency region. CH 2 rocking and CN stretching modes are associated with relatively large atomic displacements in the region between 100 and 150 meV [103,104]. Below 100 meV, CNC deformation modes and intermolecular lattice modes contribute to the intensities of the apVDOSs of hydrogen, carbon, and nitrogen [103,104]. In the case of DHMT, a similar assignment holds, albeit with different relative mode frequencies and intensities owing to isotope effects.
Neutron and x-ray diffraction We discuss the extensive benchmarking of our HLD-based predictions against their experimental counterparts by presenting the results of the refinement of the crystal structure of DHMT. The lattice constants, atomic coordinates, and atomic displacement parameters (ADPs) of all constituent atomic species present in DHMT, based on the refinement of neutron diffraction data recorded at 10 K, 50 K, 100 K, 120 K, and 300 K, are listed in Tables S1-S12 in the SI. At all temperatures, the refinement was conducted under the constraint of a bodycentred cubic system (space group I m 43 ), with a single molecule of T d 3 symmetry in the asymmetric unit [2,[11][12][13][14][15][16][17][18][19][20][21][22]. The ADPs were refined within the harmonic approximation at all temperatures, except in the case of the 300 K data, for which the ADPs were refined using an anharmonic model, as discussed below.
At 10 K, the refined cubic lattice constant a = 6.9246(10) Å, which differs from the athermal DFT value by less than 0.5%. The atomic coordinates and ADPs obtained from the analysis of the neutron data for DHMT at 10 K within the harmonic approximation are listed in table 1 and are in good agreement with the data reported in the literature [11,16,19,21,105]. Specifically, the parameters of the intermolecular C-D···N hydrogen bonds in the DHMT crystals are in good agreement with literature data [38,40,41]. A total of 24 C-D···N hydrogen bonds (12 per molecule) were detected. Each CD 2 group in one molecule forms two hydrogen bonds with the nitrogen    (14) 0.00375 (7) 0.0146 (13) atoms of the adjacent molecules, and each nitrogen atom in a molecule forms three such bonds with the deuterons in the CD 2 group of adjacent molecules. The structural parameters of the hydrogen bonds in HMT and DHMT inferred from the analysis of the ND data are listed in table 3.
As clearly visible in table 3, the differences in C-D and C-H bond lengths are smaller than the standard deviations of their measurements, while the H···N distance is always larger than the D···N distance at all temperatures by ca. 0.01 Å. Moreover, the measure of the hydrogen bond length, the C···N distance, is systematically shorter in DHMT than in HMT by ca. 0.01 Å, at all temperatures between 10 and 300 K. The observations are consistent with a reverse Ubbelohde effect in HMT [42][43][44].
Our HLD-based DFT calculation is also found to agree very well with the experimental values of the unit cell lengths (a = 6.9584 Å from the geometry optimisation) and ADPs (see table 2) at T = 10 K. The small discrepancies observed between the structural parameters and the ADPs obtained from the HLD-based geometry optimisation and the diffraction data for HMT and DHMT provide a clear signature of the interplay between the temperature, isotopic, anharmonic and nuclear quantum effects. One of the manifestations of this interplay is the reverse Ubbelohde effect, the observation that deuteration decreases hydrogen bond length in HMT. These findings are in agreement with the conclusions of the theoretical considerations present in the literature [106]. Namely, the determination of the bond length from the minimum of the total energy in ab initio theories is only an approximation since the kinetic vibration energy is neglected (it is assumed that the atoms possess an infinite mass and are thus approximated as classical, not quantum objects) [106]. In order to determine the bond length values within ab initio optimisation schemes, it is advised to perform nonperturbative minimization of the Helmholz energy rather than the minimization of the adiabatic potential energy only with the perturbative account of the zero-point motion [106]. The bond lengths at zero temperature should be determined from the minimum of the Helmholz free energy at zero temperature, i.e., from the total energy supplemented by the zero-point vibration energy [106]. In the case of HMT and DHMT, the difference in the values of the ZPE of hydrogen and deuterium shifts the minima of the Helmholtz energies of both systems to the extent that makes a marked difference in their equilibrium bond lengths. Moreover, the zero-point vibrational energy renders the equilibrium bond lengths dependent on the shape of the anharmonic correction to the interatomic potential, i.e. taking into account anharmonicity one always gets a larger bond length than  (4) 152.0 (13) that obtained from the minimization of the adiabatic potential, and this effect is isotope mass dependent [106]. Above zero temperature, two more effects render the picture even more complicated. Firstly, the HLD-obtained structural parameters can no longer be treated as a satisfactory benchmark for their experimental counterparts due to the zero-temperature nature of the HLD approximation. Secondly, Helmholtz energy involves an additional temperature-dependent term that renders the equilibrium structural parameters dependent on isotopic masses, temperature and anharmonicity [106]. Moreover, the ADPs obtained from the HLD approximation can only account for the temperature effect by introducing the Boltzmann phonon population factors for individual isotopic species. As far as the dynamic effects are concerned, in the following, we will concentrate on the analysis of TDS and the ADPs in DHMT at T = 300 K by comparing the results obtained from neutron diffraction against the HLD simulation results. In doing so, we will discuss the combined effect of nuclear kinetic energy and anharmonicity.
The clear separation of the inter-and intramolecular vibration modes in the HMT crystal led to the expectation that the contribution of the intramolecular modes to the TDS was small. Motivated by this, Duckworth and Willis conducted a neutron diffraction study on a single crystal of HMT, including TDS effects on the measured Bragg intensities [13]. The authors concluded that the thermal motion model, which takes into account deviations from the ellipsoidal atomic probability density functions brought about by libration, best fits the experimental data. Going a step further, Powell and Sandor developed a TDS model for HMT crystals in which the intermolecular mode frequencies were expressed directly in terms of axially symmetric central forces describing specific interatomic bonds [24]. The model parameters were determined from the elastic constants of HMT and the inelastic neutron scattering measurements of DHMT. In a revised study of intermolecular forces in HMT crystals by Dolling et al, an 'interatomic' lattice dynamical model was used to analyse the intermolecular mode dispersion curves of DHMT [107]. The authors found that an unambiguous choice of the best force field cannot be made for DHMT within this model, and that the rigid molecule assumption needs to be removed, allowing for coupling between the intra-and intermolecular modes. Additionally, it was suggested that, in order to achieve the best fit of the force field model, the octopole moment of the HMT molecule needs to be considered.
The characterisation of ellipsoids of thermal vibration in HMT from x-ray and neutron diffraction data has been the subject of many debates. As discussed by Stevens and Hope in the context of their XRD work on HMT [21], the geometric and thermal parameters of atoms refined from XRD and ND show significant differences in the presence of aspherical valence electron distribution, and one has to resort to using parameters refined from high-resolution XRD data to benchmark them against their counterparts obtained from the ND.
In their critical appraisal of the atomic excitations in the HMT crystal [18], Willis and Howard showed that the dissection of the vibrational modes into six intermolecular modes and 60 intramolecular vibrational modes enables the separation of the anisotropic atomic mean-square displacement tensors of hydrogen, carbon, and nitrogen in HMT into contributions from the external and internal modes. Furthermore, the authors showed that the mean-square displacements due to the external modes of motion might be described in terms of three second-rank tensors, T, L, and S, with the rigid-body translational tensor T and librational tensor L both reduced to scalar quantities ( u 2 á ñ and , 2 q á ñ respectively), and the translational-librational tensor S vanishing for symmetry reasons. The contribution of the internal modes to the anisotropic mean-square atomic displacement tensors of hydrogen, carbon, and nitrogen was then calculated using a suitable force field fitted to the observed optical frequencies of HMT and by carrying out a normal-coordinate analysis of the molecular vibrations. The agreement between the experimental and calculated values of the tensor components of the total (internal and external contributions) mean-square atomic displacement tensors was within a few standard deviations for all three atomic species [18]. Notably, a constrained version of the analysis by Willson and Howard proved that one needs to extend the TLS rigid-body model to consider the effect of internal molecular motion, especially for hydrogen atoms.
In this study, the theoretical calculation of the TDS in the HMT and DHMT crystals employed the HLDobtained values of the ADPs. The theoretically obtained ADP values of the DHMT crystal at T = 300 K are listed in table 4. These theoretical values must be contrasted with the results of the anharmonic ADP refinement obtained from the analysis of neutron data collected from DHMT at T = 300 K, as listed in table 5. As discussed above (see the Modelling subsection of the Materials and Methods section), the HLD simulations were performed in the 0 K limit, and to account for the temperature effect, Boltzmann population factors were applied. Naturally, such a treatment of the temperature effect cannot account for anharmonicity-related effects such as lattice expansion and the temperature dependence of the VDOSes and apVDOSs. Moreover, as mentioned above, the nuclear kinetic energy of individual isotopic species present in HMT and DHMT renders the equilibrium bond lengths dependent on the shape of the anharmonic correction to the interatomic potential, and this effect isotope mass dependent. This geometrical isotope effect has consequences for nuclear dynamics and, thus, for the TDS. It is usually assumed that the theoretical values of the vibrational mode frequencies depend upon the force constants determined at the minimum of the adiabatic potential energy surface [106].
However, the direct consequence of the fact that the equilibrium structure, via Helmholtz energy, depends on the isotopic masses is that the force constants entering the dynamical matrix underlying the VDOSes and apVDOSs are also rendered mass dependent. Thus, when discussing the discrepancies between the HLDpredicted and observed TDS, one should rather speak of a combined effect of nuclear kinetic energy and anharmonicity.
The combined effect of nuclear kinetic energy and anharmonicity is clearly visible in our data. To quantify this, we use an effective isotopic displacement parameter derived from the diagonal terms of the ADP tensor U U U U 1 3 .

22 33
( ) á ñ = + + It can be seen that the values of U á ñ for hydrogen, carbon, and nitrogen in DHMT at T = 300 K, obtained from the HLD calculation (table 4), are systematically lower than their counterparts obtained from the anharmonic ADP refinement based on the data obtained from neutron diffraction at T = 298 K (table 5). Drawings of the molecular structure of HMT (DHMT) as a function of temperature, with ADPs represented by ellipsoids to elucidate the meaning of the U ij reported in the tables, are presented in Supplementary figure S2 and compared graphically with their theoretical (HLD) counterparts.
The combined effect of nuclear kinetic energy and anharmonicity, clearly visible in HMT (DHMT) at 300 K, becomes much less apparent at lower temperatures, as shown in figure 4, where the experimental values of U á ñ of hydrogen, carbon, and nitrogen in HMT and DHMT are contrasted with their theoretical (HLD) counterparts. ADPs for the HMT crystals were determined from ND data measured at 15,50,80,120,160,200, and 298 K, as well as from XRD data measured at 120 and 298 K, and interpreted as non-spherical, multipole atomic scattering factors [16,19,105]. ADP analysis within the Gram-Charlier (GC) formalism included both thirdorder (C jkl ) and fourth-order (D jklm ) anharmonic thermal parameters [16,19].
A quantitative comparison of the results of an anharmonic ADP analysis in the case of HMT (DHMT) needs to be performed with caution. It should be noted that a problem is often reported in the literature regarding a   [16,19]. A similar problem is expected when the GC formalism-based analysis of ADPs is performed, including higher-order D jklm terms. This was indeed the case in our analysis of the neutron diffraction data of DHMT at 300 K. Nevertheless, fully aware of these aforementioned difficulties, we have decided to perform such an extended analysis as it is, to the best of our knowledge, the first such case in the literature. The results of the analysis are detailed in Supplementary Tables S13 and S14, but are more readily appreciated by a graphical depiction of the nuclear probability density, compared with a standard ellipsoidal representation of the harmonic U ij tensors (figure 5). We close our benchmarking of the diffraction data against the HLD-based predictions at 300K by showing figure 6. In panels A-D, we show the experimental ND and XRD images of TDS in the (h,k,0) layers of the DHMT crystal and its HLD-simulated counterparts. In panels E and F, we show the experimental ND images of TDS in the (h,0,l)-layer of the HMT and DHMT crystals, respectively. As shown in panels A -D, a good degree of agreement was reached between the data and simulations.
However, small discrepancies are visible owing to higher-order phonon terms and the combined effects of nuclear kinetic energy and anharmonicity. In panels E and F, one can see clearly that the combined effects of isotopic substitution, anharmonicity and nuclear kinetic energy produce a marked difference between the TDS images of the HMT and DHMT crystals at 300K.
As mentioned above, the results of benchmarking of the ADPs and TDS obtained from the analysis of the diffraction data against the HLD-based predictions need to be discussed in the context of the combined effects of nuclear kinetic energy and anharmonicity. We will start our discussion by highlighting the need for such an approach due to the insufficiency of traditional treatment. As mentioned by Burgi et al [35], owing to the fact that crystals commonly expand upon heating, which is the case for DHMT (see table S1 in the SI), on average, the interactions between the atoms are thereby reduced, leading to a reduction in the vibrational frequencies and mean-square amplitudes of motion that are larger than expected for a harmonic crystal (positive anharmonicity). In such a case, linear extrapolation of the high-temperature part of the anharmonic curves, describing the ADP values versus temperature in a classical limit, no longer implies zero mean-square amplitudes at 0 K, but negative values (see figure 1 in [35]). Traditionally, the effects of anharmonicity on ADPs are described in either of two ways: (a) by atomic probability density functions (PDFs) that are no longer Gaussian and require higher-order moments or (b) by Gaussian PDFs and ADPs whose increase in the classical regime is no longer linear with temperature and can be modelled using Grüneisen parameters (the so-called quasi-harmonic approach) [34][35][36]. In the case of HMT, the best results were achieved with a model in which the Grüneisen parameters for libration and translation motion were both considered and constrained to be equal [35].
However, both the quasi-harmonic approach and non-Gaussian PDFs modelling of the ADPs are insufficient, taking into account the quantum nature of the lightweight nuclear species present in both HMT and DHMT as the anharmonic effect in both systems intertwined with the nuclear kinetic energy effect. The isotopic dependence of the lattice constants and other crystal parameters is usually modelled theoretically based on the London theory, which uses, as a starting point, the expression for the Helmholtz free energy [106]. Using the first-order perturbation theory, the effect of the zero-point motion on the bond length is expressed in terms of the Grüneisen parameter (which depends on the cubic anharmonicity) and the isothermal compressibility (or the bulk modulus) [106]. Usually, both parameters are calculated at the minimum of the total energy of the crystal. From this approximation follows that two chemically identical crystals formed by different isotopes possess the same isothermal compressibility at 0 K. Generally, it is not true because of the zero-point motion of the atoms [106]. The non-perturbative minimisation of the Helmholtz energy rather than the minimization of the adiabatic potential energy only with the perturbative account of the zero-point motion avoids this approximation, and the bond lengths and lattice constants are calculated at the minimum of the free energy, i.e., at the appropriate equilibrium interatomic distance [106]. The calculation utilizing Helmholtz energy leads to the conclusion that the account of the zero-point vibration influences noticeably the value of the harmonic force constant [106]. Moreover, in the presence of anharmonicity, the influence of the nuclear kinetic energy on the harmonic, cubic, and quartic force constants is disproportionally larger than that on the equilibrium interatomic distance [106]. The nuclear quantum kinetic energy and anharmonic effects modulate the shape of the TDS and thus influence the ADP values. As our results seem to suggest, in the case of HMT and DHMT, contrary to popular belief, this combined influence is already present in the limit of the zero temperature and becomes more visible as temperature increases. The root-mean-square deviation (RMSD) value between the simulated and experimental INS curve was 2.15 cm −1 (1.98 cm −1 ). Detailed assignments of the individual vibrational modes of both HMT and DHMT have been performed elsewhere [104,[108][109][110], where only the main features relevant to our analysis of the anharmonicity in HMT (DHMT) will be considered. The assignment of the HMT (DHMT) modes that have a significant amount of projection on the hydrogen (deuteron) apVDOS is shown in table 6.
The characteristic feature of all vibrations with frequencies above 300 cm −1 listed in table 6 is that they all involve CH 2 (CD 2 ) rocking, twisting, stretching, wagging, and scissoring modes, which are strongly coupled to the motions involving the carbon and nitrogen atoms in the backbones of the HMT (DHMT) molecules. With respect to the high-frequency modes, in the case of HMT, a broad peak centred at the energy transfer value of 2956 cm −1 corresponds to the asymmetric CH 2 stretching mode. The counterpart was observed for DHMT at 2214 cm −1 . The CH 2 protons are believed to form long and weak hydrogen bonds with the nitrogen atoms of adjacent molecules that stabilise the crystal structure of HMT (DHMT) [1,[37][38][39][40][41]. The magnitude of the  vibrational isotope effect on the longitudinal C-H(D) stretching mode in HMT (DHMT) was 1.34, which is typical for a medium degree of anharmonicity in hydrogen bonds [111]. As we shall see, this result is important in the context of the analysis of the NQEs in HMT and DHMT based on the results of our NCS experiments.

Total neutron scattering cross-sections via neutron transmission
Total neutron cross-section data for HMT and DHMT are shown in figures 9 and 10, respectively. The modelling of the total cross-section data using the total (incoherent and coherent) approximation performs very well in the case of HMT at all temperatures in the range from 10 K to 300 K. Also, in the case of DHMT, the agreement between experimental data and simulation is very good at 300 K. VESUVIO transmission experiments are unique in their insight into the effects of nuclear dynamics. Namely, the entire excitation range present in the apVDOSes of all constituent nuclear species, stretching from low-energy intermolecular modes at  tens of meV towards the highest-energy intramolecular modes, underlying the theoretical simulation, is tested at an unprecedented incident neutron energy scale covering seven orders of magnitude [47,[68][69][70][71][72][73][74][75]. Such unprecedented benchmarking capability is due to the ability of VESUVIO to capture the neutron transmission in the epithermal region of incident neutron energies, where the negative logarithm of the transmission data exhibits a plateau which can be normalised to the total free neutron cross-section of each sample weighted by the sample stoichiometry, independent of variations of the effective sample density or thickness with temperature [47,[68][69][70][71][72][73][74][75]. This unique feature of VESUVIO transmission data analysis is well demonstrated in figure 9, where the epithermal plateau regions of the total neutron cross-section curves are very well reproduced, independent of the temperature conditions. The epithermal incident neutron energy region of the total scattering crosssection data of the DHMT sample is clearly visible in figure 10. It shows the isotope effect on the total neutron cross-section, whereby, upon deuteration, the total cross-section values of the DHMT sample are no longer dominated by the large value of the incoherent scattering cross-section of hydrogen. In consequence, the value of the total neutron cross-section at the incident neutron energy of 10 eV drops from 14.2 barns for HMT to a value of 5.0 barn in the case of DHMT. Another consequence of deuteration is the significantly improved visibility of the Bragg edges due to coherent neutron cross-section effects in the DHMT sample. In figure 10, two sharp Bragg edges are visible in the experimental data at incident neutron energies of 2.5 and 5.0 meV, followed by a broader feature at 9.5 meV. The modelling of this part of the total neutron cross-section curve of DHMT, based on the DFT-optimized structure of the unit cell obtained from our ab initio simulations, results in a theoretical curve that shows a good agreement with its experimental counterpart. A series of Bragg edges are clearly visible in the simulated data for the incident neutron energy values of 1.7, 2.5, 3.4, 5.0, 9.5, 11.0 and 13.5 meV, and another series of smaller and less-resolved edges at 5.9, 6.8, 7.6, 17.8, 18.7 and 19.6 meV. The latter series of Bragg edges is not reproduced in the experimental data owing to data statistics. The same series of Bragg edges are present in the simulations of the total neutron cross-section of HMT, but their visibility is much reduced by the overwhelming effect of the incoherent scattering cross-section of the hydrogen.
Importantly, the underlying assumption in all these simulations is the harmonic lattice treatment within the Born-Oppenheimer approximation (BOA). Consequently, it is not expected that simulations reproduce the effects of lattice expansion and anharmonic effects in the Debye-Waller factor because of the nonlinear dependence of the ADPs of individual nuclear species on temperature. It is also not expected that effects such as the shortening (or lengthening) of bond lengths upon deuteration due to the presence of NQEs will be reproduced by ab initio simulations and propagated into the simulations of the total neutron cross-sections. In the case of the experimental data of the HMT, no clear signs of anharmonic effects, such as lattice expansion with increasing temperature or nonlinear temperature dependence of ADPs of individual isotopic species or isotope effects on the bond lengths due to the existence of the NQEs are visible because of the combination of the dominating effect of the incoherent scattering cross-section of hydrogen, data statistics, and the finite resolution of the instrument. However, another manifestation of the NQEs is possibly visible in the low-temperature NT data of HMT (see figures 9(a) and (b)). Namely, the NT data recorded at 10 and 50 K for the incident neutron energies below 0.03 eV seem to differ from the HLD-based predictions systematically, and this discrepancy seems to vanish at higher temperatures. Moreover, the isotope effect in the experimental NT data for HMT and DHMT is most likely also visible. As a result of the thermal effect on the apVDOSs of hydrogen, carbon, and nitrogen, the minimum present in the simulated total cross-section data of HMT shifted its position from 165 meV at 10 K to 171 meV at 200 K and disappeared almost completely at 300 K. As a sign of a dynamic isotope effect, upon deuteration, this set of minima should be present in the DHMT sample at lower values of the incident neutron energy. In the case of the DHMT data recorded at 300 K, a broad minimum in the experimental values of the total cross-section was observed in the region between 100 and 150 meV, but the localisation of the minimum was made very difficult by the presence of a series of Bragg edges and data statistics.

Neutron compton scattering
Thus far, we have benchmarked the results of the HLD simulations for HMT and DHMT against the experimental data obtained using ND, INS, and NT techniques. To elevate this benchmarking exercise to another level, we compare and contrast the NCS observables with the HLD-based predictions. In doing so, we focus on the full potential of the NCS technique, whereby the comparison is conducted for all constituent atomic species, with special attention paid to the isotope effect between hydrogen and deuterium. Moreover, in the context of the combined effects of the nuclear kinetic energy and anharmonicity, we discuss a possible connection between the magnitude of the isotope effect on the widths of nuclear momentum distributions of H and D in HMT and DHMT and the isotope effect on the vibrational frequencies observed in the INS data.
We begin our discussion by presenting the mass-resolved NCS spectra for HMT and DHMT (figures 11 and 12, respectively). The corrected NCS data of HMT recorded at 300 K were fitted using stoichiometric constraints with underlying Gaussian momentum distributions of hydrogen, carbon, nitrogen, and the aluminium sample Figure 11. NCS data and fit of the HMT sample, T = 300 K. The total (focused in the TOF domain) NCS data (black points and error bars), total fitted signals (solid black lines), and the fitted recoil peaks of hydrogen (olive line and shaded area), carbon (red line and shaded area), nitrogen (blue line and shaded area), and aluminium from the sample container (grey line and shaded area). Data and fits for backscattering and forward scattering are shown in the top and middle panels, respectively. The bottom panel shows the isolated recoil peak of the hydrogen that is transformed and focused in the space of its longitudinal momentum. The black points and error bars in that plot signify the focused data, the solid red line represents the fit to the data, the green line represents the residual, and the blue line together with the shaded blue area shows the instrument resolution function for the hydrogen in its longitudinal momentum.
container. The Gaussian NMD of hydrogen is well resolved when the hydrogen recoil peak is isolated and focused in the hydrogen longitudinal momentum space (see the bottom panel of figure 11).
The corrected NCS data of DHMT recorded at 300 K, together with the underlying Gaussian momentum distribution fits of the deuteron, carbon, nitrogen, and aluminium sample container recoil peaks, are shown in figure 12. In addition, in this case, the isolation of the deuteron recoil peak and its focus in the deuteron longitudinal momentum space led to a very good fit using the underlying Gaussian NMD. It is worth mentioning that VESUVIO shows very good momentum space resolution for both hydrogen and deuteron, thereby enabling a high-quality comparison of NMDs of both nuclei, which is very important in the context of the assessment of isotope effects in nuclear quantum dynamics.
A similar situation is observed in the case of the NCS data of the HMT sample recorded at 10, 50, 120, 160, and 200 K. Gaussian NMDs of hydrogen, carbon, and nitrogen describe the nuclear quantum dynamics in the HMT down to limiting low temperatures, and the effects of anharmonicity do not appear to manifest in the entire temperature range, at least at the level of the NMD profile shapes. In the case of DHMT, data are available only at room temperature, but the lack of anharmonic effects in the NMD profile of deuterons at this temperature allows us to safely assume that no anharmonic effects can also be detected at the level of the NMD profile shapes below 300 K, down to 10 K.
A further test for the presence of anharmonic effects in NMDs in HMT and DHMT can be conducted by plotting the experimental values of the NMD widths and the values of the excess kurtosis of the NMD (denoted as c 4 ) against their theoretical counterparts based on the HLD predictions, with the temperature effect accounted for using the Boltzmann population factor. In the case of c 4 , the experimental values were all zero because no other NMD models beyond the purely Gaussian NMDs could be satisfactorily fitted to the experimental data. The theoretical predictions for NMD widths are contrasted with NMD data in figures 13 and 15, in the case of HMT and DHMT, respectively. The theoretical predictions for the values of c 4 in HMT and DHMT are plotted in figures 14 and 16, respectively. Figure 13 clearly shows that the HLD approximation accounts very well for the experimental NMD widths of carbon and nitrogen in HMT, whereas it fails to properly account for the   experimental hydrogen NMD widths in the entire temperature range. Moreover, the results of the theoretical predictions of the excess kurtosis values for HMT, based on the HLD approximation, predict Gaussian NMD shapes within experimental accuracy; thus, no signatures of the anisotropy of the local harmonic potential are expected for all constituent nuclei present in HMT. The same situation is encountered in the case of DHMT (see  plateau region in figure 8 in [111], with equilibrium C-N distances in the region of 3 Å or more, characteristic of weak hydrogen bonds and consistent with the C-H···N bond lengths of 3.777 Å inferred from our diffraction data analysis presented in this work. McKenzie et al assumed a symmetric hydrogen bond case in their model, where the donor and acceptor have the same proton affinity [111]. In this model, the X-H···Y system with X-Y separation in the region of 3 Å is characterised by effective single-dimensional potential curves for the diabatic and adiabatic states (corresponding to isolated X-H and H-Y bonds) that can be well described as Morse potentials. Using our NCS result for Ω H /Ω D of 1.28 ± 0.04 and the properties of the vibrational isotope effect for an underlying Morse potential [111], we can conclude that the anharmonicity parameter D 4 0  c = W (where D 0 is the potential well depth defined relative to the dissociated atoms) is 0.13 ± 0.03, a value typical for a moderate degree of anharmonicity. Such an anharmonic potential could be modelled by a Morse potential with the binding energy of an isolated X-H bond in the region of 1.1 eV.
In the same model, the adiabatic electronic ground state potential energy calculated within the BOA can be characterised as a double-minimum potential with a small barrier in the middle (see the top panel of figure 2 in [111]). From the NMD perspective, such an adiabatic electronic ground-state potential would produce a radial momentum distribution of protons and deuterons that are weakly bi-modal and may show weak quantum selfinterference effects of proton (deuteron) wave functions [45,53,56,[116][117][118][119][120][121][122][123][124][125][126]. However, because of the statistics of our NCS results for HMT and DHMT, one cannot infer any signatures of such interference from the analysis of the shapes of the NMD profiles of H(D) in HMT(DHMT), as both types of NMD profiles can only be fitted with Gaussian profiles.
Another reason for the experimental values of the NMD width of protons (deuterons) in HMT (DHMT), which are systematically lower than their HLD-based predictions, could be the nuclear quantum delocalisation of hydrogen (deuteron). Namely, ab initio HLD calculations provide only a harmonic framework for benchmarking the NCS experimental data. In this framework, only electrons are quantum objects, and protons and deuterons, along with other constituent nuclei, are treated as classical particles that are perfectly localised in space. Thus, the HLD-based predictions of the apVDOses and vibrational zero-point energy values may be systematically shifted from their experimental counterparts. By virtue of the Heisenberg uncertainty principle for the ground state of the QHO, both the nuclear ground-state wave function in position and its counterpart in momentum space are described as Gaussian functions, and the product of their widths is constant. Thus, a lower value of the NMD width indicates a lower value of the width of the proton (deuteron) nuclear wave function in momentum space, and, consequently, a higher width of the proton (deuteron) nuclear wave function in position space, and thus a greater degree of nuclear quantum delocalisation in space. The values of the NMD widths of a proton (deuteron) in HMT (DHMT), obtained in our NCS experiments, are 4.93 ± 0.02 and 5.59 ± 0.06 Å −1 , respectively. Using equation (2), these values correspond to the values of nuclear kinetic energies of 146.7 ± 1.2 meV and 94.4 ± 2.0 meV for proton (deuteron) in HMT (DHMT). These values are almost an order of magnitude higher than the thermal energy at 300 K, signalling a high degree of quantumness of both nuclear species, whereby the ZPE effects dominate their entire dynamics. Using the uncertainty principle for position and momentum for the ground state of the QHO, we obtained estimates of the nuclear quantum delocalisation of protons (deuterons) in HMT (DHMT) of 0.203 ± 0.001 and 0.179 ± 0.002 Å, respectively. As expected, the deuteron in DHMT, being heavier and thus closer to a classical nuclear particle description, is less delocalised in space than the proton in HMT. The increased spatial delocalisation of the protons makes them more likely to explore larger spaces around the local minima of the potential of the mean force. In doing so, protons experience more likely regions of the potential, which are anharmonic; thus, the effects of anharmonicity will be more manifested in the case of protons in HMT compared to deuterons in DHMT. However, overall, the effects of anharmonicity and nuclear quantum delocalisation will be inextricably intertwined, and only the net effect of their interplay will be observed in the form of vibrational isotope effects in the NCS and INS experiments.

Summary
In this study, an unprecedented combination of experimental techniques was applied to characterise the nuclear dynamical properties of two molecular crystals, HMT and its deuterated counterpart, DHMT. The experimental results obtained for HMT and DHMT were benchmarked against the predictions of HLD simulations.
The harmonic picture of the lattice dynamics proves to account for the total neutron cross-section curves, recorded as a function of the incident neutron energy over seven orders of magnitude at temperatures above 50 K. Below 50 K, discrepancies are visible between the HLD-based predictions and experimental data that may be due to the combined effects of anharmonicity and NQEs.
The HLD modelling also accounts satisfactorily for the diffraction results in terms of the equilibrium geometry and atomic displacement parameter values. However, the discrepancy between the experimental TDS and ADP data and their HLD-based counterparts is clearly visible. This discrepancy is attributed to the combined effects of anharmonicity and nuclear kinetic energy. Another signature of the NQEs in the diffraction data is the reverse Ubbelohde effect, clearly seen by comparing the lengths of the C-H···N and C-D···N hydrogen bonds.
As an ultimate test for the manifestation of the anharmonic and nuclear quantum effects in HMT and DHMT, we benchmarked the experimental NCS results against the predictions from the HLD-based simulations. In this case, the combined effects of the nuclear kinetic energy and anharmonicity of the local potentials of the mean force of the protons (deuterons) in the hydrogen bonds in HMT (DHMT), as well as the nuclear quantum delocalisation of both nuclear species are clearly visible in the systematic discrepancies between the HLD-based predictions for the NMD widths and the experimental values and the isotope effects on the NMD widths. Although these effects cannot be separated in the analysis of the NCS observables, it can be shown how they are intertwined, jointly dictating the nature of the nuclear quantum dynamics of protons and deuterons in both molecular crystals.
Overall, the combined use of neutron and x-ray diffraction, inelastic neutron scattering, neutron transmission, and Compton scattering, all augmented ab initio by harmonic lattice calculations, shows great potential for critical appraisal of the combined effects of anharmonicity and nuclear quantum effects of individual atomic species in molecular crystals such as HMT and DHMT.