Multiple hydrogen-bonded dimers: are only the frontier atoms relevant?

Abstract

Non-frontier atom exchanges in hydrogen-bonded aromatic dimers can induce significant interaction energy changes (up to 6.5 kcal mol⁻¹). Our quantum-chemical analyses reveal that the relative hydrogen-bond strengths of N-edited guanine–cytosine base pair isosteres, which cannot be explained from the frontier atoms, follow from the charge accumulation in the monomers.

---

Introduction

Intermolecular hydrogen bonding is one of the central interactions underlying self-assembly and molecular recognition in biochemistry, such as in protein folding and DNA duplex formation by constituting the hydrogen bonds (H-bonds) between the complementary DNA base pairs.¹ The ability to self-assemble has inspired the field of supramolecular chemistry to incorporate intermolecular H-bonding in the design of novel catalysts,² (macro)molecules,³ and materials.⁴ However, in order to rationally design new and improved molecules and materials, a profound understanding of the mechanism of intermolecular H-bonding and the prediction of the interaction strength is required.

The strength of an H-bond is often explained from the properties of the interacting frontier atoms, that is, the partially positively charged H-bond donor (i.e., H(–O) or H(–N)) and the partially negatively charged H-bond acceptor (i.e., O or N) groups.^5–7 An example of this is the secondary electrostatic interaction (SEI) model⁶ from Jorgensen and Pranata which is widely used to predict and explain trends in intermolecular H-bond strengths.⁷ This model is based on the assumption that the H-bond strength between monomers can be predicted from the position and charge of the frontier atoms. However, the SEI model is physically incorrect due to the oversimplification of the H-bonding mechanism by regarding the H-bond donor and acceptor groups as interacting point charges. Popelier and Joubert⁸ showed that the consideration of only frontier-atom electrostatic interactions, like in the SEI model, is elusive because also the electrostatic interactions between distant atoms contribute to the stability of H-bonded pairs. More importantly, it is nowadays well-established that H-bonds are not purely electrostatic in nature but also contain a significant covalent character,⁹ which arises from the donation of electronic density from the filled σ-lone pair orbital of the H-bond acceptor into the empty antibonding σ* orbital on the H-bond donor. This additional stabilizing orbital-interaction component is in fact essential to overcome the destabilizing Pauli repulsion associated with H-bond formation (see ref. 9a for an overview of all the relevant interaction components of H-bonds). Our previous work into the predictive nature of the SEI model revealed that the relative stability of H-bonded pairs follows from measuring the charge accumulation in the monomers rather than from the secondary interactions between H-bond frontier atoms.^10,11

In this work, we challenge the concept of explaining H-bond strengths by the frontier atoms further. We show that the binding strength of H-bonded pairs with identical H-bond frontier atoms can be adjusted by minimal non-frontier atom exchanges. Our dispersion-corrected density functional theory (DFT-D) based analyses at the ZORA¹²-BLYP¹³-D3(BJ)¹⁴/TZ2P¹⁵ level using the Amsterdam Modeling Suite¹⁶ (AMS2022.101) (see ESI,† Method S1 for the full computational details) show that the stability of mimics of the guanine–cytosine (GC) DNA base pair can be systematically tuned by varying the position of the non-frontier N and NH groups (see Fig. 1a). We demonstrate that these minimal non-frontier modifications have significant implications on the H-bond strength of the corresponding base pairs although the charges on the frontier atoms are barely affected in these so-called N-edited G and C isosteres (Fig. 1b). We find that the relative binding strengths can be rather explained and predicted by measuring the amount of charge accumulation in the monomers. This allows for control of the H-bond strength of these N-edited nucleobases, which can be applied in the design of innovative (bio)supramolecular building blocks.¹⁷

Fig. 1 (a) Schematic structure of the guanine–cytosine (GC) base pair and N-edited isosteres used in this work to tune the intermolecular H-bond strength by the position of the non-frontier N and NH groups highlighted in red and blue, respectively. (b) Voronoi deformation density (VDD) atomic charges Q (in milli-electrons) of the H-bond frontier atoms of the isolated guanine (G1 and G2) and cytosine (C and C1) isosteres in the geometry within the dimer with C or G1, respectively.

Results and discussion

Based on the similar charges on the H-bond frontier atoms, shown in Fig. 1b, one would expect no significant change in the H-bond strength of the GC base pair upon changing the position of the non-frontier heteroatoms in the G and C monomers.¹⁸ However, we find that the interaction energy (ΔE_int) of the GC H-bonded base pair analogues, presented in Fig. 2a, becomes more stabilizing from G1C (blue) < G2C (orange) < G1C1 (pink) < G2C1 (green), associated with a maximum stabilization of the H-bond interaction energy by ca. 6.5 kcal mol⁻¹. So, the H-bond interaction is enhanced upon (i) changing the G1 isostere to G2, thereby grouping the partially negative N and O atoms on one side of the molecule, and (ii) upon changing C to C1, in which the partially positive NH₍₂₎ groups are in closer proximity. Note that we discuss here the interaction energies because they dictate the trend in the relative H-bond strengths (see ESI,† Data S1 and see Method S2 for computational details).

Fig. 2 (a) Hydrogen-bonded GC base-pair isosteres with equilibrium hydrogen-bond distances (in Å) and interaction energy ΔE_int (in kcal mol⁻¹ in between brackets) and (b) decomposed interaction energy terms (ΔE_int = ΔV_elstat + ΔE_Pauli + ΔE_oi (+ ΔE_disp)) as a function of the middle H-bond distance r_N(H)⋯N (in Å, step size of 0.01 Å).

To understand what causes this stabilization, ΔE_int was partitioned into four physically meaningful terms using a quantitative energy decomposition analysis (EDA):¹⁹ (i) the classical electrostatic interaction (ΔV_elstat), (ii) the steric Pauli repulsion (ΔE_Pauli) arising from the repulsion between overlapping closed-shell orbitals on the interacting monomers, (iii) the orbital interaction (ΔE_oi) which accounts for charge transfer (i.e., covalency) in the σ-electronic system and polarization of the π-electronic system, and (iv) the dispersion energy (ΔE_disp) (see ESI,† Method S2 for a theoretical overview of this method).

The EDA can be performed on the equilibrium dimers (see ESI,† Data S1). However, as the H-bond distances vary for the different GC base-pair isosteres (see Fig. 2a), a more insightful picture can be obtained when performed at similar H-bond distances. This allows us to differentiate between interaction terms that are intrinsically more stabilizing, from the interaction terms that are simply enhanced by shortened H-bond distances. To this end, the EDA was performed as a function of the middle H-bond distance r_N(H)⋯N using the PyFrag 2019 software,²⁰ in which the two monomers in the geometry of the equilibrium dimer approach each other as two frozen blocks (ESI,† Method S2 for details). The results of this analysis are presented in Fig. 2b and show that the stabilization of ΔE_int from G1C (blue) < G2C (orange) < G1C1 (pink) < G2C1 (green) is preserved along the entire H-bond distance range. Furthermore, it shows that ΔE_int becomes more stabilizing upon interchanging G1 to G2 and C to C1 because both the electrostatic (ΔV_elstat) and orbital interactions (ΔE_oi) become more stabilizing along this trend. For ΔE_oi, in particular, the σ-charge transfer orbital interactions (ΔE^σ_oi) are enhanced, while the π-orbital interactions (ΔE^π_oi) do not really vary for the different isosteres. This means that the π-resonance assistance in the H-bonds, that is π-polarization, is not affected by the non-frontier modifications. Note that ΔE_Pauli does not dictate the trend, as this interaction component becomes more destabilizing from G1 to G2 and from C to C1 (ΔE_disp stays constant for all isosteres, see ESI,† Data S2).

While the similar frontier atomic charges of the isosteres cannot explain the enhancement of ΔV_elstat and ΔE^σ_oi (vide supra), we find that the origin of this effect lies in the amount of charge accumulation in the N-edited monomers, as is presented in Fig. 3. Upon going from G1 to G2 and from C to C1, the partially negative (O and N) and partially positive (NH and NH₂) functional groups, respectively, are grouped on one side of the molecule, which induces an increase of the molecular dipole moment (Fig. 3a). This grouping of the polar functional groups is associated with an increase of charge accumulation on both sides of the molecule, as is visualized by the electrostatic potential surfaces in Fig. 3b and the Voronoi deformation density (VDD)²¹ charges in Fig. 3c (see ESI,† Method S3 for details about the VDD method and Data S3 for the complete VDD charge analysis). Although the charge on the H-bond frontier atoms does not significantly change (Fig. 1b), the grouping of all polar groups on one side of the molecule by going from G1 to G2 and from C to C1, enhances the molecular charge accumulation, that is, the interacting monomers become more pronounced ‘mirrored dipoles’ of each other (see Fig. 3c), which enhances the intermolecular electrostatic interactions. This is a manifestation of the fact that the total electrostatic interaction is not exclusively determined by the frontier atoms but also by long-range electrostatic interactions between distant atoms.

Fig. 3 (a) Molecular dipole moment vectors (δ⁻ → δ⁺ direction) and magnitudes (in Debye), (b) molecular electrostatic potential surfaces (at 0.01 a.u.) from −0.1 (red) to 0.1 (blue) a.u., (c) Voronoi deformation density (VDD) atomic charges Q (in milli-electrons), and (d) isosurfaces (at 0.03 a.u.) with corresponding energies ε (in between brackets in eV) of the relevant H-bonding unoccupied (σ_LUMO) and occupied (σ_HOMO) orbitals of the isolated guanine (G1 and G2) and cytosine (C and C1) isosteres in the geometry within the dimer with C or G1, respectively.

Besides ΔV_elstat, the σ-orbital interactions ΔE^σ_oi are also enhanced upon grouping the polar functional groups on one side of the nucleobases thereby increasing the molecular charge accumulation. A larger accumulation of negative charge destabilizes occupied orbitals (i.e., they become better electron donors), while a larger accumulation of positive charge stabilizes unoccupied orbitals (i.e., they become better electron acceptors) (see Fig. 4).¹⁰ Important to note here is that molecular orbitals (MOs) involved in the H-bonds are delocalized over larger parts of the molecule and are not exclusively located on the frontier groups (Fig. 3d). So even if the charge at the frontier atoms does not change, a change in charge on other parts of the molecule causes shifts in the energies of the MOs [as long as the MO is asymmetric with respect to the molecule (see below)]. As such, we observe that upon changing from G1 to G2 and from C to C1, the larger accumulation of negative charge on the side of the H-bond acceptor group(s), that is, where the highest-occupied molecular orbitals (σ_HOMOs) involved in the intermolecular H-bonds have their coefficients, destabilizes the σ_HOMOs (see Fig. 3d). At the same time, the larger accumulation of positive charge on the side of the H-bond donor group(s) stabilizes the lowest-unoccupied molecular orbitals (σ_LUMOs) as these MOs have a coefficient on this side. Note that the σ_HOMO−1 of the cytosine isosteres is not affected upon changing C to C1 because this MO is relatively symmetric and therefore there is no significant net effect of the change in molecular charge accumulation. Nevertheless, both the destabilization of the other σ_HOMOs and the stabilization of the σ_LUMOs reduces the energy gap (Δε) between the interacting orbitals and, therefore, results in better σ-orbital interactions in the case of enhanced charge accumulation (Fig. 4), that is, upon changing from G1 to G2 and from C to C1.²²

Fig. 4 A larger accumulation of negative charge (δ⁻) destabilizes the occupied molecular orbitals (σ_HOMOs) of H-bond acceptors, while a larger accumulation of positive charge (δ⁺) stabilizes the unoccupied molecular orbitals (σ_LUMOs) of H-bond donors. This results in a smaller HOMO–LUMO energy gap (Δε) and therefore a better σ-orbital interaction in case of enhanced charge accumulation.

Conclusions

The binding strength of H-bonded aromatic dimers can be systematically tuned by varying the position of heteroatoms on the backside of the monomers while keeping the H-bond frontier atoms unchanged. These minimal non-frontier modifications can induce large interaction energy changes (ΔΔE_int up to 6.5 kcal mol⁻¹) although the charge of the H-bond frontier atoms is almost unaffected, as appears from our quantum-chemical analyses. Nevertheless, the grouping of the polar functional groups on one side of the monomer induces a larger accumulation of positive and negative charge in the molecule. This enhances the intermolecular electrostatic interactions, as well as the σ-orbital interactions by reducing the energy gaps between the interacting molecular orbitals. This is a manifestation of two effects: (i) electrostatic interactions are not only between frontier atoms but also between distant atoms (i.e., long-range); (ii) molecular orbitals are delocalized over the molecule and their energies are therefore influenced by charge changes in non-frontier parts of the molecule. Our findings challenge the concept of explaining H-bond strengths by the frontier atoms, revealing that the strength of the electrostatic interactions, as well as the orbital interactions, are determined by the whole backbone of the H-bond donor and acceptor fragments.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

The authors thank the Dutch Research Council (NWO), the Spanish MINECO (PID-2019-106830GB-I00, PID2022-138861NB-I00, and CEX2021-001202-M), and the Generalitat de Catalunya (2021SGR442) for financial support. This work was carried out on the Dutch national e-infrastructure with the support of SURF Cooperative.

Notes and references

1. (a) J. M. Berg, J. L. Tymoczko and L. Stryer, Biochemistry, W. H. Freeman and Company, New York, 5th edn, 2002; (b) G. A. Jeffrey and W. Saenger, H-bonding in Biological Structures, Springer, Berlin, 1991.
2. (a) S. Vera, A. García-Urricelqui, A. Mielgo and M. Oiarbide, Eur. J. Org. Chem., 2023, 26, e202201254; (b) G. Su, M. Formica, K. Yamazaki, T. A. Hamlin and D. J. Dixon, J. Am. Chem. Soc., 2023, 145, 12771.
3. (a) S. M. Wales, D. T. J. Morris and J. Clayden, J. Am. Chem. Soc., 2022, 144, 2841; (b) L. N. J. de Windt, Z. Fernández, M. Fernández-Míguez, F. Freire and A. R. A. Palmans, Chem. – Eur. J., 2022, 28, e202103691; (c) K. Aratsu, R. Takeya, B. R. Pauw, M. J. Hollamby, Y. Kitamoto, N. Shimizu, H. Takagi, R. Haruki, S. I. Adachi and S. Yagai, Nat. Commun., 2020, 11, 1623.
4. (a) R. B. Lin, Y. He, P. Li, H. Wang, W. Zhou and B. Chen, Chem. Soc. Rev., 2019, 48, 1362; (b) R. Sarkar, S. Majumdar, S. Kuil, J. Mallens, J. J. B. van der Tol, R. P. Sijbesma, J. P. A. Heuts and A. R. A. Palmans, J. Polym. Sci., 2023, 61, 1335.
5. (a) E. Arunan, G. R. Desiraju, R. A. Klein, J. Sadlej, S. Scheiner, I. Alkorta, D. C. Clary, R. H. Crabtree, J. J. Dannenberg, P. Hobza, H. G. Kjaergaard, A. C. Legon, B. Mennucci and D. J. Nesbitt, Pure Appl. Chem., 2011, 83, 1619; (b) IUPAC, Compendium of Chemical Terminology, ed. A. D. McNaught, A. Wilkinson, Blackwell Scientific Publications, Oxford, 2nd edn (the “GoldBook”), 1997.
6. (a) W. L. Jorgensen and J. Pranata, J. Am. Chem. Soc., 1990, 112, 2008; (b) J. Pranata, S. G. Wierschke and W. L. Jorgensen, J. Am. Chem. Soc., 1991, 113, 2810.
7. (a) M. Papmeyer, C. A. Vuilleumier, G. M. Pavan, K. O. Zhurov and K. Severin, Angew. Chem., Int. Ed., 2016, 55, 1685 ( Angew. Chem. , 2016 , 128 , 1717 ); (b) G. Bella, A. Santoro, F. Nicolò, G. Bruno and M. Cordaro, ChemPhysChem, 2021, 22, 593.
8. P. L. A. Popelier and L. Joubert, J. Am. Chem. Soc., 2002, 124, 8725.
9. (a) S. C. C. van der Lubbe and C. Fonseca Guerra, Chem. – Asian J., 2019, 14, 2760; (b) S. C. C. van der Lubbe and C. Fonseca Guerra, Chem. – Eur. J., 2017, 23, 10249; (c) L. de Azevedo Santos, T. C. Ramalho, T. A. Hamlin and F. M. Bickelhaupt, Chem. – Eur. J., 2022, 29, e202203791; (d) C. Nieuwland and C. Fonseca Guerra, Chem. – Eur. J., 2022, 28, e202200755; (e) C. Fonseca Guerra, F. M. Bickelhaupt, J. G. Snijders and E. J. Baerends, Chem. – Eur. J., 1999, 5, 3581; (f) C. Nieuwland, S. J. Lekanne Deprez, C. de Vries and C. Fonseca Guerra, Chem. – Eur. J., 2023, 29, e202300850; (g) C. Nieuwland, T. A. Hamlin, C. Fonseca Guerra, G. Barone and F. M. Bickelhaupt, ChemistryOpen, 2022, 11, e202100231; (h) F. Weinhold and R. A. Klein, Angew. Chem., Int. Ed., 2014, 53, 11214; (i) S. J. Grabowski, Chem. Rev., 2011, 111, 2597.
10. S. C. C. van der Lubbe, F. Zaccaria, X. Sun and C. Fonseca Guerra, J. Am. Chem. Soc., 2019, 141, 4878.
11. S. C. C. van der Lubbe, A. Haim, T. van Heesch and C. Fonseca Guerra, J. Phys. Chem. A, 2020, 124, 9451.
12. E. van Lenthe, A. Ehlers and E. J. Baerends, J. Chem. Phys., 1999, 110, 8943.
13. (a) A. D. Becke, Phys. Rev. A, 1988, 38, 3098; (b) C. Lee, W. Yang and R. G. Parr, Phys. Rev. B: Condens. Matter Mater. Phys., 1988, 37, 785; (c) Q. Wu and W. Yang, J. Chem. Phys., 2002, 116, 515.
14. (a) S. Grimme, J. Comput. Chem., 2004, 25, 1463; (b) S. Grimme, J. Comput. Chem., 2006, 27, 1787; (c) S. Grimme, J. Antony, S. Ehrlich and H. Krieg, J. Chem. Phys., 2010, 132, 154104; (d) S. Grimme, S. Ehrlich and L. Goerigk, J. Comput. Chem., 2011, 32, 1456.
15. E. van Lenthe and E. J. Baerends, J. Comput. Chem., 2003, 24, 1142.
16. (a) R. Rüger, M. Franchini, T. Trnka, A. Yakovlev, E. van Lenthe, P. Philipsen, T. van Vuren, B. Klumpers and T. Soini, AMS2022.101, SCM, Theoretical Chemistry, Vrije Universiteit Amsterdam, Amsterdam, 2022, https://www.scm.com; (b) C. Fonseca Guerra, J. G. Snijders, G. te Velde and E. J. Baerends, Theor. Chem. Acc., 1998, 99, 391; (c) G. te Velde, F. M. Bickelhaupt, E. J. Baerends, C. Fonseca Guerra, S. J. A. van Gisbergen, J. G. Snijders and T. Ziegler, J. Comput. Chem., 2001, 22, 931.
17. X. Li, M. Diguele Romero, S. Tcaturian, K. Kurpiewska and A. Dömling, J. Org. Chem., 2023, 88, 9823.
18. We focus in our interaction energy analysis on the smaller isostere of guanine, that is G1, which is obtained by removal of the 5-membered ring in guanine. However, this structural alteration does not change the H-bond strength of the corresponding GC base pair (see ESI† Data S4) and yields very similar frontier atomic charges (ESI† Data S3: Fig. S2). Also, there is a third possible guanine isostere (G3) which yields an average effect of G1 and G2 (see ESI† Data S1).
19. (a) F. M. Bickelhaupt and E. J. Baerends, in Reviews in Computational Chemistry, ed. K. B. Lipkowitz and D. B. Boyd, Wiley-VCH, New York, 2000, vol. 15, pp. 1–86; (b) T. A. Hamlin, P. Vermeeren, C. Fonseca Guerra and F. M. Bickelhaupt, in Complementary Bonding Analysis, ed. S. Grabowsky, De Gruyter, Berlin, 2021, vol. 8, pp. 199–212.
20. (a) X. Sun, T. M. Soini, J. Poater, T. A. Hamlin and F. M. Bickelhaupt, J. Comput. Chem., 2019, 40, 2227; (b) X. Sun, T. Soini, L. P. Wolters, W.-J. van Zeist, C. Fonseca Guerra, T. A. Hamlin and F. M. Bickelhaupt, PyFrag 2019, Vrije Universiteit, Amsterdam, 2019.
21. (a) C. Nieuwland, P. Vermeeren, F. M. Bickelhaupt and C. Fonseca Guerra, J. Comput. Chem., 2023, 44, 2108; (b) C. Fonseca Guerra, J. W. Handgraaf, E. J. Baerends and F. M. Bickelhaupt, J. Comput. Chem., 2004, 25, 189.
22. For the different GC base-pair isosteres the orbital overlap differences are small and in line with the relative H-bond distances (see ESI† Data S3: Table S4). Therefore, for the trend in the σ-orbital interactions the relative σ_HOMO–σ_LUMO energy gaps are more relevant than the overlap.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: https://doi.org/10.1039/d3cp05244c

---