2.1. Experimental data and previous results

In our previous research (Woińska et al., 2023), we performed HAR for ten metal–organic complexes featuring TM hydrogen bonds with both known X-ray structures and known neutron structures. The structures were ranked according to the quality of the neutron and X-ray datasets (based on resolution, completeness and R_int) and refinement (evaluated according to R, wR2, goodness of fit and residual density range) performed, leading to the final joint data- and refinement-quality X-ray (HAR)–neutron and X-ray (IAM)–neutron rankings. Structures were ordered from worst to best according to a given quantity and the score was equal to the position in the ranking. If a certain quantity used in the ranking procedure was unknown, the structure automatically obtained the lowest position in the ranking with respect to this quantity. The obtained scores were summed for each structure separately, leading to the data- and refinement-quality ranking, separate for the neutron and X-ray data. Next, the neutron and X-ray rankings were combined to form the final joint X-ray and neutron data and refinement quality ranking. For the top-ranked structures, HAR significantly improved agreement between the TM—H bond lengths and the neutron values compared with IAM. However, for half of the structures, particularly those lower in the ranking, HAR did not provide improvement and sometimes deteriorated the results in comparison with IAM. Consequently, one could ask whether further progression in the determination of positions of hydrogen atoms bonded to TMs is possible with HAR. Since refinement of hydrogen ADPs is often highlighted as a weakness of HAR, this study explores combining HAR with advanced models of hydrogen ADPs such as SHADE3 and NoMoRe. Six previously studied structures were chosen for this purpose. The selected structures are marked with the following REFCODEs: QOSZON (Ho et al., 2003) neutron, (Arion et al., 2001) X-ray; SITKUB (Lam et al., 2003), ZEYVAA (Bakhmutov et al., 2000) neutron, (Nikonov et al., 1995) X-ray; GOJNIF (Cammarota et al., 2019), TIWXOP (Schwamm et al., 2019), XAXMEP (Webster et al., 2005) neutron (Gross & Girolami, 2007) X-ray. Additionally, the high-resolution X-ray structure of the complex containing a Cr—H—Cr bridging bond, KCPTCR (Petersen et al., 1981) neutron, (Macchi et al., 2005) X-ray, was included to examine the role of data resolution. The ranking of structures from the previous study was supplemented with KCPTCR_std (refinement against dataset truncated to the resolution 0.59 Å⁻¹) and KCPTCR_max (refinement against the full dataset) structures. The order of structures in the joint X-ray and neutron data and refinement quality ranking, from best to worst, is as follows: (1) QOSZON, (2) KCPTCR_std, (3) GOJNIF, (4) SITKUB, (5) KCPTCR_max, (6) TIWXOP, (7) ZEYVAA, (8) XAXMEP. All the investigated structures are depicted in Fig. 1. One must bear in mind that this study is based on quite a small set of structures with varying quality of both neutron and X-ray datasets.

Figure 1 Schematic of the investigated structures with heavy metals bonded to hydrogen atoms marked.

2.2. Refinement types and theoretical calculations

The DiSCaMB library (Chodkiewicz et al., 2018) was used for HAR in various versions, all of which employed the same procedure of obtaining aspherical atomic scattering factors: the DFT method with B3LYP functional was used (Becke, 1988; Lee et al., 1988) and interactions within the crystal were included by surrounding the central molecule with a cluster of Hirshfeld partition-derived atomic charges and dipoles for all molecules with at least one atom within 8 Å from the central molecule. The basis set was chosen based on the TM atom to optimize TM—H bond lengths, as described in our previous study (Woińska et al., 2023) and detailed in Table 1. Convergence was reached when the maximum (parameter shift/sigma) was less than 0.001. Geometry optimization in the gas phase in previous work showed that the singlet state had the lowest energy for the compounds considered, so only this spin state was taken into account. In the case of KCPTCR, only geometry optimization for the singlet state was feasible. All refinements for this structure were performed against the full dataset (KCPTCR_max) and against the dataset truncated to a resolution of 0.59 Å⁻¹, which is equivalent to θ_max = 25° (KCPTCR_std) for Mo Kα radiation.

Various types of HAR were conducted with hydrogen ADPs fixed at the SHADE3/NoMoRe values. The refinements denoted HAR_SHADE3/HAR_NoMoRe involved using SHADE3/NoMoRe before each HAR cycle and refining until max(parameter shift/sigma) was below 1.0. The HAR_SHADE3 refinement for ZEYVAA was considered not converged since max(parameter shift/sigma) was oscillating and its values were much higher than 1.0. In refinements denoted HAR_SHADE3_1/HAR_NoMoRe_1, hydrogen ADPs were determined using SHADE3/NoMoRe only before the first HAR cycle and refinements continued until the maximum (parameter shift/sigma) was below 0.001. Additionally, IAM refinements using SHADE3 and NoMoRe were performed. Lastly, IAM and HAR with all hydrogen ADPs set to 0 (denoted IAM_0_HADPs and HAR_0_HADPs) were carried out. Although this is an unphysical model of thermal vibrations, it served as a boundary evaluation for the influence of the values of hydrogen ADPs on TM—H and X—H bond lengths and refinement statistics. Statistics describing all the refinements are provided in the supporting information (Figs. S1–S9 and Tables S2–S8) along with the computational details for SHADE3 and NoMoRe procedures.

2.3. Atomic thermal motions: similarity index

The similarity index S₁₂ (Whitten & Spackman, 2006) was used as a measure of similarity between thermal ellipsoids of a given atom obtained with different diffraction methods or models. S₁₂ quantifies the non-overlapping part of probability densities of a given atom (differing due to different thermal motions) and is given by

where R₁₂ is the overlap integral between the compared probability densities, calculated using only the values of the ADPs. S₁₂ = 0 indicates perfect agreement, whereas 100 denotes the highest possible difference. Notably, S₁₂ is more suitable for confirming high similarity between thermal motions, as even small differences at a level of 1 can result in noticeable visual differences in thermal ellipsoids. Moreover, depending on the estimated standard deviations (e.s.d.s) of ADPs, the e.s.d.s of S₁₂ can be very high, making S₁₂ an untrustworthy measure. Therefore, we use the mean S₁₂ calculated for all hydrogen or non-hydrogen atoms in the molecule, along with their estimated errors calculated using error propagation. The estimated error of the mean S₁₂ is lower meaning S₁₂ is a more reliable measure of thermal motion similarity. In the supporting information, the values of S₁₂ with e.s.d.s for individual atoms in all the structures are available (Figs. S10–S27), as well as the values of a mean S₁₂ with the corresponding e.s.d.s of the mean (Figs. S28–S45).

HAR with anisotropically refined hydrogen thermal motions was attainable for four of the crystal structures. However, a few hydrogen atoms in these structures had to be refined isotropically (Tables S2–S8). For TIWXOP, XAXMEP and ZEYVAA, only isotropic refinement was possible. The results of HAR for these structures have been published in a previous study (Woińska et al., 2023), with the exception of KCPTCR, which has not yet been treated with HAR, and GOJNIF, for which only isotropic HAR results were published. In the case of TIWXOP, the neutron dataset presented challenges, with large and distorted ellipsoids of anisotropically refined hydrogen atoms. One of the methyl groups in TIWXOP was significantly disordered, and only the hydrogen atoms of the largest component could be identified in the difference density map and refined isotropically, with the remaining nuclear density modelled as a ring of hydrogen nuclear density. The temperatures at which the neutron and the X-ray data for KCPTCR have been collected differ by only 8 K; therefore, we decided to compare the neutron thermal ellipsoids and those obtained from refinements against the X-ray data. Fig. 2 shows the structures of QOSZON and SITKUB obtained with HAR (anisotropic): HAR_NoMoRe and HAR_SHADE3. Fig. 3 displays the structures of XAXMEP and ZEYVAA obtained with HAR (isotropic): HAR_NoMoRe and HAR_SHADE3. Fig. 4 includes the X-ray-derived structures of TIWXOP, GOJNIF and KCPTCR along with the neutron structures.

Figure 2 Crystal structures of (a) QOSZON and (b) SITKUB obtained with HAR, HAR_NoMoRe and HAR_SHADE3. Thermal ellipsoids are depicted at the 50% probability level.

Figure 3 Crystal structures of (a) XAXMEP and (b) ZEYVAA obtained with HAR, HAR_NoMoRe and HAR_SHADE3 (HAR_SHADE3_1 for ZEYVAA). Thermal ellipsoids are depicted at the 50% probability level.

Figure 4 Crystal structures of (a) TIWXOP, (b) GOJNIF and (c) KCPTCR (high resolution) obtained based on the neutron and X-ray experiments (HAR, HAR_NoMoRe and HAR_SHADE3). Thermal ellipsoids are depicted at the 50% probability level.

2.4. The effect of data resolution on dynamic structure factors

Plotting the difference between the dynamic structure factor amplitudes [ΔF_calc = F_calc(model1) − F_calc(model2)] for the models allowed us to observe in which resolution ranges there are the highest differences between the models and how data truncation affects these differences. ΔF_calc between HAR and IAM for the full dataset [Fig. 5(a)] indicates the highest differences in low-resolution reflections (with a maximum between 0.2 and 0.3 Å⁻¹) and a sharp decrease at sin(θ)/λ = 0.6 Å⁻¹, which corresponds to the applied resolution cutoff. Applying the cutoff [Fig. 5(b)] has minimal effect on the low-resolution region, mainly reducing the concentration of diverging reflections, making HAR and IAM models slightly more similar. In both cases [Fig. 5(a) and 5(b)], ΔF_calc is quite symmetric relative to the x axis. Regarding how data trimming influences HAR and IAM individually [Figs. 5(c) and 5(d)], differences are less pronounced than between HAR and IAM, and the spread of ΔF_calc increases with higher resolution. HAR results in lower ΔF_calc compared with IAM. For both HAR and IAM, the KCPTCR_max structure tends to have slightly higher dynamic amplitudes of reflections than KCPTCR_std. The effect of thermal motions on dynamic structure factor amplitudes is very subtle; therefore, similar trends are observed for different techniques of estimating hydrogen thermal motion (Figs. S54–S59). Analogous effects for the remaining structures are presented in Figs. S60–S71.

Figure 5 The differences between the dynamic structure factors calculated in a given refinement versus resolution.

2.5. TM—H and X—H bond lengths obtained with various methods

The mean difference (MD) and mean absolute difference (MAD) between the neutron bond lengths and the X-ray bond lengths were calculated for the X—H and TM—H bonds (Fig. 6). Combined standard deviation (c.s.d.) averaged over the bond lengths was used as a measure for uncertainty of MD and MAD and was calculated according to the formula

where σ_X and σ_N are the estimated standard deviations of the X-ray and the neutron bond lengths, respectively. Aside from the plots, MD, MAD, MD/c.s.d. and MAD/c.s.d. values averaged for the X—H and TM—H bonds in all the structures are given in Tables 6 and 7. The two-tailed Welch's t-test (Welch, 1947) was employed to assess the statistical significance of these differences.

Figure 6 MAD between the X-ray/geometry-optimized and neutron values of the X—H and TM—H bond lengths with the c.s.d. averaged for the X—H/TM—H bonds in each structure as error bars. The X-ray values were obtained with various refinement methods. The structures are ordered according to their position in the X-ray (HAR)–neutron data- and refinement-quality ranking from best to worst.

The results indicate that tweaking HAR using more sophisticated methods for describing hydrogen thermal vibrations has minimal impact on the X—H and TM—H bond lengths in the investigated cases. As shown in Fig. 6, MAD between the HAR and the neutron bond lengths is fairly consistent across different hydrogen ADP estimation methods. Even setting hydrogen ADPs to 0 (IAM_0_HADPs/HAR_0_HADPs) does not significantly alter the X—H or TM—H bond lengths compared with conventional IAM/HAR, although some variability exists depending on the structure. MAD plots in Fig. 6 show that in the case of X—H bond lengths, the IAM values of MAD are consistently higher than 0.1 Å, with little variability related to different methods of estimating hydrogen ADPs. HAR generally yields lower MAD values for X—H bond lengths compared with IAM. Moreover, MAD obtained with HAR is the lowest for the structures from the top of the ranking and tends to increase with decreasing quality. For KCPTCR, pruning reflections has a minor effect on MAD for both HAR and IAM refinements.

For the TM—H bond lengths obtained with various HAR versions, it can also be observed that MAD tends to increase with lower position in the ranking and for some structures MAD is slightly higher than for the X—H bond lengths. Surprisingly, IAM yields lower MAD values for TM—H bond lengths compared with X—H bond lengths and, on average, MAD values are only slightly higher than those for HAR. For KCPTCR_std and KCPTCR_max, MAD is relatively consistent across various IAM and HAR types, although conventional HAR provides the best agreement of Cr—H bond lengths with the neutron structure.

The findings summarized in Tables 6 and 7 confirm the observations in the previous paragraph. All types of IAM, irrespective of the method of hydrogen thermal motion treatment, underestimate the X—H bond lengths in the investigated structures by on average 0.012–0.013 Å. In the case where all the variants of HAR and the geometry-optimized structure, MD is very close to 0. However, for HAR and HAR_NoMoRe(_1), MD is slightly different from 0 but still statistically significant. MAD values in these cases are significantly different from 0 but much lower compared with IAM. Geometry optimization results in the lowest MAD (0.014 Å). MAD is twice as high for HAR (0.027 Å) and HAR_NoMoRe_1 (0.028 Å). MAD increases slightly for HAR_NoMoRe (0.032 Å) and HAR_SHADE3(_1) (0.031 Å (0.035 Å), and finally for HAR_0_HADPs (0.042 Å) it is still much lower than in the case of IAM. The differences in MAD for X—H bond lengths for various types of IAM are quite small and the highest MAD is attained for IAM_0_HADPs.

Unfortunately, for TM—H bond lengths, statistical testing is impossible due to the small sample size of only 15 bonds. In this case, geometry optimization, which closely aligns with neutron values, slightly underestimates the bond lengths (MD = −0.005 Å and MAD = 0.017 Å). IAM tends to underestimate the TM—H bond lengths (MD ranging from −0.023 to −0.030 Å), resulting in MAD values of 0.038–0.046 Å, similar to what HAR achieves for X—H bond lengths. HAR, conversely, generally overestimates TM—H bond lengths (MD equal to 0.025–0.032 Å). The divergence of HAR values from neutron values, as estimated by MAD, has a similar range to IAM. However, in HAR the lowest MAD is obtained with standard HAR (0.036 Å) and the highest with HAR_0_HADPs (0.052 Å). In the case of IAM, IAM_0_HADPs results in the lowest MAD (0.038 Å) and standard IAM yields the highest MAD (0.046 Å). In summary, for the studied compounds, HAR performs slightly worse for TM—H bond lengths compared with the general case. The IAM performance for TM—H bond lengths is slightly worse than that of HAR and it appears to improve after using SHADE3 or NoMoRe, unlike the HAR performance for TM—H bond lengths. These conclusions are based on quite a small set of structures and the good agreement between the neutron and IAM TM—H bond lengths could be a coincidence. The role of high data quality, crucial for determining hydrogen positions with HAR in typical X—H bonds, seems to be even more critical for TM—H bond lengths, whereas this might not be an issue for IAM, perhaps accidentally.