Secondary Electrostatic Interaction Model Revised: Prediction Comes Mainly from Measuring Charge Accumulation in Hydrogen-Bonded Monomers.

The secondary electrostatic interaction (SEI) model is often used to predict and explain relative hydrogen bond strengths of self-assembled systems. The SEI model oversimplifies the hydrogen-bonding mechanisms by viewing them as interacting point charges, but nevertheless experimental binding strengths are often in line with the model's predictions. To understand how this rudimentary model can be predictive, we computationally studied two tautomeric quadruple hydrogen-bonded systems, DDAA-AADD and DADA-ADAD. Our results reveal that when the proton donors D (which are electron-donating) and the proton acceptors A (which are electron-withdrawing) are grouped together as in DDAA, there is a larger accumulation of charge around the frontier atoms than when the proton donor and acceptor groups are alternating as in DADA. This accumulation of charge makes the proton donors more positive and the proton acceptors more negative, which enhances both the electrostatic and covalent interactions in the DDAA dimer. The SEI model is thus predictive because it provides a measure for the charge accumulation in hydrogen-bonded monomers. Our findings can be understood from simple physical organic chemistry principles and provide supramolecular chemists with meaningful understanding for tuning hydrogen bond strengths and thus for controlling the properties of self-assembled systems.

1. Introduction

The self-assembly characteristics of hydrogen bonds in biochemical processes have inspired supramolecular chemists to employ them effectively in bottom-up synthesis for applications in materials science,[1,2] rational drug design,[3,4] and nanotechnology.[5] A profound understanding of the hydrogen bond mechanism and the prediction of its strength are therefore essential in the process of designing new materials because improper models lead to a waste of research time and resources.

Synthetic hydrogen-bonded dimers often resemble DNA base pairs because they are aromatic monomers connected through multiple hydrogen bonds. Experimentalists have measured different association constants for these dimers with the same number of hydrogen bonds connecting the monomers.[6−8] Because hydrogen bonds were believed to be primarily an electrostatic interaction of spherical (point) charges between a positively charged hydrogen atom of a proton donor (D) and a negatively charged proton acceptor (A) of the other monomer (Figure a), these experimental observations were not understood.

Figure 1

Hydrogen bond energy components that can influence the bonding strength. (a) Electrostatic interaction, (b) charge-transfer interactions, (c) Pauli (steric) repulsion, (d) dispersion, (e) hydrogen bond cooperativity, (f) π-resonance assistance, and (g) secondary electrostatic interactions.

Jorgensen and Pranata attempted to explain these puzzling phenomena by introducing the concept of the secondary electrostatic interaction (SEI) based on the hydrogen bond energy differences between guanine-cytosine (GC) and 2,6-diaminopyridine-uracil (PU).[9] SEIs are defined as the diagonal interactions between adjacent hydrogen bonds (Figure g). They can be either attractive between A and D atoms (green arrows in Figures , 2, and 5) or repulsive between two diagonally opposed A or D atoms (red arrows in Figures , 2, and 5). On the basis of these assumptions, the strongest hydrogen-bonded pairs are formed between monomers in which all of the D atoms are aligned on one monomer and all of the A atoms are aligned on the other monomer (e.g., AAA-DDD), while the weakest pairs are formed between pairs with alternating A and D atoms (e.g., ADA-DAD). The model accounts for 2 to 3 kcal mol–1 per SEI, so a DD-AA dimer (−2 SEIs) is predicted to be 8–12 kcal mol–1 more stable than a DA-AD dimer (+2 SEIs).

Figure 2

(a) Chemical formulas with repulsive (red arrows) and attractive (green arrows) SEIs, and their optimized structures with hydrogen bond distances (in Å) and interaction energies between brackets (in kcal mol–1) and (b) decomposed energy terms (in kcal mol–1) as a function of the middle hydrogen bond distance r (in Å) for DDAA-AADD (blue) and DADA-ADAD (red). All data was obtained at the BLYP-D3(BJ)/TZ2P level of theory.

Figure 5

(a) Chemical formulas with repulsive (red arrows) and attractive (green arrows) SEIs and their optimized structures with hydrogen bond distances (in Å) and interaction energies in brackets (in kcal mol–1), (b) decomposed energy terms (in kcal mol–1) as a function of the average hydrogen bond distance r (in Å) for GC (blue) and PU (red), and (c) electrostatic potential surfaces (at 0.01 au) from −0.1 (red) to 0.1 (blue) a.u. All data was obtained at the BLYP-D3(BJ)/TZ2P level of theory.

The difference in binding strength between GC and PU can also be explained by considering the polarity of the monomers. Monomers with larger molecular dipole moments (such as G and C) are generally expected to form stronger hydrogen-bonded complexes than monomers with smaller molecular dipole moments (such as P and U, see Supporting Information Figure 1). This monomeric polarity model has been used successfully to explain, for example, hydrogen bond strengths of DNA base pairs.[10−12]

Nowadays, the SEI model is teaching material in organic chemistry textbooks[13,14] and is widely used in supramolecular chemistry to predict and explain the relative strengths of self-assembled systems.[15−22] However, the validity of the model has been questioned, and the predicted trends in hydrogen bond stabilities are not always in line with experimental trends.[23−27] Popelier and Joubert,[26] Tiwari and Vanka,[24] and most recently Hernández-Rodríguez, Rocha-Rinza, and co-workers[28,29] have demonstrated that all possible atom–atom interactions should be considered, not just the interactions between the frontier atoms as is done in the SEI model. Clark, Murray, and Politzer state that hydrogen bonds can be described in terms of Coulombic interactions, but only when the exact electronic density (or a reasonable approximation) is known and the polarization of each molecule is taken into account.[30,31]

Another concern is that the SEI model gives an incomplete picture of the bonding mechanism, as theoretical and experimental studies have shown that hydrogen bonds are not purely electrostatic but also partially covalent in nature.[32−36] These charge-transfer interactions enhance the hydrogen bonding via donor–acceptor interactions between the σ-lone pair orbital on the hydrogen acceptor atom and the antibonding σ* empty orbital on the opposing H-A bond (Figure b).[37−39] Besides this covalent component, there are also other components that can play a decisive role in relative hydrogen bond strengths and lengths, including Pauli repulsive interactions[40] (Figure c), dispersion interactions[41] (Figure d), hydrogen bond cooperativity,[38,42] (Figure e) and resonance assistance by the π electrons[43−46] (Figure f). The SEI model is entirely based on point charges and does not account for the long-range electrostatic interaction or any of the other components as in Figure b–f. Given this rudimentary description of hydrogen bonding, it is remarkable that experimental binding strengths are often in line with the model’s predictions.

To understand how this simple model can be predictive, we studied two tautomeric quadruple hydrogen-bonded dimers with the DDAA and DADA motif. We found that hydrogen bonds are strengthened when proton donors (which are electron-donating) and proton acceptors (which are electron-withdrawing) are grouped because they give rise to favorable charge accumulation around the frontier atoms. This monomeric charge accumulation is also responsible for the differences in binding strength between the GC and PU pairs on which the model is originally based.

Methodology

Computational Settings

All calculations were performed using the density functional theory (DFT)-based program Amsterdam Density Functional (ADF) 2017.208.[47−49] We used the dispersion-corrected BLYP-D3(BJ) functional in combination with a TZ2P basis set for geometry optimizations and energies,[50−53] which accurately reproduces the structural and energy properties of hydrogen-bonded systems.[54−56] The molecular figures were illustrated using CYLview.[57] Full computational details are available in Supporting Methods 1.

Energy Decomposition Analysis (EDA)

The hydrogen bond energy can be decomposed into the preparation energy ΔEprep and interaction energy ΔEint:The preparation energy ΔEprep is the energy needed to deform the monomers from their optimal geometry to the geometry that they acquire in the interacting dimer. The interaction energy ΔEint accounts for the actual chemical interaction between the prepared monomers and can be further decomposed into physically meaningful terms:The term ΔVelstat corresponds to the classical electrostatic interactions between the charge distributions of the prepared monomers and is usually attractive. The Pauli repulsion ΔEPauli comprises the destabilizing interactions between the overlapping, occupied orbitals of the two monomers and is responsible for any steric repulsion. The orbital interaction ΔEoi accounts for charge transfer (i.e., donor–acceptor interactions between the hydrogen-bonded monomers) and polarization (i.e., electron density redistribution on one monomer due to the presence of another monomer). In planar systems, the orbital interaction energy can be further decomposed into the contributions from the σ- and π-electron system. Finally, the term ΔEdisp is added to account for the dispersion corrections. A theoretical overview of this energy decomposition energy (EDA) scheme is given in Supporting Methods 2 and ref (58).

The hydrogen bond energies were analyzed as a function of the hydrogen bond distance r. In this approach, the hydrogen bond distances were varied over a certain interval while keeping the monomers frozen in the same geometry that they have in the fully optimized dimer (schematic representation given in Supporting Figure 2). The advantage of this approach is that we can compare the dimers with similar hydrogen bond lengths, which allows us to differentiate between the energy terms that are effectively stronger from the energy terms that are simply enhanced by the shortened bond distances. In other words, comparing the dimers with similar hydrogen bond distances r allows us to address the question of whether the relative hydrogen bond strengths and lengths are determined by the electrostatic interaction ΔVelstat, Pauli repulsion ΔEPauli, orbital interaction ΔEoi, or a combination between those terms (Figure in ref (40)). We also used a second approach in which we reoptimized the dimers with constrained hydrogen bond lengths while keeping the O···H–N and N···H–N angles linear (Supporting Figures 3 and 4). This approach gave us the same results as the one in which the monomers approach each other as blocks (i.e., without geometry reoptimization).

Voronoi Deformation Density (VDD) Charge

The charge distribution was analyzed by using the Voronoi deformation density (VDD) method. The VDD charges Q represent the flow of electron density when going from a fictitious promolecule (summation of atomic densities) to the final molecular density of the interacting system. A positive VDD charge Q thus corresponds to the loss of electrons, whereas a negative charge Q is associated with the gain of electrons. This VDD scheme can be extended to the analysis of hydrogen bonding by computing the change in electron density ΔQ that is associated with hydrogen bond formation. In this approach, the sum of densities of the prepared monomers is taken as the initial density, which offers direct insight into the redistribution of the electronic density caused by the formation of the interacting dimer. As a further analysis tool, the ΔQ charge can be decomposed into σ and π components. A theoretical overview is given in Supporting Methods 3 and ref (59).

Results and Discussion

Quadruple-Hydrogen-Bonded Systems

We have studied two quadruple hydrogen-bonded dimers with the DDAA (7-amino-1,8-naphthyridin-2(8H)-one) and DADA (7-amino-1,8-naphthyridin-2(1H)-one) motifs (Figure a). These tautomeric systems vary only by a simple proton transfer, which makes the direct comparison between them relatively straightforward. The dimer with the DDAA motif, which has two repulsive and four attractive SEIs, has an interaction energy of −71.8 kcal mol–1. Its tautomeric counterpart with the DADA motif has six repulsive SEIs and has an interaction energy of −32.9 kcal mol–1. The DDAA dimer is thus 38.9 kcal mol–1 more stable than DADA-ADAD. This energy difference might be attributed to the SEIs but can also be explained by using the monomeric polarity model because the DDAA monomer has a larger molecular dipole moment (11.5 D) than the DADA monomer (7.1 D, see Figure a). Therefore, to what extent can we actually attribute the difference in hydrogen bond stability to the polarities of the monomers?

Figure 3

(a) Molecular dipole moments (in Debye) with their arrows pointing from negative to positive charge and (b) optimized structures with hydrogen bond distances (in Å) and interaction energies between brackets (in kcal mol–1). All data was obtained at the BLYP-D3(BJ)/TZ2P level of theory.

To address this question, we altered the monomeric dipole moments by substituting a nitrogen atom in the aromatic rings, resulting in two new structures, DDAA* and DADA*. The DDAA* monomer has a quenched molecular dipole moment of 9.1 D (Figure a), and as a result, its interaction energy becomes 4.4 kcal mol–1 less stable than that of the original DDAA dimer (Figure b). On the other hand, the DADA* monomer has an enhanced molecular dipole moment of 9.5 D. The interaction energy of the resulting DADA* dimer is therefore 2.2 kcal mol–1 more stable than that of the original DADA dimer. These results demonstrate that the polarities of the monomers do indeed influence the hydrogen bond strength, but only to a small degree; the DDAA* dimer is still 32.3 kcal mol–1 more stable than the DADA* dimer, even though the monomeric dipole moments are slightly larger for DADA*. Clearly, there are other factors that are responsible for the differences in hydrogen bond strengths and lengths.

To understand the actual origin of these energy differences, we analyzed the hydrogen bond energy of each dimer by computing the interaction energies around their point of equilibrium. In this approach, the monomers approached each other as blocks by varying the middle hydrogen bond distance over an interval from 2.80 to 3.25 Å with 0.01 Å per step. The essential results are graphed in Figure b; the complete data set can be found in Supporting Tables 1 and 2.

The DDAA-AADD dimer is still more stable than the DADA-ADAD pair by up to 42.1 kcal mol–1, even though both dimers have now similar bond distances (Figure b). A large contribution of this enhanced interaction energy comes from the electrostatic interaction ΔVelstat, which is up to 20.1 kcal mol–1 more attractive for the DDAA-AADD pair. The remaining part of the enhanced interaction energy comes from the orbital interaction ΔEoi, which is up to 14.5 kcal mol–1 more favorable for the σ component and 11.3 kcal mol–1 more favorable for the π component in DDAA-AADD. The Pauli repulsion ΔEPauli is approximately the same for both dimers. We will now demonstrate that the differences in electrostatic interaction ΔVelstat and orbital interaction ΔEσ can be easily understood from the charge accumulation around the frontier atoms in the monomers.

Because the hydrogen bond donor groups are relatively electron-donating in nature (Supporting Figure 5), there is an accumulation of positive charge around their frontier atoms. As can be seen in the electrostatic potential surfaces in Figure a, this accumulation of positive charge is more pronounced when the proton donors are grouped together as in the DDAA monomer. On the other hand, the hydrogen bond acceptor groups are more electron-withdrawing in nature (Supporting Figure 5), which leads to an accumulation of negative charge around their frontier atoms (Figure a). Again, this accumulation of charge is more significant when both acceptor atoms are grouped together as in DDAA. We have also analyzed the electrostatic potential surfaces of the DDAA* and DADA* monomers (Supporting Figure 6). Even though these systems share similar molecular dipole moments, the charge accumulation around the frontier atoms is still more pronounced in DDAA* than in DADA*. This strongly suggests that this larger accumulation of charge is really a consequence of the grouping of proton donor and acceptor atoms.

Figure 4

(a) Electrostatic potential surfaces (at 0.01 au) from −0.1 (red) to 0.1 (blue) a.u., (b) isosurfaces (at 0.03 au) and energies (in eV) of the most important orbitals of the prepared fragments for DDAA (up) and DADA (down), and (c) Voronoi deformation density (VDD) π charges ΔQπ (in millielectrons) associated with the formation of the dimer. Computed at the BLYP-D3(BJ)/TZ2P level of theory.

The stronger accumulation of charge within the DDAA monomer enhances the hydrogen bond strength in two ways. First, the DDAA-AADD dimer is strengthened because it can participate in more favorable electrostatic interactions. That is, the proton atoms are more positively charged and the proton acceptor atoms are more negatively charged, resulting in enhanced primary electrostatic interactions. Of course, the attractive secondary interactions as identified by the SEI model are also augmented by the larger accumulation of charge, but to a smaller degree than the primary interactions because of their larger distances.

Second, the orbital interactions ΔEσ in DDAA-AADD are enhanced because the buildup of positive charge stabilizes the N–H antibonding acceptor orbitals σ*N–H while the buildup of negative charge destabilizes the σ lone pair orbitals (Figure b). For example, the σHOMO orbital in DDAA is 1.0 eV higher in energy while its σLUMO orbital is 0.8 eV lower in energy than in DADA, resulting in a 38% decrease in the HOMO–LUMO gap. The accumulation of charge in DDAA is thus responsible for a smaller occupied virtual energy gap and therefore a better orbital interaction.

Both dimers are further stabilized by polarization in the π-electron system. This so-called π-resonance assistance (Figure f) reinforces the hydrogen bonds by making the proton acceptors more negative and proton donors more positive upon formation of the dimer. As can be seen in Figure c, this favorable rearrangement of π density is almost twice as large when the proton donor and acceptors atoms are grouped together as in DDAA. This explains why the DDAA dimer has a stronger π orbital interaction ΔEπ than its tautomeric counterpart DADA-ADAD.

Recently, Hernández-Rodríguez, Rocha-Rinza, and co-workers were able to rationalize trends in hydrogen bond stabilities by considering the Brønsted-Lowry acid/base properties of the proton donor and acceptor groups.[28,29] In this acidity–basicity interplay (ABI) model, the hydrogen strength increases with the acidity of the proton donor groups and the basicity of the proton acceptor groups. Because a larger accumulation of charge is associated with stronger acidities and basicities, our results explain their findings from a molecular orbital point of view.

Finally, we comment on dimers with all hydrogen bond donor groups on one monomer and all hydrogen bond acceptor atoms on the other monomer (i.e., A–D where n ≥ 2). The SEI model considers these systems to be the most stable because the number of attractive SEIs has been maximized. Indeed, AAA-DDD and AAAA-DDDD systems with exceptional strong binding strengths have been reported.[16−19,22] However, because these monomers contain either all hydrogen bond donor atoms or all hydrogen bond acceptor atoms, their hydrogen bond properties are uniquely tunable. Unsurprisingly, the strongest AAA-DDD[18] and AAAA-DDDD[17] dimers reported in literature have a +1 charge on their donor monomer D, which enhances both the electrostatic and orbital interactions tremendously (Supporting Discussion 1). The holy grail in the quest for finding extremely strong hydrogen-bonded arrays is therefore not the maximization of the attractive SEIs but the maximization of a favorable charge accumulation that enhances the binding strength. This tunability of A-D complexes can be exploited to find new polymer building blocks.

In summary, the grouping of electron-donating and electron-withdrawing groups in DDAA systems results in a larger accumulation of charge than for DADA systems. The monomeric charge accumulation results in (1) an enhanced electrostatic interaction and (2) an enhanced orbital interaction due to the decreased σ-HOMO–LUMO gap. Upon formation of the dimer, the hydrogen bonds are further enhanced by π resonance assistance. These charge accumulation effects also explain the energy difference between GC and PU (i.e., the dimers on which the SEI model was originally based).

GC and PU Dimers

The GC pair has an interaction energy of −33.7 kcal mol–1, while the interaction energy of PU is −21.5 kcal mol–1 (Figure a). GC is thus 12.2 kcal mol–1 more stable than PU. To understand the origin of this difference, we again applied the energy decomposition analysis (EDA) scheme to the dimers with constrained hydrogen bond distances. This was done by varying the middle hydrogen bond distance over an interval from 3.1 to 2.7 Å with 0.01 Å per step while keeping the monomers frozen in the same geometry that they have in the fully optimized dimer. The essential results are graphed in Figure b; the complete data set is given in Supporting Table 3.

The GC dimer is still more stable than the PU pair by up to 13.0 kcal mol–1, even with similar hydrogen bond distances (Figure b). This enhanced stabilization comes from both the electrostatic interaction ΔVelstat and orbital interaction ΔEoi, which are more attractive for GC by up to 6.8 and 6.1 kcal mol–1, respectively. Both the ΔEσ and ΔEπ orbital interactions contribute around 3.0 kcal mol–1 to the enhanced ΔEoi. The Pauli repulsion ΔEPauli is approximately the same for both dimers.

Because the GC dimer has two attractive SEIs, it has two proton donors and acceptors grouped together. Just as in the DDAA monomer, this leads to a stronger accumulation of charge around the frontier atoms in G and C than in the P and U monomers with alternating donors and acceptors (Figure c). The monomeric charge accumulation in GC enhances the primary hydrogen bonds by making the proton donors more positive and the proton acceptors more negative and furthermore allows for favorable secondary electrostatic interactions as identified by the SEI model.

The stronger accumulation of charge in G and C is also responsible for the stronger σ-orbital interaction in GC, which becomes evident from their orbital energies (Supporting Figure 7). The occupied virtual energy gaps for the upper, middle, and lower hydrogen bonds in GC are 4.8, 4.9, and 4.3 eV, respectively, while they are 5.8, 5.7, and 6.6 eV for PU. The smaller gaps result in a stronger orbital interaction and thus a stronger overall interaction strength in GC.

Again, both dimers are further stabilized by polarization in the π-electron system by making the proton donors more positive and proton acceptors more negative upon formation of the dimer. Because this favorable rearrangement of π charge is more pronounced in systems with grouped proton donors and acceptors (Supporting Figure 8), the π-orbital interaction is stronger for GC than for PU.

Just as for the quadruple hydrogen-bonded DDAA and DADA dimers, the enhanced binding strength for GC is thus determined by a favorable charge accumulation around the frontier atoms. Because counting the number of attractive SEIs gives an indication of the charge accumulation in the monomers, the trends in hydrogen bond strengths are often in line with the SEI model’s predictions.

Finally, we emphasize that hydrogen bonds are a complex interplay of many bonding components, including electrostatics, covalency, steric repulsion, and π-resonance assistance. The importance of each component is system-dependent, which is very difficult (if not impossible) to capture by an easy-to-use predictive model. For example, doubly hydrogen-bonded mismatched DNA base pairs GG and CC[40] have the same number of SEIs but still vary significantly in hydrogen bond strength, which follows entirely from the difference in steric repulsion. These subtle effects might be captured by using state-of-the-art quantum chemical software and can assist supramolecular chemists in understanding and controlling the properties of self-assembled systems.

Conclusions

The secondary electrostatic interaction (SEI) model oversimplifies the hydrogen-bonding mechanism by viewing it as interacting point charges. Nevertheless, experimental binding strengths are often in line with the model’s predictions. We have therefore proposed a new view of the rationalization of hydrogen bond stabilities that is both chemically intuitive and grounded on quantum chemical insight.

Our dispersion-corrected DFT computations on tautomeric quadruple hydrogen-bonded DDAA and DADA homodimers show that the SEI model is often predictive because it provides a measure of favorable charge accumulation in the monomers. When both hydrogen bond donors (which are electron-donating) are grouped on one side while both hydrogen bond acceptors (which are electron-withdrawing) are grouped on the other side, there is a stronger accumulation of charge around the frontier atoms. This monomeric charge accumulation results in (1) a more favorable electrostatic interaction and (2) an enhanced σ-orbital interaction due to reduced occupied virtual orbital energy gaps. Upon formation of the dimer, the hydrogen bond energy is further enhanced by π-resonance assistance, which is stronger in systems with the proton donor and acceptors grouped. These insights also enabled us to explain the hydrogen bond strengths of the GC and PU dimers on which the model is originally based.

Counting the number of attractive and repulsive SEIs, which involves analyzing the alternation of hydrogen bond donors and acceptors on a monomer, gives an indication of the charge accumulation in the monomer and explains why relative hydrogen bond stabilities are often in line with the model’s predictions. Exceptions to the predicted trends remain because the hydrogen bond strength is always determined by an interplay of electrostatic, covalent, steric, and other interaction components (Figure ), and the relative importance of these components is system-dependent.