A combined crystallography and DFT study on ring-shaped Cucurbit[n]urils: structures, surface character, and host-guest recognition.

A combined crystallography and DFT study of cucurbit[n]urils (n = 5-8, 10) was carried out, and PBE0 was certified to be the most rational density functional method for optimization task. Steric hindrance and electronic effect of the hindered lone pair electrons in cucurbit[n]urils were qualitatively measured by bond order analysis, lone pair electron (LP) visualization and electrostatic potential (ESP) study. Together with energy decomposition analysis of some selected host-guest systems, we quantitatively verified the effect of size/cavity and noncovalent interaction in host-guest recognition. This solid study revealed that lone pairs electrons affect not only on host-guest identification mode but also on geometry stability, which pave the avenue for further sophisticated applications. This journal is © The Royal Society of Chemistry.

Introduction

The first preparation of Cucurbit[n]uril (CB[n]) can be dated back to 1904 by Behrend and co-workers,[1] who reported the condensation reaction of formaldehyde and glycoluril. However, limited by the poor solubility and difficulty in crystallization, the single crystal structure of the condensed product was not obtained until 1981,[1] which hindered the understanding and application of CB[n]s. In the past twenty years, 7 more homologues were synthesized and CB[n] family now contains skeletons with 5 to 15 glycoluril units.[1] As shown in Fig. 1, CB[n]s (n ≤ 10) has pumpkin-like shapes. Crystallography study shows that all the single-ring CB[n]s have a congenial height approximate to 6.08 Å, and their inner/outer diameters range from 2.64 Å/5.41 Å to 6.46 Å/8.84 Å. Thanks to the hydrophobic cavities with adjustable volumes (from 12.38 Å3 to 416.49 Å3, more details in Fig. S1†), CB[n]s have been widely applied in catalysis,[2] photochemistry,[2] drug delivery,[2] biosensing,[2] gas storage,[2] and separation science,[2]etc.[2]

Fig. 1

Structural parameters of the reported ring-shaped CB[n]s.

With the prosperity of porous materials in the past decades, scientists have gained developing understandings of the recognition mode between CB[n]s and guest molecules. And some more sophisticated applications were developed as a result. Liu group[3] prepared the syn/anti-excimer complexes by assembling dipolar dyes styrylpyridiniums into CB[n] (n = 7 or 8), which performed reversibly tunable white-light emissions capability when opposite stacking orientations embedded. Agasti and co-worker[4] exploited dynamic supramolecular recognition between fluorescently labelled guest molecules and complementary CB[7] hosts, based on which a super-resolution imaging method for live-cell microscopy was established. Among these applications, cavity size and noncovalent interactions provided the vital contribution, and it's believed that theoretical calculation would supply more detailed information to the mechanism.

However, to the best of our knowledge, little theoretical studies on CB[n]s were conducted. Most reported theoretical studies focused on explaining the solubility or host–guest bonding strength of CB[n]s by referring to energies, while a general comparison about the similarity/difference among different CB[n]s was rarely discussed. What's more, we deemed that crystal structures should be paid more attention when theoretical study was involved, for which supplied the most solid information about CB[n] species. In this paper, we carried out our investigation on CB[n]s (n ≤ 10) through combining crystallography and ab initio calculation. By focusing on the structure, surface character, and host–guest binding properties, a deeper understanding of the host–guest recognition was achieved, which provided fundamental ideas for designing of such supramolecular self-assembly systems.

Computational details

Gaussian 16 package[5] and xtb program[5] were applied in the calculation. Original structures were set up based on single crystals structures obtained from CCDC database.[6] Semi-empirical method of PM7 (ref. 7) and GFN2-xTB,[5] together with DFT methods including B3LYP,[7] PBE0,[7] wB97-D3,[7] M06-2X,[7] and B3LYP combined D3 version of Grimme's dispersion[7] were applied in geometry optimization, using 6-31G(d) basis set.[7] All the geometry optimization tasks were performed in vacuum to reproduce the crystal environment for CB[n]s. Single point energies were calculated at M06-2X/def-TZVP[7] (SDD for I) level for the optimized structures and host–guest systems extracted from crystal structures. Wavefunction analysis including Mayer bond orders,[8] Mulliken charge population,[8] LP distribution,[8] ESP, energy decomposition analysis (EDA) and non-covalent interaction (NCI) analysis were performed by Multiwfn program,[5] using the fchk files generated in M06-2X/def-TZVP level. RDSM calculation and the visualization of structures and isosurface were carried out by VMD program.[5]

Results and discussion

Structural information

Firstly, some commonly used density functional methods were screened for the geometry optimization task, and the single point energies of the optimized structures were summarized. As shown in Table 1, the relative energies of the optimized structures were much lower than those of crystal structures, which implied that the monomer structure extract from periodic crystal structure was not the minimal point and geometry optimization was prerequisite for ab initio study. As shown in the first raw, energies of CB[n]s optimized by semi-empirical method PM7 were much higher (about 40 kcal mol−1) than those optimized by DFT methods, implying insufficient precision of PM7 for CB[n]s optimization task. Comparatively, GFN2-xTB provided a much better result with accuracy approaching to B3LYP. Considering of the calculation cost and hardware dependencies, GFN2-xTB is suitable and recommended for pre-optimized calculation. Among the selected DFT methods, B3LYP-D3, PBE0, and M06-2X showed similar capabilities for geometry optimization, while wB97-D3 was much inferior (about 8 kcal mol−1 higher). The B3LYP functional did not combine D3 version of Grimme's dispersion raised energies about 2 kcal mol−1 higher than those calculated from B3LYP-D3, inferring dispersion correction had a certain impact on the optimization of CB[n]s structures. By combining with the RMSD data (see Table S1†), it can be seen that PBE0 changed the crystal structure in the minimal level. In order to fit the crystal structures to the greatest extent, structures optimized by PBE0 were chosen for the following wavefunction analysis.

Relative energies of CB[n] optimized by selected theoretical methods. Basis set 6-31G(d) was applied for geometry optimizations; energies calculated at M06-2X/def-TZVP level, and setting the calculated energy of crystal as zero point

ΔG (kcal mol−1)	CB[5]	CB[6]	CB[7]	CB[8]	CB[10]
PM7	−196.2	−234.9	−179.6	−301.2	−403.4
GFN2-xTB	−236.4	−288.1	−242.6	−376.7	−502.9
B3LYP	−236.1	−287.1	−241.5	−375.1	−502.0
B3LYP-d3(BJ)	−238.0	−288.8	−243.9	−377.6	−504.9
PBE0	−239.0	−290.1	−245.2	−379.3	−507.2
wB97-D	−233.1	−282.9	−236.7	−369.6	−494.9
M06-2X	−239.7	−290.9	−246.1	−380.5	−508.6

With the optimized structures in hand, bond order index analysis was then performed. Bond connection in CB[n] skeletons could be summarized in 5 types, as shown in Scheme 1, which are C1–C2, C2–N3, C4–N3, C4–O5, and C6–N3. With n enlarging, bond order indexes of C4–O5 significantly decreased (from 1.9521 to 1.9161), while those of C1–C2 and C4–N4 tinnily increased (from 0.9474 to 0.9561 and from 1.0446 to 1.0603). Specially, bond order indexes of C2–N3 and C6–N3 held a tendency to increase first and then decrease (0.8999 → 0.9158 → 0.9103 and 0.900 3 → 0.9207 → 0.9103). As the chemical structures shown on the right, eight-membered rings (blue ring) constituted of C2–N3 and C6–N3, and five-membered rings (red ring) constituted of C1–C2, C2–N3, and C4–N3, hence the rigidity of the ring could be described by the involved bond order indexes. In other word, the growth of the bond order indexes of C1–C2 and C4–N3 implied the stability raised in the five-membered rings with n increasing. And eight-membered rings were most stable when n equalled to 7. As a result, CB[7] was less deformed than other CB[n]s, and CB[7] should be the most stable and rigid compound among the selected structures.

Scheme 1

Mayer bond order index variation of CB[n] skeletons with n enlarging.

Apparently, the bond order index changes of C1–C2, C2–N3, C4–N3, and C6–N3 were closely related to those of C4–O5, and we assumed that the bond order index decreases of C4–O5 as n increasing were associated with the LP and ring tension. As shown in Fig. 2, the LP were visualized by Suresh's method,[8] where isosurface value equalled to the minimal add 10 kcal mol−1 were drawn. Apparently, owing to the limitation of the ring size, LP overlapped in CB[5] and CB[6]. The LP overlapping reserved electrons in high density and made hydrogen bond prone to form with those doners contributing hydrogen. Similar to steric effect, the overlap of LP brought repulsive force, which pressed on the C4–O5 bond and shorten it. In addition, the overlap of LP might squeeze each other and force the electrons to reflow to the C4–O5 bond. With the ring size enlarging, the overlap of LP gradually decreased in CB[7] and CB[8], and finally became independent in CB[10]. Ring enlargement decreased LP overlap and made room for keeping LP, which determined the Mulliken charges population of oxygen as CB[5] owned the minimal negative value (−0.3522) while CB[10] owned the maximum negative one (−0.3976). Ring tension decreased with n enlarging (more evidences in Fig. S2†), leading to enhancement of coplanarity of the urea functional groups, thus improved the conjugation of CN bonds. Hence, those electrons formerly restricted to CO bonds could be more shared by the urea functional groups, resulting in the decrease and increase of bond order indexes of C4–O5 and C4–N3 respectively, and Mulliken charge population of nitrogen became more negative (−0.0858 to −0.0981).

Fig. 2

Lone pair electron distribution in CB[n]s.

Surface character

ESP plays an important role in the noncovalent interaction for polarized structures, which determines both the interaction area and strength of the electrostatics and inductive force. Thus, ESP study brings insight into more in-depth understanding of host–guest recognition. The van der Waals surface map colored by ESP of CB[n]s (n = 5, 6, 7, 8, 10) were summarized in Fig. 3. It can be clearly identified that the negative ESP value distributed in the both ends, above/below the carbonyl oxygen, and the positive ESP value distributed in the middle of the structures, on the outside of hydrogen. The up-down distributed negative ESP area were electron-rich, where cation or electron-deficient species were prone to be trapped. Owing to the enlargement of the pumpkin-like shape from CB[5] to CB[10], the surface outside of hydrogens expanded, resulting in gradually increasing ratio of the positive ESP surface area (the red area increased in histogram). However, the same change brought opposite effect to the negative ESP surface. The enlargement of the CB[n]s departed the carbonyl oxygen, which decreased their LP overlap and made negative ESP value less negative (the blue area shifted to 0 in histogram). Combined with the LP distribution analysis above, we can infer that electron density increased because of ring tensions in compact rings. Therefore, electrostatic attraction held a greater strength among the smaller CB[n]s. By influence of the distances and angles of the LP, CB[5] and CB[6] usually preferred to co-crystallize with cations. On the other hand, the electron-neutral ESP value mostly distributed inside the structures. The inner surface gradually became larger with the enlargement of CB[n]s, leading those low-polarized structures tended to be trapped into the cavity. Lots of crystallographic examples support these conclusions, which will be further discussed below.

Fig. 3

Electrostatic potential (ESP) character of CB[n]s.

Host–guest recognition

Independent gradient model (IGM)[9] takes advantages in analysing host–guest interaction by dividing the integral complex into fragments and visualizes the noncovalent interaction between fragments optionally, which have been widely applied in pioneer's work.[9] As shown in Fig. 4, IGM isosurface of some typical ion/neutral molecule-CB[n] complex structures were presented with attachment of the result of EDA based on AMBER force field[10] and Merz–Kollmann ESP fitting atomic charge.[11] The selected host–guest systems can be treated as typical representatives of several types, which were described as metal ion (cluster) outside of CB[n] (Mg2+(H2O)4@CB[6]), organic ion outside of CB[n] (NH4+@CB[5]*), polarized molecule outside of CB[n] (H2O@CB[5]*), charged organic compound inside of CB[n] (Az+@CB[8]), polarized molecule inside of CB[n] (CO2@CB[5]*), and non-polarized molecule inside of CB[n] (I2@CB[6]).

Fig. 4

Independent gradient model (IGM) analysis and energy decomposition analysis (EDA) of the selected host–guest systems.

Consistent with the above conclusions, ions and hydrogen bond doners, such as Mg2+, NH4+, and H2O were caught by carbonyl oxygens. Hydrogen bonds formed in these complexes as the IGM isosurface showed (colored in pale bule and blue-green). According to the hydrogen bond energies (noted in red), hydrogen bond strength was slightly stronger in charged complexes than that in neutral ones (−17.82 vs. −11.75 kJ mol−1). Strong attraction (isosurface colored in deep blue) and repulsion (isosurface colored in deep red) existed between Mg2+(H2O)4 and CB[6] of Mg2+(H2O)4@CB[6], and we inferred it should be a manifestation of strong coordination bonds. Neutral molecules and those compounds which cannot donor hydrogen to form hydrogen bonds preferred to be trapped within the cavities and form complex structures, for example, CO2@CB[5]*, I2@CB[6], and Az+@CB[8]. IGM pointed out the repulsion and dispersion betwixt the host and guest in CO2@CB[5]* and I2@CB[6] were in medium intensity (isosurface mostly colored in bright green and bright grey-green), whose isosurface were tattered and in small area. Relatively, isosurface in Az+@CB[8] was much more diffuse, and coloured method showed that the host–guest interactions were stronger (isosurface colored in green and dark grey-green).

The result of EDA qualitatively matched with the basis set superposition error (BSSE) corrected complexation energies. The value of EDA/corrected complexation energies were −859.31/−609.56, −419.11/−340.59, −34.98/−37.66, −346.02/−408.79, −20.71/−40.05, and −67.32/−84.04 (kJ mol−1) for Mg2+(H2O)4@CB[6], NH4+@CB[5]*, H2O@CB[5]*, Az+@CB[8], CO2@CB[5]*, and I2@CB[6], respectively. Generally, ions held a 5 to 10-times of electrostatic intensity than neutral molecules did, and electrostatic force was a vital factor to keep polarized guest outside the host's cavity. On the other hand, dispersion was related to molecular volume and interaction surface area, and dispersion was the main contributor for maintaining guest trapped in cavities.

EDA histogram indicated that electrostatic force played a negative role in stabilizing the complex structures for CO2@CB[5]*, which could be ascribed to LP repulsion. Carbonyl oxygens of CB[5]* in both ends were in negative ESP value (Fig. 3), the same as oxygens of CO2. Groups with same charge repels each other and repels and raises up energy. Therefore, the most stable conformation was formed when CO2 was stuck inside the cavity, where electrostatic repulsion of the both ends offset each other out and dispersion can be improved to the maximum level as the largest contact surface was achieved.

The interaction modes found in these complex structures consistent with the conclusions drawn by ESP analysis, that is high-polarized structures tended to interact at both end of CB[n]s while low-polarized structures preferred to be trapped in the central cavities. The example of I2@CB[6] implied that non-polarized molecule trapped inside the cavity could be stabilized without electrostatic, and maintained only by van der Waals forces. And the example of Az+@CB[8] disclosed that dispersion intensity could be expanded remarkably when fully-trapped guests owned large VDW surface. Hence, more attention should be paid for the contribution of dispersion, which was even more effective in reducing the energy of the system than electrostatics when large-sized molecular were trapped in CB[n]s. The effect of dispersion was a key role in the applications such as drug delivery, molecular recognition and separation science for CB[n]s, and it was the first time discussed by visualization method in qualitative/quantitative.

Conclusions

In summary, we confirmed that the PBE0 was suitable for geometry optimization task of CB[n]s, which got the balance between least structural changes and lowest relative energies. We disclosed that the eight-membered ring became more flexible while five-membered glycoluril ring turned to be firmer with ring size enlarging, which was controlled by the conjugation effect, ring tension, and LP repulsion. To address, LP repulsion in CB[n] was carefully discussed for the first time, which not only effect on the geometry structures, but also on the electrostatics potential distribution. Combined with IGM and EDA analysis, we revealed how electrostatic dominated polar host–guest recognition in both ends of CB[n]s and dispersion played an important role in stabilizing guests trapped in cavities, in qualitative/semi-quantitative. These studies provided solid evidences and clear understandings of the CB[n]-guest recognition mode/mechanism, bringing insight to the further research and sophisticated applications.

Conflicts of interest

There are no conflicts to declare.