**4. Conclusions**

We have presented computational results of two naphthazarin derivatives substituted with methyl and methoxy groups in diverse manner. We have examined various factors influencing the molecular features exhibited by the aforementioned derivatives in relation to the properties of the substituents and symmetry breaking by their introduction. The presence of the substituents and changes in the chemical composition have led to changes in the bridged proton dynamics and intermolecular interactions in comparison to the parent compound, naphthtazarin. The computations were performed in the electronic ground state, both in the gas phase and solid state. In order to shed light on the

intermolecular interactions, the dimers of compounds **1** and **2** were investigated. Our computational findings were compared with the experimental data available (structural and spectroscopic). The application of the DFT method with three different functionals, each using a 6-311++G(2d,2p) basis set, complemented with the single-point MP2 and CCSD calculations with the def2-TZVP basis set, provided information of the proton reaction path and the energy barrier for the proton transfer. The highest DFT energy barrier equals ca. 10 kcal/mol, while MP2 and CCSD provided the barrier heights of ca. 6.8 and 9.6 kcal/mol, respectively. Moreover, two energy minima were located in both molecules and in both examined hydrogen bridges. The application of the AIM theory gave a quantitative picture of the electron density distribution in the molecular and proton transferred forms of the studied compounds. The topological analysis confirmed the presence of the intramolecular hydrogen bonds (in agreemen<sup>t</sup> with experimental X-ray findings in the literature). Additionally, it was shown, on the basis of electron density and its Laplacian values, that the hydrogen bonds are stronger in the tautomeric PT forms. The SAPT analysis gave an insight into energy partitioning and provided information on the primary factors responsible for dimer stabilization. It was found that the primary factors are the dispersive forces. Using the SAPT method, we could identify and describe quantitatively external forces influencing the molecular features of compounds **1** and **2**. The CPMD results showed that protons in the hydrogen bridges are very labile. Proton transfer phenomena were observed in the gas phase as well as in the solid state. In compound **2**, there is a clearly visible influence of environmental factors on the hydrogen bridge dynamics. The vibrational analysis confirmed, by the broad absorption regions observed in the computed power spectra, a strong anharmonicity of the studied hydrogen bonds as well as their dynamics. It is especially visible in compound **2**, where in the solid state only one very broad absorption (700 cm<sup>−</sup>1–3100 cm<sup>−</sup>1) region was found. The incorporation of nuclear quantum effects to the hydrogen bridges showed a stronger delocalization of the bridged protons, especially at shorter, but not the shortest, distances between the donor and acceptor heavy atoms.

**Supplementary Materials:** The following are available online. Figure S1. The structures of the investigated naphthazarin derivatives: 2,3-dimethylnaphthazarin (**1**) and 2,3-dimethoxy-6-methylnaphthazarin (**2**), with atom numbering scheme for hydrogen bridges. Coloring scheme: oxygen atom—red, carbon atom—grey and hydrogen atom—white. Figure S2. The models for gas phase and solid state CPMD simulations. Left—the isolated molecule model of 2,3-dimethylnaphthazarin (**1**); right—the model used for solid state simulations of 2,3-dimethoxy-6-methylnaphthazarin (**2**). Table S1. Energy for compounds **1** and **2** with different proton positions in the hydrogen bridges computed using DFT method. Electronic as well as vibrational zero point-corrected values are given. Table S2. Selected geometric parameters related to the intramolecular hydrogen bonds of 2,3-dimethylnaphthazarin (**1**) and 2,3-dimethoxy-6-methylnaphthazarin (**2**). Comparison of experimental and computed data. Metric parameters are given in Å and degrees. CPMD results are presented as average ± standard deviation. Sets of coordinates for the minima and transition state estimates from the DFT scans (XYZ format).

**Author Contributions:** Conceptualization, A.J.; methodology, A.J. and J.J.P.; computations, K.K., M.P., A.J. and J.J.P.; data analysis, K.K., M.P., A.J. and J.J.P.; writing—original draft preparation, K.K., M.P., A.J. and J.J.P.; writing—review and editing, K.K., M.P., A.J. and J.J.P.; supervision, A.J. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was supported by the National Science Centre gran<sup>t</sup> number UMO-2015/17/B /ST4/03568.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** All the relevant processed data (energy values, AIM data, SAPT energy terms, structural information, time evolution of distances in the hydrogen bridges, vibrational signatures) are reported within the manuscript.

**Acknowledgments:** The authors thank the Wrocław Center for Networking and Supercomputing (WCSS) in Wrocław and Academic Computer Centre Cyfronet AGH (PL-Grid infrastructure Prometheus) in Kraków for generous allocation of computing time and technical support. In addition, A.J. and J.J.P. thank the National Science Centre (Poland) for supporting this work under the gran<sup>t</sup> no. UMO-2015/17/B/ST4/03568.

**Conflicts of Interest:** The authors declare no conflict of interest.
