Assessing Parameter Suitability for the Strength Evaluation of Intramolecular Resonance Assisted Hydrogen Bonding in o-Carbonyl Hydroquinones

Intramolecular hydrogen bond (IMHB) interactions have attracted considerable attention due to their central role in molecular structure, chemical reactivity, and interactions of biologically active molecules. Precise correlations of the strength of IMHB’s with experimental parameters are a key goal in order to model compounds for drug discovery. In this work, we carry out an experimental (NMR) and theoretical (DFT) study of the IMHB in a series of structurally similar o-carbonyl hydroquinones. Geometrical parameters, as well as Natural Bond Orbital (NBO) and Quantum Theory of Atoms in Molecules (QTAIM) parameters for IMHB were compared with experimental NMR data. Three DFT functionals were employed to calculated theoretical parameters: B3LYP, M06-2X, and ωB97XD. O…H distance is the most suitable geometrical parameter to distinguish among similar IMHBs. Second order stabilization energies ΔEij(2) from NBO analysis and hydrogen bond energy (EHB) obtained from QTAIM analysis also properly distinguishes the order in strength of the studied IMHB. ΔEij(2) from NBO give values for the IMHB below 30 kcal/mol, while EHB from QTAIM analysis give values above 30 kcal/mol. In all cases, the calculated parameters using ωB97XD give the best correlations with experimental 1H-NMR chemical shifts for the IMHB, with R2 values around 0.89. Although the results show that these parameters correctly reflect the strength of the IMHB, when the weakest one is removed from the analysis, arguing experimental considerations, correlations improve significantly to values around 0.95 for R2.


Introduction
Hydrogen bonding (HB) is generally described as a stabilizing electrostatic interaction with a partly covalent character. Some of the most direct experimental evidence of the HB formation is the deshielding of the hydrogen atom involved in the interaction, which is easily observed by NMR spectroscopy [1,2]. It is well known that intramolecular hydrogen bonding (IMBH) has a key role in determining the structure and properties of biologically-active molecules, and their study in model
Molecules 2019, 24, x 3 of 13 rearrangement of these benzofurans, yielding the corresponding bicyclic hydroquinones (18,20,(22)(23)(24). Monomethyl ethers of hydroquinone 19 and 21 were obtained by methylation of hydroquinones 18 and 20, respectively; they were synthesized to assess the effect of blocking the hydroxyl group not involved in IMHB. Three hybrid functional were employed in this work: B3LYP, composed of Becke three parameters hybrid functional B3 and the non-local correlation functional of Lee, Yang, and Parr (LYP) [34][35][36]. The second corresponds to a hybrid meta-generalized gradient-approximation (hybrid meta-GGAs) functional, M06-2x, which have been described as one the of the best Truhlar group functionals for non-covalent interactions [37]. The last functional corresponds to the long-range corrected functional ωB97XD, which is designed to avoid rapid die-off of the non-coulomb part of the exchange functionals, which affects the modelling of long distance processes [38].
In this first part, we study main geometrical parameters from the IMHBs and their correlation with 1 H-NMR chemical shifts. The results are presented in Table 1. Calculations were carried out at DFT level using the B3LYP, M06-2X, and ωB97XD functionals with the 6-311++g(d,p) basis set. Among the geometrical parameters, only the O15 … H16 distances displayed enough differentiation to be used as comparison parameters for correlation with the chemical shifts of chelated protons. Compound 21 displayed the largest chemical shift and a more de-shielded proton, indicating the strongest IMHB, followed in strength by 20. The later indicates that methyl substituents on C1 and C2 increase the IMHB strength. Additionally, the methylation at O12 and also a methyl group at C9 increase the strength of the IHB in 21 and 22, respectively. The three functionals give the shortest O15 … H16 distance for 21, indicating the strongest IHB, which is in agreement with the largest downfield chemical shift. The smallest chemical shift for 18 indicates the weakest IHB in the series. Nevertheless none of the three functional gives the largest O15 … H16 distance for 18. These discrepancies can reflect the presence of intermolecular interactions competing with the IMHB, contributing in a non-negligible degree to the 1 H-NMR spectra at the probe temperature, especially Three hybrid functional were employed in this work: B3LYP, composed of Becke three parameters hybrid functional B3 and the non-local correlation functional of Lee, Yang, and Parr (LYP) [34][35][36]. The second corresponds to a hybrid meta-generalized gradient-approximation (hybrid meta-GGAs) functional, M06-2x, which have been described as one the of the best Truhlar group functionals for non-covalent interactions [37]. The last functional corresponds to the long-range corrected functional ωB97XD, which is designed to avoid rapid die-off of the non-coulomb part of the exchange functionals, which affects the modelling of long distance processes [38].
In this first part, we study main geometrical parameters from the IMHBs and their correlation with 1 H-NMR chemical shifts. The results are presented in Table 1. Calculations were carried out at DFT level using the B3LYP, M06-2X, and ωB97XD functionals with the 6-311++g(d,p) basis set. Among the geometrical parameters, only the O15 . . . H16 distances displayed enough differentiation to be used as comparison parameters for correlation with the chemical shifts of chelated protons. Compound 21 displayed the largest chemical shift and a more de-shielded proton, indicating the strongest IMHB, followed in strength by 20. The later indicates that methyl substituents on C1 and C2 increase the IMHB strength. Additionally, the methylation at O12 and also a methyl group at C9 increase the strength of the IHB in 21 and 22, respectively. The three functionals give the shortest O15 . . . H16 distance for 21, indicating the strongest IHB, which is in agreement with the largest downfield chemical shift. The smallest chemical shift for 18 indicates the weakest IHB in the series. Nevertheless none of the three functional gives the largest O15 . . . H16 distance for 18. These discrepancies can reflect the presence of intermolecular interactions competing with the IMHB, contributing in a non-negligible degree to the 1 H-NMR spectra at the probe temperature, especially between pairs of molecules where one phenolic hydroxyl is replaced by a methoxy group. The same considerations can be applied to correlations with the next parameters. Global correlation among the chemical shift and O15 . . . H16 distance for all molecules were carried out ( Figure 1). Results indicate that the ωB97XD functional gives the best correlation with an R 2 = 0.89, followed by M06-2X with a R 2 = 0.85 and B3LYP, and the worst correlation with R 2 = 0.80. between pairs of molecules where one phenolic hydroxyl is replaced by a methoxy group. The same considerations can be applied to correlations with the next parameters. Global correlation among the chemical shift and O15 … H16 distance for all molecules were carried out ( Figure 1). Results indicate that the ωB97XD functional gives the best correlation with an R 2 = 0.89, followed by M06-2X with a R 2 = 0.85 and B3LYP, and the worst correlation with R 2 = 0.80. In order to estimate the energy of IMHB, the open-close method was employed, which defines the H-bond energy as the difference between the open and close conformers [39,40]. Calculations were carried out with the functional ωB97XD, which presents the best correlation between experimental chemical shift and O16 … H15 distance. Table 2 shows that energy for IMHB in the studied hydroquinones are around 15 kcal/mol. Hydroquinone 20 presents the highest value, 16.37 kcal/mol, while 19 presents the lowest IMHB energy, 14.22 kcal/mol. These energy differences are large enough to discard the contribution of the open form to the experimental NMR shieldings. Figure 2 shows the correlation between Hydrogen bond energy (EHB) and 1 H-NMR chemical shifts (δ), which presents a R 2 = 0.71, a value lower than the correlation with the O15 … H16 distance (R 2 = 0.89), calculated with the same functional. In order to estimate the energy of IMHB, the open-close method was employed, which defines the H-bond energy as the difference between the open and close conformers [39,40]. Calculations were carried out with the functional ωB97XD, which presents the best correlation between experimental chemical shift and O16 . . . H15 distance. Table 2 shows that energy for IMHB in the studied hydroquinones are around 15 kcal/mol. Hydroquinone 20 presents the highest value, 16.37 kcal/mol, while 19 presents the lowest IMHB energy, 14.22 kcal/mol. These energy differences are large enough to discard the contribution of the open form to the experimental NMR shieldings. Figure 2 shows the correlation between Hydrogen bond energy (E HB ) and 1 H-NMR chemical shifts (δ), which presents a R 2 = 0.71, a value lower than the correlation with the O 15 . . . H 16 distance (R 2 = 0.89), calculated with the same functional.   The results from NBO analysis are presented in Table 3. Second order stabilization energies ΔEij (2) between donor (φi) and acceptor orbital (φj) allow one to study orbital interactions responsible for the IMHB. For all studied compounds, the HB donor atom is always an oxygen (two lone-pairs) and the acceptor is the antibonding sigma orbital of an O-H groups (σ * OH). As we found in the geometrical analysis, strong IMHB presence in 21 (the largest chemical shift) is in agreement with the largest ΔEij (2) obtained from the calculations with the three functionals. Nevertheless, ΔEij (2) did not predict well the weakest IMHB, corresponding to 18, and the lowest stabilization energy was not predicted by any of the three calculations using the studied functionals. Global correlation among the chemical shift and ΔEij (2) for all HQs were carried out ( Figure 3). As in the geometrical analysis, correlation results indicate that the best correlation is obtained with the ωB97XD functional, with R 2 = 0.89. However, in this case it was followed by B3LYP with R 2 = 0.87 and then M06-2X, the one that gave the worst correlation with R 2 = 0.85. Table 3. Second order stabilization energies (ΔEij (2) ) between donor LP O15 and acceptor σ* O11-H16 orbitals involved in the IMHB. Energies in kcal/mol.

Molecule
ΔEij (  The results from NBO analysis are presented in Table 3. Second order stabilization energies ∆E ij (2) between donor (ϕ i ) and acceptor orbital (ϕ j ) allow one to study orbital interactions responsible for the IMHB. For all studied compounds, the HB donor atom is always an oxygen (two lone-pairs) and the acceptor is the antibonding sigma orbital of an O-H groups (σ * OH ). As we found in the geometrical analysis, strong IMHB presence in 21 (the largest chemical shift) is in agreement with the largest ∆E ij (2) obtained from the calculations with the three functionals. Nevertheless, ∆E ij (2) did not predict well the weakest IMHB, corresponding to 18, and the lowest stabilization energy was not predicted by any of the three calculations using the studied functionals. Global correlation among the chemical shift and ∆E ij (2) for all HQs were carried out ( Figure 3). As in the geometrical analysis, correlation results indicate that the best correlation is obtained with the ωB97XD functional, with R 2 = 0.89. However, in this case it was followed by B3LYP with R 2 = 0.87 and then M06-2X, the one that gave the worst correlation with R 2 = 0.85.  The QTAIM analysis is presented in Table 4. It shows negative values for  2 ρ and H for all studied compounds, indicating the covalent nature of their IHBs. |V|/G ratio is also a sensitive parameter for evaluating the covalency [41,42]. V can be interpreted as the pressure exerted on the electrons by the HB system, while G can be interpreted as the pressure exerted (as a reaction) by the electrons on the HB system, at the critical point. For negative values of  2 ρ, as occurred in all our cases, values of |V|/G > 2 were indicative of covalent interaction [2]. Hydrogen bond energy (EHB) was obtained from the potential energy density (V) [13]. All values were slightly above 30 kcal/mol, while NBO analysis gave all energies slightly below 30 kcal/mol for the same interaction. These values were approximately twice those calculated by the open-close method, which indicates that NBO and AIM overestimate the interaction energy for the IMHB. EHB values calculated by different methods have shown to be generally inconsistent [43]. As we found in both geometrical and NBO analysis, the strongest IHB is obtained for compound 21 and is well represented by QTAIM analysis for the three functionals, giving the highest EHB for 21. On the other hand, similarly to NBO analysis, here QTAIM EHB also did not adequately reflect the weakest IMHB present in 18. Global correlations between QTAIM EHB and chemical shifts are presented in Figure 4. A better performance is shown by ωB97XD and B3LYP with a R 2 = 0.88. As in the above cases, M06-2X gave the worst correlation with R 2 = 0.85. Table 4. Atoms-in-molecule parameters for IHB of 18-24. Electron density at the critical point ρ BCP (a.u), its Laplacian  2 ρ(a.u.), electron kinetic energy density G (a.u.), potential energy density V (a.u.), total electron energy density H (a.u.), and hydrogen bond energy EHB (kcal/mol).

Functional
Molecule  The QTAIM analysis is presented in Table 4. It shows negative values for ∇ 2 ρ and H for all studied compounds, indicating the covalent nature of their IHBs. |V|/G ratio is also a sensitive parameter for evaluating the covalency [41,42]. V can be interpreted as the pressure exerted on the electrons by the HB system, while G can be interpreted as the pressure exerted (as a reaction) by the electrons on the HB system, at the critical point. For negative values of ∇ 2 ρ, as occurred in all our cases, values of |V|/G > 2 were indicative of covalent interaction [2]. Hydrogen bond energy (E HB ) was obtained from the potential energy density (V) [13]. All values were slightly above 30 kcal/mol, while NBO analysis gave all energies slightly below 30 kcal/mol for the same interaction. These values were approximately twice those calculated by the open-close method, which indicates that NBO and AIM overestimate the interaction energy for the IMHB. E HB values calculated by different methods have shown to be generally inconsistent [43]. As we found in both geometrical and NBO analysis, the strongest IHB is obtained for compound 21 and is well represented by QTAIM analysis for the three functionals, giving the highest E HB for 21. On the other hand, similarly to NBO analysis, here QTAIM E HB also did not adequately reflect the weakest IMHB present in 18. Global correlations between QTAIM E HB and chemical shifts are presented in Figure 4. A better performance is shown by ωB97XD and B3LYP with a R 2 = 0.88. As in the above cases, M06-2X gave the worst correlation with R 2 = 0.85.
The results shown in Figure 3 corroborates that correlations with the weakest IMHB of this series are not well represented by the computational calculations employed. In order to test what is the impact of the weakest IMHB on the correlations studied previously, compound 18 was eliminated from the data and correlations were obtained again. The results showed a significant improvement, without change in the relative position among different functionals and parameters studied. Table 4. Atoms-in-molecule parameters for IHB of 18-24. Electron density at the critical point ρ BCP (a.u), its Laplacian ∇ 2 ρ (a.u.), electron kinetic energy density G (a.u.), potential energy density V (a.u.), total electron energy density H (a.u.), and hydrogen bond energy E HB (kcal/mol).  The results shown in Figure 3 corroborates that correlations with the weakest IMHB of this series are not well represented by the computational calculations employed. In order to test what is the impact of the weakest IMHB on the correlations studied previously, compound 18 was eliminated from the data and correlations were obtained again. The results showed a significant improvement, without change in the relative position among different functionals and parameters studied. Figure 5 shows correlations between 1 H-NMR chemical shifts and O11-H16 … O15 distance for the IHBs. R 2 for B3LYP changed from 0.80 to 0.88, from 0.85 to 0.91 for M06-2x, and from 0.89 to 0.95 for ωB97XD.  Similarly, correlations between 1 H-NMR chemical shifts and NBO energy for the IHBs are presented in Figure 6. R 2 changed from 0.88 to 0.95 for B3LYP, from 0.85 to 0.92 for M06-2X, and from 0.89 to 0.96 for ωB97XD.  Similarly, correlations between 1 H-NMR chemical shifts and NBO energy for the IHBs are presented in Figure 6. R 2 changed from 0.88 to 0.95 for B3LYP, from 0.85 to 0.92 for M06-2X, and from 0.89 to 0.96 for ωB97XD. Similarly, correlations between 1 H-NMR chemical shifts and NBO energy for the IHBs are presented in Figure 6. R 2 changed from 0.88 to 0.95 for B3LYP, from 0.85 to 0.92 for M06-2X, and from 0.89 to 0.96 for ωB97XD.  Finally, correlations between 1 H-NMR chemical shifts and QTAIM energy for the IHBs, presented in Figure 7, showed a change in R 2 from 0.88 to 0.94 for B3LYP, from 0.85 to 0.91 to M06-2X, and from 0.88 to 0.94 for ωB97XD.

Functional
Finally, correlations between 1 H-NMR chemical shifts and QTAIM energy for the IHBs, presented in Figure 7, showed a change in R 2 from 0.88 to 0.94 for B3LYP, from 0.85 to 0.91 to M06-2X, and from 0.88 to 0.94 for ωB97XD.

Synthetic Procedures
Melting points are uncorrected. All NMR spectra were acquired using a Bruker AVANCE DRX 300 spectrometer (Bruker BioSpin GmbH, Rheinstetten, Germany) operating at 300.13 MHz ( 1 H) or 75.47 MHz ( 13 C). Measurements were carried out at a probe temperature of 300 K.

Synthetic Procedures
Melting points are uncorrected. All NMR spectra were acquired using a Bruker AVANCE DRX 300 spectrometer (Bruker BioSpin GmbH, Rheinstetten, Germany) operating at 300.13 MHz ( 1 H) or 75.47 MHz ( 13 C). Measurements were carried out at a probe temperature of 300 K.

Computational Calculations
The calculations were carried out using the Gaussian 09 [44] program package (Revision a.01; Gaussian, Inc.: Wallingford, CT, USA). No symmetry constraints were imposed on the optimizations, which were performed at the DFT B3LYP/6-311++G (d,p) level. No imaginary vibrational frequencies were found at the optimized geometries, indicating that they are true minima of the potential energy surface. NBO analysis was performed using the NBOPro 6.0 [45] program package (NBO 6.0; University of Wisconsin: Madison, WI, USA). QTAIM analysis was performed using the AIM2000 [46] program package (Büro für innovative software, Bielefeld, Germany).

Conclusions
In this work, we found that O . . . H distance is the most suitable geometrical parameter for distinguishing among IMHBs of similar strengths. Second order stabilization energies ∆E ij (2) from NBO analysis and hydrogen bond energy (E HB ) obtained from QTAIM analysis were also able to properly distinguish among IMHBs in the studied molecule series. In all cases, ∆E ij (2) from NBO give values for the IMHB lower than E HB obtained from QTAIM analysis. ∆E ij (2) are slightly below 30 kcal/mol, while E HB from QTAIM gave values somewhat above 30 kcal/mol, which are approximately twice of those calculated by the open-close method. In all cases, parameters calculated using ωB97XD gave the best correlations with experimental 1 H-NMR chemical shifts for the IMHBs, with values of R 2 around 0.89. Despite this, the results showed that the employed parameters correctly described the strength of the IMHB, when the weakest one was removed from the analysis, correlations improved significantly to values around 0.95 for R 2 . These results can help to select the most suitable parameters for modeling molecules with IMHBs of similar strengths and achieving an adequate distinction between them through computational calculations.