A quantum crystallographic approach to short hydrogen bonds

Abstract

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.

---

Introduction

Hydrogen bonds are important in biology and chemistry for their roles in maintaining structure,^1,2 for molecular recognition³ as well as in facilitating reaction pathways.^4,5 The formation of hydrogen bonds occurs within a set of conditions, as defined by Etter's rules,⁶ where available hydrogen bond donors interact with available acceptors in an order determined by the strength of the resulting interaction. Hydrogen bonds involving oxygen or nitrogen donor/acceptor atoms such as O—H⋯O, N—H⋯O or N—H⋯N tend to be stronger. There is particular interest in the shorter O—H⋯O hydrogen bonds whose donor–acceptor (D–A) distances approach 2.45 Å as they are found to be on the border between strong and very strong interactions.⁷

Specific features of these short, strong hydrogen bonds include a sharing of electron density distribution across the hydrogen bond,⁸ in some cases this may equate to covalency in a 3-centre 2-electron system.^9,10 The interactions can have hydrogen potentials approaching that of a single well with a diminished barrier to proton transfer¹¹ with large values of hydrogen bond energy (>100 kJ mol⁻¹).^12,13 Varying proton positions are often found, either bonded to the original donor (a neutral hydrogen bond), transferred to the acceptor (salt formation) or, in some cases, the donor and acceptor atoms compete for the hydrogen atom such that it may be centred¹⁴ (salt-cocrystal continuum) or appear on the edge of proton transfer.¹⁵ Where proton transfer occurs, these are charge assisted hydrogen bonds (CAHBs) and are even stronger due to the additional charges increasing the electrostatic component of the interaction.¹⁶ Protonation state is particularly important in the regulation of pharmaceutical cocrystals, determining whether the final form is an approved drug or a new chemical species.¹⁷ ΔpK_a of the interacting components (pK_base − pK_acid) can be used to predict salt formation by proton transfer, expected if ΔpK_a > 2 or 3.¹⁸

Temperature dependent proton migration behaviour can be observed across short O—H⋯O hydrogen bonds^19–21 when formed in organic molecular crystals between acidic and basic molecular components.²² The number of systems identified with this behaviour are still relatively few but the proton hopping can have useful applications in ferroelectrics²³ or colour changing²⁴ materials. This behaviour is only observed for systems where the D–A distance is less than 2.45 Å whilst a less negative ΔpK_a (i.e. approaching zero) between components may result in larger proton shifts.¹⁹ Combined inelastic neutron scattering and molecular dynamics simulations also indicate migration across short O—H⋯O hydrogen bonds to be a result of changes in hydrogen atom (H-atom) potentials induced by thermal fluctuations of the molecular environment.²⁵

In this work, we explore short O—H⋯O/O⁺—H⋯O⁻ hydrogen bonds formed in a selection of known substituted urea organic acid cocrystals and salts with O⋯O D–A distances in the region of 2.45 Å, where their character and temperature dependent behaviour are varied (Fig. 1).^19,20 The short hydrogen bonds in 1 (N,N′-dimethylurea oxalic acid 1 : 1) and 3 (N,N-dimethylurea 3,5-dinitrobenzoic acid 2 : 2) are neutral whilst in 2 (N,N-dimethylurea 2,4-dinitrobenzoate 1 : 1) the hydrogen bond is charge assisted. The proton position is known to be sensitive to temperature in the short hydrogen bonds of 1 and 2, gradually shifting across the hydrogen bond as a function of temperature (called proton migration).

Fig. 1 The neutral O—H⋯O and charge assisted O⁺—H⋯O⁻ short hydrogen bonds formed in 1N,N′-dimethylurea oxalic acid 2 : 1, 2N,N-dimethlyurea 2,4-dinitrobenzoate 1 : 1 and 3N,N-dimethlyurea 3,5-dinitrobenzoic acid 2 : 2 (two symmetry independent hydrogen bonded dimers form in the asymmetric unit). The short hydrogen bond proton positions in 1 and 2 are known to be sensitive to temperature, exhibiting temperature-dependent proton migration (highlighted by purple circle).

Through charge density analysis,²⁶ we examine the experimental electron density distribution (EDD) interpreted using Bader's quantum theory of atoms in molecules (QTAIM)^27,28 from single crystal samples of 1–3. The experimental EDDs are obtained following a multipolar crystal structure refinement²⁹ against high-resolution single crystal synchrotron X-ray diffraction data. Synchrotron X-ray facilities are less routinely used to obtain the EDD despite offering high X-ray flux and tuneable wavelengths.³⁰ The experimental EDD analysis is combined with first principle calculations of properties including the reduced density gradient from non-covalent interaction analysis,³¹ the electrostatic potential,^32–34 PIXEL interaction energies^35,36 and H-atom potentials.^37–39 A combined analysis approach is commonly used to assess interaction characteristics in terms of material properties.^40–42 Here it allows a complete characterisation of the short hydrogen bonds (Fig. 2) so that specific features that determine their properties, such as their differing migration behaviour, can be established.

Fig. 2 The key components of the combined approach used to study the short O—H⋯O hydrogen bonds of systems 1–3 including topological analysis and first principle calculations. The start point is the aspherical atom model obtained from the multipolar refinement.

As this study focuses on hydrogen bonds, additional methods (experimental and computational) are used to improve on the H-atom parameters determined from the X-ray data; typically the X—H distances are underestimated whilst only isotropic displacement parameters for the H-atoms can be refined. Obtaining the best models for the H-atoms is important in charge density analysis where any deficiencies will significantly affect the derived properties involving these atoms (and thus hydrogen bond characterisation).^43–45 The technique of neutron diffraction allows the most accurate determination of H-atom parameters and has been performed for 2. In the case of 1 and 3, the larger single crystals required for neutron diffraction could not be grown and so calculation-based alternatives were implemented on the X-ray structures. Hirshfeld atom refinement (HAR)⁴⁶ with NoSpherA2,⁴⁷ was performed to obtain accurate X—H distances. The SHADE3 server⁴⁸ is a recommended method to obtain ADPs for charge density analysis.⁴⁹ However, in this study, HAR was implemented to obtain the H-atom ADP model, giving the better fit to the experimental data.

Experimental

The details of the experiments are described in the ESI† and include sample preparation (ESI† 1), diffraction data collection and processing (ESI† 2), crystal structure solution and refinement (ESI† 3) and procedures for the charge density analysis (ESI† 3.4, 4). The implementation of the ab initio first principles computational methods is also described (ESI† 5).

Structure optimisation

Within this study, crystal structures have been optimised either by multipolar refinement for topological analysis of the experimental electron density (ESI† 3.4) or for first principles calculations to obtain complementary properties (ESI† 5). These complementary properties include reduced density gradients from non-covalent interaction (NCI) plotting and analysis³¹ (ESI† 5.3), electrostatic potentials (ESP)^32–34 (ESI† 5.4), CLP-PIXEL dimer interaction energies^35,36 (ESI† 5.5) and potentials for H-atom motion (ESI† 5.6). The first principle energetics calculations have been performed on the hydrogen bonded dimers of 1–3 as extracted from the crystal structure following multipolar refinement, ‘as in crystal’, and following geometry optimisation using Gaussian09 (ref. 50) code, B3LYP functional⁵¹ and 6-31G+(d,p) basis set (ESI† section 5.1).⁵² Calculations are also extended to account for the effects of crystal packing, an important consideration in studies of proton transfer,^53,54 evaluated using a pseudo-Ewald embedding calculation⁵⁵ and giving a ‘cluster’ structure (ESI† section 5.2). These are discussed in the section describing crystal packing effects.

A gas-phase calculation does not fully account for the deformation of the electron density of the molecule due to the crystal environment, and therefore the conclusions from the modelling must be treated with caution. On the other hand, the use of molecular calculations allows us to use a hybrid functional, a level of theory which is not routinely accessible to solid-state calculations, and it has recently been shown that such hybrid functionals provide a more accurate description of hydrogen bonding.⁵³ A particular issue is that B3LYP only accounts partially for dispersive interactions. This would be a particular problem for clusters with π-stacking. We partly mitigated this by keeping the molecule in the crystal structure geometry, as described by our ‘cluster’ structures. We also note that we may expect some errors due to the small 6-31G+(d,p) basis set used. However, this was selected as a balance between cost and accuracy, particularly for the larger cluster models.

Results

Structural properties

The crystal structures of 1 (N,N′-dimethylurea oxalic acid 2 : 1), 2 (N,N-dimethylurea 2,4-dinitrobenzoate 1 : 1) and 3 (N,N-dimethylurea 3,5-dinitrobenzoic acid 2 : 2) have a common structural motif of a R²₂(8) carboxylic acid : amide hydrogen bonded dimer formed between substituted urea and organic acid molecular units (Fig. 3).^19,20 The focus of this study is on the short O—H⋯O hydrogen bonds formed within this motif. In 1–3 the donor–acceptor (D–A) distances of the short hydrogen bonds are ca. 2.45 Å (Table 1); the hydrogen bond in 3d1 is slightly longer at 2.46 Å. The O—H⋯O hydrogen bonds are neutral in the case of 1 and 3; the proton resides on the acid group. The O⁺—H⋯O⁻ hydrogen bond formed in 2 is the only charge assisted interaction of this set; the acidic proton has fully transferred to the N,N-dimethylurea molecule. Two symmetry independent hydrogen bond dimers form the asymmetric unit of 3.

Fig. 3 The carboxylic acid : amide R²₂(8) hydrogen bonded dimers in the multipole refined crystal structures with respective H-atom models implemented of (a) 1, (b) 2, (c) 3 dimer 1 (3d1) and (d) 3 dimer 2 (3d2) at 100 K. The oxalic acid molecule in 1 occupies an inversion centre through the central C—C bond and thus the symmetry related portion is shown for completeness.

Table 1 Short O—H⋯O hydrogen bond parameters from the final multipolar refined crystal structure with corresponding HAR or neutron H-atom positions for 1–3

System	dO⋯O (Å)	dO—H (Å)	dH⋯A (Å)	⁁DHA (°)
1	2.4466(3)	1.179(9)	1.280(9)	168(1)
2	2.4442(7)	1.16(1)	1.31(1)	170(2)
3	d1 2.4602(6)	1.11(1)	1.36(1)	170(1)
	d2 2.4436(6)	1.15(1)	1.31(1)	169(1)

---

Within the short hydrogen bonds of 1–3, the O—H bond distances are elongated (from 1.015 Å (ref. 56)) and this feature is more pronounced where the donor–acceptor distances are shorter, in agreement with that found elsewhere.⁵⁷ The short hydrogen bond H-atom (Fig. 3) also has an elongated anisotropic thermal parameter in the case of 1 and 3d1 indicating increased thermal motion of this atom along the hydrogen bond.

Electron density distribution ρ(r)

The electron density distributions (EDD) of the short hydrogen bonds are examined in the topological parameters of the experimentally obtained EDD together with isosurfaces of the calculated reduced density gradient following non covalent interaction (NCI) analysis³¹ (Fig. 4a). The NCI analysis³¹ can be used to directly compare the reduced density gradient (Laplacian ∇²(r_BCP)) observed during the multipole reduction of experimental charge densities to those of a comparable molecular adduct in vacuum.⁵⁸ Both analyses are performed at bond critical points of the short hydrogen bonds (BCP; a saddle point in the electron density ρ(r) where its gradient is ∼0 being a minimum in the direction of the bond and a maximum perpendicular to it).

Fig. 4 (a) NCI isosurface plot for the hydrogen bonded dimer in 1. The figure was obtained using a reduced density gradient of 2 a.u. and the blue–green–red values ranging from (−0.05 to +0.05 a.u). Green/red arrows point to regions of vdW interactions, while yellow arrows point to mild repulsions and blue arrows indicate hydrogen bonds. (b) The electrostatic potential surface in 1 calculated for the hydrogen bonded dimer from the optimised or ‘as in crystal’ structure. Red indicates positive regions, blue negative regions and green van der Waals regions. Equivalent plots (a) and (b) for 2 and 3 are located in the ESI,† Fig. S7 and S9.

The topological parameters (Table 2) show that the value of ρ(r)_BCP at the O—H/H⋯O BCP decreases with increasing atomic separation (Fig. 5) whilst only a negative value of the Laplacian ∇²(r_BCP) for the H⋯O interaction is found in 2 (−1.71(11) e Å⁻⁵) where full proton transfer has occurred across the hydrogen bond. In agreement with previous findings,⁹ as the electron density reflects proton position in 1–3, it therefore appears sensitive to the extent of static transfer of the H-atom across the hydrogen bond. The large values of the electron density in the H⋯O region of the short hydrogen bonds given by the topological parameter ρ(r)_BCP (0.79 to 0.9 e Å⁻³) combined with the ring shaped deep blue/covalent bond of the respective calculated NCI isosurface (Fig. 4b and S7†), which exceeds the mapping limits, suggest the short O—H⋯O hydrogen bonding interactions to be strong⁵⁹ and approaching a 3-centre, 2-electron system (ρ(r)_BCP ≈ 1 e ų).¹⁰ These observations agree with the short hydrogen bond O⋯O separation distance of ca. 2.45 Å often being associated with stronger hydrogen bonds. As an indicator of strength, the values of ρ(r_BCP) correspond to ca. 30% of that determined in the covalent bond of O—H in H₂O (ρ(r_BCP) 2.63 e Å⁻³).^27,60 In weak to moderate strength hydrogen bonds, typical values are significantly lower on the order of 0.05 to 0.23 e Å⁻³ (ref. 57) and have lighter blue to green calculated NCI isosurfaces, as can be seen for the complementary N—H⋯O interaction in the dimer motifs (Fig. 4a). This analysis shows that both the NCI isosurface and ρ(r)_BCP are useful in assessing hydrogen bond strength whilst ρ(r)_BCP gives a better indicator of position (and can be used as a scaling tool) within the salt–cocrystal continuum.

Table 2 Selected topological parameters of the short O—H⋯O/O⁺—H⋯O⁻ hydrogen bond H⋯O BCPs for each system. The energy densities (kinetic G_BCP, potential V_BCP and total H_BCP) are in units of Hartrees Å⁻³. See ESI† (Table S5) for full list of topological parameters including those of the O—H bond

System	Interaction	R ij (Å)	ρ(rBCP) (e Å^−3)	∇^2ρ(rBCP) (e Å^−5)	G(rBCP)	V(rBCP)	H(rBCP)	|V|/G
1	O1—H1⋯O3	1.2806	0.91(4)	3.0(2)	0.82	−1.43	−0.61	1.74
2	O1—H1⋯O2	1.2953	0.88(3)	−1.7(1)	0.58	−1.27	−0.7	2.19
3d1	O1—H1⋯O13	1.3609	0.78(4)	0.95(1)	0.58	−1.09	−0.51	1.87
3d2	O9—H9⋯O14	1.3049	0.88(4)	1.89(2)	0.74	−1.35	−0.61	1.82

---

Fig. 5 ρ(r_BCP) for the O—H and H⋯O BCPs of the short hydrogen bonds in 1–3 as a function of atom–atom separation.

A comparison of the values of the topological parameters of the Laplacian ∇²(r_BCP) and energy densities (G_BCP, V_BCP, H_BCP) across the systems (Table 2) identify a difference in 2. For this short hydrogen bond, the ∇²(r_BCP) is negative (−1.71(11) e Å⁻⁵) indicating a concentration of charge in the H⋯O region and, when found alongside a negative total energy density (H_BCP), is characteristic of a ‘very strong’ hydrogen bond with a shared-shell covalent character. In contrast, the values of the Laplacian are positive (∇²(r_BCP) > 0) for 1 and 3 corresponding to a depletion of charge within the H⋯O atom–atom region and is characteristic of a closed shell electrostatic interaction or ionic type bonding.⁶¹ The ratio of the Lagrangian kinetic and potential energy densities (|V_BCP|/G_BCP),⁶² which can be used to further define the hydrogen bonding interactions at the short hydrogen bond BCPs (Table 2), additionally suggest a covalent character for 2, where the ratio is greater than 2 and in the region attributed to covalent bonds. In the case of systems 1 and 3, a partial covalency is suggested by the |V_BCP|/G_BCP ratio, where it lies between 1 and 2 in the intermediate region between electrostatic and covalent. Systems 1–3 therefore add to the cases of short hydrogen bonds that are known to exist in this in-between region, with both electrostatic and covalent contributions.⁶³ Short O—H⋯O strong hydrogen bonds formed in related substituted urea organic acid co-crystals and salts are found to have either covalent (urea phosphoric acid 1 : 1 (UPA))⁶⁴ or partial covalent character (urea oxalic acid)⁶⁵ from similar experimental electron density analysis. Of these, UPA⁶⁴ is a system which also exhibits temperature dependent proton migration.²¹ Across the systems exhibiting this type of behaviour (UPA, systems 1 and 2), there does not appear to be any correlation to the electron density distribution across the short hydrogen bond and the migration behaviour observed. Therefore, despite the electron density being sensitive to the extent of static proton transfer, it does not appear to be useful in predicting the likelihood of further proton transfer events as a function of temperature.

When analysing the topological parameters at the H⋯O BCP, the O—H bond directed quadrupole (Q₀) was found to be important in influencing the value of the Laplacian ∇²(r_BCP). A larger Q₀ elongates the charge concentration of the H-atom along the interaction increasing likelihood of overlap with the H⋯O BCP and therefore affecting its value. This effect is made clear in the 2D Laplacian maps (Fig. 6a–d). When trialling different ADP models for 1, they were found to ‘tune’ Q₀ (a less elongated ADP increased the value of Q₀). The importance of obtaining the ‘best’ ADP model for H-atoms in short hydrogen bonds is therefore highlighted, as extracted parameters at the BCPs such as the Laplacian can vary significantly.

Fig. 6 The short hydrogen bonds of 1–3 are visualised in the plane of the O—H⋯O SSHBs in each system. In (a–d), in Laplacian maps showing regions of charge concentration (red contours) or depletion (blue dashed lines) over the interacting groups. The (3, −1) bond critical points are shown as red spheres and importantly show the character of the H⋯O BCP whether located in a charge concentrated or depleted area. Contour lines are at 0.05 e Å⁻⁵. The refined value of the bond directed quadrupole (Q₀) of the hydrogen bond H-atom is additionally included. In (e–h) static deformation densities. Blue and red colours indicate positive and negative regions of electron density, respectively and reveal lone pairs and bonding density. Contours are drawn at intervals of 0.05 e Å⁻³. The populations for the H-atom, modelled and refined as a monopole, are additionally included in e–h.

Atomic charges

2D deformation density maps obtained from the experimental electron density distribution (EDD) and showing bonding and lone pair features (Fig. 6e–h) indicate that there is no evidence of charge transfer across the short hydrogen bonds. A connection of the O—H deformation lobes with the acceptor oxygen would be observed in such a case.⁶⁶ This is found to occur across an O—H⋯O hydrogen bond in an identical dimer motif formed in a 1 : 1 urea oxalic acid system.⁶⁵ The nature of the short hydrogen bond is again shown to be different for 2 where, in contrast to 1 and 3, the H-atom (H1) occupies a negative (red) contour region with a correspondingly low population (0.68(1)e), modelled and refined as a monopole, and has a large, positive net charge (+0.32(1)e) (Table S7†). This evidence points to a deshielded H-atom⁵⁷ significantly depleted of electrons. The large, net positive charge is as found for H-atoms in other charge assisted [O⋯H⋯O]⁻ hydrogen bonds⁶⁷ and is evidence of being involved in the formation of a stronger hydrogen bond, as corroborated by the topological parameters. In 1 and 3 there is no evidence of deshielding of the hydrogen bond H-atom in the deformation density maps, the H-atom position is overlapped by the bonding density, whilst its monopole populations (0.8e) and net atomic charges (ca. +0.2e) are correspondingly larger and less positive, respectively.

Across all the short hydrogen bonds, the low monopole population for the H-atom (<1.0e) indicates an extent of charge sharing to the oxygen atom in the O—H bond. The electrostatic potentials for 1–3, calculated using Gaussian09 (ref. 50) (ESI† 5.4), for the ‘as in crystal’ and optimised hydrogen bonded dimer (Fig. 4b and S9†), visually reflect this reduced population where the H-atom has a positive (red region) ESP in a highly polarised O—H bond (the oxygen opposite the hydrogen in the bond is correspondingly dark blue).

Gaussian09 (ref. 50) calculates the electron density as present in Gaussian orbitals and hence, the ESP diagrams enable the location and contribution of the lone pair (lp) of electrons to the different interactions. As expected, in each case, one of the lp of electrons on the carboxylic acid O(H) is delocalised into the O—H bond while it is easy to spot the location of the other O(H) lp as an increase in electron density, i.e. deeper blue in the ESP (Fig. 4b). For the optimised dimer ESP, both oxygen acceptor lps can be located and the lp involved in the hydrogen bond formation sits opposite the positive local potential (red region) on the H-atom. This complementarity generates the idea of a purely electrostatic type interaction for the optimised dimer. However, in the case of the ‘as in crystal’ ESP only one of the lps on the urea O can be located, while the other appears to be diffused over the short hydrogen bond. In contrast to the deformation density maps, the ESP therefore suggest an extent of charge transfer across the short O—H⋯O hydrogen bonds. This effect is reduced in 3d1, potentially the longer donor–acceptor distance for the short hydrogen bond inhibits the extent of charge redistribution/sharing. This remains the case also in ‘cluster’ structures (Fig. S10 and S11†), but where there is a general delocalisation of electron density over the molecules. The ‘as in crystal’ dimer is as extracted from the crystal structure and so it appears that the sharing of ESP across the components is an element of the crystal packing. The comparison of the optimised versus ‘as in crystal’ ESP therefore gives a clear representation of the crystal packing effects on molecular charge transfer.

Pixel interaction energies

Total (Pixel) interaction energies⁶⁸ have been calculated using the CLP-PIXEL method^36,69,70 (ESI† 5.5) for each ‘as in crystal’ hydrogen bonded dimer including their decomposition into contributing terms (Table 3). They define the energy of the dimer interaction allowing further distinction in terms of hydrogen bond characteristics.^71,72

Table 3 Decomposition of the total dimer interaction energy (E_TOTAL) into contributing components from CLP-PIXEL calculations for the hydrogen bond dimers in 1–3. Energies are in kJ mol⁻¹

System	E COULOMBIC	E POLARISATION	E DISPERSION	E REPULSION	E TOTAL
1	−227.3	−158.2	−25.7	350.0	−61.1
2	−545.8	−209.1	−24.7	287.0	−492.7
3d1	−188.8	−123.3	−23.7	258.1	−77.8
3d2	−214.2	−150.6	−24.3	319.8	−69.3

---

As by the other analyses, the dimer in system 2 is highlighted as being different having a total dimer interaction energy (E_TOTAL) of −493 kJ mol⁻¹versus those for 1 and 3 that are in the region of −60 to −80 kJ mol⁻¹. The E_TOTAL is significantly larger for 2; to compare, an O—H covalent bond has a bond energy of ca. 465 kJ mol⁻¹.⁷³ However, this is not unusual for interactions between charged species, thought to occupy the strongly stabilising, ionic coulombic contact zone range.⁷⁴ Using the Hibbert and Emsley classification,^75,76 the hydrogen bonding interactions forming the dimers in 1 and 3 are considered strong (E_TOT of −60 to −80 kJ mol⁻¹) whilst they are very strong in 2 (E_TOTAL is −493 kJ mol⁻¹). The co-crystal/salt nature of 1–3 is further reflected in the dimer interaction energies where the salt has a significantly larger total interaction energy.

The decomposition of E_TOTAL into its contributing terms (attractive E_COULOMBIC, E_POLARISATION, E_DISPERSION and repulsive E_REPULSION) shows that, for all dimers, the coloumbic component contributes most to the stabilisation of the dimer, where E_COULOMBIC is the largest negative energy. E_COULOMBIC defines the electrostatic interaction in terms of Coulomb interactions between charges at points ‘pixels’ in space and is not always found to be the dominating term for hydrogen bonds.⁷⁷ The large repulsive contribution (E_REPULSION) could be due to unfavourable secondary diagonal interactions⁷⁸ in the dimer motif between the diagonally opposite carbonyl groups (in the urea and acid) oxygen lone pairs or H-atoms. These repulsive interactions are shown to exist by the yellow NCI isosurfaces (Fig. 4a and S7†) in the centre of the hydrogen bonded ring.

H-atom potentials

Potential energy surface calculations were performed for the hydrogen bonded dimer with respect to the proton position between the interacting O⋯O atoms in the short hydrogen bonds for 1–3. This allowed H-atom potentials to be generated. The potentials were produced using Gaussian09 (ref. 50) code and B3LYP functional and 6-31G+(d,p) basis set (ESI† section 5.6) for the ‘as in crystal’ dimer and following its geometry optimisation.

A clear difference is observable between the potentials calculated for the ‘optimised’ versus the ‘as in crystal’ dimers (Fig. S12†). Both have the shape of an asymmetric single well potential however, a very small barrier is present in the ‘optimised’ potential which is absent for the ‘as in crystal’ dimer. Optimisation of the dimer fragment also results in a shift in the minimum position, for 1–3, towards the acid group whilst the O⋯O distance is altered by approximately 0.005 Å. In general, DFT calculations are known to underestimate the strength of short hydrogen bonds resulting in shorter calculated O—H distances and longer calculated O⋯O distances.²⁵ This may explain the presence of the ‘shoulder’ like barrier in the optimised potential, where the donor–acceptor atoms are slightly further apart such that well overlap is reduced and barrier height is increased.

The effects of crystal packing are also absent in the optimised in vacuo dimer. While the ‘as in crystal’ model does not fully reproduce the bulk crystal environment, it is well known that molecules undergo conformational adjustments on crystal packing⁷⁹ and therefore the difference between the gas-phase ‘as in crystal’ and ‘optimised’ models is not surprising. It is important, however, to consider these effects as they can significantly influence structure,⁴⁰ especially proton transfer state.^80–82 Furthermore, NCI analysis shows that there are significant weak and strong interactions in the hydrogen bond dimer region (Fig. 7 and S8†) and so should be considered. Therefore, for a truer picture of the H-atom potentials, only the ‘as in crystal’ potentials are considered in the following analysis.

Fig. 7 NCI isosurface plot for the cluster models of (a) 2 and (b) 3d1. The figure was obtained using a reduced density gradient of 2 a.u. and the blue–green–red values ranging from (−0.05 to +0.05 a.u). Green arrows point to regions of vdW interactions, while yellow and blue arrows point to mild repulsions in the middle of hydrogen bonded rings and strong hydrogen bonds, respectively. Equivalent plots for 1 and 3d2 are located in the ESI,† Fig. S8.

‘as in crystal’ potentials. The asymmetric ‘as in crystal’ potentials (Fig. S12†) indicate unsymmetrical short hydrogen bonds and therefore an asymmetric molecular environment,¹⁰ as might be expected for a hetero hydrogen bonded dimer. In the absence of a barrier to a second site, the potentials are ‘slide like’ with no clear second minimum. Such a behaviour is also observed from experimental XPS³⁷ and molecular dynamics studies⁸³ of 3,5-pyridinedicarboxylic acid and urea phosphoric acid²⁵ reported in the literature. We also see this for related materials containing the same R²₂(8) carboxylic acid : amide hydrogen bonded dimer (Fig. S13†).

Following on from the other analyses, differences are seen in the case of 2 where the ‘slide’ is significantly shallower such that the potential is closest to resembling a single well flat-bottomed potential; one site is only just favoured energetically over the other. This flatter potential is characteristic of a hydrogen bond with added strength,¹¹ as seen from the electron density analysis and Pixel interaction energies. Furthermore, the flatter potential may explain the greater extent of migration observed for this system, as a function of temperature, of the H-atom across the hydrogen bond from the donor to acceptor site. The flatter potential favours the occupation of more H-atom sites along the hydrogen bond. In contrast, for other migration material 1, the H-atom potential is less flat and a reduced extent of migration is observed; the H-atom undergoes a small shift migrating back towards the acid as a function of temperature rather than across the hydrogen bond.

The donor–acceptor distances for the short hydrogen bonds are similar (ca. 2.45 Å) and therefore not likely to be the cause for the flatter potential in 2. Instead, it could be the added charge assisted component of this interaction, seen in this study to affect other properties of the hydrogen bond. pK_a values are known to determine the energy barrier height in potentials for moving a H-atom between donor–acceptor wells.^19,84 For systems 1 to 3, this does not seem to be the case where the ΔpK_a of components are similar (ΔpK_a1: −1.94, 2: −1.74 and 3: −1.81) whilst the ‘slide’ in the potentials occupy different regions for 1 and 3 (30–40 kJ mol⁻¹) versus in 2 (ca. 10 kJ mol⁻¹). One study also formed a link between the extent of temperature dependent migration and ΔpK_a.¹⁹ As was found by Jones et al. (2012), the least negative (closer to zero) ΔpK_a in 2 results in the greatest extent of migration observed. However, the trend is not reflected in 3, where the short hydrogen bonds show no evidence of temperature dependent proton migration yet have a less negative ΔpK_a than 1, which is a migration material. If taken from the literature, ΔpK_a values should be used with caution as they are often determined for the solution state and are affected by intermolecular interactions in a crystal.¹⁴ Here they are obtained from the literature^19,85 and as such are used only as a guide.

Crystal packing effects on potentials. The ‘as in crystal’ H-atom potentials are extended to consider nearest neighbours giving ‘cluster’ calculated potentials. Specifically, nearest neighbours are included in differing intermolecular pairs to determine their individual effects on the hydrogen bond dynamics.

Inspection of the extended H-atom potentials (Fig. 8) indicates that, of the intermolecular interactions, the expanded hydrogen bonding contacts cause the most significant changes to the potentials. This might be expected where hydrogen bonding contacts are stronger than π-type interactions and more likely to perturb the local environment. In the most extreme case of 2, the presence of the hydrogen bonding to the neighbouring dimer unit leads to the inversion of the minima from the acceptor back to the donor, thereby identifying the proton transfer as a feature of these interactions in the crystal environment. In the case of 1 and 3d1, the extended hydrogen bonding leads to a considerable decrease in the energies of the acceptor and hence in the slope of the PES curve. The local interactions are shown to be weaker in 3d2 by the NCI analysis (ESI† Fig. S8), and there is correspondingly little change between the potential of the dimer alone (asu) versus the ‘cluster’ stack. The π-stacking effect leads to an increase in the slope of the potential but in most cases is off-set by the strength of the extended hydrogen bonding with the neighbouring molecules.

Fig. 8 H-Atom potentials in vacuo of various ‘cluster’ models of (a) 1, (b) 2 and (c) 3, where the proton on the central dimer highlighted by an ellipse (yellow) is moved by 0.02 Å increments between the donor and acceptor oxygens. The green box highlights the presence of extra H-bonding in the expanded H-bonding model. The H-atom position as found in the crystal structure is indicated by black arrows (H_POS).

As is shown for the extreme case of 2, by calculating the extended hydrogen bonding cluster potentials for charge assisted hydrogen bonds, if the crystal packing significantly shifts the minima, then this approach may be able to predict whether a material will undergo a complete donor–acceptor temperature dependent proton migration. System 2 could be used as a benchmark for this type of behaviour. This is a hypothesis and other related materials would need to be studied to confirm this.

Temperature effects on potentials. The H-atom potential calculations for the ‘as in crystal’ dimers are extended to consider temperature effects for 1 and 2 in order to shed light on their known temperature dependent migration behaviour.^19,20 Small changes start to occur in the PES at 200 K in 1 and from 250 K in 2, the ‘slide’ becomes slightly shallower whilst the potential broadens (Fig. S14†). The changes also coincide with a lengthening of the D–A distances of the short hydrogen bonds at these temperatures (Table S9†) and could be the cause of the broadening effect. However, a barrier might be more likely to develop as the donor–acceptor wells separate. Instead, a broadening in the potentials with temperature can be due to thermal fluctuations in the molecular environment, as seen in related migration material of urea phosphoric acid (UPA) 1 : 1.²⁵ Furthermore, from molecular dynamics simulations and inelastic neutron scattering^83,86,87 the mechanism for migration in UPA is attributed to these thermal fluctuations, in particular, low frequency vibrations modes are known to contribute most to the migration occurring.^25,83

Due to the similarity in the behaviours of the H-atom potentials as a function of temperature for 1, 2versus UPA, and because crystal environment appears to have a significant effect on the H-atom minima positions in the ‘cluster’ structures of 1 and 2, we propose that their temperature dependent migration behaviour could be due to the same mechanism as UPA, caused by thermal fluctuations in the molecular environment. Furthermore, Boltzmann distributions for 1 and 2, calculated as a function of temperature, favour this explanation showing a broadening of the population distribution with temperature (greatest in 2) and not an increasingly populated first excited state (Fig. S15†). The increased population of higher energy levels is the alternative explanation advanced for migration behaviour, as found in other N⋯H⋯O short hydrogen bonded systems.^86,88 To confirm the nature of the temperature dependent migration behaviour in 1 and 2, these systems would benefit from further molecular dynamics simulations and inelastic neutron scattering studies. This should be the focus of future work for these systems.

Conclusions

In this work we have explored short O—H⋯O and O⁺—H⋯O⁻ hydrogen bonds with donor–acceptor distances in the region of 2.45 Å in substituted dimethylurea organic acid salts and co-crystals. The combined approach of experimental charge density analysis and first principles calculations has allowed a complete characterisation of hydrogen bonds within this donor–acceptor distance region.

In the first part we analysed properties of the experimentally determined charge density alongside calculated properties from non-covalent interaction analysis and electrostatic potentials. These properties were also compared to the calculated total Pixel dimer interaction energies. From this, the short hydrogen bonds in N,N′-dimethylurea oxalic acid 2 : 1 (1) and N,N-dimethylurea 3,5-dinitrobenzoic acid 2 : 2 (3) are characterised as strong with a more electrostatic contribution whilst in N,N-dimethylurea 2,4-dinitrobenzoate 1 : 1 (2) the short hydrogen bond is characterised as very strong with a significant covalent contribution. This donor–acceptor distance region therefore remains a zone of fluctuation in terms of the nature of the hydrogen bond. The identification of O—H⋯O hydrogen bonds with a covalent nature is still uncommon. However, a charge assisted component is seen to favour covalency in this interaction and proves to be the determining factor for this set of substituted urea organic acid molecular complexes, in terms of the strength and character of the short hydrogen bonds. Therefore, if salt formation is targeted, this presents a route to access hydrogen bonds with a covalent component.

The relation found between the atom-atom separation or protonation state and the topological properties of the charge density highlights how the electron density in these types of interactions remains highly sensitive to the static proton transfer process. It is therefore a useful parameter in studies of the salt-cocrystal continuum. This has implications for both pharmaceutics and the design of crystalline functional materials. However, for these systems, the electron density distribution does not appear to be useful to predict further proton transfer events as a function of temperature.

In the second part, hydrogen bond characteristics were examined in H-atom potentials, as calculated, and as a function of the crystal packing and temperature. The extensive studies of the H-atom potentials suggest asymmetric hydrogen bonds whose energetics are different when a charge assisted component is included. Extending the potentials to consider crystal packing neighbours, identified by the non-covalent interaction analysis, confirms that proton transfer is an effect of the local environment. Here we further identify that, for these systems, the interactions forming the extended hydrogen bonding have the biggest effect on the H-atom energetics. The similarities in the changes of the potentials as a function of temperature to well-studied migration material urea phosphoric acid suggest that the reported migration of 1 and 2 is likely to also be due to fluctuations in the molecular environment and presents these systems as suitable cases for future inelastic neutron scattering and molecular dynamics studies. System 2 (N,N-dimethylurea 2,4-dinitrobenzoate 1 : 1) is highlighted as a model system, to add to the well-studied urea phosphoric acid 1 : 1 adduct, for learning about migration materials due to the significant perturbations exhibited by the short hydrogen bond.

Overall, the use of high-resolution synchrotron X-ray diffraction combined with neutron diffraction or Hirshfeld atom refinement has allowed good models of the experimental electron density to be produced in which trends and conclusions about the short hydrogen bonds could be made. Obtaining reliable H-atom models was shown to be important in these studies of short hydrogen bonds, where the H-atom exhibits atypical behaviour as a result of the strength of the hydrogen bond it is involved in. The Hirshfeld atom refinement method, in the absence of neutron data, enabled the study to be extended to additional systems for a better picture of hydrogen bonding characteristics.

Author contributions

Contributions are defined using CRediT for standardised contribution descriptions. We highlight joint first co-authorship of L. K. Saunders (LKS) and A. R. Pallipurath (ARP) having significant shared data curation and investigation roles with LKS taking lead in conceptualisation of the study and writing the original draft and ARP for writing – review and editing. Credit LKS and ARP for formal analysis, methodology, validation and visualisation. Credit M. J. Gutmann for neutron diffraction data curation and formal analysis. Credit H. Nowell for X-ray diffraction data curation, N. Zhang for software assistance with Gaussian09 and CLP-Pixel calculations and D. R. Allan for providing resources and supervision.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

We thank Florian Kleemiss and Horst Puschmann from the HAR/Olex2 team for useful discussions on running the NoSphereA2 refinement software. We thank Emil Darmsgaard and Lennard Krause of the Overgaard Group at Aarhus University for sharing practical knowledge on performing charge density refinements in XD2006. We thank Diamond Light Source, funded by the Wellcome Trust and the Science and Technology Facilities Council, for access to beamline I19-1⁸⁹ (NR18193) for data for samples 1–3. We thank the ISIS Pulsed Neutron and Muon Source, funded by UKRI and the Science and Technology Facilities Council, for access to beamline SXD⁹⁰ for data on sample 2. The computational work was undertaken on ARC3, part of the high performance computing facilities at the University of Leeds, UK. ARP acknowledges the financial support from the Future Continuous Manufacturing and Advanced Crystallization (CMAC) Hub (EPSRC Grant EP/P006965/1). The authors also thank Prof Sven Schroeder and Dr Jonathan Skelton for fruitful discussions about the potential energy surface calculations.

References

1. R. E. Hubbard and M. Kamran Haider, in Encyclopedia of Life Sciences (ELS), John Wiley & Sons, Ltd, Chichester, 2010,  DOI:10.1002/9780470015902.a0003011.pub2.
2. Q. Huang, W. Li, Z. Mao, L. Qu, Y. Li, H. Zhang, T. Yu, Z. Yang, J. Zhao, Y. Zhang, M. P. Aldred and Z. Chi, Nat. Commun., 2019, 10, 3074.
3. L. Schaeffer, in The Practice of Medicinal Chemistry, ed. C. G. Wermuth, D. Aldous, P. Raboisson and D. Rognan, Academic Press, San Diego, 4th edn, 2008, pp. 359–378,  DOI:10.1016/B978-0-12-417205-0.00014-6.
4. J. Trylska, P. Grochowski and J. A. McCammon, Protein Sci., 2004, 13, 513–528.
5. F. Ibba, G. Pupo, A. L. Thompson, J. M. Brown, T. D. W. Claridge and V. Gouverneur, J. Am. Chem. Soc., 2020, 142, 19731–19744.
6. M. C. Etter, Acc. Chem. Res., 1990, 23, 120–126.
7. M. V. Vener, E. O. Levina, A. A. Astakhov and V. G. Tsirelson, Chem. Phys. Lett., 2015, 638, 233–236.
8. V. R. Hathwar, R. Pal and T. N. Guru Row, Cryst. Growth Des., 2010, 10, 3306–3310.
9. I. Majerz, Mol. Phys., 2007, 105, 2305–2314.
10. K. Molčanov, C. Jelsch, E. Wenger, J. Stare, A. Ø. Madsen and B. Kojić-Prodić, CrystEngComm, 2017, 19, 3898–3901.
11. C. L. Perrin, Acc. Chem. Res., 2010, 43, 1550–1557.
12. G. Gilli and P. Gilli, J. Mol. Struct., 2000, 552, 1–15.
13. J. Overgaard and B. B. Iversen, Struct. Bonding, 2012, 146, 53.
14. M. Schmidtmann and C. C. Wilson, CrystEngComm, 2008, 10, 177–183.
15. M. Schmidtmann, L. J. Farrugia, D. S. Middlemiss, M. J. Gutmann, G. J. McIntyre and C. C. Wilson, J. Phys. Chem. A, 2009, 113, 13985–13997.
16. A. J. L. Jesus and J. S. Redinha, J. Phys. Chem. A, 2011, 115, 14069–14077.
17. Regulatory Classification of Pharmaceutical Co-Crystals, Food and Drug Administration Guidance for Industry, Silver Spring, MD, USA, 2013.
18. P. H. Stahl and C. G. Wermuth, Handbook of Pharmaceutical Salts: Properties, Selection, and Use, Wiley-VCH, 2002.
19. A. O. F. Jones, M.-H. Lemee-Cailleau, D. M. S. Martins, G. J. McIntyre, I. D. H. Oswald, C. R. Pulham, C. K. Spanswick, L. H. Thomas and C. C. Wilson, Phys. Chem. Chem. Phys., 2012, 14, 13273–13283.
20. L. K. Saunders, H. Nowell, L. E. Hatcher, H. J. Shepherd, S. J. Teat, D. R. Allan, P. R. Raithby and C. C. Wilson, CrystEngComm, 2019, 21, 5249–5260.
21. C. C. Wilson, K. Shankland and N. Shankland, Z. Kristallogr., 2001, 216, 303–306.
22. A. O. F. Jones, C. K. Leech, G. J. McIntyre, C. C. Wilson and L. H. Thomas, CrystEngComm, 2014, 16, 8177–8184.
23. S. Horiuchi, K. Kobayashi, R. Kumai and S. Ishibashi, Nat. Commun., 2017, 8, 14426.
24. C. L. Jones, C. C. Wilson and L. H. Thomas, CrystEngComm, 2014, 16, 5849–5858.
25. F. Fontaine-Vive, M. R. Johnson, G. J. Kearley, J. A. K. Howard and S. F. Parker, J. Am. Chem. Soc., 2006, 128, 2963–2969.
26. A. Krawczuk and P. Macchi, Chem. Cent. J., 2014, 8, 68.
27. R. F. W. Bader, Atoms in Molecules-A Quantum Theory, Oxford University Press, Oxford, 1994.
28. C. Matta and R. J. Boyd, in The Quantum Theory of Atoms in Molecules, ed. C. Matta and R. J Boyd, Wiley-VCH, New York, 2007, pp. 1–34,  DOI:10.1002/9783527610709.ch1.
29. N. K. Hansen and P. Coppens, Acta Crystallogr., Sect. A: Cryst. Phys., Diffr., Theor. Gen. Crystallogr., 1978, 34, 909–921.
30. A. Genoni, L. Bucinsky, N. Claiser, J. Contreras-Garcia, B. Dittrich, P. M. Dominiak, E. Espinosa, C. Gatti, P. Giannozzi, J.-M. Gillet, D. Jayatilaka, P. Macchi, A. Ø. Madsen, L. J. Massa, C. F. Matta, K. M. Merz, P. N. H. Nakashima, H. Ott, U. Ryde, K. Schwarz, M. Sierka and S. Grabowsky, Chem. – Eur. J., 2018, 24, 10881–10905.
31. J. Contreras-García, E. R. Johnson, S. Keinan, R. Chaudret, J.-P. Piquemal, D. N. Beratan and W. Yang, J. Chem. Theory Comput., 2011, 7, 625–632.
32. J. S. Murray and P. Politzer, WIREs Comput. Mol. Sci., 2017, 7, e1326.
33. M. A. Spackman, J. J. McKinnon and D. Jayatilaka, CrystEngComm, 2008, 10, 377–388.
34. I. Otero-Navas and J. M. Seminario, J. Mol. Model., 2012, 18, 91–101.
35. N. Walshe, M. Crushell, J. Karpinska, A. Erxleben and P. McArdle, Cryst. Growth Des., 2015, 15, 3235–3248.
36. A. Gavezzotti, New J. Chem., 2011, 35, 1360–1368.
37. J. S. Stevens, S. Coultas, C. Jaye, D. A. Fischer and S. L. M. Schroeder, Phys. Chem. Chem. Phys., 2020, 22, 4916–4923.
38. P. Gilli, V. Bertolasi, L. Pretto, V. Ferretti and G. Gilli, J. Am. Chem. Soc., 2004, 126, 3845–3855.
39. P. Durlak, K. Mierzwicki and Z. Latajka, J. Phys. Chem. B, 2013, 117, 5430–5440.
40. A. A. Melekhova, A. S. Novikov, K. V. Luzyanin, N. A. Bokach, G. L. Starova, V. V. Gurzhiy and V. Y. Kukushkin, Inorg. Chim. Acta, 2015, 434, 31–36.
41. A. S. Mikherdov, A. S. Novikov, M. A. Kinzhalov, V. P. Boyarskiy, G. L. Starova, A. Y. Ivanov and V. Y. Kukushkin, Inorg. Chem., 2018, 57, 3420–3433.
42. C. Kalaiarasi, M. S. Pavan and P. Kumaradhas, Acta Crystallogr., Sect. B: Struct. Sci., Cryst. Eng. Mater., 2016, 72, 775–786.
43. A. O. Madsen, H. O. Sorensen, C. Flensburg, R. F. Stewart and S. Larsen, Acta Crystallogr., Sect. A: Found. Crystallogr., 2004, 60, 550–561.
44. A. A. Hoser, P. M. Dominiak and K. Wozniak, Acta Crystallogr., Sect. A: Found. Crystallogr., 2009, 65, 300–311.
45. P. Roversi and R. Destro, Chem. Phys. Lett., 2004, 386, 472–478.
46. S. C. Capelli, H.-B. Burgi, B. Dittrich, S. Grabowsky and D. Jayatilaka, IUCrJ, 2014, 1, 361–379.
47. F. Kleemiss, O. V. Dolomanov, M. Bodensteiner, N. Peyerimhoff, L. Midgley, L. J. Bourhis, A. Genoni, L. A. Malaspina, D. Jayatilaka, J. L. Spencer, F. White, B. Grundkötter-Stock, S. Steinhauer, D. Lentz, H. Puschmann and S. Grabowsky, Chem. Sci., 2021, 12, 1675–1692.
48. A. O. Madsen and A. A. Hoser, J. Appl. Crystallogr., 2014, 47, 2100–2104.
49. L. A. Malaspina, A. A. Hoser, A. J. Edwards, M. Woińska, M. J. Turner, J. R. Price, K. Sugimoto, E. Nishibori, H.-B. Bürgi, D. Jayatilaka and S. Grabowsky, CrystEngComm, 2020, 22, 4778–4789.
50. M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, G. Scalmani, V. Barone, B. Mennucci, G. A. Petersson, H. Nakatsuji, M. Caricato, X. Li, H. P. Hratchian, A. F. Izmaylov, J. Bloino, G. Zheng, J. L. Sonnenberg, M. Hada, M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, T. Vreven, J. A. Montgomery Jr., J. E. Peralta, F. Ogliaro, M. Bearpark, J. J. Heyd, E. Brothers, K. N. Kudin, V. N. Staroverov, R. Kobayashi, J. Normand, K. Raghavachari, A. Rendell, J. C. Burant, S. S. Iyengar, J. Tomasi, M. Cossi, N. Rega, J. M. Millam, M. Klene, J. E. Knox, J. B. Cross, V. Bakken, C. Adamo, J. Jaramillo, R. Gomperts, R. E. Stratmann, O. Yazyev, A. J. Austin, R. Cammi, C. Pomelli, J. W. Ochterski, R. L. Martin, K. Morokuma, V. G. Zakrzewski, G. A. Voth, P. Salvador, J. J. Dannenberg, S. Dapprich, A. D. Daniels, Ö. Farkas, J. B. Foresman, J. V. Ortiz, J. Cioslowski and D. J. Fox, Gaussian09, Gaussian Inc., Wallingford, USA, 2009.
51. S. Grimme, J. Comput. Chem., 2006, 27, 1787–1799.
52. R. Ditchfield, W. J. Hehre and J. A. Pople, J. Chem. Phys., 1971, 54, 724–728.
53. L. M. LeBlanc, S. G. Dale, C. R. Taylor, A. D. Becke, G. M. Day and E. R. Johnson, Angew. Chem., Int. Ed., 2018, 57, 14906–14910.
54. S. L. Childs, G. P. Stahly and A. Park, Mol. Pharmaceutics, 2007, 4, 323–338.
55. M. Rivera, M. Dommett and R. Crespo-Otero, J. Chem. Theory Comput., 2019, 15, 2504–2516.
56. F. H. Allen, O. Kennard, D. G. Watson, L. Brammer, A. G. Orpen and R. Taylor, J. Chem. Soc., Perkin Trans. 2, 1987, S1–S19,  10.1039/P298700000S1.
57. G. K. H. Madsen, B. B. Iversen, F. K. Larsen, M. Kapon, G. M. Reisner and F. H. Herbstein, J. Am. Chem. Soc., 1998, 120, 10040–10045.
58. S. Scheiner, J. Comput. Chem., 2018, 39, 500–510.
59. J. S. Lomas, ChemistrySelect, 2019, 4, 4238–4243.
60. D. Madsen, C. Flensburg and S. Larsen, J. Phys. Chem. A, 1998, 102, 2177–2188.
61. C. Flensburg, S. Larsen and R. F. Stewart, J. Phys. Chem., 1995, 99, 10130–10141.
62. E. Espinosa, I. Alkorta, J. Elguero and E. Molins, J. Chem. Phys., 2002, 117, 5529–5542.
63. B. Dereka, Q. Yu, N. H. C. Lewis, W. B. Carpenter, J. M. Bowman and A. Tokmakoff, Science, 2021, 371, 160–164.
64. B. L. Rodrigues, R. Tellgren and N. G. Fernandes, Acta Crystallogr., Sect. B: Struct. Sci., 2001, 57, 353–358.
65. A. Krawczuk, M. Gryl, M. B. Pitak and K. Stadnicka, Cryst. Growth Des., 2015, 15, 5578–5592.
66. F. Weinhold and R. A. Klein, Angew. Chem., Int. Ed., 2014, 53, 11214–11217.
67. P. R. Mallinson, G. T. Smith, C. C. Wilson, E. Grech and K. Wozniak, J. Am. Chem. Soc., 2003, 125, 4259–4270.
68. A. Gavezzotti, New J. Chem., 2016, 40, 6848–6853.
69. A. Gavezzotti, J. Phys. Chem. B, 2003, 107, 2344–2353.
70. A. Gavezzotti, Mol. Phys., 2008, 106, 1473–1485.
71. S. J. Grabowski, Chem. Rev., 2011, 111, 2597–2625.
72. A. G. P. Maloney, P. A. Wood and S. Parsons, CrystEngComm, 2015, 17, 9300–9310.
73. L. Pauling, General Chemistry, Dover Publications, E-book edn, 2014.
74. A. Gavezzotti, Acta Crystallogr., Sect. B: Struct. Sci., 2010, 66, 396–406.
75. I. Rozas, Phys. Chem. Chem. Phys., 2007, 9, 2782–2790.
76. F. Hibbert and J. Emsley, in Adv. Phys. Org. Chem., ed. D. Bethell, Academic Press, 1990, vol. 26, pp. 255–379.
77. J. Hoja, A. F. Sax and K. Szalewicz, Chem. – Eur. J., 2014, 20, 2292–2300.
78. S. C. C. van der Lubbe, F. Zaccaria, X. Sun and C. Fonseca Guerra, J. Am. Chem. Soc., 2019, 141, 4878–4885.
79. A. J. Cruz-Cabeza and J. Bernstein, Chem. Rev., 2014, 114, 2170–2191.
80. I. Majerz and M. J. Gutmann, RSC Adv., 2011, 1, 219–228.
81. A. O. F. Jones, N. Blagden, G. J. McIntyre, A. Parkin, C. C. Seaton, L. H. Thomas and C. C. Wilson, Cryst. Growth Des., 2013, 13, 497–509.
82. I. O. Majerz and I. Olovsson, Acta Chim. Slov., 2011, 58, 379–384.
83. F. Fontaine-Vive, M. R. Johnson, G. J. Kearley, J. A. Cowan, J. A. K. Howard and S. F. Parker, J. Chem. Phys., 2006, 124, 234503.
84. K. Saito and H. Ishikita, Proc. Natl. Acad. Sci. U. S. A., 2012, 109, 167–172.
85. W. M. Haynes, CRC Handbook of Chemistry and Physics, CRC Press, Boca Raton, 96th edn, 2016.
86. J. A. Cowan, J. A. K. Howard, G. J. McIntyre, S. M. F. Lo and I. D. Williams, Acta Crystallogr., Sect. B: Struct. Sci., 2005, 61, 724–730.
87. C. A. Morrison, M. M. Siddick, P. J. Camp and C. C. Wilson, J. Am. Chem. Soc., 2005, 127, 4042–4048.
88. J. A. Cowan, J. A. K. Howard, G. J. McIntyre, S. M.-F. Lo and I. D. Williams, Acta Crystallogr., Sect. B: Struct. Sci., 2003, 59, 794–801.
89. D. Allan, H. Nowell, S. Barnett, M. Warren, A. Wilcox, J. Christensen, L. Saunders, A. Peach, M. Hooper, L. Zaja, S. Patel, L. Cahill, R. Marshall, S. Trimnell, A. Foster, T. Bates, S. Lay, M. Williams, P. Hathaway, G. Winter, M. Gerstel and R. Wooley, Crystals, 2017, 7, 336.
90. D. A. Keen, M. J. Gutmann and C. C. Wilson, J. Appl. Crystallogr., 2006, 39, 714–722.

---

Footnote

† Electronic supplementary information (ESI) available. Crystallographic data CCDC 2068354–2068359; Experimental section and additional figures and tables are included in the ESI. For ESI and crystallographic data in CIF or other electronic format see DOI: 10.1039/d1ce00355k

---