**3. Computational Methodology**

*3.1. Static Models on the Basis of Density Functional Theory (DFT)*

The models of monomers and dimers were constructed on the basis of X-ray structures of 2,3-dimethylnaphthazarin (**1**) (CCDC deposition number—1125030) and 2,3-dimethoxy6-methylnaphthazarin (**2**) (CCDC deposition number—1161869) [64,65,82]. The geometry optimization for the molecular forms of monomers was performed using Density Functional Theory (DFT) [66,67] and three functionals: B3LYP [83], PBE [84,85] and *ω*B97XD [86] with valence-split triple-zeta Pople's style basis set denoted as 6-311++G(2d,2p) [87,88]. The choice of functionals was devised to represent the current spectrum of the most widely used approaches: the PBE functional is of the Generalized Gradient Approximation (GGA) type used frequently in the context of plane-wave calculations (including Car–Parrinello MD), and does not use the exact exchange. On the other hand, B3LYP is a hybrid functional, and so is the *ω*B97XD, but the latter includes empirical dispersion correction. Following the geometry optimization, harmonic frequencies were computed to confirm that the obtained structures correspond with the minimum on the Potential Energy Surface (PES). Additionally, models with diverse proton positions were constructed and optimized as well using the DFT method (for details, see Table S1 of the Supplementary Information). In the next step, the single-point simulations at the MP2 [89] and CCSD [90,91] levels with def2-TZVP basis set [92] were carried out for the structures of the minima and transition state on the PT pathway. Next, the structures with OH groups on the proton-donor side were taken to investigate the proton potential paths using the scan method with geometry optimization (the O-H increment was set to 0.05 Å, the O8H*BP*1O1 and O5H*BP*2O4 valence angles were frozen while the remaining parts of the molecules were optimized). The results of the scans formed a discrete set of points, from which a proton potential function was derived. Thus, the barrier height is determined with accuracy depending on the discrete steps of energy in the vicinity of the transition state; the error estimate is the internal property of the procedure based on the discrete series of points, not the absolute uncertainty of a particular DFT functional. The zero-point vibrational correction is not included in the reported values. Finally, the wavefunctions for the Atoms In Molecules (AIM) theory [69] analysis were prepared with assistance of the B3LYP functional and 6-311++G(2d,2p) basis set for molecular and proton transferred forms of monomers. The theory was applied for the electronic structure as well as molecular topology investigations. Special attention was paid to the electron density and its Laplacian values at Bond and Ring Critical Points (BCPs and RCPs) related to the intramolecular hydrogen bonding. Next, for the dimeric structures extracted from the crystal data of compounds **1** and **2** [64,65], the energy minimization was performed using the *ω*B97XD functional [86] and 6-311++G(2d,2p) basis set. The simulations were carried out in the gas phase with the Gaussian 09 rev. D.01 [93] and Gaussian 16 rev. C.01. suite of programs [94]. The single-point MP2 and CCSD calculations were conducted with the Turbomole 6.5 program [95]. The AIM analysis was performed using the AIMAll program [96]. In addition, sets of coordinates are provided in the Supplementary Information for the current study.

#### *3.2. An Application of Symmetry-Adapted Perturbation Theory (SAPT) to Dimers*

The Symmetry-Adapted Perturbation Theory (SAPT) [70] enables energy decomposition between interacting molecules, in our case dimers. The method divides an exact Hamiltonian into Hartree–Fock contribution of monomers, *F*ˆ *A* and *F*ˆ *B*, correlation components interacting inside the monomers, *W*ˆ *A* and *W*ˆ *B*, and the contribution covering interaction between monomers, *V*ˆ :

$$
\hat{H} = \hat{F}\_A + \hat{F}\_B + \hat{W}\_A + \hat{W}\_B + \hat{V} \tag{1}
$$

An important advantage of the SAPT scheme is the fact the individual components could be grouped into four principal groups with precisely defined physical interpretation: (i) electrostatic (*Eelst*)—approximate Coulombic interactions of electron density decomposition of isolated monomers (without the effect of polarization by the neighboring molecule); (ii) exchange (*Eexch*—which is the short-range Pauli repulsion; (iii) Induction (*Eind*) and exchange-induction (*Eex*−*ind*—which is based on mutual polarization of the monomers; (iv) dispersion (*Edisp*)—consideration of short-lived instantaneous multipoles. Depending on the considered energy components, the SAPT hierarchy of interactions is obtained. The

SAPT levels most commonly used are SAPT0 (in agreemen<sup>t</sup> with Hartree–Fock method) and SAPT2 (with accuracy approximate to the MP2 method):

$$E\_{SAPT0} = E\_{elst}^{10} + E\_{exch}^{10} + E\_{ind,r}^{20} + E\_{ex-ind,r}^{20} + \\\delta E^{HF} + E\_{disp}^{20} + E\_{ex-disp}^{20} \tag{2}$$

$$E\_{SAPT2} = E\_{SAPT0} + E\_{elst,r}^{12} + E\_{exch}^{11} + E\_{exch}^{12} + {}^tE\_{ind}^{22} + {}^tE\_{ex-ind}^{22} \tag{3}$$

These equations show the fundamental difference between the SAPT0 and SAPT2 approximations: the SAPT0 components never use intramonomer electron correlation, so—generally speaking—the resulting components of interaction energy are based on the non-correlated Hartree–Fock wavefunctions of the monomers. SAPT2, on the other hand, includes intramonomer correlation up to the second perturbative order, which is especially important for very weak interactions. In our experience with hydrogen-bonded systems, SAPT0 results are overestimated in comparison to the more accurate SAPT2 approach, but the general trends are reproduced with quite a high degree of correlation between the methods. Regarding the computational efficiency and memory requirements, SAPT2 can be prohibitively demanding for systems of ca. 60 atoms. However, due to the electron density expansion on specially fitted basis functions (density fitting technique), the SAPT0 computational cost is comparable to the MP2 method.

The energy decomposition of the naphthazarin derivative dimers (see Figure 5) was performed for: (i) data extracted from the X-ray structures of the investigated compounds [64,65] in order to reproduce the intermolecular forces in the crystal structure responsible for the crystal unit cell arrangement; (ii) the data obtained as a result of gas phase DFT simulations at the *ω*B97XD/6-311++G(2d,2p) level of theory. The interaction energy was calculated at the SAPT2/jun-cc-pVDZ level of theory (truncation of the diffuse functions in the jun-cc-pVDZ basis is derived in [97]). The basis set superposition error (BSSE) correction [98] was included in the simulations of the dimers (the studied dimers were divided into "monomers" in order to fulfil the requirements of the Boys–Bernardi method). The SAPT calculations were carried out using the Psi4 1.2.1 [99] program.

#### *3.3. Car-Parrinello Molecular Dynamics in the Gas Phase and Solid State*

The dynamical nature of the studied naphthazarin derivatives (compounds denoted as **1** and **2**, see Figures 1 and S2) [64,65] were examined in the light of First-Principle Molecular Dynamics (FPMD) method. The simulations were performed for the isolated molecules as well as for the molecular crystals. The gas phase simulations results were further used for the comparative study of differences introduced by the interatomic forces present in the solid state. Our attention was placed on the intramolecular hydrogen bonds' dynamics and properties. We have analyzed the hydrogen bridges dynamics as a function of simulation time. For this purpose, detailed analysis of metric parameters was performed for O1...O8/O5...O4 interatomic distance, O1-H*BP*1/O2-H*BP*<sup>2</sup> covalent bonds and H*BP*1...O8/H*BP*2...O4 intramolecular hydrogen bonds in compound **1**. Compound **1** is symmetric; therefore, we could expect that the bridged proton dynamics will be similar. However, we placed emphasis on a detailed view of protons motion in the hydrogen bridges. Compound **2** has a broken symmetry due to the presence of the CH3 substituent. Both hydrogen bridges were taken into consideration in the analysis of metric parameters. We were looking for any correlations in the hydrogen bridge dynamics. Another aspect related to the data analyses were vibrational signatures provided by the OH groups. The Fourier transformation of the autocorrelation function of atomic velocity was employed to develop power spectra. The models used for Car–Parrinello molecular dynamics (CPMD) in the gas phase are presented in Figure 1. The initial geometries for the isolated molecules were extracted from the X-ray data [64,65] and placed in cubic boxes with a = 15 Å for compound **1** and a= 16 Å for compound **2**. The models for CPMD in the solid state were prepared on the basis of crystallographic unit cells [64,65]. The unit cell dimensions for compound **1** are as follows: a = 16.429 Å, b = 6.524 Å, c = 9.136 Å and *β* = 90.19◦ with Z = 4, while for compound **2**, a = 3.873 Å, b = 20.21 Å, c = 15.00 Å and *β* = 96.05◦ with Z = 4. The computational setup for the simulations in both studied phases was prepared bearing in mind the fact that intramolecular hydrogen bond dynamics were being studied. The simulations were divided into geometry optimization of the studied compounds, **1** and **2**, and subsequent CPMD runs in the gas phase and solid state. The exchange correlation functional by Perdew, Burke and Ernzerhof (PBE) [84,85] and Troullier–Martins [100] pseudopotentials were applied. The fictitious electron mass (EMASS) was equal to 400 a.u. and the time-step was set to 3 a.u. The kinetic energy cutoff for the plane-wave basis set was 80 Ry. The CPMD calculations were performed at 295 K, controlled by Nosé–Hoover thermostat chain assigned to ions [101,102]; the electronic system was thermostatted at the orbital kinetic energy values determined in separate short non-thermostatted runs for each system. Hockney's scheme [103] was applied to remove interactions with periodic images of the cubic cell during the gas phase dynamics. The translational and rotational movements were removed from the CPMD data collection as well. The crystalline phase CPMD was carried out with Γ point approximation [104] and Periodic Boundary Conditions (PBCs). The real-space electrostatic summation was set to TESR = 8 nearest neighbours in each direction. The CPMD simulations were divided into two parts: (i) equilibration (the initial part of the trajectory—ca. 10,000 steps—was removed from further analyses); (ii) production run, which lasted for 21 ps.

The CPMD simulations were performed using the CPMD 3.17.1 program [105]. The post-processing was carried out using home-made scripts and the VMD 1.9.3 [106] program. The graphical presentation of the obtained results in the current study was conducted with assistance of the VMD 1.9.3 [106] and Gnuplot [107] programs.

#### *3.4. Estimation of the Nuclear Quantum Effects on the Structural Properties in the Gas Phase and Solid State*

The nuclear quantum effects for the bridge proton motion were studied using an a posteriori approach based on the CPMD trajectory [79,80]. In short, the method consists of selecting several snapshots from the CPMD trajectory, calculating proton potential functions for each snapshots, and then, finally, solving the vibrational Schrödinger equation (see, e.g., [79,81]). The particular details for the current study are as follows. Four cases were considered: compounds **1** and **2** in the gas phase and solid state. For each case, five snapshots were extracted from the CPMD trajectory with constant time intervals. For each snapshot, a set of 16 to 20 bridge proton positions (depending on the donor–acceptor distance) was generated for the scan using the donor, proton and acceptor coordinates to define a fragment of an arc. The generated proton positions were then used to calculate single-point energies for the studied systems using the corresponding computational setup of the CPMD code—see the section above. Then, each of the generated proton potential profiles was fitted with a 9th degree polynomial, and a one-dimensional vibrational Schrödinger equation was solved using a grid basis set of 400 points spanning the O8-H*BP*<sup>1</sup> region from 0.7 Å to 2.0 Å. Finally, the expectation value of the O8-H*BP*<sup>1</sup> distance operator at 295 K was calculated taking into account the three lowest-lying vibrational levels. The electronic structure calculations were carried out with the CPMD 3.17.1 program [105], while the quantum vibrational effects were studied with the software developed by Stare and Mavri [81].
