On the temperature dependence of H-U iso in the riding hydrogen model

The temperature dependence of hydrogen U iso and parent U eq in the riding hydrogen model is investigated by neutron diffraction, aspherical-atom refinements and QM/MM and MO/MO cluster calculations. Fixed values of 1.2 or 1.5 appear to be underestimated, especially at temperatures below 100 K.


Introduction
The riding hydrogen model is widely used in refining smallmolecule X-ray diffraction data. Three positional and one isotropic displacement parameter can be constrained to a 'parent atom' that the H atom is 'riding' on, improving the data-to-parameter ratio and ensuring a chemically meaningful geometry. Alternatively, a single isotropic displacement parameter per riding H atom can be included in the leastsquares refinement model while still constraining hydrogen positional parameters.
Predicted H-atom positions usually lead to comparable figures of merit to a free refinement of H-atom positional parameters. This holds even for high-quality X-ray data, extending far into reciprocal space, since the scattering contribution of hydrogen is small and limited in resolution. Therefore predicted positions, e.g. by SHELXL (Sheldrick, 2008), have also been used for 'invariom' (Dittrich et al., 2004) aspherical-atom refinements (Schü rmann et al., 2012;Prö pper et al., 2013). Such H-atom treatment, in combination with elongating X-H vectors to bond distances computed by quantum chemical optimizations of model compounds, provides structures of high quality from conventional diffraction data.
As stated above, the riding hydrogen model can include constraints for isotropic hydrogen displacement parameters. Ratios of 1.2 and 1.5 of H-U iso with respect to the U eq of the parent atom are being used in most refinement programs today. These ratios had been empirically derived for use with room-temperature data. However, most of today's data sets are collected at temperatures of 100 K or lower, making full use of reduced thermal motion, e.g. to reduce the bias of anisotropic displacement parameters on bond distances (Busing & Levy, 1964). We will show that the ratio of H-U iso / X-U eq is temperature dependent, which indirectly follows from Bü rgi & Capelli (2000). Therefore constant H-U iso multipliers are inaccurate; the simple remedy of using temperature-dependent multipliers is proposed herein.
Taking into account the temperature dependence of riding hydrogen treatments of H-U iso is a detail of increasing importance in X-ray diffraction, as experimental data quality is improving with modern detectors and X-ray sources. Taking the effect into account allows removal of a resolutiondependent systematic error that would otherwise only affect low-resolution data, which is where the hydrogen scattering contributes. While the effect of underestimating H-U iso might seem unimportant when only looking at the R factor (which is practically unchanged), the effect can be frequently detected when aspherical scattering factors, 1 which take into account bonding and lone-pair electron-density distribution, are used for least-squares refinement of positional and atomic displacement parameters (ADPs). 2 For charge-density studies, where the aim is to adjust the scattering factor via multipole parameters to the X-ray data, the anisotropic description of atomic displacements should be used. Hydrogen ADPs are usually estimated in such studies. Munshi et al. (2008) have compared competing approaches for such estimates, and the SHADE (simple hydrogen anisotropic displacement estimator) server (Madsen, 2006) is the approach most frequently used for that purpose today. Since the focus of this work is the most frequently used isotropic treatment of hydrogen displacements, we will not discuss the anisotropic description here.

Experimental
Single crystals of the compound N-acetyl-l-4-hydroxyproline monohydrate (NACÁH 2 O) were grown by slow evaporation of saturated solutions prepared in hot acetone. Crystals grow to sizes suitable for neutron diffraction. A series of multitemperature X-ray diffraction data collections 3 at 9, 30, 50 and 75 K on the same specimen with dimensions of 0.34 Â 0.28 Â 0.28 mm (0.5 mm pinhole) were collected at the DORIS beamline D3 at the HASYLAB/DESY synchrotron in Hamburg. The experimental setup consisted of an Oxford Diffraction open-flow helium gas-stream cooling device, a Huber type 512 four-circle diffractometer and a 165 mm MAR CCD detector. A wavelength of 0.5166 Å and a detector distance of 40.3 mm were chosen, allowing a high resolution of d ¼ 0:50 Å or sin = of 1.0 Å À1 to be reached with a single detector setting. The XDS program (Kabsch, 2010) was used for data integration and scaling. Standard deviations of the unit-cell parameters ( Fig. 1) were obtained by calculating the variance of intermediate cells during integration.
A detector correction (Johnas et al., 2006) was applied to properly correct for the effect of oblique incidence (Wu et al., 2002) on the measured intensities. An empirical absorption correction was not performed at this short wavelength; Friedel opposites were merged. The structural model, cell settings but not the atom notation of the original structure determination by Hospital et al. (1979) as given in the CIF file of the Cambridge Structural Database refcode NAHYPL were used as input. Preliminary least-squares refinements were initialized with this model and performed with the program SHELXL (Sheldrick, 2008).
Data sets at 100, 150 and 200 and 250 K were collected on an Xcalibur S diffractometer equipped with an Mo K sealed tube. Here an analytical absorption correction was performed following the method of Clark & Reid (1995) as implemented in the program CRYSALIS RED (Oxford Diffraction Ltd, 2006) employed for data reduction; Friedel mates were not merged. A second specimen was used for these four higher temperatures. High-resolution data (sin = ! 1) were again measured with the exception of the data set at 250 K.
Neutron diffraction data were collected at the OPAL reactor on the Koala beamline at ANSTO, the Australian Nuclear Science and Technology Organization, in Lucas Heights, Australia. Laue neutron data were collected for a single specimen 1.8 Â 1.4 Â 0.5 mm using an unmonochromated thermal neutron beam on KOALA (Edwards, 2011) at 9, 150, 200 and 250 K. Each data set comprises 16, 12, 12 and ten images (each exposure = 42 min) accumulated on the image plate and from which intensities were extracted using LaueG (Piltz, 2011) employing unit-cell dimensions from the corresponding X-ray determination. The CRYSTALS program (Betteridge et al., 2003) was used for the refinement of positions and ADPs for all atoms. An isotropic extinction parameter was required at 9 K due to good crystal quality and comparably large specimen size for the neutron data. CCDC 977814-977817 contain the supplementary crystallographic information for the neutron data. These files can be obtained free of charge from the Cambridge Crystallographic Data  Temperature dependence of the lattice constants of the X-ray data of Nacetyl-l-hydroxyproline monohydrate. Unit-cell parameters and volume are normalized to the lowest data point at 9 K. E.s.d's are also plotted (but may not be visible when small). Connecting lines are only guidelines for the eye.
1 This was only tested for scattering factors from the generalized invariom database , not for those from the UBDB2011 (Jarzembska & Dominiak, 2012), the ELMAM2 (Domagala et al., 2012) nor the SBFA (Hathwar et al., 2011) libraries. For modelling hydrogen scattering, theoretically derived databases have the advantage of higher precision, since experimental scattering factors for hydrogen can only be reliably determined to the dipolar level of the multipole expansion. 2 Since displacement parameters used in this work are either isotropic or anisotropic, we use the abbreviation ADPs, which was recommended to be used only for anisotropic displacement parameters (Trueblood et al., 1996), in a different manner here.
Centre via http://www.ccdc.cam.ac.uk/data_request/cif. A depiction of the molecule with its atomic numbering scheme and anisotropic ADPs at 9 K from neutron diffraction is shown in Fig. 2.

Methods
H-U iso /X-U eq ratios reported here were derived from four independent methods. Benchmark results for NACÁH 2 O were obtained from multi-temperature neutron diffraction. Values for the hydrogen ADPs from multi-temperature single-crystal X-ray diffraction evaluated with the independent-atom model (IAM) cannot reach the accuracy achievable by neutron diffraction. To improve the physical significance of ADPs and their accuracy from X-ray diffraction (Jelsch et al., 1998;Dittrich et al., 2008), we therefore performed aspherical-atom refinements [either Hirshfeld-atom (Jayatilaka & Dittrich, 2008) or invariom refinement (Dittrich et al., 2004), see below]. QM/MM and MO/MO quantum mechanical cluster calculations (for details of how to run such computations see Dittrich et al., 2012) yield normal modes within the 'molecular Einstein approximation'. These were combined with a TLS fit in the TLS+ONIOM approach (Whitten & Spackman, 2006) and were subsequently converted to give anisotropic ADPs for H atoms. Such computations were performed to complement the experimental results (see x3.3).

Aspherical-atom refinements
Two types of aspherical-atom refinements were performed: in invariom refinement (Dittrich et al., 2005(Dittrich et al., , 2006, 2013) the molecular electron-density distribution is reconstructed from Hansen/Coppens' multipole-model parameter values (Hansen & Coppens, 1978) as tabulated in the generalized invariom database . In Hirshfeld-atom refinement (HAR) (Jayatilaka & Dittrich, 2008) the electron density of the asymmetric unit is obtained by a single-point energy calculation.
Invariom refinements. For structure models based on invariom refinements a least-squares refinement of positions and displacement parameters was carried out using the program XDLSM of the XD2006 package (Volkov et al., 2006). The program INVARIOMTOOL (Hü bschle et al., 2007) was used to set up XD system files for that purpose. Refinement was against F 2 with a SHELXL-type weighting scheme, and the R 1 factor was calculated for all reflections with F > 4ðFÞ. Crystallographic details are given in Table 1. CCDC 990102-990109 contain the supplementary crystallographic data for the X-ray structures. CIF files including intensities are only provided for the invariom refinements, since the same intensities were also used for HAR.
Scattering factors, their local atomic site symmetry and invariom names as well as the model compounds these were derived from are given in Table 2. H-atom positions were initially calculated with SHELXL. In invariom refinement the X-H bond distances were then elongated during initial scalefactor refinement to optimized bond distances of the respective model compound for the invariom assigned to the H atom. This new H-atom position was then constrained to have the same shift as the parent X atom. Only U iso values were freely refined. This procedure was followed because it is also feasible when conventional data of lower resolution than the data studied here are available. Moreover, idealized H-atom positions provided better input for the MO/MO cluster computations (see x3.3), since idealized positions facilitate reaching convergence. Ratios of hydrogen U iso to U eq of the parent atom were then averaged for H atoms sharing the same invariom name using the program APD-TOOLKIT. 4 For direct comparison with HAR, free refinement of H atoms was also performed and the results obtained (not shown) are very similar.
Hirshfeld-atom refinement. In Hirshfeld refinement the electron density from single-point energy calculations is used and partitioned into atomic contributions using Hirshfeld's fuzzy boundary partitioning scheme (Hirshfeld, 1977). Fourier transform (Jayatilaka, 1994) then gives aspherical atomic scattering factors. Atomic positions and ADPs are adjusted to best fit the experimental data using these scattering factors. In an improved implementation of HAR in the quantum crystallography program TONTO (Jayatilaka & Grimwood, 2003), cycles of molecular electron-density calculations, asphericalatom partitioning and least-squares refinement are now iterated to convergence in an automatic manner (Capelli et al., 2014). The Hartree-Fock method was used in combination with the basis set cc-pVTZ (Dunning, 1989). A supermolecule cluster approach was chosen to calculate a wavefunction for both molecules of the asymmetric unit for use in HAR (Woinska et al., 2014). The structural model used in HAR included individual positional parameters and isotropic ADPs for H atoms. Ratios of hydrogen U iso and U eq of the parent atom were again averaged for H atoms sharing the same invariom name. As expected and shown before for three urea ADPs of N-acetyl-l-hydroxyproline monohydrate from neutron diffraction at T = 9 K. Ellipsoids at 50% probability (Burnett & Johnson, 1996). derivatives (Checiń ska et al., 2013), both types of asphericalatom models, the Hansen & Coppens multipole model and HAR, give similar figures of merit and anisotropic ADPs of the non-H atoms with experimental X-ray data.

Neutron diffraction
As mentioned before, the H-U eq /X-U eq ratios from neutron diffraction provide benchmark values for this study. One of the advantages of neutron diffraction is that the scattering lengths of the elements that correspond to atomic scattering factors in X-ray diffraction are constant. Stewart (1976) demonstrated that U iso and U eq from single-crystal X-ray and neutron diffraction differ, and that U eq will be in between the arithmetic and geometric mean of the diagonal elements of the mean-square displacement matrix. Since we are interested in the ratio of hydrogen U iso and U eq of the parent atom, conventional least-squares adjustment can nevertheless provide relative reliable experimental estimates of atomic motion at a particular temperature. Equivalent isotropic displacements H-U eq 5 (orthorhombic system) were obtained both by geometric and by arithmetric averaging the diagonal elements of the matrix of the anisotropic displacements of H atoms (Fischer & Tillmanns, 1988), and both give the same ratios within the estimated uncertainty. In contrast to the deposited structural model, refinements were evaluated without using split-atom sites to model rotational disorder in the methyl group above 150 K. Structural models are given in the supporting information.  Table 1 Crystal data of N-acetyl-l-4-hydroxyproline monohydrate from invariom refinements. GoF = goodness of fit; GoFW = goodness of fit (weighted).  (7) 9.884 (6) 9.9026 (2) 9.9408 (2) 9.9748 (2) 10.0123 (2) b (Å ) 9.249 (3) 9.251 (5) 9.250 (7) 9.253 (6) 9.2485 (2) 9.2479 (2) 9.2492 (2) 9.2556 (2) c (Å ) 10.144 (2) 10.145 (2) 10.149 (6)

Theoretical computations
A quantum mechanical cluster computation was performed to complement the experimental results. The computation was initiated using the experimental geometry from invariom refinement at the lowest temperature of 9 K with idealized hydrogen positions and elongated X-H distances. The method/basis set for optimizing these model compounds was B3LYP/D95++(3df,3pd). The utility program BAERLAUCH  was used to generate a cluster of 17 asymmetric units packed around a central unit. The water solvent molecule was optimized together with the main molecule. Preliminary QM/MM calculations [HF/6-31G(d,p): UFF] ensured that this cluster size leads to convergence and is suitable to reproduce experimental ADPs at low temperature. Calculations to obtain final results employed the MO/MO basis-set combination B3LYP/cc-pVTZ:B3LYP/3-21G. Only the central molecule was optimized, whereas the surrounding 16 asymmetric units were kept at fixed positions. Normal modes were calculated and transformed to Cartesian atomic displacements after optimization.
On the basis of the discussion by Dunitz et al. (1988), the temperature dependence of atomic motion can be described in analogy to a Boltzmann-type distribution of the harmonic oscillator. Atomic motion at higher temperatures can therefore be estimated by the formula given by Blessing (1995): Although the molecular Einstein approach underlying the MO/MO calculations is not able to take into account lattice vibrations with acceptable accuracy, such a cluster calculation can provide a H/parent-atom U iso =U eq ratio, which is however dominated by internal atomic motion. Estimates so derived predict a higher ratio than the experimentally observed ratios from neutron and X-ray diffraction and require a U eq scale factor. To reach agreement between theory and experiment, and to take the temperature dependence into account, it was therefore necessary to go back to the TLS+ONIOM approach (Whitten & Spackman, 2006) and to include the experimental TLS contribution, treating the whole asymmetric unit as a rigid body. In this process the internal atomic relative displacement predicted by the MO/MO computation was subtracted from the experimental ADP data at a given temperature prior to the TLS fit (Schomaker & Trueblood, 1968). Both TLS fit and subtraction were performed by the program APD-TOOLKIT. A more sophisticated (and computationally more demanding) theoretical method based on periodic computations of different-sized unit-cell assemblies was studied by Madsen et al. (2013); for reproducing temperature dependence the TLS+ONIOM approach was sufficient.

Results and discussion
4.1. Temperature dependence of U iso /U eq of riding hydrogen and parent atom from X-ray data Prior to further analysis, a way to distinguish H atoms and their chemical environment is required. One choice would be the well established SHELXL AFIX groups. However, this would not distinguish H atoms exhibiting a distinct vibrational behaviour in theoretical calculations, e.g. an OH group in ethyl alcohol and one in phenol. The invariom formalism  allows a finer distinction. Here H atoms that share the same invariom name are in the same covalent bonding environment and have the same number of next-nearest non-H neighbours, so it was used for classification throughout. 7 Vibrational modes of individual invarioms (as derived from their model compounds) in other molecules will be investigated in a forthcoming study.
We can now consider the ratio of hydrogen U iso and parentatom U eq from X-ray diffraction at different temperatures. Initial observations with invariom refinements on d,l-serine (Dittrich et al., 2005) indicated a temperature dependence at very low temperatures. Subsequent tests using the IAM showed that the IAM does not provide the model precision required to obtain significant results (Thorn, 2012). Our first question was therefore whether the U iso =U eq ratios from aspherical-atom refinements on NACÁH 2 O are able to reproduce the temperature dependence seen for d,l-serine. Fig. 3 shows such ratios for several hydrogen invarioms. Values were either obtained from invariom refinements with constrained riding hydrogen positions, but adjusted hydrogen U iso (Fig.  3a), or by free least-squares refinement of positional and isotropic displacement parameters with HAR (Jayatilaka & Dittrich, 2008) (Fig. 3b). As one can see, the programs/ methods used, XD (Fig. 3a) (Volkov et al., 2006) Table 2 Scattering factor assigned during invariom refinement with atom names, invariom names, local atomic site symmetry and model compounds they were derived from.  (Jayatilaka & Grimwood, 2003) (Fig. 3b), give comparable results. Both refinements do indeed show the expected temperature dependence and even distinguish different hydrogen invarioms from X-ray data, although the standard deviation associated with each value is non-negligible (not shown for clarity). 8 At very low temperatures the relative motion of H atoms relative to their parent atoms is clearly appreciably extended compared to higher temperatures. This temperature dependence can be understood by looking at the low-and high-temperature limits of equation (1) as well as the transition temperature between both limits. Such dependence can be understood by considering the evolution of ADPs of atoms of different mass with temperature (Bü rgi & Capelli, 2000). Division of the mass-and temperaturedependent functions f 1 ð!; T; m 1 Þ and f 2 ð!; T; m 2 Þ with masses m 1 < m 2 yields a function with the observed shape. Vibrations of atoms with larger contributions from higher internal frequencies are more prominent at lower temperatures, while at higher temperatures the atomic mass independent external modes dominate the overall amplitudes.
Interestingly, the data set that was an outlier in the expansion of the unit-cell volume at 67-75 K shows further deviations in atomic displacements: hydrogen invarioms of the type H1c[1c1h1h] mainly show a deviating behaviour with respect to the other atoms at higher temperatures. This is due to rotational disorder of the methyl group, and it is easily conceivable that the differences in the lattice constants seen at 67-75 K are due to the rotation becoming more frequent, either starting from this temperature or due to temperature fluctuations at this data point. We have previously studied similar disorder in methylaminobutyric acid hydrochloride by difference electron-density plots and molecular-dynamics simulations (Dittrich et al., 2009). The abnormal temperature dependence of the three H1c[1c1h1h] methyl-hydrogen invarioms is proof that rotational disorder is also present in the acetyl group of NACÁH 2 O. We will study rotational disorder in this molecule and its anhydrous form in more detail in a subsequent study.
Since positional and displacement parameters are correlated, limiting model flexibility (constrained hydrogen positions) seems to make the onset of additional rotational motion more apparent in invariom refinement, whereas free refinement of positions and U iso seems to lead to an overparameterized model in HAR.

QM/MM and MO/MO calculations
We were interested in reproducing the temperature dependence of U iso by using the TLS+ONIOM approach (Whitten & Spackman, 2006). For this purpose the abovementioned two-layer ONIOM computation using the coordinates from the 9 K X-ray diffraction experiment was combined with a TLS fit at each temperature to provide another set of results that includes information independent from experiment.
Details of the procedure need to be highlighted prior to a discussion of the results. Before performing the TLS fit the computed internal modes were subtracted from the experimental non-H-atom ADPs. Contrary to expectation, this is not accompanied by an improvement of the TLS R factors (not shown) with temperature, since the internal ADPs of heavy atoms are mostly spherical in shape and get almost completely absorbed in the TLS ADPs. Nevertheless, such a correction is physically reasonable and we recommend that it is performed. Furthermore, the agreement of the ratio of U iso =U eq seen in the aspherical-atom refinements of the X-ray data improves when this internal TLS contribution is taken into account.
Another detail concerns the low-frequency modes describing the movement of the asymmetric unit in the crystal framework. Low-frequency modes have a very large impact on the overall displacements, which can be derived directly from equation (1). Since the approximations present in the theoretical method do not allow the estimation of these frequencies with sufficient accuracy, low-frequency modes were omitted in the calculation of ADPs. A frequency threshold of  Temperature dependence of the U iso /U eq ratio from the X-ray data of Nacetyl-l-hydroxyproline monohydrate. (a) Invariom refinement with constrained hydrogen positions and refined U iso . (b) Hirshfeld-atom refinement with freely refined hydrogen positions/U iso , both using the same X-ray diffraction data.
200 cm À1 was found to be adequate (Madsen et al., 2013). The required information on the overall displacement is instead taken from the TLS fit, which yields more reliable values. The TLS+ONIOM approach is hence an attractive computational method to understand the ratio of H-U iso /parent U eq when experimental TLS contributions are available. It reproduces the temperature dependence nicely, although rotational disorder cannot be predicted. More work is required for ab initio prediction of the temperature dependence by theoretical computations without any experimental input. So far it can only provide an independent source of information for the internal modes and requires the application of the TLS fit; theoretical methods are nevertheless best suited to provide H-U iso /X-U eq ratios since experimental errors are limited to ADPs used in the TLS fit.
We now compare these TLS+ONIOM results to those from neutron diffraction, our experimental benchmark (Fig. 4). Since neutron diffraction data sets were not collected at temperatures of 30, 50, 75 and 100 K, these data points are absent in the comparison with high-resolution X-ray and TLS+ONIOM results. The TLS+ONIOM approach confirms that individual displacements of hydrogen invarioms are distinguishable mainly at temperatures below 100 K. Rotational disorder cannot be predicted from this approach, whereas it is also detectable in the neutron data at higher temperatures. Good agreement of neutron diffraction and the TLS+ONIOM approach is found for the temperaturedependent ratio of hydrogen U iso or U eq and the parent atom, which also agrees rather well with the X-ray results in Fig. 3. However, at higher temperatures HAR and neutron diffraction show a trend that deviates from the other methods, with the ratio being higher than 1.5, whereas the ratio is smaller in invariom refinement. This is probably due to the predicted/ constrained hydrogen positions in invariom refinement, which impose a beneficial limit on the flexibility of the structural model here. Since the rigid-body fit in the TLS+ONIOM approach is based on the invariom results, a similar temperature dependence to that in invariom refinement is observed.
A comparison of Figs. 3 and 4 indicates an overall surprisingly good agreement for each curve of H-U iso /X-U eq versus temperature over the whole range, especially when taking into account that different methods/experiments were used. All four methods consistently indicate that at very low temperatures the ratio H-U iso /X-U eq can be as high as four, e.g. for H atoms attached to an sp 3 C atom with three non-H-atom neighbours (corresponding to AFIX 13 in SHELXL). Moreover, all four methods consistently confirm (or reproduce in the case of TLS+ONIOM) the temperature dependence that is predicted from equation (1). Conventional IAM structure determinations employing riding hydrogen constraints -and likewise models with aspherical scattering factorsshould therefore take the temperature-dependent ratio into account.

Conclusion and outlook
Four different methods providing the temperature-dependent ratio of H-U iso to X-U eq in the riding hydrogen treatment have been evaluated and compared. Neutron diffraction experiments provide benchmark values. 'Invariom' and 'Hirshfeldatom' aspherical-atom refinements with high-resolution X-ray diffraction data yield very similar results, with the invariom model using constrained hydrogen positions giving a more consistent result, but the Hirshfeld-atom model being closer to neutron diffraction at higher temperatures. Implementing restraints in the TONTO program would therefore be useful. Furthermore, experimental findings can be well reproduced by the TLS+ONIOM approach. Here a single quantum chemical MO/MO cluster calculation is combined with a temperaturedependent rigid-body fit of the non-hydrogen ADPs from aspherical-atom X-ray refinements. All methods show that the ratio of H-U iso=eq /X-U eq , which is usually assumed to be 1.2 or 1.5 independent of temperature, is frequently more than twice as high at lower temperatures. Fixed values of 1.2 or 1.5, as usually used in conventional spherical-atom 'IAM' refinements, are therefore underestimating the relative displacement of H atoms at cryogenic temperatures. Since all methods used here consistently show or reproduce that the H-U iso=eq / X-U eq ratio is temperature dependent, the effect should be taken into account in low-temperature structure determina-Acta Cryst. Temperature dependence of the U iso /U eq ratios obtained by TLS+ ONIOM and neutron diffraction. (a) Internal vibrations from MO/MO ONIOM calculations and TLS fit against ADPs obtained by invariom refinement. (b) Refinement against neutron diffraction data.