Molecular Tailoring Approach for the Estimation of Intramolecular Hydrogen Bond Energy.

Hydrogen bonds (HBs) play a crucial role in many physicochemical and biological processes. Theoretical methods can reliably estimate the intermolecular HB energies. However, the methods for the quantification of intramolecular HB (IHB) energy available in the literature are mostly empirical or indirect and limited only to evaluating the energy of a single HB. During the past decade, the authors have developed a direct procedure for the IHB energy estimation based on the molecular tailoring approach (MTA), a fragmentation method. This MTA-based method can yield a reliable estimate of individual IHB energy in a system containing multiple H-bonds. After explaining and illustrating the methodology of MTA, we present its use for the IHB energy estimation in molecules and clusters. We also discuss the use of this method by other researchers as a standard, state-of-the-art method for estimating IHB energy as well as those of other noncovalent interactions.

1. Introduction

The hydrogen bond (HB) is a dominant noncovalent interaction found in chemical and biological systems [1,2,3,4]. The term “hydrogen bond” seems to have emerged around 1930, from the works of Pauling [1] and Huggins [5,6]. However, the mention of weak, yet specific interactions involving the hydrogen atom is much older. The dimeric association of molecules with hydroxyl groups was suggested by Nernst in 1892 [7]. The term “Nebenvalenz” by Werner [8] and “weak union” by Moore and Winmill [9] are other early stipulations of this noncovalent interaction. In 1920, Latimer and Rodebush [10] suggested that the hydrogen nucleus in an aqueous solution of amines is held jointly by two octets, constituting a weak bond. Barnes, while studying the structure of ice [11], suggested that the hydrogen atoms were midway between the two oxygen atoms, though he did not explicitly mention the hydrogen bond. Huggins [12] claimed that he was the first to propose the term “H-bond” in 1919. His later usage of the term “hydrogen-bridge” may have led to the German word “Wasserstoffbrücke.” The concept of HB gained popularity after Pauling published his classic book, The Nature of the Chemical Bond in 1939 [1]. Pimentel and McClellan [13] suggested that an HB exists when (i) there is evidence of a bond and (ii) there is evidence that this bond involves a hydrogen atom already bonded to another atom. The recent definition of HB by IUPAC [14] is similar in spirit to that in Ref. [13]. The former [14] states that “the hydrogen bond is an attractive interaction between a hydrogen atom from a molecule or a molecular fragment X–H in which X is more electronegative than H, and an atom or a group of atoms in the same or a different molecule, in which there is evidence of bond formation.”

An HB may be generally repren class="Gene">sented as X–H···Y, where X–H is the proton donor and Y is a proton acceptor. The X–H···Y interactions such as O–H···O, N–H···O, N–H···N, S–H···O, etc., in neutral molecular systems, exhibit interaction energies lying between ~1 to 20 kcal/mol. Typical H···Y distances suggested in the literature fall in the range of ~1.2–3.0 Å and X–H···Y angles lie between 100 and 180° [2,3,4]. The HB’s in liquid water are central to water’s life-providing properties [1,15]. It is stipulated in the literature [16] that if HBs in water were 7% stronger or 29% weaker, water would not be a liquid at room temperature. HBs provide a significant driving force for the native structures and functions of biomolecules [17,18]. Hence, it is of great importance to reliably estimate these X−H···Y HB strengths for shedding light on several physicochemical phenomena and life processes.

The theoretical estimation on class="Chemical">f intermolecular X–H···Y HB strength in a complex A···B is routinely performed using a supermolecular approach, in which the HB energy (EHB) is estimated as EHB = EA···B − (EA + EB). Several methods for estimating the intermolecular interaction energies are reported in the literature (see, e.g., Refs. [19,20,21,22]). On the other hand, quantifying an intramolecular hydrogen bond (IHB) strength is not as straightforward as the intermolecular one. The main difficulty lies in isolating the X–H···Y interaction present within a molecule than in a dimer or a complex. Many studies for gauging the strength of the IHB in the literature are based on spectroscopic- [23,24,25,26] and electron density topological approaches [27,28,29,30]. Nevertheless, some empirical, semiempirical, and ab initio procedures [31,32,33,34,35,36,37,38,39] have also been reported in the literature for estimating the IHB energy. These have been nicely summarized by Jablonski [40] in his article in this Special Issue, and we shall discuss only the aspects of these methods (e.g., their merits and demerits) that are not explicitly covered in Ref. [40].

One of the early approaches is the conformational analysis (CA), in which two different conformers of the reference molecule are considered. These conformers are chosen such that the HB is kept intact in one of the conformers and is broken in another. The energy difference between these two conformers is then taken as the measure of IHB energy [32,40,41]. A significant disadvantage of this method is that the estimated IHB energy is erroneous due to the incorporation of attractive (syn-anti) or repulsive (anti-anti) additional interaction in one of the conformers [42]. In another similar procedure, viz., the ortho-para method [43], the IHB energy of the X-H···Y bond formed by two substituents, which are ortho to each other, is taken as the energy difference between the ortho and para forms of the reference molecule. However, this method applies only to aromatic systems in which an HB is present in two substituents, which are ortho to each other. The main drawback of this method is that it assumes that the electronic effects caused by the substituent at different positions are similar between the ortho- and para-conformers [44,45,46].

Yet, another indirect approach is the isodesmic/homodesmic reactions. In the former, the IHB making/breaking reaction is written so that the number and type of bonds on either side of the reaction are equal, except the HB, which is retained in one of the reactants [33,47,48]. A further assumption is that the atomic hybridization is conserved on both sides of the reaction. In that case, the method is called homodesmic reaction, which is supposed to give more reliable energy estimates than the isodesmic reaction [49,50,51]. The main drawback of the isodesmic/homodesmic reaction approach is that it does not give HB energy but includes strain energy due to the formation of a ring structure [52]. Another major disadvantage of these indirect methods is that they are applicable only to the evaluation of energy of a single HB present in the system and cannot be employed in a system containing multiple HBs.

Another popular but indirect method is based on the quantum theory of atoms in molecules (QTAIM) [53]. In this method, the presence of a (3, −1) bond critical point (BCP) of the molecular electron density (MED) between H···Y, is considered as the signature of an HB. The large/small value of the MED at the BCP is seen to correlate with strong/weak X–H···Y interaction [54,55,56,57,58]. Espinosa et al. [59] proposed an empirical relation, EHB = 0.5 V(rcp), where V(rcp) is the potential energy density at BCP. Interacting quantum atoms (IQA) [60] framework, leading to QTAIM-compatible energy partition, is another indirect approach wherein the HB energy is calculated as the sum of the classical Coulombic interaction between groups involved in the HB and the exchange-correlation energy. It has been pointed out that there is no check on the reliability of HB energy provided by both of these methods [61,62]. Further, these empirical equations are applicable only when a (3, −1) BCP is present. For instance, a (3, −1) BCP at O–H···O bond is conspicuous by its absence in all the polyols having an O–H···O interactions between the vicinal -OH groups. However, the weak O–H···O HB was confirmed in ethylene glycol-based on the vapor phase OH-stretching overtone spectroscopy [63,64,65]. Some other variants and indirect empirical equations proposed in the literature for estimating the strength of IHBs are summarized in Ref. [40]. Hence, we skip the discussion of these methods.

A brief review on class="Chemical">f the above literature suggests that these empirical and indirect methods for estimating IHB energy are generally limited to singly H-bonded systems. These cannot be readily extended to a system containing multiple HBs and hence also to the estimation of HB cooperativity due to an interconnected network of HBs. Most importantly, there is no check on the reliability of the estimated HB energies. Thus, it was felt necessary to come up with a direct theoretical method for a reliable estimation of IHB energy. Deshmukh and Gadre [66,67] proposed such a procedure, based on the in-house developed molecular tailoring approach (MTA) for the ab initio treatment of large molecular systems [68,69,70,71,72,73,74,75,76,77,78,79,80,81,82]. The MTA currently enables the calculation of one-electron properties, geometry optimization, and the calculation of vibrational infrared and Raman spectra of large molecules/clusters using DFT or correlated methods. Before discussing the application of MTA for the IHB energy estimation, we explain the working principle of MTA, with an illustrative example, in Section 2.

2. Molecular Tailoring Approach

The molecular tailoring approach (MTA) is a fragmentation-based technique developed by Gadre et al. [68,69,70,71,72,73,74,75,76,77,78,79,80,81,82]. Within MTA, a spatially extended molecular system under consideration is partitioned into a set of overlapping fragments (called the “main fragments”) on which ab initio calculations for one-electron property or the energy are carried out. The fragmentation may be carried out automatically or manually. The quality of the fragmentation scheme is gauged by a parameter called R-Goodness (Rg), which may be estimated as follows: Put a sphere of radius R centered on a reference atom i so that all the atoms of the parent system lying within this sphere belong to at least one of the main fragments. The maximum value of R obeying this condition is called the Rg value of atom i in the given fragmentation scheme. The minimum of such atomic Rg values is called the Rg value of the scheme. In general, the larger the Rg value of a fragmentation scheme, the better is the chemical environment of each atom mimicked [69]. After choosing an appropriate fragmentation scheme, molecular energy E, of the spatially extensive parent system is estimated approximately by patching those of the individual fragment energies [69] using Equation (1). where the energy E of the parent molecule is estimated as the sum of energies of primary fragments { minus the sum of energies of binary overlap fragments { plus the sum of energies of ternary overlap fragments, etc. Here, k stands for the degree of overlap., e.g., for binary overlap, k = 1, for ternary overlap, k = 2, etc. Equation (1) is generalized for estimating an electronic property of the molecule, such as the energy gradients, the Hessian matrix elements, etc.

We now illustrate the n class="Chemical">fragmentation procedure with a test example, viz., the α-tocopherol molecule (shown as M) in Scheme 1, fragmented into three primary fragments F1, F2, and F3 (shown by appropriate circles). The fragments F4 and F5 are the binary overlaps of fragments F1, F2, and F2, F3, respectively. Here, the term binary overlap means the common structural part of two primary fragments. In fragmentation Scheme 1, the ternary fragment (overlap of three fragments F1, F2, and F3) is absent. The valencies of the cut regions are satisfied by placing the H-atoms at the appropriate C–H distance of 1 Å along the cut C–C bond. The calculations for the single point energy (or the property) of these fragments F1 to F5 are performed. For this test case, we report the MTA energy calculation of M at HF/6-31+G(d,p) level theory using Equation (1). It should be noted that the HF method employed here is only for illustrative purposes. MTA method works at any correlated level of theory. All the calculations are performed with the Gaussian 16 package [83].

Scheme 1

Illustration of the molecular tailoring approach (MTA)-based fragmentation procedure for α-tocopherol, shown as the parent molecule, M. See text for details.

In the present case, the energy (EMTA) of the parent molecule (M) is obtained by Equation (2) as E

In the present case of α-tocopherol this energy is calculated utilizing the energies of the fragments as EMTA = {(−769.56869) + (−469.57421) + (−391.51579)} − {(−118.26152) + (−235.36444)} = −1277.03273 hartrees (a.u.). The actual energy (EFC) at the HF/6-31+G(d,p) level of theory is EFC = −1277.03288 a.u. Thus, the error (EError) in the estimation of molecular energy is given EError = EFC − EM = −0.00016 a.u. This error can be further reduced using the so-called grafting procedure embedded in the current version of MTA [77,78].

Gadre et al. n class="Chemical">first proposed the MTA methodology for estimating the electrostatic properties of large, closed-shell molecules [68]. However, in the last 24 years, the scope of MTA was extended for the estimation of molecular energy [69], geometry optimization [70,71], the estimation of the Hessian matrix [72], the computation of vibrational infrared [73] and Raman spectra [74], and binding energies of large molecular clusters and complexes [75,76,77,78,79,80,81,82]. Since the fragment computations are independent of each other, MTA has the advantage that the energy computation of the parent molecule is intrinsically parallel. Further, MTA can currently work with Gaussian, GAMESS, and NWChem at the back-end, thereby becoming a powerful tool for ab initio treatment of large molecules and clusters, when used in conjunction with a high level of theory, such as MP2 or CCSD(T) employing a large basis set. One important application of MTA is in estimating the IHB energy [52,66,67,84,85,86,87,88,89,90,91,92]. This will be discussed in the next section.

3. Intramolecular Hydrogen Bond Energy Estimation by Molecular Tailoring Approach

With the above brief introduction to MTA, we now discuss its application for estimating the IHB energy. As discussed in the introduction section, intermolecular X–H···Y HB energy in a complex A···B is estimated as EHB = EA···B − (EA + EB). This estimation is possible because the energies of the two monomers A and B can be separately calculated. In the case of intramolecular X–H···Y interaction, such a separation is, in general, difficult. However, the MTA procedure allows the generation of fragments so that the atoms/functional groups involved in the HB formation are parts of two different fragments.

The fragmentation procedure is illustrated in Scheme 2 for the test molecule of 1,2,4-butanetriol. In Scheme 2, the parent test molecule, denoted as M, is shown at the center. The geometry of 1,2,4-butanetriol was optimized at the MP2/6-31+G(d,p) (default option: frozen core) level of theory using the Gaussian package [83]. The energy of the optimized structure is −383.01926 a.u. at MP2/6-31+G(d,p) level. The three oxygen atoms are shown as O1, O2, and O3 (see Scheme 2), with the two HB’s, viz., HB1 (O2-H···O1) and HB2 (O3-H···O2) whose energy is to be estimated. For this purpose, the parent molecule is cut into three primary fragments F1, F2, and F3, obtained by replacing −O1H, −O2H, and −O3H groups, respectively, with an H atom each. Dotted circles show these cut regions on the original molecule. The H-atoms are added along the respective C–O bonds (which are cut to form these primary fragments) so that the C–H distance is 1 Å. Hydrogen is the simplest monovalent atom that can be used for satisfying the valencies of cut regions. It is emphasized here that H-atoms placed at slightly different distances (say at 0.9 or 1.1 Å) from the C-atom do not change the results appreciably. This is because of the cancellation of errors in estimating the molecular energy using these fragments. Fragments F4, F5, and F6 are obtained by taking the binary intersection of these primary fragments, i.e., (F1∩F2), (F2∩F3), and (F1∩F3), respectively. Here, intersection means the common structural parts between two primary fragments apart from added H-atoms. For instance, in fragment F4, C1(H2)–C2(H2)–C3(H2)–C4(H2)O3(H) is the common structural part that is also present in fragments F1 and F2. Similarly, fragment F7 is the common intersection of three primary fragments F1, F2, and F3, i.e., (F1∩F2∩F3). A single point energy evaluation, at MP2/6-31+G(d,p) level of theory, is carried out on all seven fragments obtained by the above fragmentation procedure. The fragment geometries are not optimized to avoid the conformational changes in them so that they lead to reliable estimates of IHB energies. It is necessary first to provide a check on MTA application to the parent molecule, M. As discussed above, the actual energy of the original molecule (M) is EM = −383.01926 a.u. at the MP2 level of theory. Using the MP2 single point energies of these fragments, the estimated molecular energy of M is: EM = {EF1 + EF2 + EF3} − {EF4 + EF5 + EF6} + EF7 = {−307.97304 + (−307.9642) + (−307.97765)} − {−232.92410 + (−232.92514) + (−232.93273)} + (−157.88552) = −383.01844 a.u. The error, ΔE = |MTA energy - actual energy| in molecular energy indeed turns out to be very small, viz., 0.00082 a.u. This excellent agreement between the MTA-estimated and actual energy suggests that the present fragmentation scheme is reliable for evaluating HB energies.

Scheme 2

Fragmentation procedure for estimating the energies of the H-bonds, HB1, and HB2 in 1,2,4-butanetriol (Parent M) molecule. See text for details.

Now we estimate the energy of two n class="Chemical">hydrogen bonds HB1 and HB2, in the parent 1,2,4-butanetriol. Recall that the hydroxyl groups involved in the formation of hydrogen bond HB1 are O1–H and O2–H. These hydroxyl groups were replaced in fragments F1 and F2, respectively, by H-atoms. Putting the geometry of fragment F1 over F2, we regenerate the parent molecule except following two things: (i) the O–H···O H-bond, i.e., the HB1 interaction between O1–H and O2–H present in the parent molecule is missed out and (ii) there is double counting of common structural part between F1 and F2 (viz., the secondary fragment, F4). Upon addition of single-point energies of fragments F1 and F2, followed by subtraction of the energy of fragment F4 would give the energy of the parent molecule except that the energy of the HB, viz., HB1 is missed out. If the energy of the parent 1,2,4-butanetriol EM is subtracted from (EF1 + EF2 – EF4), the HB energy E is obtained as E = (EF1 + EF2 − EF4) − EM = 3.84 kcal/mol. In a similar fashion, E is obtained as E = (EF2 + EF3 − EF5) − EM = 1.60 kcal/mol. It should be noted here that these estimated HB energies are in the gas phase. However, the MTA-based method in principle can provide HB energies in the solvent phase, wherein the energies of the fragments in solvent (using continuum solvation model) could be employed.

We note that the two HBs, HB1 and n class="Gene">HB2, are interconnected, forming an H-bond network. Such networking of H-bonds leads to a phenomenon called cooperativity [67]. In general, it is anticipated that the strengths of HB1 and HB2 are enhanced because of this networking effect. To estimate the contribution of cooperativity toward each of these two H-bonds, we reestimated the HB energy of these two HBs by isolating them from each other. The difference between the HB energy estimated earlier (in the presence of network) and the one when they are isolated (in the absence of a network) is the cooperativity contribution toward this HB. For example, consider fragment F3 in which only HB1 is present and fragment F1 in which HB2 is present. To estimate the energy of HB1 in the absence of the networking effect of HB2, we consider fragment F3 as our parent molecule. In the present case, fragments F5 and F6 are the two primary fragments that, when placed over each other, would give us the parent fragment F3 except HB1, and fragment F7 is the binary overlap of F5 and F6. Therefore, utilizing these fragments’ energies, the energy of HB1 is obtained as EHB1 = (EF5 + EF6 − EF7) − EF3 = 3.32 kcal/mol. Similarly, the energy of HB2 in the absence of the networking effect of HB1 is obtained as EHB2 = (EF4 + EF6 − EF7) − EF1 = 1.09 kcal/mol. These reestimated HB energies are indeed smaller than those estimated in the presence of the networking effect. The difference in the energy is cooperativity contribution. The cooperativity contribution to HB1 is = 3.84 − 3.32 = 0.52 kcal/mol and that for HB2 is = 1.60 − 1.09 = 0.51 kcal/mol. In the present test case, the estimated cooperativity contributions are not large because only two HBs are present. The later sections will show that cooperativity values in some molecules can indeed be as large as a typical HB energy.

The HB energies obtained by applying the above procedure to some n class="Chemical">alkanetriol molecules are shown in Table 1 [66]. The estimated HB energies fall in a range between 1.50 and 4.97 kcal/mol (see Table 1). This is the expected energy range from chemical intuition. Further, these HB energies are in a qualitative agreement with those expected from the corresponding HB distances. For instance, the strongest HB in 1,2,5-pentanetriol has an energy of 4.97 kcal/mol, with the corresponding HB distance being the shortest (1.80 Å) among all the alkanetriols reported in Table 1. One of the noteworthy results in Table 1 is that the error in estimating molecular energies of all the alkanetriols is quite small, viz., between 0.40 to 0.65 kcal/mol. By considering this accuracy, we estimate the maximum error associated with our calculated HB energies to be 0.3 kcal/mol. The present method is thus capable of calculating accurately the IHB energies and cooperativity values of multiply H-bonded systems. This is a significant advantage of the current method over the other indirect approaches reported in the literature.

Table 1

The H-bond (HB) distances (in Å), HB energies (in kcal/mol), and the error in the molecular energy estimation for alkanetriols using similar fragments, ΔE=|EM − Ee|. The corresponding O–H stretching frequencies (cm−1) and the molecular electron density (MED) value at the (3, −1) bond critical point (BCP) (a.u.) are also shown. The calculations are performed at MP2 (full)/6-311++G(2d,2p) level theory.

Molecule	HBLabel	HBDistances a(in Å)	HB Energy(kcal/mol)	ΔE(kcal/mol)	O–H Stretch Frequency(cm−1)	MED at(3, −1) BCP (a.u.)
1,2,3-propanetriol b	HB1HB2HB3	2.162.082.58	1.902.471.63	0.50	378437653845	0.0201
1,2,3-butanetriol b	HB1HB2HB3	2.132.052.58	2.132.721.60	0.50	376837453844	0.0211
1,2,4-butanetriol	HB1HB2	1.982.22	2.901.75	0.40	378938283875	0.0219
1,2,5-pentanetriol	HB1HB2	1.802.25	4.971.78	0.55	366938253865	0.0334
1,3,5-pentanetriol	HB1HB2	1.941.96	2.912.90	0.58	376337923875	0.02250.0239
2,3,4-pentanetriol b	HB1HB2HB3	2.122.022.56	2.182.941.50	0.52	375937313820	0.0223
2,4,6-heptanetriol	HB1HB2	1.921.93	3.022.94	0.65	375337733857	0.02500.0242

a The MP2 (FC)/6-311++G(2d,2p) optimized geometries were employed. Table 1 is partially reproduced from our earlier study reported in Ref. [66]; Copyright (2006) The American Chemical Society. b The triols wherein three OH groups are present on the successive C-atoms show three H-bonds. See text for details.

4. Critical Comparison of MTA with Other Methods

We now compare the estimated HB strengths in then class="Gene">se molecules with those qualitatively estimated by other indirect measures. These measures include the O–H stretching frequency, molecular electron density (MED) value at (3, −1) BCP, and the HB energy estimated using the isodesmic reaction approach (IDRA) and that by using Espinosa’s equation. Both the MED value at (3, −1) BCP and the shift in the stretching frequency of the O–H involved in the HB show a good qualitative correlation with the estimated HB energies (see Table 1). For instance, the strongest HB found in 1,2,5-pentanetriol (4.97 kcal/mol) corresponds to the highest MED value (0.0334 au) at the (3, −1) BCP. The stretching frequency of the O–H bond involved in the HB is 3669 cm−1, showing a redshift. The weakest HB is seen in 2,3,4-pentanetriol (1.50 kcal/mol), with the (3, −1) BCP being absent and a large O–H stretch frequency of 3820 cm−1. In general, the calculated HB energies show a good qualitative rank–order relationship with the corresponding O–H stretching frequencies and the MED value at the respective (3, −1) BCP.

We now critically compare our MTA-based results with those obtained by yet another indirect method, viz., the isodesmic reaction approach (IDRA) [33,48]. In IRDA, the IHB-making/breaking reaction is written such that, except for the HB under consideration, the number and type of other bonds are conserved on both sides of the reaction. Within IDRA, all the reactant and product geometries are optimized at the appropriate level of theory. The energy change for such a reaction is taken as the HB energy. The main disadvantage of IDRA is that there is no unique way of writing an isodesmic reaction for estimating the HB energy. For example, in Scheme 3, four possible reactions are shown for evaluating the HB energy E in 1,2,5-pentanetriol, which retain HB1 on both sides of these isodesmic reactions. See Ref. [52] for details of several other isodesmic reactions for estimating the HB energies of HB1 and HB2.

Scheme 3

Some possible isodesmic reactions for the estimation of H-bond energy, E in 1,2,5-Pentanetriol. See text and Ref. [52] for details.

Figure 1 shows a histogn class="Chemical">ram of the estimated HB energy E by these four reactions and also by the MTA-based method. The HB energies estimated by the isodesmic reactions vary significantly across the levels of theory and from each other. These estimated energy values are much smaller (30 to 40%) than the respective MTA-based ones. The reasons for these smaller HB energy values by IDRA can be understood as follows: In IDRA reactions presented in Scheme 3, on the reactant side, the HB2 bond formation between the -OH groups at C2 and C5 positions leads to a seven-membered ring-like structure involving one O-H···O, two C-O, and three C-C bonds. This ring formation has the ring strain effect in the parent 1,2,5-pentanetriol molecule. Since the reactant and product geometries are optimized, this ring strain effect is not preserved on the product side of these reactions. Therefore, the molecules on the product side are more relaxed, losing most of their ring strain due to the loss of the HB2 bond. It may be noted here that the ring strain may be small or canceled out when one estimates the energy of HB1 using IDRA. This is because the formation of HB1 leads to the formation of a five-membered ring-like structure involving one O-H···O, two C-O, and only one C-C bonds. This ring (five-membered) backbone is expected to be maintained (to some extent if not fully) on both reactant and product sides as it involves only one C-C bond. For more details about HB1 energy by IDRA, see Ref. [52].

Figure 1

The H-bond energy, E, in 1,2,5-pentanetriol calculated at different levels of theory using isodesmic reactions (see Scheme 3) and also by molecular tailoring approach (MTA). See text for details.

Thus, the finer bonding en class="Chemical">ffects are not mimicked evenly on both sides of the reaction, resulting in poor estimates of HB energy by IDRA. On the contrary, in the MTA method, the geometry of the fragments is not optimized, and the fragment backbone structure is preserved, leading to accurate estimates of the HB energies. For instance, the evaluation of HB energy E as (EF1 + EF2 − EF4) − EM in 1,2,4-butanetiols involves the energies of F1, F2, and F4 (see Scheme 2). Thus, the MTA-based method leads to unique isodesmic/homodesmic reaction, viz., M + F4 → F1 + F2. As seen in Scheme 2, the carbon backbone structure is retained on either side of this unique reaction. In other words, as advocated in the literature [49,50,51,52], IDRA does not yield the true HB energies. Another drawback of the IDRA is for estimating the HB strengths in multiply hydrogen-bonded systems. Here, one has to write different reactions for different HBs. In summary, the MTA-based method for the estimation of IHB energy is accurate and can be applied to the evaluation of multiple IHB energies in a given molecule. In the following sections, we present and discuss the results of IHB energies obtained by applying the MTA-based method to a variety of systems.

We now present a comparison of the HB energy estimated by the MTA-based method with those by Espinosa’s empirical relation [59]. For this purpose, we consider two molecular systems having O–H···O=C HB, viz., 2′-Hydroxyacetophenone, and methyl 2-hydroxybenzoate. The geometries and HB energies by Espinosa’s method in these two molecules were taken from the Ref. [93]. The reported HB energies at B3LYP/6-311++G(d,p) level, by Espinosa’s method, in 2′-Hydroxyacetophenone and methyl 2-hydroxybenzoate are 15.3 and 11.4 kcal/mol, respectively. The HB energies for these molecules were calculated by us at the given geometries, at the same level of theory, by MTA-based method are 9.7 and 7.8 kcal/mol, respectively. As can be observed, the HB energies estimated by Espinosa’s method are significantly overestimated, compared to the HB energies by the MTA method. These results are in agreement with an earlier report that the use of Espinosa’s method significantly overestimated the energy of HBs [62,93]. Although our MTA-based method requires extensive calculations and more computational time, reliable direct estimates of the HB energies are thereby obtained. Further, the internal benchmarking of the total energy is possible for the MTA-based method. Such benchmarks are not available in Espinosa’s method.

5. Application to Large Molecules and Clusters with Multiple Hydrogen Bonds

We now discuss the application of the MTA-based method for IHB energy estimation in several large systems. One such interesting system is a class of carbohydrates, viz., sugar molecules. These molecules play an important role in biological processes, which are mainly governed by weak interactions such as hydrogen bonding, hydrophobic effects, etc. Hence, understanding the interactions such as H-bonding is of utmost importance. Further, the O–H groups in carbohydrates form a network of interconnected O–H ···O HBs. It is suggested that the strengths of these individual O–H···O HBs are enhanced due to such networking of HBs.

In our earlier work [67], the IHB energies and the contribution to cooperativity were investigated in eight aldopyranose monosaccharides, which vary in the position of hydroxyl groups [axial (ax) or equatorial (eq)], as shown in Figure 2. The estimated HB energies are in the range of 1.2 to 4.1 kcal/mol at the MP2(full)/6-311++G(2d,2p) level of theory. It is found that the OH···O eq-eq interaction energies are between 1.8 and 2.5 kcal/mol, with axial-equatorial ones being stronger (2.0 to 3.5 kcal/mol). The strongest bonds involve nonvicinal ax-ax O–H groups (3.0 to 4.1 kcal/mol). The cooperativity contribution to the HBs is seen to fall between 0.1 and 0.6 kcal/mol for eq-eq HBs and is seen to be higher (0.5 to 1.1 kcal/mol) for ax–ax HBs. This work [67] was one of the first attempts for estimating the IHB energies and the respective cooperativity contributions in sugar molecules.

Figure 2

General structure of the aldopyranose sugar. In Table, ax represents axial, and eq represents equatorial orientations of the hydroxyl group at carbons C2-C4. Figure partially reproduced from Ref. [67] with the permission from American Chemical Society (ACS). Copyright (2008) The American Chemical Society.

A similar class on class="Chemical">f compounds having no ring oxygen atom is hexahydroxy–cyclohexane, called inositol. Inositol derivatives function as intracellular signal transduction molecules, playing an important role in biological processes. It is hence important to understand the structure and stability of various inositol conformers [84]. On applying the MTA-based method, the estimated HB energies in the presence of a cooperative HB network are seen to be in the range of 2.2 to 3.8 kcal/mol at the MP2/6-311+G(d,p) level of theory. The sum of all the HB energies in these conformers of inositol falls between 7.2 to 18.1 kcal/mol. The total cooperativity contribution in these conformers is rather large, between 2 to 5 kcal/mol. It increases on going from isomers with more eq O-H groups to those with more ax ones. Importantly, the highest stability of the scyllo isomer in the solvent was attributed [84] to weaker intramolecular OH···O HBs between eq hydroxyl groups. It is suggested that these weak OH···O HBs in scyllo isomer may facilitate favorable intermolecular interactions with solvent molecules. In contrast, the inositol isomers with ax O-H groups are involved in the formation of relatively strong HBs. Therefore, they are less stable due to large steric factors in the gas phase and unfavorable intermolecular interaction in the solvent phase.

The MTA-ban class="Gene">sed method was also employed for understanding the conformational stability of fructose [85]. The experimental rotational spectroscopic studies suggest that the molecules of fructose and ribose preferentially adopt the β-pyranose structure in the gas phase. It was noted that [85] the relative stability of different conformers of fructose in the gas phase could be explained in terms of three collective effects: (i) the sum of HB energies in a given conformer, (ii) the strain energy of a bare fructose ring, and (iii) the sum of anomeric stabilization (endo + exo) energies. It was concluded [85] that the small ring strain, sufficiently large sum of the IHB energies, and the higher stabilization due to anomeric interactions in β-fructo-pyranose makes it a conformationally locked predominant structure in the gas phase.

Another molecular system wherein the D-n class="Chemical">glucose units are joined to each other by 1–4 linkage is cyclodextrins (CDs), which are macrocyclic oligosaccharides. They possess a unique ability to entrap guest molecules in their cavities owing to their bucket/bowl shape. Such inclusion complexes are used in the pharmaceutical industry for a variety of formulations. The most commonly known CDs are α-, β- and γ-CDs which have six, seven, and eight glucose units, respectively (see Figure 3). The IHB’s between the secondary O–H groups of glucose units result in the formation of a smaller rim of the CD bowls, whereas the primary O–H groups form the larger ones. The strength of the IHBs is suggested to be an important factor governing the stability of the CDs [94,95,96]. Moreover, one would expect a larger cooperativity contribution due to the formation of a more extended network of HBs between different types of O–H groups. The estimated IHB energies obtained by using the MTA-based method [86] belonged to a wide range of 1.1 to 8.3 kcal/mol at the B3LYP/6-311++G(d,p) level, suggesting that strong HBs are seen in CDs. For the O6H···O6′ HBs, HB energies fall between 6.7 to 8.3 kcal/mol, and for O3H···O2′ HBs, they are between 3.3 to 5.5 kcal/mol and 1.9 to 2.8 kcal/mol for O2H···O3 HBs. The cooperativity contribution by the O6H···O6′ HB is larger (1.3 to 2.7 kcal/mol) than that for O3H···O2 (0.3 to 1.0 kcal/mol) and O2H···O3 (0.25 to 1.10 kcal/mol) HBs. Note that the cooperativity contribution in glucopyranose (0.1 to 1.1 kcal/mol) is much smaller than that in CDs. The higher HB strength in CDs could be one of the possible reasons for the much lower aqueous solubility of the natural CDs than that of comparable acyclic polysaccharides [97,98,99].

Figure 3

(a) A schematic structure of β-cyclodextrin molecule (containing seven glucose units) indicating different types of hydroxyl groups and (b) a CD bowl, secondary hydroxyl groups at smaller and primary ones at the larger rim, respectively. See text for details.

We now discuss the HB energies and cooperativity contributions calculated by using MTA in yet another macrocyclic system, viz., p-substituted Calix[n]arenes CX[n] (n = 4, 5), for exploring substitution effect on the strength of the HBs [87]. The estimated HB energies were between 4.6 to 8.2 kcal/mol, with cooperativity contribution being 0.2 to 2.6 kcal/mol, both calculated at B3LYP/6-311G(d,p) level of theory. The estimated HB energies showed [87] the following order: CX[n]-t-Bu ≈ CX[n]-NH2 > CX[n]-CH2Cl > CX[n] > CX[n]-SO3H > CX[n]-NO2 > S-CX[n]-t-Bu > S-CX[n]. It was also observed that the HB energies in CX[5] derivatives are larger than those in CX[4]. Additionally, as expected, large HB energy values and the respective cooperativity contributions were found in CX[n] hosts with electron-donating substituents.

Recently, we have widened the application of the MTA-based method for the estimation of individual O–H···O HB energy and their cooperativity contribution in water clusters [88]. Many works reported in the literature focus on the global minimum structure of water clusters of variable sizes [73,75,76,100,101,102,103]. In spite of a large number of studies, the nature of water clusters at the molecular level is not fully understood. For instance, the strengths of individual O-H···O HB interactions reported in the literature are mostly limited to the water–dimer. Larger clusters are expected to have a cooperativity contribution due to the networking of HBs therein. We recently applied the MTA-based method for the reliable estimation of individual O–H···O HB energies and their cooperativity contribution in small water clusters Wn, n = 3 to 8 (see Figure 4).

Figure 4

The MP2-optimized geometries of various water clusters (Wn). Figure reproduced from Ref. [88] with the permission from American Chemical Society (ACS). Copyright (2020) The American Chemical Society.

The fragmentation procedure for the estimation of O–H···O HB energy in Wn is similar to the one for molecules, discussed in the previous section. The difference is that no dummy atoms are needed here because no covalent bond is cut. In the present case, two water molecules involved in the formation of an O-H···O HB were removed for generating the two primary fragments keeping the other water molecules within the cluster intact (see Ref. [88]). The calculated HB energies in Wn, for n = 3 to 8, are in the wide range of 0.3 to 10.7 kcal/mol at the CCSD(T)/aug-cc-pVDZ level, with the respective cooperativity contribution being 1.2 to 7.0 kcal/mol. To check the reliability of the results, the sum of all the HB energies for a given cluster was added to the sum of monomers energies. The molecular energy of a water cluster thus estimated agreed well with the actual energy (typical error less than 8.3 mH), suggesting HB energies obtained by our MTA-based method [88] are reliable.

6. Applications to Biomolecules

We now discuss the application of the MTA-based method for the energy estimation of the N–H···O=C, N–H···N, N–H···S, N–H···O, C–H···N, and resonance-assisted O–H···O=C HBs in biomolecules. Quantitative estimates of the strengths of such individual HB interactions are scarcely available, although they play a vital role in the determination of the structure of biomolecules such as polypeptides/proteins. In order to gauge the strengths of the N–H···O=C HBs, a model tetrapeptide was taken as a test case [89]. It has two different N–H···O=C HBs, the fragmentation scheme being similar to that discussed earlier. In the present case, we remove a –(H)N–(O=C)- functional group of amino acids, involved in the formation of N–H···O=C HB to generate the primary fragments (see Ref. [89]). The MTA-based methodology is applied for the estimation of these IHB energies in five different substituted tetrapeptides of polyglycine abbreviated as GGGGG apart from the capped acetyl and NH2 groups. The five substituted tetrapeptides (viz., GAGGG, GVGGG, GLGGG, and GIGGG) in which the second amino acid residue is replaced by alanine (A), valine (V), leucine (L), and isoleucine (I), respectively, were considered (see Figure 5). The corresponding completely substituted tetrapeptides, AAAAA, VVVVV, LLLLL, and IIIII, were also employed with a view to addressing the effect of substituents on the strengths of the different types of N–H···O=C HBs (Chain 1 and Chain 2 in Figure 5). The estimated HB energies at B3LYP/D95(d,p) level were in the range of 4 to 6 kcal/mol in the partially and fully substituted polypeptides. These values were in concurrence with the geometric parameters and reflected the subtle substituents effects for the substituted polypeptides. The MTA-based procedure was thus considered to be applicable for the IHB energy estimation in polypeptides and it would be fascinating to apply it to an actual protein.

Figure 5

Hydrogen-bonded chains in the helical peptide structures. “X” represents the substituent at the second position. A, B, C and D are four hydrogen bonds. See text for details. Reprinted from Ref. [89] with permission from American Chemical Society. Copyright (2009) The American Chemical Society.

Recently, the n class="Chemical">MTA-based procedure was applied to another biologically important class of molecules, viz., porphyrin analogs called meta-benziporphodimethenes (1) and N-confused isomers containing γ-lactam ring (2 and 3) [90,91]. In γ-lactam-containing isomers 3 (not shown in Figure 6) the ring O-atom of lactam ring is down, whereas, in isomer 2, O-atom is up (see Figure 6). The substitution at meso sp (R′) and on the benzene at sp carbons (R1 to R5) may affect the strengths of N–H···N, N–H···S, N-H···O, C-H···N HBs, which, in turn, decide the stability of these conformers. Hence, the IHB strengths in substituted 1, 2, and 3 were determined by the MTA-based methodology. The estimated N–H···N HB energies were between 11.0 to 15.6 kcal/mol at B3LYP/6-311+G(d,p) level, the individual bifurcated N–H···N interactions being ~6.0 kcal/mol. The N-H···O (~10.0 kcal/mol) and N–H···S (11.2 kcal/mol) HBs were found to be weaker than the N–H···N HBs, with the C-H···N HBs being the weakest (0.1 to 4.4 kcal/mol) [90]. The substituent effect on IHB strength was also investigated [90,91].

Figure 6

Structures of meta-benziporphodimethene 1 and N-confused meta-benziporphodimethene containing γ-lactam ring isomer (O-up, 2). See text for details.

It has been suggested that intramolecular O–H···O and resonance-assisted (RA) O–H···O=C HBs may have a pronounced effect on the antioxidant activity of the natural products [104,105,106,107]. Hence, an attempt was made to obtain reliable estimates of IHB energies in these molecules [92]. The presence of an RA IHB was seen in all the antioxidant molecules studied. In some molecules, this type of HB was found to be in the nature of two bifurcated HBs along with IHBs such as O–H···OCH3, O–H···OH, and O–H···O(ring) also being present. The energies of O–H···OCH3 HBs were found to lie between 2.2 to 2.4 kcal/mol, with the O–H···O (ring) HB energies being much smaller (~0.5 to 0.6 kcal/mol) at the MP2/6-311G(d,p) level. The energy of RA O–H···O=C IHB being largest (10.0 to 17.9 kcal/mol.) Recently, Restropo et al. [108] suggested that the weak HBs are directly related to the antioxidant properties of these molecules. The reliable estimates of the strengths of these weak HBs by the MTA-based method could be useful for quantitative structure–activity relationship (QSAR) studies on a larger set of antioxidant molecules.

7. Use of the MTA-Based Method by Other Researchers

We summarize here the use on class="Chemical">f our MTA-based method for the estimation of IHB energy by other researchers [109,110,111,112,113,114,115,116,117]. In an early use of the method, Lopes Jesus et al. [109] studied the 65 local minima conformers of 1,4-butanediol molecule. The H···O IHB energies in the two lowest energy conformers (c1, c2) were estimated using the MTA method to be 5.5 and 4.1 kcal/mol, respectively. The IHB energies were also obtained by two other methods, viz., conformational analysis (CA) [33,41,42] and use of the empirical Iogansen equation (IE) to spectroscopic data [110]. The IHB energy values turned out to be CA: (4.5 and 3.1 kcal/mol) and IE: (4.1 and 2.2 kcal/mol), in good qualitative agreement with their MTA counterparts.

Rusinska-Roszak et al. extensively used the n class="Chemical">MTA-based method [111,112,113]. In an early work [111], they estimated the O–H···O=C IHB energies in several aliphatic systems. However, their fragmentation scheme consisted of two appropriately made primary fragments and an overlapping fragment. The IHB energy was estimated, by using the energies of these three fragments and the actual energy of the parent system, ignoring the ternary and quaternary overlap fragments. The estimated IHB energies at MP2/6-311++(2d,2p) level in the saturated hydroxyl compounds were between 1.4 to 7.0 kcal/mol and 13.7 kcal/mol in the unsaturated ones [111].

In another work [112], they estimated the O−H···O=C IHB energies in 299 structures of hydroxycarbonyl aliphatic compounds involving resonance-assisted HB. The estimated IHB energies showed a wide range, from 8.2 to 26.3 kcal/mol. The HB energies showed a good qualitative correlation with other indirect measures of HB strength, viz., geometrical parameters, the O–H stretching frequencies, and the MED values at the BCP [112].

These authors [113] also applied the n class="Chemical">MTA-based method for the estimation of Ar–O–H···O=C HBs in mono-, di-, and triphenols substituted (by electron-donating/withdrawing groups) at the ortho- position by carbonyl-containing functional groups. The estimated HB energies for phenolic O–H···O=C (six-membered ring) fall in the range of 5.4 to 15.4 kcal/mol at MP2/6-311++G(2d,2p) level. These HBs energies are smaller (4.6 to 9.6 kcal/mol) when the O=C group is a part of seven or eight-membered rings. These HBs’ energy range is similar to corresponding HBs involving saturated carbonyl substituted alcohols.

Very recently, Afonin et al. [114] applied our method for estimating the energy of push–pull effect in β-diketones. In this push–pull system, the intramolecular charge transfer (ICT) occurs as a result of interaction between the π-donor and -acceptor parts, joined by a π linker. Further, the IHB is also present in the Z-conformation of β-diketones. Here, basic idea is to estimate the π component of conjugation energy in these systems. For this purpose, an E conformation of the parent β-diketones was considered wherein this HB is not present. The molecule in this configuration was fragmented using the MTA-based method. The donor and acceptor groups involved in ICT interaction were separated into two primary fragments, F1 and F2, and the overlapping fragment, F3, is the conjugation unit connecting these groups. The energy expression for the π conjugation energy is similar to that discussed in Section 3, viz., E(MTA) = [(EF1 + EF2 − EF3) − EM], EM being the energy of the parent molecule in E configuration. The effect of electron-donating and -withdrawing groups on the π component of conjugation energy was also estimated. The authors stated [114] that “although the choice of the fragmentation scheme and the computational protocol used did not play a decisive role in estimating the energies of the same IMHB, this issue also needs to be investigated when estimating the conjugation energy.”

Afonin et al. [115] applied the methodology for the quantitative decomposition of RAHB energy in β-diketones into resonance and hydrogen bonding components. They also compared the estimated HB energy by the MTA method with that obtained by the functional-based approach (FBA). In FBA, the HB energy is written as an empirical function of HB descriptors, EHB = f(D), the parameter D is one of the H-bond descriptors (i.e., geometrical, topological, and/or spectral characteristics of the H-bond). For details of these empirical FBA equations, see Ref. [115]. It has been shown by these authors that the FBA method evaluates “the component of RAHB interaction corresponding to the energy of the pure H-bond without resonance component.” However, the MTA method implicitly takes into account the π component in the RAHB interaction and hence “the difference in the energy of the IMHB as evaluated by means of the MTA and FBA yields a quantitative estimate of the resonance component in the case of the resonance-assisted hydrogen bonds” [115]. The resonance component energy was reported to be 6 to 7 kcal/mol for the weak to strong RAHB, respectively. The HB component varied in the wide range from 2 to 20 kcal/mol in a series of the β-diketone molecules [115]. In summary, the energy of IHB by the MTA-based method provides reliable values for the RAHB, including both the resonance and HB components.

Recently the IHB energies estimated by the n class="Chemical">MTA-based method were compared with those obtained by other indirect HB descriptors for the wide range of malonaldehyde derivatives by Nowroozi et al. [116]. The latter included structural, spectroscopic, topological, and molecular orbital parameters in the intramolecular RAHBs. The substituent effect of electron-donating and -accepting groups on the HB energy was also estimated by the MTA-based method [116]. Further, the significance of π-electron delocalization (π-ED) of RAHB rings was evaluated by the geometrical factor and the harmonic oscillator model of aromaticity (HOMA). The authors [116] stated that “the excellent linear correlations with MTA energies, which may be implied on the validity of RAHB theory”. Further, it was emphasized by authors that “On the basis of these results, one can claim that the MTA method is reliable for estimation of IMHB energies of RAHB systems.”

A significant application of our MTA-method was in the determination of the push–pull π+/π- (PPππ) effect in the Henry reaction [117], an organocatalytic reaction catalyzed by squaramide [118]. In this reaction, several noncovalent interactions such as HBs, π···H, π···O, etc. were observed between the squaramide and benzaldehyde. It was suggested that the reaction proceeds via an unprecedented mode of activation, modulated by π···H (π···δ+) and π···O (π···δ-) interactions formed with the two rings of a naphthyl group and benzaldehyde. These authors employed the MTA-based approach for determining the energies of these two π interactions [118]. Here, the naphthyl group involved in the PPππ interactions, observed in the intermediates (INTs) and transition state (TS) structures, along the most favorable pathway (P1), was replaced with an H atom. Thus, generated INTs and TS structures along the modeled pathway (P1–H) have the same noncovalent interactions as P1, except the two π interactions whose energy is to be determined. The interaction energy of these two π interactions in the INTs and TS structures were estimated as the energy difference in the total interaction energies between the catalyst and the substrates in P1 and that in modeled P1−H pathways. The estimated sum of two π interactions in INTs and TS were found to be between 2.7 to 4.7 kcal/mol at the ωB97X-D level of theory. In summary, the MTA-based method proposed by us has been successfully employed by several other active research groups for exploring the strengths of IHBs and other intramolecular noncovalent interactions in a variety of systems.

8. Summary and Concluding Remarks

This review article has summarized a direct and simple procedure for the estimation of X–H···Y IHB energy employing a fragmentation method, viz., the molecular tailoring approach (MTA). It has been applied to a variety of systems having multiple HBs over the last one and half decades. A plus point of the method is that it provides the reliable energy of every individual HB and can also estimate the corresponding cooperativity contribution as a result of interconnected networks of HBs. In this present review article, we have discussed the application of the MTA-based method for the estimation of X–H···Y (X–H = O–H, N–H, C–H, etc. and Y = N, O, S, OH, OCH3, O=C, etc.) HB energies in a variety of systems. The systems covered ranged from small alkanediols to large systems such as cyclodextrins and biomolecules such as polypeptides, meta-benzoporphodimethenes, and antioxidant molecules. Further, being of general nature, our method has been utilized as a standard method for a reliable estimation of HB energies by other research groups. It may be emphasized here that the MTA-based method is applicable for the estimation of other noncovalent interactions as well. In recent years, this approach has been applied for the estimation of O–H···O HB energies in water clusters, the π component of the conjugation energy in resonance-assisted, hydrogen-bonded push–pull systems, etc. In a recent impressive application, the method is employed for determining the favorable organocatalytic reaction pathways [118].

It may be noted that the present method can, in principle, be applied to the estimation of IHB energy in larger molecular systems. However, with the increase in the size of a molecule, the size of the fragments would also increase, making the evaluation of HB energies computationally rather demanding. This difficulty can, in principle, be overcome by the use of the MTA methodology discussed in Section 2. For instance, the energy of the parent and fragment molecules can be reliably calculated using MTA at a correlated level of theory which may be further used for estimating the HB energies in larger molecular systems.

Being general in nature, the n class="Chemical">MTA-based method can be employed for exploring other intramolecular interactions such as π···π [119], C-H···π [120], the so-called halogen bonds [121], dihydrogen bonds [122], sulfur bonds [123], metal-H···S and metal-H···Se bonds [124], etc. Although identified in the recent literature by these labels, they are cut from the same cloth called the non-covalent interactions. The prowess of the MTA-based method is that it can be applied to all such inter- and intramolecular interactions.