**1. Introduction**

Hydrogen bond (HB) can be considered to be the most important type of non-covalent interaction [1]. It is important to emphasize that hydrogen bonds are of great importance for the properties of water, the binding of drugs to receptors or the stability of macromolecules [2,3]. Furthermore, HBs are some of the most important factors that affect the packing in crystals [4]. A hydrogen bond is usually defined as X-H. . . :Y, the interaction of a bridging hydrogen attached to an electronegative donor atom (denoted as X) with another electron-rich species (denoted as Y) [5]. A common feature of most types of hydrogen bonds is the elongation of the X-H covalent bond with the co-existing redshift in the X-H stretching vibrations and a decrease in H...Y distance. However, this classical definition does not encompass the whole diversity of hydrogen bonds. As shown, HBs can also be formed with carbon as a hydrogen donor or acceptor [1]. Theoretical studies of C-H proton donors with benzene and ethylene oxide species as acceptors have shown that redshift is not a definitive characteristic of all HBs—in the cited studies, the authors have demonstrated the existence of so-called anti-hydrogen or blue-shifting hydrogen bonds, in which hydrogen bond formation is accompanied by the C-H covalent bond contraction and an increase in its stretch (blue shift) [6,7]. A closely related group of non-covalent interactions, which are based on the same physical principles, is the so-called *σ*-hole bond family [8–13]. The *σ*-hole concept is related to the depletion of electron density on the bridge atom at the extension of its covalent bond, X-A (A denotes the bridge atom and X is any electronegative species).

**Citation:** Wojtkowiak, K.; Jezierska, A. Exploring the Dynamical Nature of Intermolecular Hydrogen Bonds in Benzamide, Quinoline and Benzoic Acid Derivatives. *Molecules* **2022**, *27*, 8847. https://doi.org/10.3390/ molecules27248847

Academic Editor: Miroslaw Jablonski

Received: 14 October 2022 Accepted: 25 November 2022 Published: 13 December 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

Due to a local increase in molecular electrostatic potential (MEP), an atom can participate in a highly directional interaction with Lewis bases [14,15]. These electrostatically driven interactions are named after the family from which the bridge atom is derived [16]. Among the most well known are halogen, chalcogen and pnicogen bonds [17–23]. In general, the strength of these interactions increases with the increasing electronegativity, polarizability and basicity of the donor atom, bridge atom and Lewis base, respectively [24,25]. Redshift is usually a reliable measure of their strength, however, the change in their X-A stretching frequency is less correlated with the interaction strength compared to HBs [18,26]. Moreover, the heavier and more polarizable the bridge atom is, the smaller the co-occurring redshift [25].

Theoretical methods such as the quantum theory of atoms in molecules (QTAIM) [27], Reduced Density Gradient (RDG) [28], Electron Localization Function (ELF) [29] or various perturbational or variational energy-decomposition schemes [30,31] have proven to be invaluable in the detailed characterization of HBs and other secondary bonds [24,32]. From a theoretical perspective, there is also a way to characterize hydrogen bonding and it is usually performed using the QTAIM conceptual apparatus—the indicators of the presence of the HB in the examined system are: (i) the bond path between the proton donor and acceptor atoms with the Bond Critical Point (BCP) located on it; and (ii) relatively small electron density at the BCP between atoms (about an order to even two orders smaller than *ρBCP* for typical covalent bonds) with the Laplacian values close to zero. Considering the energy-decomposition of the hydrogen bonds it is usually described as mainly covalent in nature, especially in the case of strong HBs with an interaction energy above 24 kcal \* mol−<sup>1</sup> [33,34]. In this context, the term "covalent" is associated with a shortening of the H. . . Y distance and an increase in the importance of the induction and dispersion terms (in terms of the nomenclature used in the SAPT energy-partitioning scheme [30]). On the contrary, weak hydrogen bonds are described as interactions in which the electrostatics play the most important role [35].

However, it should not be forgotten that biological or chemical systems are inherently dynamic, hence the above-mentioned, so-called static approaches which, despite being very useful, can only provide information about one particular arrangement of atoms in the complex. Thus, Car–Parrinello molecular dynamics (CPMD) [36] and path integral molecular dynamics (PIMD) [37] are among the most often employed dynamical approaches to study H/D isotope effects [38], non-covalent interactions and hydrogen bridges in particular [39–41]. Furthermore, the quantization of nuclei using PIMD allows one to take into account the quantum nature of the examined system—the hydrogen atom is particularly sensitive to quantum effects, even at a standard temperature due to its small mass, and thus has a relatively large value of thermal de Broglie wavelength compared to other atoms. Noteworthy is the fact that the quantum-tunneling phenomena (which can be efficiently studied using path integration techniques [42,43]) are important in biochemistry and are associated with enzyme-mediated electron or hydrogen transfer [44,45]. Nowadays, attempts are being made to study the quantum effects at the ligand–receptor binding site. An example of this can be found in the work where binding affinities for histamine receptor ligands were studied [46]. Another important aspect is associated with the spectroscopic features' investigation on the basis of the CPMD method. The vibrational properties could be studied using standard approaches, that is, the Fourier transformation of the autocorrelation function of atomic velocity or dipole moment. However, it is also possible to apply a method, which enabled the a posteriori inclusion of quantum effects to the O-H, N-H stretching etc. The method was successfully used in studies where strongly anharmonic systems were investigated [47–50].

Due to the aforementioned reasons, we decided to take a hybrid approach in order to characterize in detail the HBs present in the studied compounds—the dynamical features as well as detailed static characteristics were taken into account. Interactions for dimers extracted from the crystal were quantified and assessed using the SAPT and QTAIM approaches. The metric parameters and spectroscopic signatures of the investigated compounds were obtained and thoroughly analyzed. Obtaining quantum statistics for the nuclear degrees of freedom via PIMD application allowed to estimate the importance of quantum effects in the description of intermolecular HBs, when compared to the classical-quantum CPMD approach.

In the current study, we investigated three aromatic compounds from the benzamide, quinoline and benzoic acid groups. The choice was dictated by the network of hydrogen bonds present in the crystal structures [51]. Benzamide is a derivative of benzoic acid. Some substituted benzamides are well-known commercial drugs, e.g., procainamide, imatinib and veralipride [52–54]. Benzamides are still an attractive group of compounds, especially in drug design, where very often their derivatives are taken into consideration as compounds that show a specific type of biological activity. Therefore, they are studied both experimentally and by molecular modeling methods, e.g., [55–59]. We chose 2,6-difluorobenzamide, which is a metabolite of pesticide diflubenzuron [60], to theoretically study its hydrogen bonding network. It was found that benzamide as well as 2,6-difluorobenzamide can form mutual intermolecular hydrogen bonds. However, the extended amide. . . amide dimer synthon in benzamide can form a network of HBs via NH2 group. Concerning the 2,6-difluorobenzamide, the presence of fluorine atoms allows the formation of other intermolecular interactions [51]. The next compound chosen for our theoretical investigations is 5-hydroxyquinoline, where the O-H. . . N intermolecular hydrogen bond is present as the strongest intermolecular interaction, however, C-H. . . O interaction was noted in the crystalline phase [35]. Taking into account the fact that quinoline and its derivatives have diverse applications, e.g., in medicine as drugs, compounds exhibiting various biological activity, as dyes, and as solvents [61–67], it is of interest to investigate the properties of hydrogen bonds in the class of compounds. The last compound taken into account was 4-hydroxybenzoic acid, which is the simplest aromatic carboxylic acid. Benzoic acid occurs naturally in many plants [68]. Its salts are used as food and cosmetics preservatives [69]. Generally speaking, it is an important precursor for the industrial synthesis of many other organic substances [70,71]. 4-hydroxybenzoic acid is primarily known as the basis for the preparation of its esters (parabens), which are used as preservatives in, e.g., cosmetics [72,73]. In our case, it was interesting to explore the dynamical nature of intermolecular hydrogen bonds, wherein the carboxylic as well as hydroxy groups were involved [51].

We hope that our research contributes to the knowledge of hydrogen bond dynamics, and will help with the rational design of new derivatives with specific properties. Therefore, the main objective of this research is to perform multi-factor studies of non-covalent interactions in the examined compounds. In order to reproduce the dynamical nature of a hydrogen bonds network, we employed the Car–Parrinello molecular dynamics [36]. To be able to make comparisons the time evolution simulations were performed in vacuo and in the crystalline phase. The nuclear quantum effects (NQEs) was taken into account and path integral molecular dynamics (PIMD) simulations were performed for this purpose [37,74]. The quantum theory of atoms in molecules (QTAIM) [27] was applied for electron density topological studies enabling the estimation of the interaction strength and the detection of weaker interactions. An application of symmetry-adapted perturbation theory (SAPT) [30] method allowed the energy decomposition in the studied dimers.
