**1. Introduction**

The nature of hydrogen bonds is complex and still presents open questions. In addition to conventional hydrogen bonds, during recent decades so called unconventional hydrogen bonds have appeared as important scientific topics [1–8]. Hydrogen bonds are located on

 Jezierska, A.; Panek, J.J. Naphthazarin Derivatives in the Light of Intra- and Intermolecular Forces. *Molecules* **2021**,*26*, 5642. https://doi.org/10.3390/ molecules26185642

K.;

Poche´c,

M.;

Academic Editors: Mirosław Jabło ´nski and Benedito José Costa Cabral

**Citation:** Kułacz,

Received: 12 July 2021 Accepted: 14 September 2021 Published: 17 September 2021

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

**Copyright:** © 2021 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/).

the interaction strength ladder in the middle position, between covalent, ionic, and van der Waals interactions [9]. They are much weaker than chemical covalent bonds, but their presence is of grea<sup>t</sup> significance in nature. The strength of most hydrogen bonds is between 10 kJ/mol and 40 kJ/mol [10]; however, they are ubiquitous and cannot be neglected in the discussion of factors decisive for the structure of bulk materials: liquids and solid states. They are key elements of many processes at the molecular level, as well as influencing molecular and macroscopic properties of various systems [11–15]. They were found to be, e.g., present in enzymatic reactions [16,17], responsible for structure stabilization [18–22], arrangemen<sup>t</sup> of molecules in crystals [23,24] and supporting molecular engineering [25,26]. Therefore, it is evident that hydrogen bonds are non-covalent interactions relevant in various branches of contemporary science [27].

Hydrogen bonds are generally divided into intra- and intermolecular ones. While the latter provide the possibility of supramolecular assembly and molecular recognition, the former are relevant to the molecular structure, stabilizing the hydrogen-bonded conformations and giving rise to the tautomeric forms. The intramolecular hydrogen bonds are more open to types of enhancement such as resonance-assisted phenomena [28,29], but on the other hand the charge-assisted bonds can be easily formed between two separate molecules [30–32]. Another aspect related to hydrogen bonds is the proton transfer phenomenon and associated changes in geometric and electronic structure parameters [33–35]. This phenomenon has been studied using experimental and theoretical approaches, because its role cannot be neglected in biomolecules and other compounds where tautomeric forms occur, e.g., see refs. [36–42].

The current study is a continuation of our long-term research effort to shed light on the hydrogen bridge dynamics in molecular crystals with diverse chemical compositions— cyclic compounds in particular. In order to make our study comprehensive, we focus not only on intramolecular interactions, but also on intermolecular forces. Our studies covered, e.g., description of molecular properties in monocyclic aromatic o-hydroxy Schiff and Mannich bases [43,44], bicyclic N-oxides [45–47] and proton sponges [41,48]. Following the line of our research devoted to intramolecular hydrogen bonds investigations in compounds possessing fused rings, we have focused our attention on naphthazarin and its derivatives. The 5,6-dihydroxy-1,4-naphtoquinone (commonly referred to as naphthazarin) and its derivatives are a part of the wider naphthoquinone group. Members of this class of compounds are widely distributed in natural sources [49] and have been proven to possess many interesting biological properties. Recent years have brought reports on their antibacterial, antifungal and biostatic activities [50,51]. Naphthazarin, alone or in conjunction with other compounds, can be used as a potent biopesticide or insecticide [52]. Medicine also has many potential ways of utilizing naphthazarin and its derivatives. They have been proven to possess high anti-inflammatory potential and the positive effect on the healing of wounds [53]. More elaborate derivatives can also be applied in oncological treatment—they have been reported to have an inhibitory effect on DNA Topoisomerase-I [54], heat-shock factor and glutathione status in the aftermath of hypoxia [55]. All of the aforementioned properties show that there is still need to study both the possible ways of synthesis and design of derivatives with desired properties [56,57]. One of the most pronounced characteristics of naphthazarin is the presence of two hydrogen bonds between the hydroxyl groups and their neighboring carbonyl oxygen atoms. This not only forms two distinct quasi-rings in the structure, but also allows the opportunity to study the effects of substitution in the fused rings and double hydrogen bonding properties [58]. In case of the naphthazarin and its selected derivatives, investigations into physico-chemical properties are reported by [59–61]. The computational studies have been also used to assess compatibility of experimental and theoretical data in the IR and Raman spectra of naphthazarin, which in turn allowed precise assignment of bands [62].

Here, we present our theoretical results obtained for two naphthazarin derivatives: 2,3-Dimethyl-5,8-dihydroxy-1,4-naphthoquinone (**1**) and 5,8-Dihydroxy-2,3-dimethoxy-6-methylnaphtho-1,4-quinone (**2**), presented in Figures 1 and S1. The motivation for the

current work and the choice of the compounds was the comparison of symmetric and asymmetric substitution with diverse, but not very strong, substituent properties represented by the classical physico-chemical parameters, Hammett constants. The –Me and –OMe substituents possess different properties regarding their electro- and nucleophilicity. The methyl groups are relatively mild on the Hammett scale (their classical Hammett constants are only −0.07 and −0.17 for *σm* and *<sup>σ</sup>p*, respectively [63]. The –OMe substituents are reported to reach the values of +0.12 and −0.27 for *σm* and *<sup>σ</sup>p* [63]. Such values indicate the possibility of an interesting interplay between local properties of the substituents and an environmental influence in the solid state (with the presence of other molecules). The physico-chemical properties of both compounds were studied on the basis of X-ray methods as well as, e.g., NMR spectroscopy for **1** by Rodríguez et al. [64] and for **2** by Cannon et al. [65]. Concerning compound **1** the crystal structure is built of molecules, which are stacked up the c axis and the molecules overlap forming the charge-transfer complex [64]. Compound **2** contains two methoxyl groups, which are slightly different with regard to geometry: one of the groups lies in the plane of the ring, but the methyl group of the methoxyl deviates from the ring plane by 1.08 Å [65]. This experimentally observed difference in geometry, which may result from packing forces, would not necessarily be detected in solution. According to the authors, electronic effects were responsible for observed non-equivalence of the methoxyl resonances in the NMR spectrum in solution. The electron-donating properties of the methyl group are well known and characterized [63]. Inductive effects usually cover short distances; however, in the case of compound **2**, they could influence the methoxy groups as well [65]. Therefore, the main aim of the study was further examination of internal and external forces responsible for the molecular features of the compounds.

**Figure 1.** Molecular structures of the studied 2,3-dimethylnaphthazarin (**1**) and 2,3-dimethoxy-6- methylnaphthazarin (**2**) with the atom numbering scheme applied in the current study.

In order to achieve the goals delineated in the previous paragraphs, the fundamental issues of comparisons of intra- and intermolecular phenomena, influence of substituents, and correlation of bridge proton motions, diverse theoretical approaches were considered and employed. Static and dynamical models were developed on the basis of Density Functional Theory (DFT) [66,67] and Car–Parrinello molecular dynamics [68]. The simulations were carried out in the gas phase and in the solid state. Our particular attention was focused on: (i) The proton reaction path and related energy changes in the monomeric forms of the studied compounds. We put emphasis on the substituent effects on the hydrogen bond properties; (ii) The electronic structure and topology changes—the comparison of molecular and proton transferred (PT) forms on the basis of the Atoms in Molecules (AIM) theory [69]; (iii) The energy partitioning in the dimers extracted from the X-ray data as well as obtained theoretically on the basis of Symmetry-Adapted Perturbation Theory (SAPT) [70]; (iv) The hydrogen bridges dynamics in the gas and crystalline phases, which enabled us to detect differences derived from environmental effects—gas phase vs. solid state comparison

and the structural impact of the nuclear quantum effects (NQE) for the bridge protons; (v) Vibrational signatures present in the studied compounds, but particularly associated with the intramolecular hydrogen bond—gas phase vs. solid state comparison. To the best of our knowledge this is the first study that examines these particular naphthazarin derivatives considering intra- and intermolecular forces.

## **2. Results and Discussion**

*2.1. Geometric and Electronic Structure Description of Naphthazarin Derivatives Monomers with Special Emphasis on Intramolecular Hydrogen Bonds*

This section contains results from diverse theoretical approaches, ranging from gradientcorrected DFT to post-Hartree–Fock schemes. However, the common trait is that the computational models represent the gas phase molecules using atom-centered Gaussian basis sets providing wavefunctions, Kohn–Sham orbitals and electron density. The chosen DFT functionals belong to the most widely used approaches: the PBE functional is of the Generalized Gradient Approximation (GGA) type and does not use the exact (Hartree– Fock) exchange. B3LYP is a hybrid functional with Hartree–Fock admixture. The *ω*B97XD is also a hybrid functional, additionally including empirical dispersion correction. The chosen post-Hartree–Fock schemes are Møller–Plesset second order perturbative calculus, MP2, and Coupled Cluster theory with single and double excitations, CCSD. The Car– Parrinello scheme, relevant to the further sections, is based on the delocalized plane-wave basis sets with inherent periodicity. This makes CPMD calculations technically easier for solids and liquids than for the gas phase, where periodicity has to be artificially removed. Moreover, exact (Hartree–Fock) exchange is not easily implementable on the plane-wave basis. Therefore, CPMD calculations utilize almost solely the GGA functionals, such as PBE, due to their efficiency and well tested performance.

The molecular structures of the monomeric forms of studied 2,3-dimethylnaphthazarin (**1**) and 2,3-dimethoxy-6-methylnaphthazarin (**2**) are presented in Figure 1. The intramolecular non-covalent interactions present in the molecules are classified as Resonance-Assisted Hydrogen Bonds (RAHBs) [28]. The geometry of the studied compounds was modified to analyze the energy of various isomers (see Table S1 for details; this table includes not only the electronic energy, but also the energy values corrected for the vibrational zero point energy (ZPE) contribution). The lowest energy was found for the molecular forms. The double-PT forms are slightly higher in energy (we will use the convention: electronic/ZPE-corrected value in kcal/mol): for **1**, the difference is 1.94/1.91 kcal/mol (B3LYP), 1.61/1.38 kcal/mol (PBE) and 2.02/1.54 kcal/mol ( *ω*B97XD). For **2**, the corresponding increases in energy of the double-PT form with respect to the molecular form are: 2.94/2.74 kcal/mol (B3LYP), 2.78/2.29 kcal/mol (PBE) and 2.91/2.58 kcal/mol ( *ω*B97XD). This shows that from a thermodynamic point of view the molecular forms are preferred, but not strongly, over the PT forms. The trans forms through which the double-proton transfer proceeds [61] are even higher in energy, e.g., for **1** the differences with respect to the molecular forms are: 5.30/4.65 kcal/mol (B3LYP), 2.82/1.26 kcal/mol (PBE) and 6.81/6.23 kcal/mol ( *ω*B97XD). However, these values are all well within the range thermally accessible to the extended molecules at room temperature, especially when the effects of donor–acceptor distance modulation are accounted for (see the Car–Parrinello study below).

The geometric details of the intramolecular hydrogen bonds (comparison of experimental and computational data) are summarized in Table S2. It is shown that the DFT method was able to reproduce the metric parameters related to the hydrogen bonding with a good agreement. The impact of the substituents, as already noted in the Introduction, corresponds to their inductive and resonance properties. The Hammett constants are, respectively: *σm*<sup>=</sup> −0.07 and *<sup>σ</sup>p*<sup>=</sup> −0.17 for the -Me, while for –OMe they are *σm*<sup>=</sup> +0.12 and *<sup>σ</sup>p*<sup>=</sup> −0.27 [63]. The bridge protons are located at the non-substituted ring of **1** and methyl-substituted ring of **2**, which shows that substitution affects aromaticity and the modified rings prefer participation in the quinoid-like structure.

The proton potential functions for the proton motion are presented in Figure 2. The proton reaction energy paths were investigated using B3LYP, *ω*B97XD and PBE functionals with the 6-311++G(2d,2p) basis set. As shown in the figure, two energy minima were obtained in the case of both compounds and both studied bridges. In the case of compound **1**, only one hydrogen bridge was analyzed due to the symmetry exhibited by the compound. Concerning compound **2**, both intramolecular hydrogen bridges were analyzed. The compound is not symmetrical due to the presence of the methyl group. Before discussing the details of the proton potential functions, it is necessary to consider whether singleor double-proton transfer should be pursued. Our DFT and CPMD results for naphthazarin [61] indicate that the simultaneous double-proton transfer is less probable and leads to higher barriers, which is assumed to be the effect of deeper modification of aromaticity than in the case of the single-proton transfer event. However, the single-proton transfer enables the second proton transfer (PT) event to happen very fast but not simultaneously, in the order of several O-H stretching periods.

The deeper energy minimum is localized at ca. 1 Å of the O8-H*BP*1/O5-H*BP*<sup>2</sup> covalent bond length in both compounds. The elongation of the bond towards the acceptor atom (O1/O4) provided information of the energy barrier, which was found—depending on the applied functional—to be very similar in both studied compounds. The highest energy barrier was obtained for the *ω*B97XD functional (10.1 ± 0.05) kcal/mol for Bridge 1 of both compounds and (9.55 ± 0.05) kcal/mol for Bridge 2 of compound **2**. The results obtained with assistance of the B3LYP functional are (7.95 ± 0.05) kcal/mol for Bridge 1 and (7.45 ± 0.05) kcal/mol for Bridge 2 in the case of compound **2**. The lowest energy barrier was noticed as a result of PBE functional application (3.55 ± 0.05) kcal/mol for all cases. The DFT results were validated with the single-point energy calculations at the post-Hartree–Fock MP2 and CCSD level for the PBE geometries. Both of these approaches yielded the same ordering of relative energies than the DFT functionals: the barriers for Bridge 1 are almost equal for **1** and **2**, and the barrier for Bridge 2 of compound **2** is slightly lower. The MP2 perturbative calculus provided a barrier height of 6.72 kcal/mol for **1**, 6.86 kcal/mol for Bridge 1 of **2**, and 6.59 kcal/mol for Bridge 2 of **2**. The corresponding CCSD values of barrier height estimate are: 9.59 kcal/mol for **1**, 9.37 kcal/mol for Bridge 1 of **2**, and 9.13 kcal/mol for Bridge 2 of **2**. Our previous results on naphthazarin [61] have shown that MP2 and CCSD methods provide PT barrier heights correspondingly lower and higher from the accurate CCSD(T) barrier height, and we expect similar performances from these methods in the current study. This indicates that the barrier height estimates from the post-Hartree–Fock methods and DFT functionals are in agreement.

**Figure 2.** The potential energy profiles for the proton motion in the hydrogen bridges of compounds **1** (**a**) and **2** (**b**), respectively. The hydrogen bridges denoted as O8-H*BP*1...O1 (Bridge 1) for compound **1** and O8-H*BP*1...O1 (Bridge 1) and O5-H*BP*2...O4 (Bridge 2) for compound **2** are presented. In the case of **2**, the solid line denotes Bridge 1 while the dotted line denotes Bridge 2.

Returning to the discussion of structure–energy relations, we note that the second energy minimum is shallow; therefore, we could expect that the bridged proton is mostly localized on the donor (O8/O5) atom. However, the presence of the second energy minimum indicates that the bridge protons are labile and they could approach the proton-acceptor atom domain. Almost identical energy barriers showed that the substituent effects as well as the lack of the symmetry (in the case of compound **2**) have not significantly influenced the proton transfer reaction path in the investigated napththazarin derivatives.

The electronic structure analysis was carried out based on AIM theory. The selected results of the analysis are presented in Tables 1 and 2. The partial atomic charges are reported for atoms forming quasi-rings in the studied compounds (see Table 1). We have analyzed molecular and tautomeric (proto transferred (PT)) forms of both compounds. The electron density of the donor atom (O8/O5) is smaller when the bridged proton is attached to it. A decrease in the electron density at the acceptor atom (O1/O4) is observed for the tautomeric (PT) form. It can be seen that the hydrogen atoms (H*BP*<sup>1</sup> and H*BP*2) are more positively charged when they are transferred to the acceptor atom side. Next, the sum of partial atomic charges in the quasi-rings was computed. It was found that for compound **1**, Bridge 1, the sum decreased from −0.1532 [e] in the molecular form to −0.1591 [e] in the proton transferred form. A similar observation was made for compound **2**—the sum of the quasi-ring (Bridge 1) atomic charges decreased from −0.1042 [e] to −0.1089 [e]. Concerning the hydrogen bridge denoted as Bridge 2 (see Table 1), in compound **1**, there was a decrease in the sum of the partial atomic charges in the quasi-ring from −0.1554 [e] to −0.1583. However, an opposite situation was found in the case of compound **2**— there was an increase in the values of the sum of the atomic charges from −0.1707 to −0.1663 [e]. This could be associated with the presence of the methyl group in the vicinity of the quasi-ring and asymmetry introduced by it to compound **2**. It is also known from the crystal structure of the compound [65] that the methoxy groups are sterically not equivalent; moreover, these groups are not chemically equivalent due to their having different relative positions with regard to the methyl group. There was also an interaction between the methoxy group and O4 proton-acceptor atom (for details, see the text below). The interatomic O8...O1 and O5...O4 distances determined experimentally are equal to 2.551 Å and 2.589 Å, respectively. This could also be the reason why an opposite tendency concerning the electron density distribution was observed for compound **2** (Bridge 2). The values of electron density and its Laplacian at Bond Critical Points (BCPs) of intramolecular hydrogen bonds for both compounds are shown in Table 2. The electron density values at the hydrogen bridge BCPs are consistent with our previous calculations performed for 2,3-dichloronaphthazarin [71]. The covalent O8-H*BP*1/O5-H*BP*<sup>2</sup> bonds are stronger than those formed after proton transfer (O1-H*BP*1/O4-H*BP*2). This observation was made after the electron density and its Laplacian examination at BCPs. The electron density values at BCPs are higher for the OH covalent bonds in the molecular forms of compounds **1** and **2**. However, the intramolecular hydrogen bonds are stronger (higher electron density values at BCPs) for the proton transfer (PT) forms. The values at BCPs obtained based on AIM theory are rather similar for molecular and PT forms. They do not much differ, even comparing compound **1** with compound **2**. This could sugges<sup>t</sup> that the proton transferred form is best described not as O<sup>−</sup>...<sup>+</sup>H-O, but as simply O...H-O, in parallel with the molecular form O-H...O. The topology maps of electron density are presented in Figure 3. They contain molecular properties common for the AIM description of the electronic structure: critical points (BCPs and RCPs), which are stationary points of the electron density field (i.e., the density gradient is zero at the critical point). In the graphical presentation of Figure 3, these critical points are recognizable as maxima (nuclear positions), saddle points, and minima. Due to the presence of intramolecular hydrogen bonds, the typical quasi-rings were formed and recognized by the BCPs of covalent bonds and the indicated bond paths of the hydrogen bridges. In addition, in the case of compound **2**, some intramolecular contacts were detected between the methoxy groups as well as between the hydrogen of the methoxy group with the O4 atom from the second hydrogen bridge. The presence of an intramolecular hydrogen

bonds stabilizes the conformation of the molecules. However, the topology maps showed (in the case of compound **2**), that the electron density distribution in the hydrogen bridge (Bridge 2) could be affected by the competitive interactions introduced by the methoxy substituents. As is shown, two additional quasi-rings were found, indicating that the C-H...O intramolecular hydrogen bonds were formed. They are characterized by the presence of the BCPs and RCPs. However, the presence of such intramolecular interactions was not identified experimentally in the crystal structure [65]. Therefore, the presence of the interactions could be driven by steric effects and degrees of freedom introduced to the isolated molecule model.


**Table 1.** Atoms In Molecules (AIM) atomic charges calculated for selected atoms of the studied compounds, compounds **1** and **2**, and their proton transferred (PT) forms at the B3LYP/6-311++G(2d,2p) level of theory.

**Table 2.** Atoms In Molecules (AIM) Bond Critical Point properties calculated for selected bonds of the studied compounds, compounds **1** and **2**, and their proton transferred forms (PT) at the B3LYP/6-311++G(2d,2p) level of theory. Electron density *ρBCP* is given in *e* · *a*<sup>−</sup><sup>3</sup> 0 atomic units, and its Laplacian <sup>∇</sup>2*ρBCP* is given in *e* · *a*<sup>−</sup><sup>5</sup> 0 units.


**Figure 3.** Topology maps of electron density obtained on the basis of AIM theory at the B3LYP/6- 311++G(2d,2p) level of theory for compounds **1** (**a**) and **2** (**b**). The molecular (left) and proton transferred (right) forms are presented. The black solid and dashed lines indicate the intramolecular interaction paths. The green and red dots mark the presence of BCPs and RCPs, respectively.

#### *2.2. Intermolecular Forces in Naphthazarin Derivatives Dimers Based on Symmetry-Adapted Perturbation Theory (SAPT)*

The presence of two distinct types of stacked dimers in the crystals of **1** or **2** (antiparallel vs. parallel arrangemen<sup>t</sup> of molecules, respectively; see Figure 4) indicates that even the relatively mild substitution can influence crystal packing forces. It is therefore necessary to investigate the molecules of **1** and **2** on the basis of interaction energy partitioning schemes. Symmetry-Adapted Perturbation Theory (SAPT, see Ref. [70]) has become a de facto standard for such investigations, although many other approaches exist, some of them capable of tackling covalent bonding, for example, a DFT-based energy decomposition analysis [72] or localized orbital energy decomposition, LMOEDA, useful for non-covalent forces such as beryllium bonds [73]. While analyzing the results presented in this section, it is necessary to remember that SAPT is a perturbative approach, in which intra- and intermonomer correlations are treated separately. SAPT0 omits the intramonomer correlation, while SAPT2 includes this effect up to the second order of perturbation. Both levels account for the intermonomer electron correlation, which is the source of polarization and dispersion contributions. The SAPT partitioning divides the interaction energy into "static" contributions (electrostatic interaction of frozen electron densities, and Pauli exchange repulsion) and "correlated" terms (induction–mutual polarization of monomers—and dispersion).

The crystal structures of both compounds contain the following basic types of dimeric structures, depicted in Figure 5: d1—the head-to-head planar arrangements, d2—molecular planes tilted at a shallow angle, d3—stacking, d4—planar arrangemen<sup>t</sup> with a side-to-side skew (present only for **2**). Direct use of the experimental structures in the SAPT calculations leads to the results gathered in Table 3, while the DFT-optimized structures are described in Table 4. It must be stressed that the DFT optimization leads to the collapse of the d2-type dimers into the stacking arrangement, underlining the role of the confinement of molecules leading to the formation of diverse structural motifs.

**Figure 4.** Two distinct types of stacking in the crystal structures of **1**, 2,3-dimethylnaphthazarin (antiparallel stacking), and **2**, 2,3-dimethoxy-6-methylnaphthazarin (parallel stacking). The interplanar distances are 3.614 Å and 3.442 Å, respectively.

**Figure 5.** Dimers extracted from the crystal structures of compounds **1** (upper part) and **2** (lower part), used in the SAPT study.

**Table 3.** SAPT2/jun-cc-pVDZ results of energy partitioning for the dimers of compounds **1** and **2** (see Figure 5) with structures taken directly from the diffraction experiments. All energy terms in kcal/mol: Elst—electrostatics; Exch—exchange (Pauli) repulsion; Ind—induction (polarization); Disp—dispersion; SAPT0 and SAPT2 are defined according to Ref. [74].


**Table 4.** SAPT2/jun-cc-pVDZ results of energy partitioning for the dimers of compounds **1** and **2** (see Figure 5) with structures taken from the DFT structural optimization. All energy terms in kcal/mol: Elst—electrostatics; Exch—exchange (Pauli) repulsion; Ind—induction (polarization); Disp—dispersion; SAPT0 and SAPT2 are defined according to Ref. [74].


The results gathered in Tables 3 and 4 show that the energetically most important structural motif (stacking dimer d3) is formed with the dominant role of dispersion. The role of dispersion is visible especially when the d3 dimers of experimental solid state structure are compared with their DFT-optimized analogues. Surprisingly, the latter are more strongly bound. This is an outcome of two competing factors. On the one hand, the crystal electrostatic and steric field tends to squeeze the molecules together, so that no empty voids remain in the structure. This promotes smaller intermolecular separations and stronger stacking forces. On the other hand, the presence of neighbouring molecules means that the capacity of the molecule to interact with its neighbours must split between much more interactions than in the dimer. The latter factor prevails, and the DFT-optimized stacking dimers are bound stronger by ca. 4–5 kcal/mol than their crystal structure equivalents. It is interesting to note that the DFT-optimized d3 structures exhibit not only stronger dispersion, but also electrostatic and induction contributions.

It seems paradoxical that the d1, d2 and d4 dimers, relying mostly on electrostatic forces including hydrogen bonds, present more equalized distribution of the interaction energy terms than the stacked d3 dimers (both in the gas phase and in the arrangemen<sup>t</sup> from the crystal structure). However, there is another factor which is closely related to the type of force dominating the interactions. Two levels of theory, SAPT0 and SAPT2, are provided in Tables 3 and 4 to explain this factor. We note that the total SAPT0 and SAPT2 interaction energies are very close to each other when the studied molecules do not engage in hydrogen bonding, highlighting the role of intramonomer electron correlation in the formation of hydrogen bonds. For example, the d3 dimer of **1** has the interaction energy of −13.681 kcal/mol at the SAPT0 level and −13.196 kcal/mol at the SAPT2 level. This means that the hydrogen bonds and electrostatic forces, displaying larger differences between the SAPT0 and SAPT2 energies, contain significant contributions of higher-order corrections connected with electron correlation, not present at the SAPT0 level. On the other hand, the presence of hydrogen bonds in the dimers is associated in this case with a relatively weak interaction (ca. 4–5 kcal/mol). The weakest dimers (d2 type) are rather just multipolar, electrostatic contacts and their particular shape is governed by steric hindrance of the substituents (especially for **2**).

Summarizing the SAPT study, we stress that the stacked arrangemen<sup>t</sup> is the principal structural motif of the crystal from the geometric point of view. This fact agrees with the role of dispersion forces as the most important factor from the energetic point of view. However, the details of the solid state structure are modified by the substituents and the polar nature of the compounds introduced not only by the intramolecular hydrogen bonds, but also by the substituents, even relatively mild on the Hammett scale (the methyl groups in **1**, with classical Hammett constants of only −0.07 and −0.17 for *σm* and *<sup>σ</sup>p*, respectively [63]).

#### *2.3. First-Principle Molecular Dynamics (FPMD) in the Gas and Crystalline Phases*

The applied Car–Parrinello molecular dynamics (CPMD) enables the investigation of molecular and spectroscopic features of the naphthazarin derivatives based on ab initio Potential Energy Surface (PES), which is of grea<sup>t</sup> importance when we are expecting to register proton transfer phenomena events. The time-evolution study provides an insight into the dynamical nature of the hydrogen bonding present in the studied systems. Therefore, special attention was paid to the intramolecular hydrogen bridges present in both compounds. The CPMD simulations were performed in the gas phase and in the solid state. The two phase study enabled detection of differences related to the environmental effects' influence on the hydrogen bond dynamics, e.g., the crystal field and the presence of neighbouring molecules. The details of the hydrogen bonds' average metric parameters are presented in Table S2. The reported values are in good agreemen<sup>t</sup> with the experimental data available [64,65], as well as the static DFT results.

Figures 6 and 7, showing the time evolution of the distances related to the hydrogen bridges, use the same color coding to aid data interpretation. The black line corresponds to the O...O donor–acceptor distance, and it simply oscillates around an equilibrium value throughout the simulation time. The red line is the donor-proton bond length, of rather small amplitude, while the green line is the proton-acceptor distance, oscillating in a wide range. The objects of our interest, proton transfer events, are accompanied by a sudden increase in the donor-proton bond length, accompanied by the lowering of the proton-acceptor separation (the red and green lines cross over).

In Figure 6, the hydrogen bridge dynamics is presented for compound **1**. The upper part shows data obtained in the gas phase. There are many proton-sharing events during the CPMD simulation run. The bridged protons exhibit strong mobility, which results in proton transfer phenomena registered during the 20 ps of the CPMD run. The protons moved to the acceptor-atom side, stayed there for a short time and kept moving again to the proton-donor side. In the solid state (lower part of the Figure 6), the proton transfer events were noticed as well. However, there were less proton-sharing events comparing to the gas phase results (see also the discussion of proton possession statistics two paragraphs below). This could be explained by the presence of intermolecular hydrogen bonds and molecular overlapping present in the crystal structure [64]. The presence of neighbouring molecules forming intermolecular hydrogen bonds (O-H...O), where the O-H group is involved in the intramolecular hydrogen bond, as well as interacts with an oxygen atom (proton-acceptor from another molecule) and introduces competition in the interactions. This shows a significant difference between the isolated molecule and the crystalline phase dynamics, where many factors are included during the CPMD simulations. There is a visible correlation in the bridged protons dynamics in both phases. Compound **1** exhibits symmetry; therefore, one could expect that the dynamical nature of the bridged protons will be similar, but it will depend on the phase discussed—the crystal packing lowers the effective symmetry perceived by the analyzed molecule.

**Figure 6.** Hydrogen bridge structural parameters during the CPMD simulation of 2,3- dimethylnaphthazarin (**1**). The graphs show gas phase results for (**a**) Bridge 1 and (**b**) Bridge 2, and solid state results for (**c**) Bridge 1 and (**d**) Bridge 2. For atom numbering scheme, see Figure 1.

**Figure 7.** Hydrogen bridge structural parameters during the CPMD simulation of 2,3-dimethoxy-6- methylnaphthazarin (**2**). The graphs show gas phase results for (**a**) Bridge 1 and (**b**) Bridge 2, and solid state results for (**c**) Bridge 1 and (**d**) Bridge 2. For atom numbering scheme, see Figure 1.

The CPMD results concerning the intramolecular hydrogen bond dynamics of compound **2** are presented in Figure 7. The bridged protons exhibit strong mobility in both studied phases. In the gas phase (upper part of Figure 7), there are frequent proton-sharing events and proton transfer phenomena were noticed as well. During the 20 ps run, there were 3 ps long proton transfer events, and after this time, the bridged protons moved back again to the proton-donor atom. There is also a correlation in the bridged protons dynamics—the 3 ps long PT events happened at the same time for both bridges. A solid state study provided a different picture of the proton mobility in the hydrogen bridges (lower part of Figure 7). The bridged protons were strongly delocalized between the donor and acceptor atoms in both hydrogen bonds. The compound did not exhibit symmetry due to the presence of the methyl group in the sixth position as well as methoxy groups, which are not equivalent [65]. There were also intermolecular hydrogen bonds, involving (similarly to compound **1**) OH groups from the intramolecular hydrogen bonds and protonacceptor atoms from the neighbouring molecules. There was also molecular overlapping according to the X-ray measurements [65]. Comparing gas phase results with the solid state of compound **2**, it is visible from the structural data analysis that external forces influence the bridged protons dynamics. In the case of compound **2**, we can draw the conclusion that the presence of methoxy groups and the lack of symmetry introduced inductive and steric effects, which provided us with a different dynamical nature of the intramolecular hydrogen bonds present in compound **2** with respect to **1**.

The diverse properties of hydrogen bonds were further analyzed using statisticsbased approaches. First, we calculated the proton possession statistics, i.e., percentages of the time spent by the given bridge proton at the donor or acceptor site. The proper association of the proton with its site at a given time was determined by the Voronoi geometric criterion—donor-proton vs. proton-acceptor distance comparison. The results, gathered in Table 5, indicate that in case of compound **1**, gas phase and solid state statistics are very similar. The degree of convergence of the dynamics trajectory can be estimated

by comparison of the two equivalent bridges, O8...O1 vs. O5...O4—the differences are no more than 0.5%. The differences between the gas phase and solid state results are 0.2% for the O8...O1 bridge and 1.1% for the O5...O4 bridge, which is close to the difference between the two bridges. This means that the solid state environment does not seem to change the overall partitioning of the proton residence time between the donor and the acceptor sites, and it promotes slower dynamics (fewer proton sharing events). The results for the gas phase CPMD simulation of **2** are also similar to the case of **1**: the two bridges, which are not equivalent, are still similar enough to provide the same statistics of proton possession. The solid state case is more interesting: the packing forces (the presence of neighbours and their electrostatic field) lead to almost equally shared protons; however, the H*BP*<sup>1</sup> tends to reside more at the acceptor site than the H*BP*<sup>2</sup> proton.


**Table 5.** Proton possession statistics for the CPMD runs. Percentages of the time spent by the bridge proton at the donor or acceptor site, determined by Voronoi geometric criterion–distance comparison.

The second part of the CPMD trajectory statistical analysis is provided by the histograms for the donor-proton positions in the two hydrogen bridges—see Figure 8. The histograms (probability density plots) show how the proton positions are correlated in the sense of averaging over the CPMD run. It is visible that for compound **1**, regardless of the simulation conditions—the gas phase or the solid state—the two protons H*BP*<sup>1</sup> and H*BP*<sup>2</sup> are strongly correlated and located mostly at the donor site. This is also true for the gas phase trajectory of **2**. These results are in agreemen<sup>t</sup> with the data for naphthazarin [61]: it was shown that when an asynchronous proton jump occurs, it is very probable that a second proton transfer will follow within a few O-H oscillation periods. From this point of view, it is interesting to note that the solid state simulation of **2**, where the protons are more delocalized, also exhibits important motion correlations. The histogram shown in panel (d) of Figure 8 consists of four more populated regions forming a square shape. These regions correspond to the molecular form, the PT form, and the two less stable forms with single-proton transfer. There are no indications of a synchronous double-proton transfer, which would result in formation of a populated region in the center of the square shape.

While the current study was carried out within the Newtonian classical nuclear dynamics, corresponding to the Born–Oppenheimer picture, it is recognized that in some instances the nuclear quantum effects (NQE) are important for qualitative and quantitative agreemen<sup>t</sup> with experiments. An excess proton in water migrates due to a complicated mechanism in which quantum fluctuations, rather than tunneling, play a crucial role [75,76]. Quantum disorder in the hydrogen bonds is required to explain the X-ray absorption spectra of water and ice [77]. Intramolecular hydrogen bonds can be diversely affected by nuclear quantization. Picolinic acid N-oxide with a very strong O-H...O bond requires anharmonic, quantum treatment of the proton motion to rationalize enormous red shifts of the ν*OH* mode [78]. Weaker hydrogen bonds, such as those in o-hydroxy Mannich bases, exhibit a single-well potential with the minimum clearly at the donor side [79], while in the N-oxides of Mannich bases the potentials are very flat and broad, allowing the proton

to move almost freely within the bridge [45]. The current study contains an assessment of the importance of nuclear quantum effects for the H*BP*<sup>1</sup> proton. The results, shown in Figure 9, are obtained with the snapshot-based a posteriori approach [79,80] involving numerical solution of a vibrational Schrödinger equation [81]. Our attention was focused on the impact of the NQE phenomena on the O8-H*BP*<sup>1</sup> distance. It is visible in Figure 9 that the NQE tend to increase the proton delocalization between the donor and acceptor sites, making the H*BP*<sup>1</sup> atom shift towards the center of the O8...O1 bridge (the red crosses, indicating the NQE-corrected positions, are located closer to the half of the actual O8...O1 distance than are the green circles—classical positions). For each of the four investigated cases, one of the snapshots presents the PT structure, where the O8-H*BP*<sup>1</sup> distance is larger than 1.5 Å. In such cases, the NQE shift the proton position towards lower O8-H*BP*<sup>1</sup> values. The impact of NQE is not decisive for the proton localization in the studied compounds **1** and **2**, with a very interesting exception of the crystalline phase of **2**. The proton at the O8...O1 distances above 2.5 Å behaves in a way similar to the other cases, but at 2.48 Å the impact of NQE is particularly large. The same distance for **1** and gas phase **2** does not lead to such large NQE; therefore, it seems that this is the precise region of bridge length at which the combination of the molecular structure of **2** and the crystal environment make the NQE (including tunneling) very effective. However, when the bridge is compressed even further—to 2.37 Å—the impact of NQE is again very small. The explanation is as follows: at such a short bridge length, the proton potential is already of the flat single-well type, making this bridge temporarily a "low-barrier hydrogen bond" for which the tunneling effects are negligible [75]. As a final remark to the study of NQE, we note that the classical CPMD trajectory is able to sample this region of the molecular phase space, as seen in Figure 8. This fact indicates that the NQE should not have a qualitative impact on the properties of the investigated systems.

**Figure 8.** Histograms for the donor-proton distances in the two hydrogen bridges of the studied compounds—results of the CPMD simulation for (**a**) **1** in the gas phase, (**b**) **1** in the solid state, (**c**) **2** in the gas phase, (**d**) **2** in the solid state. Color scale represents probability density in Å−2.

**Figure 9.** Impact of nuclear quantum effects for the H*BP*<sup>1</sup> bridge proton on the O8-H*BP*<sup>1</sup> distance. Green circles—classical value of the distance; red crosses—quantum expectation value of the O8-H*BP*<sup>1</sup> distance operator. Results of a posteriori quantum treatment of CPMD trajectory for (**a**) **1** in the gas phase, (**b**) **1** in the solid state, (**c**) **2** in the gas phase, (**d**) **2** in the solid state.

Vibrational signatures of the bridged protons, corresponding to the *ν*OH at the highwavenumber region, are presented in Figure 10. Since the IR spectra of these compounds are not available in the literature, the most natural source of comparison is the parent compound, naphthazarin. Investigation of the bridge proton features has two main goals. First, it is possible to trace the presence of strong interactions in the crystal. On the other hand, the bridges in **2** are not symmetrical, and their asymmetry can lead to slightly diverse positions of the normal modes. This is an interesting issue in relation to the parent naphthazarin itself, where the skeleton, devoid of the substituents, does not prefer any of the proton positions. The broad absorptions of the *ν*OH/OD and *γ*OH/OD stretching modes were experimentally identified in naphthazarin at 3060/2200 cm<sup>−</sup><sup>1</sup> and 793/560 cm<sup>−</sup>1, respectively; the upper wavenumber region is the most relevant for the fast proton dynamics corresponding to the stretching mode [62].

The first goal, detection of strong interactions, can be accomplished by comparison of the *ν*OH band positions. Compound **1** exhibits similar positions of this band in the gas phase (from 2300 to 3400 cm<sup>−</sup>1) and crystal (from 2200 to 3400 cm<sup>−</sup>1). The band center at ca. 2800 cm<sup>−</sup>1–2900 cm<sup>−</sup><sup>1</sup> is at a slightly lower wavenumber than the experimental value of 3060 cm<sup>−</sup><sup>1</sup> for naphthazarin [62]. The lower wavenumber absorptions, 700–1700 cm<sup>−</sup><sup>1</sup> in the gas phase and 600–1800 cm<sup>−</sup><sup>1</sup> in the solid state, should be attributed to the mechanical influence of the heavy-atom motions. These values indicate on the one hand a middlestrong O-H...O hydrogen bonding, and on the other a relatively small impact of the crystal packing effects on the vibrational features of **1**. These facts agree well with the not too frequent proton transfer events in this compound (see Figure 6, which also confirms that the PT occurrence in the gas phase and in the crystal is very similar). We have noted already that the PT events are not strictly synchronous, but they are strongly correlated. This makes

the vibrational signatures of H*BP*<sup>1</sup> and H*BP*<sup>2</sup> virtually identical. This is not strictly true for compound **2**, where the chemical nature of the substituents is different in the vicinity of H*BP*<sup>1</sup> than in the vicinity of H*BP*2. The difference is almost not visible in the results of the gas phase simulation of **1**—the *ν*OH vibrational features of both protons fall into the 2300 cm<sup>−</sup><sup>1</sup> to 3300 cm<sup>−</sup><sup>1</sup> range (the 700 cm<sup>−</sup>1–1700 cm<sup>−</sup><sup>1</sup> region is associated with heavy atom motions, as already noted for compound **1**), and the signature of H*BP*<sup>1</sup> is centered at ca. 2800 cm<sup>−</sup>1, while the signature of the H*BP*<sup>2</sup> proton peaks at ca. 100 cm<sup>−</sup><sup>1</sup> has a lower wavenumber. The difference is small, and it is also in agreemen<sup>t</sup> with the time evolution of the distance parameters (see Figure 7). The lowering of the band center position with respect to naphthazarin (3060 cm<sup>−</sup><sup>1</sup> in the experimental spectrum) is also not large. Quite unexpectedly (if one has not ye<sup>t</sup> appreciated the solid state distance parameters shown in Figure 7), the crystal field makes the bridge protons very strongly delocalized. The resulting vibrational signature is extremely broad and forms a continuous background feature from ca. 700 to 3100 cm<sup>−</sup>1. This feature does not differentiate the two bridge protons. The reason for such a profound change in the bridge proton dynamics should be sought after for a particular arrangemen<sup>t</sup> of molecules in crystal; thus, the competition between inter- and intramolecular contacts turns out to be cooperation in the case of the solid state compound, compound **2**.

**Figure 10.** Vibrational signatures (atomic velocity power spectra) of the bridge protons calculated from the CPMD simulation of **1** and **2**. In the case of **2**, signatures of the non-equivalent bridge protons are presented as separate curves placed back to back. For atom numbering scheme, see Figure 1.
