A quantum crystallographic approach to short hydrogen bonds

In this work we use high-resolution synchrotron X-ray diffraction for electron density mapping, in conjunction with ab initio modelling, to study short O—H⋯O and O+—H⋯O− hydrogen bonds whose behaviour is known to be tuneable by temperature. The short hydrogen bonds have donor–acceptor distances in the region of 2.45 Å and are formed in substituted urea and organic acid molecular complexes of N,N′-dimethylurea oxalic acid 2 : 1 (1), N,N-dimethylurea 2,4-dinitrobenzoate 1 : 1 (2) and N,N-dimethylurea 3,5-dinitrobenzoic acid 2 : 2 (3). From the combined analyses, these complexes are found to fall within the salt-cocrystal continuum and exhibit short hydrogen bonds that can be characterised as both strong and electrostatic (1, 3) or very strong with a significant covalent contribution (2). An additional charge assisted component is found to be important in distinguishing the relatively uncommon O—H⋯O pseudo-covalent interaction from a typical strong hydrogen bond. The electron density is found to be sensitive to the extent of static proton transfer, presenting it as a useful parameter in the study of the salt–cocrystal continuum. From complementary calculated hydrogen atom potentials, we attribute changes in proton position to the molecular environment. Calculated potentials also show zero barrier to proton migration, forming an ‘energy slide’ between the donor and acceptor atoms. The better fundamental understanding of the short hydrogen bond in the ‘zone of fluctuation’ presented in a salt-cocrystal continuum, enabled by studies like this, provide greater insight into their related properties and can have implications in the regulation of pharmaceutical materials.

1. Single crystal sample preparation All samples were prepared by a method of slow evaporation from solution. Single crystals of N,N'-dimethylurea oxalic acid 2:1 (1) were grown from an ethanolic solution of the two components (2:1 stoichiometry) at room temperature. Single crystals of N,N-dimethylurea 2,4-dinitrobenzoate 1:1 (2) were grown from ethanolic or acetonitrile solutions of the two components (1:1 stoichiometry) at room temperature. Single crystals of N,N-dimethylurea 3,5-dinitrobenzoic 2:2 (3) were grown from an ethanolic solution of the two components (1:1 stoichiometry) in the fridge at 4 °C. Crystals of this 2:2 C2/c form are identified by needle crystals whilst the larger plates indicate the 1:1 P2 1 form of this system. 1 2. Single crystal diffraction data 2.1 High-resolution X-ray diffraction Single crystal synchrotron X-ray diffraction data were collected on systems 1 -3 on beamline I19-1 2 at Diamond Light Source, U.K (NR18193) using a Fluid Film Devices Ltd diffractometer equipped with a PILATUS 2M detector. The sample temperature was controlled at 100 K using an Oxford Cryosystems Cryostream Plus. Data were collected at a wavelength of λ = 0.6889 Å. This wavelength was selected to keep the number of 2θ detector positions to a minimum whilst operating at an energy where the efficiency of the beamline Pilatus 2M detector is improved.
To collect the high-resolution data, the detector arm was positioned up to 2θ 80 ° (0.43 Å). At the lower angles, and to ensure good overlap of reciprocal space up to the highest resolution, data were also collected at a 2θ of 0°, 30° and 55°. At each 2θ position, a complete sphere of data were collected (3 x ω, 1 x φ scans). At the very low angle, the first set of four runs comprised a single φ scan at 2θ = 0° combined with three ω scans at 2θ = 30° and a 0.2°/0.2s oscillation speed. A screening run was first performed at 2θ 0° to determine the optimum level of beam transmission required for each sample. This method implements the inhouse screen19 tool, a programme that provides the user with an assessment of the screen-run X-ray intensities based on maximum resolution of photons and their incidence rate on the detector. For the higher 2θ positions, to ensure accurate high-angle data were obtained, images were collected for longer exposure times using a 0.2°/0.4s oscillation speed whilst the level of beam transmission was increased from the screen19 suggested start value, with this being doubled at 2θ = 55° and quadrupled at 2θ = 80°.
Data collection was performed using the in-house General Data Acquisition (GDA) 3 software and processed using the in-house DIALS 4 for Small Molecule software. The data were additionally scaled and merged using SORTAV 5, 6 software within the WinGX 7 GUI. In SORTAV, an empirical absorption correction was applied based on crystal dimensions and the linear absorption coefficient obtained from the WinGX 7 gui. This method derives an empirical absorption correction for absorption anisotropy by fitting real spherical harmonic functions to the transmission surface. 8 For 2 and 3, data were truncated at 0.48 Å at the limit of observed diffraction (I/σ<2) ( Figure S1, Figure S2). Crystal data are included in Table S1.

Neutron diffraction
Only crystals of 2 were suitable for neutron diffraction measurements. A single crystal of N,Ndimethylurea 2,4-dinitrobenzoate 1:1 (2) was mounted on an Al pin using adhesive Al tape and cooled to 100 K inside a closed-cycle helium refrigerator under He exchange gas. Data were collected at 6 fixed orientations of the crystal on the neutron time-of-flight Laue single crystal diffractometer SXD installed at the ISIS spallation neutron source. 9 The data were processed using the locally available SXD2001 software. 10 The unit cells were taken from X-ray determinations as starting points in the data processing. Crystal data are included in Table S1. An initial starting model for use in the subsequent refinement techniques ( Figure S3) was obtained for 1 -3 by structure solution using SHELXS 11 and refinement using spherical scattering factors in the independent atom model (IAM). A full-matrix least-squares refinement on F 2 was performed implemented by SHELXL 12 in the WinGX GUI. 7 Non-hydrogen atoms (C N O) were refined with anisotropic thermal displacement parameters. H-atoms were nearly all located from Fourier difference maps and their positions and isotropic thermal parameters allowed to refine freely. Exceptions to this are in 3 where some disorder is indicated in the dimethylurea methyl groups and so the methyl H-atoms for C17 were placed geometrically using the HFIX 137 constraint. Refinement data are given in Table S2.

Hirshfeld atom refinement
For 1 and 3 Hirshfeld atom refinement (HAR) 14 was performed on the initial IAM X-ray structures to obtain X-H distances and models for the hydrogen atom anisotropic thermal parameters. HAR was first implemented in the Olex2 15 gui using the NoSpherA2 16 interface under the olex2.refine G-N option. NoSphera2 extras included using ORCA 17 for wavefunction calculation, the def2-TZVPP basis set alongside the PBE method. The HAR process was iterated at a maximum of 10 cycles to reach convergence. Anisotropic refinement of H-atoms was also selected. The asymmetric unit was used as the initial crystal fragment in 3 whilst in 1, the oxalic acid molecule occupies an inversion centre in the crystal structure and so was completed for the calculation.

Multipolar refinement
The experimental electron density distribution of 1 -3 was obtained from the X-ray diffraction data by implementing aspherical scattering factors in a multipolar refinement model. In this method, crystal structure refinement is expanded from considering atoms as point scatterers to considering their aspherical contribution to the atomic scattering (the multipolar model). This model uses the multipole formalism of Hansen and Coppens, 18 implemented in the XD2006 19 program package using the XDLSM refinement routine.
The aspherical electron density is obtained in the multipolar model: Where the first two terms correspond to the spherically averaged core and valence ( ) ( ) electron densities whilst the third term corresponds to the non-spherical valence density. and are the population parameters, κ and κ' are the expansion/contraction parameters, ± ± includes the spherical harmonic angular functions and is the radial function.
The IAM model was used as the start point for 1 -3. Refinement was performed on F 2 and using the following steps. Spherical atom refinement was performed initially including (1) the refinement of a scale factor against all data. This was followed by (2) a high order refinement (sin θ/λ>0.6 Å -1 ) of non-Hatom positional and thermal parameters and (3) a low order refinement (sin θ/λ<0.6 Å -1 ) of H-atom positional and thermal parameters. This allowed optimal atomic coordinates and thermal parameters to be achieved for all atoms. To derive a good starting model for the multipole refinement, (4) starting values for the multipole and κ parameters were imported from the Invariom database 20 using the programme InvarionTool. 21 Each imported Invariom contains theoretically predicted multipole population parameters with relative orientation defined (including local site symmetry and coordinate system) to fit with the starting molecule. These first steps were performed in Linux using a bash script.
In step (5), the H-atom X-H bonds were elongated to distances obtained from either the HAR (1, 3) or neutron refinements (2). In the case of 3, the methyl positional parameters were reset to those from HAR. (6) Anisotropic models for the H-atom thermal parameters were imported into each structure. In the case of 2, these were from the 100 K neutron structure whilst they were used from HAR in 1 and 3. The aspherical atom refinement was carried out by (7) a stepwise multipole expansion (M, κ, D, Q, O, H) for all atoms and truncated at the quadrupole level for H-atoms (refining Q0 only) and the hexadecapole level for non H-atoms. All parameters were refined together in the last step (8).
Only for 2 were the non-H κ parameters refined. κ' parameters were kept fixed for all atoms at Invariom database values. For the H-atoms, κ and κ' were kept fixed at the data base values (ca 1.11 for κ and 1.2 for κ'). The X-ray data used in the refinement were truncated at an appropriate sinθ/λ, determined by considering the agreement between F obs and F calc at the high resolution (1 1.0245 Å -1 , 2 1.0245 Å -1 and 3 1.0233 Å -1 ). The difference mean-square displacement amplitudes (DMSDAs) for all bonds involving non-H atoms were within Hirshfeld limits (< 10 -4 ). Exceptions to this were in 3, involving the methyl groups having values at 10 -14 x 10 -4 . This may indicate a slight rotational motion of this group, which is not adequately modelled by the ADP tensor; it was not however large enough to be included in the model.
Residuals around several of the oxygen atoms in 2 and 3 indicated an extent of anharmonic motion for these atoms. 22,23 This was found mainly for nitro group oxygen atoms but was also present on the acid carbonyl oxygen atom in 2. The anharmonic motion was modelled up to the third order of Gram-Charlier expansion and this improved the residual density around these atoms.
The validity of the multipolar refinement was established by the assessment of the residual density. The residuals in 1 -3, visualised in 3D maps using MoleCoolQT, 24 were found to be low whilst systematic features indicative of an incorrect model, such as un-modelled bonding density, or error in the diffraction data, were absent. The fractal dimension plots resulting from the Meindl and Henn 25 analysis ( Figure S4a to Figure S6a) gave Gaussian distributions and do not contain any significant features that would suggest an error with either the model or in the diffraction data. Plots of the variation in ratio of versus (DRK-plot) 26,27 ( Figure S4b to Figure S6b) remain / within the 5% tolerance required for a satisfactory model, free of significant data problems. 22 Hirshfeld atom refinement (HAR) 14 ADP models were used over (SHADE3 28 ) alternatives in the final multipolar refinements where they resulted in fewer residuals in the vicinity of the short hydrogen bond H-atom position whilst the versus DRK-plot and fractal dimension / plots showed fewer deviations from unity or a gaussian shape, respectively.

Topological properties of the electron density distribution
The electron density distribution (EDD) is then examined by extracting its topological properties (the relative arrangement in space) 44 at bond critical points (BCPs) 45-47 between atoms where a bond path exists and where the gradient of the electron density is equal to zero ( ). Important extracted ∇ ( ) = 0 parameters at BCPs include the electron density ρ(r) and the Laplacian revealing where ∇ 2 ( ) electronic charge is concentrated or depleted. Other parameters include energy densities, quantifying the energetics of the electron density at the BCP. 36 Net atomic charges and populations of the atoms involved in hydrogen bond formation can also be derived using the QTAIM 33, 34 analysis.
Topological parameters at bond critical points (BCPs) between atoms were extracted from the experimental electron density using the XDPROP module in XD2006. 19 These include , and ( ) ∇ 2 ( ) the kinetic ( ) and potential energy densities, calculated from the electron density using electrons. 30 The total energy density (the balance between the two) is: Using XDPROP, BCPs between atoms were searched in rho (ρ) using the CPSEARCH instruction with atom distance limits of 0.7 to 1.7 Å (these limits cover the range of atomic separations from bonds to intermolecular interactions). The topological parameters at each identified BCP were then extracted from the electron density using the molecular graph instruction (MOLGR) in XD.mas instruction file. Monopole populations (P VAL ) and net atomic charges were obtained from the XDGEOM module.
All Laplacian and deformation density maps were generated using the XDPROP module and were plotted in MAPVIEW of the XD2006 19 suite.     5. Ab initio first principle computational methods

As in crystal and optimisation
The molecular adducts used for the first principle calculations are the hydrogen bonded dimers in 1 -3 extracted from the crystal structure following multipolar refinement. The molecular adduct was used either 'as in crystal' or optimised using Gaussian09 32 code, B3LYP functional 33 and 6-31G+(d,p) basis set. 34 Optimisation of the Kohn-Sham orbitals was carried out to the default tolerances of changes in the density matrix, and the geometry to associated default tolerances of the change in forces and displacements of atoms.

Influence of neighbouring molecules
Calculations were also carried out to identify the effect of interactions in the crystal packing and local to the hydrogen bonded dimer molecular adducts, such as the effect of π-stacking and/or extended hydrogen bonding interactions. Here a pseudo-Ewald embedding type calculation was carried out where the additional molecules and interactions were fixed in place and the proton position of the central dimer was moved manually to generate 51 starting structures for the H-atom potential (potential energy surface) studies, as explained in section 5.5 for a single dimer.
Single point energy calculations, electrostatic potential (ESP) and non-covalent interaction (NCI) analysis for the various expansions of the dimer interactions in the gas phase, to resemble those close interactions in the crystal structure, were carried out. This is in the spirit of the 'Ewald embedding' or 'electrostatic embedding' 35 and these models will be referred to as the 'cluster' structures in the text henceforth. Unlike in 'Ewald embedding' we do not convert these interactions into point charges but merely include them as is, in the structure models used for the following studies. Similar to the 'as in crystal' calculations, the structures were not optimised to prevent changes in both the O···O distances and the relative orientations of groups. A comparison of these with optimised periodic structures could prove useful in studying these interactions holistically, but we consider the current model to be sufficient for the present study. The xyz files with the starting structures are available as SI.

NCI analysis and plotting
Non-covalent interaction (NCI) analysis 36 of molecular adducts was carried out for 1 -3 using the Multiwfn code 37 on the .fchk files obtained from optimisation / single point energy calculations using Gaussian09 32 code (B3LYP functional and 6-31G+(d,p) basis set). A high quality grid with 1728000 mesh points was used for these calculations to obtain the reduced density gradient (first derivative of the electron density) and the sign of the Laplacian (second derivative of the electron density). The plots were created using the VMD software along with the visualisation code supplied along with the Multiwfn code. 37 An isosurface value of 0.05 is used for the representation, where the isosurface indicates reduced density gradient (RDG), colour coded using the sign of the Laplacian, with deep blue (-ve Laplacian) indicating strong H-bonding, green (near 0 Laplacian) indicating vdW interactions and red (+ve Laplacian) indicating steric hindrance or an absence of electron density.

Electrostatic potentials
Molecular adducts of the hydrogen bonded dimers were optimised in the gas phase using density functional theory (DFT) using Gaussian09 32 code, B3LYP functional and 6-31G+(d,p) basis set. 34 Electrostatic potentials (ESP) and electron densities were also calculated using Gaussian09. 32 Images were generated using the VESTA software package 38 where the electron density is represented and colour based on the ESP. An isosurface of 0.6 was used. Blue regions indicate regions of high electron density and red indicates regions with low electron density.    43

Boltzmann distributions
The probability (p) of finding a proton at a given position in the PES is calculated from the partition function Z. By varying T we can predict the change in the probability of finding the proton with temperature. This model however, does not account for the change in the O…O distances with temperature. This could further widen the probability distribution of proton positions.