Structure-Property Relationship in Selected Naphtho- and Anthra-Quinone Derivatives on the Basis of Density Functional Theory and Car–Parrinello Molecular Dynamics

: Intra- and inter-molecular interactions were studied in 2,3-dichloro-5,8-dihydroxy-1,4-naphthoquinone and 1,4-dihydroxy-anthraquinone to shed more light on the molecular assembly phenomena. The electronic ground and excited states features of the compounds were investigated to ﬁnd structure-property dependencies. The theoretical study was carried out on the basis of Density Functional Theory (DFT), its Time-Dependent (TD-DFT) extension, and using Car–Parrinello Molecular Dynamics (CPMD). In order to show how the environmental effects modulate the physico-chemical properties, the simulations were performed in vacuo, with the solvent reaction ﬁeld (Polar-izable Continuum Model (PCM) and water as a solvent) and crystalline phase. The intramolecular hydrogen bonds and the bridged proton dynamics were analyzed in detail. The aromatic rings and electronic structure changes were estimated using the Harmonic Oscillator Model of Aromaticity (HOMA) and Atoms in Molecules (AIM) theory. The Symmetry-Adapted Perturbation Theory (SAPT) was employed for interaction energy decomposition in the studied dimers and trimers. It was found that the presence of a polar solvent decreased the energy barrier for the bridged proton transfer. However, it did not signiﬁcantly affect the aromaticity and electronic structure. The SAPT results showed that the mutual polarization of the monomers in the dimer was weak and that the dispersion was responsible for most of the intermolecular attraction. The intermolecular hydrogen bonds seem to be much weaker than the intramolecular bridges. The TD-DFT results conﬁrmed that the electronic excitations do not play any signiﬁcant role in the intramolecular proton transfer. The CPMD results indicated that the protons are very labile in the hydrogen bridges. Short proton transfer and proton-sharing events were observed, and a correlation between them in the twin bridges was noticed, especially for the ﬁrst investigated compound.


Introduction
Naphtho-and anthra-quinone are parent compounds for 2,3-dichloro-5,8-dihydroxy-1,4-naphthoquinone 1 [1] and 1,4-dihydroxy-anthraquinone 2 [2], which are the objects of the current study; see Figure 1. Naphtho-and anthra-quinone occur in many natural products, possessing diverse biological activity, being building blocks of many dyes; see for example [3][4][5][6][7][8][9][10][11][12][13]. They are also interesting classes of compounds from the materials science point of view. Their derivatives have been mentioned as, e.g., potential electrode-active material candidates in an acidic aqueous electrolyte [14]. They are able to improve the performance of flow batteries [15], and also, long-life lithium-organic batteries were developed based on the core structure of naphthoquinone [16]. The anthraquinone oligomers were reported as anode-active materials [17] in rechargeable nickel/oligomer batteries. In addition, naphthoquinone derivatives were found as colorimetric and fluorescence chemosensors for the selective recognition of Sn 2+ in an aqueous medium [18]. The examples mentioned above show how important and interesting from diverse points of view naphtho-and anthra-quinone derivatives are. The main aim of the study covers theoretical investigations of intra-and inter-molecular interactions present in Compounds 1 [1] and 2 [2]. The detailed analysis of such interactions sheds more light onto molecular level processes involving the molecules, e.g., proton transfer phenomenon, aggregation processes, etc., which are important to design new molecules with desired properties. The crystal structure of Compound 1 is built up of stacked molecules, as shown in Figure 2 (see the CCDC Database Identifier DCDHNQ01, Deposition Number 1137406) [1]. The structure organization of the stacked dimers is as follows: due to the symmetry centers, the quinoid part of one molecule partially overlaps the benzenoid part of the neighboring molecules. According to the X-ray measurements, the distance between the overlapping planes of the dimers is 3.40(3) Å, and a charge-transfer complex with itself is formed in the solid-state [1]. Concerning Compound 2 (see Figure 3 for its unit cell), the X-ray measurements detected the presence of remarkably short intermolecular contacts with a O...O distance of 2.644 Å and 2.778 Å (see the CCDC Database Identifier DHXANT10, Deposition Number 1140045) [2]. This was explained by the fact that each O-H group is involved in the intra-, as well as the inter-molecular hydrogen bonding [2]; see Figure S4 of the Supplementary Materials. Both compounds possess strong, Resonance-Assisted Hydrogen Bonds (RAHBs) [19] forming quasi-rings and additionally stabilizing the structures. The presence of the intra-and inter-molecular hydrogen bonds made the compounds interesting and suitable for the discussion of the role of interactions as the building blocks of larger systems. The compounds are well-suited for our Density Functional Theory (DFT) study [20,21] to discuss the hydrogen bonds' properties and the substituents' influence on molecular features. Compounds 1 and 2 were investigated in the electronic ground and excited states. In order to study the Excited State Intramolecular Proton Transfer (ESIPT), the time-dependent extension of the classical DFT method was applied [22]. The electronic ground state computations were performed in vacuo , using the continuum solvation model (Polarizable Continuum Model (PCM)) [23] to reproduce the polar environment impact on the molecules, as well as in the crystalline phase. Special attention was paid to the theoretical description of intramolecular hydrogen bonds. We employed static and first principles molecular dynamics models to get a deeper insight into the hydrogen bridges' dynamics and related phenomena-internal electronic and geometric parameters' reorganization. The Harmonic Oscillator Model of Aromaticity (HOMA) [24] index was estimated to follow the aromaticity changes upon the O-H distance elongation in the hydrogen bridge. The Atoms In Molecules (AIM) theory [25] was applied to study the electronic structure evolution in the same manner. The intermolecular interactions and the energy decomposition in the dimers and trimers of the compounds were investigated using the Symmetry-Adapted Perturbation Theory (SAPT) method [26]. On the basis of the method, the stability of the formed complexes could be estimated and discussed in detail showing which energy components play the most important role. In order to reproduce the molecular properties of both hydrogen bridges, Car-Parrinello Molecular Dynamics (CPMD) [27] was employed. The isolated molecule simulations provided us with the details of the hydrogen bridges' dynamics when the molecule exhibits degrees of freedom not limited by the presence of the neighboring molecules and the crystal field. However, the crystalline phase simulations gave us an unrelieved description of the intramolecular hydrogen bonds' molecular features. The correlations of the bridged protons during the CPMD runs were analyzed to give a deeper insight into the bridged protons dynamics. The molecules' arrangement in the unit cell is shown at the left side, while the intramolecular hydrogen bonds at the right side. The dotted lines indicate the presence of intramolecular hydrogen bonds. The data were used to construct models for CPMD simulations in the crystalline phase [1]. Figure 3. The crystallographic unit cell for Compound 2. The molecules' arrangement in the unit cell is shown at the left side, while the intramolecular hydrogen bonds at the right side. The dotted lines indicate the presence of intramolecular hydrogen bonds. The data were used to construct models for CPMD simulations in the crystalline phase [2].

Static Models on the Basis of Density Functional Theory
The models for Density Functional Theory (DFT) [20,21] simulations were prepared based on the experimentally determined crystal structures of 2,3-dichloro-5,8-dihydroxy-1,4-naphthoquinone 1 and 1,4-dihydroxy-anthraquinone 2 [1,2]; see Figure 1. The dimers and trimers of the compounds were also extracted from the crystallographic data; see Figures 4 and 5 for the dimers and Figures S3 and S4 of the Supplementary Materials for the trimers, respectively. The geometry optimization was performed in vacuo and with the solvent reaction field reproduced by the Polarizable Continuum Model (PCM) [23] (in its IEF-PCM formulation) using water as a solvent. The B3LYP [28], PBE [29,30], and ωB97XD [31] functionals were applied to calculate metric and electronic structure parameters. The Pople-style triple zeta split-valence 6-311+G(2d,2p) orbital basis set [32][33][34] was employed in the static DFT. Next, the harmonic frequencies were computed to check if the obtained geometries corresponded to the minimum on the Potential Energy Surface (PES). The proton reaction path in the intramolecular hydrogen bond was investigated by means of the scan with the optimization method (the OHO valence angle was fixed; the O-H distance was elongated by a 0.05 Å increment; while the remaining part of the molecules was optimized) using the ωB97XD functional. The same functional was applied to prepare wavefunctions for the calculations of aromaticity index and Atoms In Molecules (AIM) theory [25]. The geometry-based Harmonic Oscillator Model of Aromaticity (HOMA) [24] was used to estimate the aromaticity changes upon the O-H distance elongation in the hydrogen bridge. The HOMA index was computed in vacuo and in water. The initial geometries of the dimers and trimers for further Symmetry-Adapted Perturbation Theory (SAPT) analysis [26] were calculated at the B3LYP/6-311+G(2d,2p) level of theory. This part of the quantum-chemical simulations was performed with the Gaussian 16, Rev. C.01 program [35]. The aromaticity index was calculated with the Multiwfn program [36].  (1,3) and experimental X-ray data (2,4,5) [1] of the dimers taken into account.

Atoms in Molecules Theory
Atoms In Molecules (AIM) [25] analysis was performed with the AIMAll package [37] in the gas phase. The analysis was carried out to investigate the electronic structure changes in the monomeric forms upon the intramolecular hydrogen bridge deformation by the O-H distance changes. The analysis was performed only for one hydrogen bridge due to the symmetry exhibited by the studied compounds. AIM theory served for the atomic charges' calculations, as well as to estimate the changes of the molecular properties at the Bond Critical Points (BCPs): the electron density and its Laplacian. It is well known that the AIM formalism should be applied at the equilibrium geometry, and the standard method nomenclature related to the Bond and Ring Critical Points (BCPs and RCPs) could be used. Here, we investigated non-equilibrium structures obtained as a result of the O-H distance changes in the studied molecules. Therefore, the electronic structure and topological features cannot be strictly related to the nomenclature and terminology developed in the classical AIM theory. However, we labeled the located BCPs following the convention used in the literature in the studies of non-equilibrium structures where the AIM theory was applied [38,39].

Symmetry-Adapted Perturbation Theory
The energy decomposition of the dimers and trimers of 1 and 2 (see Figures 4 and 5 and Figures S3 and S4) was carried out using Symmetry-Adapted Perturbation Theory (SAPT) [26]. The energy was decomposed for two sets of structures: (i) extracted from the experimental X-ray data [1,2] (ii) obtained as a result of quantum-chemical simulations at the B3LYP/6-311+G(2d,2p) level of theory. The SAPT analysis was done at the SAPT2 level of theory [40] for both sets of investigated dimers and trimers. The SAPT2 calculations for the DFT-optimized structures were performed with the aug-cc-pVDZ basis set [41] applied for the atomic orbital expansion. The interaction energy computation at the SAPT2/aug-cc-pVDZ level of theory involved also the counterpoise correction for the Basis Set Superposition Error (BSSE) [42], dividing the investigated dimers and trimers into "monomers". The SAPT calculations were performed using the Psi4 1.2.1 [43] program.

Time-Dependent Density Functional Theory
The time-dependent DFT calculations were carried out using the ORCA program (Version 4.2.1) [44,45]. All equilibrium and non-equilibrium structures were investigated in order to evaluate the differences in the excitation energies with respect to the geometry. We used the hybrid, three-parameter B3LYP functional [28] and range-separated hybrid functional, with dispersion corrections: ωB97XD [31]. To maintain consistency with the static, ground state calculations, the Pople-style triple zeta split-valence 6-311+G(2d,2p) orbital basis set [32][33][34] was used.
The solvent effects were accounted for using the Conductor-like Polarizable Continuum Model (CPCM) as implemented in ORCA, with the dielectric constant of = 80.4 and refractive index of 1.33.

Car-Parrinello Molecular Dynamics In Vacuo and the Crystalline Phase
Car-Parrinello molecular dynamics simulations were performed in vacuo and in the crystalline phase for Compounds 1 and 2 [1,2]. The models used for CPMD in the gas phase are presented in Figure 1. The initial geometries were extracted from the X-ray data and placed into cubic boxes with a = 14 Å for Compound 1 and a = 16 Å for Compound 2. The models for CPMD in the crystalline phase are presented in Figures 2 and 3, respectively. They were prepared on the basis of experimental monoclinic unit cells with a = 14.110 Å, b = 7.070 Å, c = 9.810 Å, and β = 102.09 • with Z = 4 for 1 and a = 20.339 Å, b = 6.017 Å, c = 10.391 Å, and β = 125.46 • with Z = 4 for 2. The exchange correlation functional by Perdew, Berke, and Ernzerhof (PBE) [29,30] and Troullier-Martins [46] pseudopotentials were applied. The fictitious Electron Mass (EMASS) was equal to 400 a.u. for both compounds and both phases' simulations. The time-step was 3 a.u., and the kinetic energy cutoff for the plane-wave basis set was 80 Ry. The CPMD calculations were performed at 295 K controlled by a Nosé-Hoover thermostat [47,48]. Hockney's scheme [49] was applied to remove interactions with periodic images of the cubic cell. The translational and rotational movements were removed from the CPMD data collection in the gas phase. The crystalline phase CPMD was performed with Γ point approximation [50] and Periodic Boundary Conditions (PBCs). The real-space electrostatic summations were set to TESR= 8 nearest neighbors in each direction. The initial part of the production run of the CPMD was taken as an equilibration (ca. 1.8 ps), and it was excluded during the post-processing procedures. The CPMD trajectories were collected further for 20 ps for both compounds in both phases. The data analysis focused mostly on the intramolecular hydrogen bonds dynamics. The CPMD simulations were carried out with the CPMD 3.17.1 program [51]. The graphics was prepared with the assistance of the VMD 1.9.3 [52] and Vesta (Version 3) [53] programs.

Results and Discussion
We focused on metric and electronic structure analyses of 2,3-dichloro-5,8-dihydroxy-1,4-naphthoquinone 1 [1] and 1,4-dihydroxy-anthraquinone 2 [2] in the electronic ground and excited states. We investigated the molecular features of the compounds using various theoretical approaches reproducing forces involved in the assembly of the compounds at the molecular level. We applied static DFT models, as well as first principles molecular dynamics using the Car-Parrinello formulation for the comprehensive theoretical description of molecular features in the studied naphtho-and anthra-quinone derivatives. The calculated geometric parameters of intramolecular hydrogen bonds of Compounds 1 and 2 are presented in Tables S1 and S2, respectively. The molecular structures of the compounds with the atoms numbering scheme are shown in Figure 1 and Figure S1. The literature X-ray data for the intramolecular hydrogen bonds of both compounds [1,2] are presented in Figure S2. Concerning Compound 1, the comparison of experimental and computed metric parameters showed that the O1...O2 interatomic distances were reproduced correctly by the B3LYP and ωB97XD functionals, while the PBE functional shortened the distance ca. 0.06 Å in both phases (results of the gas phase simulations and with the PCM). The O1 − H BP1 bond length was reproduced correctly as well by the three applied functionals in both discussed phases. The computed values of intramolecular hydrogen bond differed from the experimental data. The biggest discrepancies compared to the X-ray structure were obtained for the PBE functional, as shown in Table S1. The computed values of the O1 − H BP1 − O2 valence angle were larger compared to the experiment. The B3LYP and ωB97XD functionals were able to reproduce the geometric parameters of the intramolecular hydrogen bond more accurately than the PBE functional. The DFT simulations of the intramolecular hydrogen bond in Compound 2 provided us with a similar conclusion; see Table S2. The biggest discrepancies were noticed for the PBE functional compared to the X-ray data and the results obtained using the B3LYP and ωB97XD functionals. As we can see, the polar environment influence on the metric parameters is rather subtle concerning both compounds. This is valuable information while designing new naphtho-and anthra-quinone derivatives, similar to the discussed compounds, with specific molecular properties. The proton potential functions are presented in Figure 6. The proton reaction path was simulated in the gas phase and using the PCM (water as a solvent) with the assistance of the B3LYP and ωB97XD functionals. On the basis of the obtained results, we can draw the conclusion that the spontaneous proton-transfer phenomenon does not occur in Compound 1. However, the energy barriers to move the bridged proton to the proton-acceptor atom side are not extremely high. They are ca. 8.2 kcal/mol and ca. 10.5 kcal/mol according to the applied B3LYP and ωB97XD functionals and the basis set in the gas phase. The presence of the polar solvent slightly decreased the energy barriers up to 8 kcal/mol (B3LYP) and 10.2 kcal/mol (ωB97XD). As is shown in Figure 6, the second shallow energy minimum is present. Similar results were obtained for Compound 2. The energy necessary to move the bridged proton to the proton-acceptor side is not very high, and it is ca. 8 kcal/mol and 10.1 kcal/mol depending on the functional for the gas phase simulations. The application of the continuum solvation model and water to reproduce the polar environment influence on the molecular properties provided us with slightly lower energy barriers, which were equal 7.5 kcal/mol and 9.5 kcal/mol, respectively. The second energy minimum was detected as well. The changes of the aromaticity in both compounds were estimated based on HOMA index. The HOMA values were calculated as a function of the O1 − H BP1 distance changes in the gas phase and within the application of the PCM. The obtained results are summarized in Tables S3 and S4. In addition, they are presented in Figure S5 graphically to illustrate the observed dependencies. The structure in the equilibrium of Compound 1 consists of two rings, which are aromatic and anti-aromatic according to the HOMA index; see Table S3. While the bridged proton was moving to the proton-acceptor O2 atom, we noticed a consistent loss of aromaticity in the A ring; see Figures S1 and S2, respectively. Contrary results were obtained for the 1,4-benzoquinone ring (denoted as B). When the bridged proton was moving to the proton-acceptor O2 atom, the loss of anti-aromaticity was noticed. The aromaticity of the B ring increased, as shown in Figure S5. The introduction of the polar environment did not qualitatively affect the observed tendencies. This is in agreement with our geometric parameters' analysis, which showed that the metric parameters did not change remarkably upon the polar solvent. Compound 2 consists of three fused rings: two aromatic and one anti-aromatic. The loss of aromaticity was noticed for the A ring when the bridged proton was moving to the acceptor atom side (in our case, the O2 atom) along the hydrogen bridge. In the B ring, the increased of aromaticity was observed. The aromaticity in the C ring was slightly affected due to the internal geometric and electronic structure reorganization; see Figure S5. The analysis of the HOMA index values in the polar solution provided us with the same dependencies as for the gas phase results. The aromaticity of the rings was only slightly affected quantitatively, as shown in Table S4.

Electronic Structure and Topological Analysis of the Intramolecular Hydrogen Bonds Using Atoms in Molecules Theory
In Figure 7, the changes of the AIM net charges of atoms involved in the intramolecular hydrogen bond (O1 − H BP1 ...O2) are presented. The values of the net atomic charges computed at the ωB97XD/6-311+G(2d,2p) level of theory are listed in Tables S5 and S6, respectively. The elongation of the O1 − H BP1 distance influenced the electronic structure of the studied Compounds 1 and 2, as discussed below. In Compound 1, the value of the net charge of the proton-donor O1 atom was increasing up to 1.1298 Å of the O-H distance, then the value decreased again, but at 1.3298 Å, it was constantly increasing, while the bridged proton was moving to the proton-acceptor side. The value of the net charge of H BP1 was decreasing, while the O-H distance was elongated up to 1.1298 Å. Then, the net charge value was increasing, and finally, when the O-H distance was equal to 1.6298 Å, it was decreasing again. The value of the net charge of the proton-acceptor O2 atom was sharply decreasing, while the bridged proton was approaching the acceptor atom side. However, when the O-H distance was equal to 1.5798 Å, the increase of the net charge value was noticed. Other important AIM parameters of the covalent O1 − H BP1 bond, namely electron density and its Laplacian at the Bond Critical Points (BCPs), are provided in Figure S6. The electron density decreased, while the O-H distance elongated and the bridged proton was moving to the acceptor side. The electron density Laplacian values increased due to the O-H distance elongation. According to the AIM theory, the negative values of the electron density Laplacian indicate the covalent character of the interaction. The loss of the covalency was noticed while the bridged proton was moving closer to the protonacceptor atom. In our case, the changes of the electron density Laplacian values indicate that the intramolecular hydrogen bond has a partially covalent character depending on the bridged proton position. In the case of Compound 2, the net atomic charge value changes of the atom O1 as a function of the O-H distance elongation showed similar dependencies compared to Compound 1. First, the net atomic charge value increased upon the O-H distance elongation to 1.1312 Å. The decrease of the net charge value was noticed for the O-H distance between 1.1812 Å and 1.2812 Å. When the bridged proton was moving closer to the proton-acceptor atom, the increase in the net charge value was noticed again. The changes in the net charge values of the bridged proton H BP1 were as follows: the decrease of the values was noticed while the O-H distance was elongated up to 1.1812 Å. Then, the value increased, and when the O-H distance reached 1.5812 Å, it decreased again. Concerning the O2 proton-acceptor atom, the value of the net atomic charge decreased until the O-H distance was elongated up to 1.5812 Å. When the bridged proton was closer to the acceptor side, the value was increasing, as shown in Figure 7. The electron density and its Laplacian showed very similar dependencies compared to Compound 1; see Figure S6. The obtained results showed that for compounds with fused rings and strong, short intramolecular hydrogen bonds, we obtained very similar results. This kind of knowledge is valuable in the design of new molecules on the basis of naphthoand anthra-quinone core parts. The substituents do not significantly affect the properties of the intramolecular hydrogen bond. Moreover, we can conclude that the influence of substituents on the AIM properties of the atoms involved in the hydrogen bridge is rather subtle. Compound 1 contains two electronegative chlorine substituents, while Compound 2 contains an additional benzene C ring extending the zone of aromaticity, but the overall picture of the changes during the proton scan is the same. The role of the C ring of Compound 2 as the aromaticity buffer was noted already in the HOMA study, where this ring kept a constant, very high HOMA value even when the A and B rings underwent significant aromaticity modifications.
An interesting point is the location of the transition point when the bridge proton is crossing from the donor to the acceptor side. The curves of the net atomic charges of the donor and acceptor oxygen atoms cross at ca. 1.35 Å for Compound 1 and 1.3 Å for Compound 2. The corresponding crossings between the electron densities at the O1 − H BP1 and H BP1 − O2 bond BCPs, as well as their Laplacians, occur however at a lower O1 − H BP1 distance, ca. 1.24 Å for both compounds. This inconsistency shows that there are other factors influencing the properties of the hydrogen bridge. Obviously, the hydrogen bridge is a part of a quasi-ring, and the A and B rings strongly influence the substituents. This inequality in a case that could be expected to exhibit more symmetry shows that the levels of aromaticity in the A and B rings strongly interact with the donor and acceptor atoms (resonance effect). Additionally, the second hydrogen bridge O3 − H BP2 − O4 influences the studied bridge through the resonance.

The Intermolecular Interaction Energy Decomposition in the Dimers and Trimers on the Basis of Symmetry-Adapted Perturbation Theory Method
The dimers and trimers of Compounds 1 and 2 are presented in Figures 4 and 5 and in Figures S3 and S4. The SAPT calculations were performed for the optimized gas phase structures, as well as for the dimers and trimers extracted from the crystal structure. The analysis of interactions was carried out to determine the relative importance of the factors responsible for the crystal structure integrity.
The network of intermolecular hydrogen bonds is shown in Figure S2, prepared on the basis of the source article [2]. The summary of the SAPT energy partitioning is shown in Table 1. The comparison of the two levels of theory, the very fast SAPT0 and the more demanding SAPT2, shows a constant tendency of SAPT0 to bind the monomers stronger than SAPT2, with the only exception of the stacked dimer 5 of Compound 1, where the two SAPT levels are virtually identical. This structure is exceptional among the studied ones because of its domination by dispersion. The stacked monomers are aligned in an anti-parallel (centrosymmetric) arrangement; thus, their dipole moments are also aligned in the most advantageous manner. Indeed, among the electrostatic terms for Compound 1, the largest in magnitude is the one for the stacked dimer (−5.30 kcal/mol). The ratio of electrostatics and exchange terms is close to one in most cases, but again, in the stacked dimer, the strong dispersion contribution forces the monomers to become very close, while the Pauli exchange repulsion is exceptionally large.
The relatively low importance of the induction term indicates another feature of the studied naphtho-and anthra-quinone derivatives: they are molecules of low polarity (the symmetry of the -OH substituents' arrangement is one of the factors responsible), and their polarizability (response to an electrostatic stimulus) cannot be used to the advantage of the dimer.
Finally, let us note that the intermolecular forces corresponding to the hydrogen bonds are rather weak, especially in comparison with the intramolecular contact. The external interactions compete rather unsuccessfully with the intramolecular resonanceassisted hydrogen bonds. Thus, the driving force for the arrangement of the molecules in the crystal is dispersion, but the electrostatics can determine the orientation of the dimer components. When the dimer is planar and cannot be stabilized by dispersion as strongly as the stacked dimer, the molecules reorient from the X-ray-determined positions to maximize the electrostatic attraction. This is true for Compound 1, but the planar dimer of Compound 2 minimizes the exchange contribution instead of maximizing the electrostatics term. All in all, the data for the dimer of Compound 2 are not very different from the dimers of Compound 1, even if the total interaction energy is lower (i.e., better) by ca. 1 kcal/mol. This lowering can be explained by stronger dispersion forces due to an additional benzene ring. Thus, the bifurcated hydrogen bond in the dimer of Compound 2, postulated on the basis of the experimental crystal structure [2] and mentioned in the Introduction, cannot be strong in its intermolecular part.
The SAPT partitioning of the interaction energy for the trimers is summarized in Table 2. Each trimer was divided into three dimers, AB, AC, and BC; these fragments were then treated separately. The trimer selected for Compound 1 was composed of three interactions of similar magnitude. Note that the optimized structures displayed an interesting relation between the interaction energy terms: electrostatics and exchange energies canceled out almost completely. This was already seen for the dimers; see above. The induction was small, and the stability of the trimers relied on dispersion. The structures taken directly from the crystal coordinates behaved in a similar way; however, the presence of an extended interaction network in the crystal forced the molecules closer, and the repulsive exchange term was not balanced by the electrostatics and induction. The triangular arrangement of the investigated trimer of Compound 1 was very different from that of Compound 2. In this case (see Figure S4), an infinite ladder of molecules was formed, and the molecules were aligned in an anti-parallel way along the ladder. The formation of such a ladder was facilitated, but not driven, by the interplay between intra-and inter-molecular hydrogen bonds, of which the intermolecular bridges were significantly weaker, as already discussed above. Please note also that the interaction of the two molecules set farthest apart (the dimers AC for Compound 2) was virtually non-existent. The molecules were too far away to interact strongly by dispersion forces, and their parallel arrangement made the electrostatics and induction terms equally negligible.

An Application of Time-Dependent Density Functional Theory
The TD-DFT results were in agreement with our intuitive expectations. The intramolecular proton transfer did not strongly influence the distribution of the electron density within the system, and as such, it did not lead to the significant differences in the excitation energies. The calculated values for the non-equilibrium geometries along the proton transfer are shown in Figure 8.
All investigated systems were characterized by the ability to absorb a photon within the UV-Vis range between 2.7 eV and 3.3 eV, which corresponds to the wavelengths of 460 nm and 375 nm, respectively. In all cases, the first excited state was determined to be the bright state, with an oscillator strength of approximately 0.2-0.4. The calculated wavelengths were blue-shifted with respect to the available experimental data on Compound 1: its UV-Vis absorption bands in chloroform solution were located at 570, 528, 496, and 468 nm [1]. A similar problem was noted for the parent compound, naphthazarin, which exhibited vibronic couplings in its first excited state (absorption maxima at 564, 547, 524, and 490 nm [54]). On the other hand, theoretical predictions placed the absorption of naphthazarin at 436 nm (CAM-B3LYP range-separated hybrid) and 490 nm (PBE0 hybrid functional). Such a discrepancy might be considered small, but it is believed that TD-DFT should be rather accurate in the case of organic dyes [55]. This issue is further discussed at the end of this section. At this point, it is important to stress that the results were similar for both Compounds 1 and 2, indicating that the substituent effects were rather limited.
Interestingly, taking into account the solvent effect did not bring significant differences in the excitation energies, and in all investigated cases, the curves for the gas phase and continuous solvent calculations overlapped for the most part. The biggest differences were observed for the optimized (equilibrium) geometries, with the presence of the water slightly decreasing the activation energies; however, the differences were still minor.
As far as the role of the functional used is concerned, the values obtained from the three-parameter B3LYP functional were smaller by approximately 0.02 a.u. (0.5 eV) compared to ωB97XD. A similar situation was found for naphthazarin, where the PBE0 hybrid gave better absolute agreement with the experiment, but the CAM-B3LYP functional was able to reproduce vibronic coupling pattern with the accuracy much improved over PBE0 [55]. In the case of Compound 1, the experimental bands were located at 570, 528, 496, and 468 nm, and the lowest energy value (570 nm) was still far away from the B3LYP estimate of ca. 460 nm. The explanation could be given in a manner similar to the case of naphthazarin: the symmetric double C=O stretching mode modulated the zwitterionic character of the enedione moiety (the quinoid ring), strongly affecting the HOMO-LUMO transition [55]. The visualization of the HOMO/LUMO orbitals (see Figure S7 in the SI) confirmed the π/π * character of these orbitals.
In general, despite the discrepancies between the exact values of the computational and experimental results, the trends still hold. There was a negligible effect of the substituent on the photoactive properties of the investigated compounds. In addition, the proton transfer did not significantly influence the electron configuration, and therefore, its role in the photoactivity was also negligible.

Car-Parrinello Molecular Dynamics Results in In Vacuo and in the Crystalline Phase-the Analysis of the Intra-and Intermolecular Hydrogen Bonds
Following our gas phase and solvent reaction field findings based on static DFT models, we applied Car-Parrinello molecular dynamics (CPMD) to study the proton dynamics in both intramolecular hydrogen bridges (O1 − H BP1 ...O2 and O3 − H BP2 ...O4) in Compounds 1 and 2. The time-evolution of the metric parameters of atoms involved in the intramolecular hydrogen bonds in Compound 1 is presented in Figure 9. The q parameter is defined as the difference between the O-H and the hydrogen bond (H...O) distance in our case. When the q-parameter value oscillates around 0.5, the bridged proton is localized on the donor side. The q-parameter values around zero indicate that the proton is in the middle of the hydrogen bridge. Finally, the q-parameter value around −0.5 means that we observed proton transfer phenomena. The CPMD simulations of the isolated molecule showed that in both hydrogen bridges, the proton was labile, but localized on the donor side. However, as is shown in Figure 9a, there were very frequent events of the proton-sharing revealed by the CPMD simulation. The average interatomic distance of O...O in both hydrogen bridges was equal to 2.549 Å (due to the symmetry exhibited by the molecule). In the crystalline phase, the bridged proton dynamics was affected by several factors, e.g., crystal field and the presence of the neighboring molecules. The latter one provided additional intermolecular interactions, e.g., stacking type weak interactions induced by the arrangement of the aromatic rings in Compound 1. In Figure 9b, the results of bridged protons dynamics are presented. Similar to the gas phase results, the bridged protons were located on the donor side. However, the proton-sharing events and short proton transfer phenomena were noticed during the collection of the CPMD trajectory. The proton moved to the acceptor side and immediately returned to the donor side. This showed even larger proton mobility compared to the isolated molecule results. The average O...O interatomic distance was equal 2.543 Å, and it was shortened slightly by external forces associated with the crystalline phase. The hydrogen bridges' correlation was examined as well. It was found that, similar to the gas phase, the cis arrangement of the bridged protons was preferable. Compound 2 possessed strong, Resonance-Assisted Hydrogen Bonds (RAHBs). The CPMD results obtained for the isolated molecule indicated that the bridged proton was located mostly on the donor side. However, as is shown in Figure 10, the proton exhibited strong mobility. The proton-sharing events were noticed during the simulation time, as well as proton transfer phenomena were observed (very rare events). The bridged proton moved to the acceptor side for a very short period of time and returned again to the donor side. The average value of the O...O distance was 2.542 Å, and it was shorter compared to Compound 1. A very similar proton dynamics was obtained for the crystalline phase simulations. This could be related to the competition between the intra-and inter-molecular interactions. For the isolated molecule, only intramolecular interactions could be analyzed; therefore, it is a good reference to discuss the intermolecular interactions in the crystalline phase. Figure 10b presents results describing the intramolecular hydrogen bridges ' dynamics. The bridged protons were located on the donor side, but they were strongly labile. The proton-sharing events were noticed, and even proton transfer phenomena occurred, but rare and very short. The bridged protons always returned to the protondonor side. Interestingly, the average O...O interatomic distance was equal to 2.550 Å, and it was longer in the crystalline phase comparing to the gas phase results. This result could be explain by the arrangements of the molecules in the crystal unit cell. According to the experimental X-ray results [2] the 1,4-dihydroxy-anthraquinone molecules were connected by remarkably short contacts between oxygen atoms; see Figure S4. Therefore, in Compound 2, the external forces affected the intramolecular hydrogen bridge metric parameters, but did not significantly change the proton dynamics. This is in agreement with the SAPT discussion showing that the intermolecular contacts were much weaker than the intramolecular bridges.
The last, but very important point of the structural time evolution study is the question of the correlation between the motions of the protons in the twin hydrogen bridges. The investigated molecules exhibited C 2v symmetry in their minima on the potential energy surface. The aromaticity analysis in the structural analysis section revealed that the quinoid ring did not become significantly more aromatic when a single proton transfer occurred; see Tables S3 and S4 and Figure S5 in the Supplementary Information. Even worse, when such a transfer occurred, there was also a large loss of aromaticity in the benzenoid ring. We can then expect that the proton-sharing events should be correlated between the twin bridges. This is indeed the case: most of the rare, short proton-sharing or transfer events of one proton were accompanied by a similar action of the other bridged proton. Figure 9 shows such events for Compound 1 in the gas phase at 1.7/1.9 ps (the convention we will use, t1/t2, is: t1-time for the first event at one of the protons, t2-the accompanying event for the other bridge proton), 9.1/9.6 ps, and 15.5/15.9 ps. There were however many cases where a short proton-sharing or transfer event of one proton was not accompanied by a respective motion in the other bridge, such as at 1.0 ps in the gas phase simulation of Compound 1. The solid state trajectory for this compound was more interesting, because it contained the sole observed instance of a longer sustained double-proton transfer lasting from 3.75 to 4.4 ps. Figure S8 in the Supplementary Materials shows the details of this event. The double-proton transfer was not ideally synchronous, but the response of one of the protons was delayed by just two-three periods of the bridge O-H stretching. The event happened because fortunately just before 3.75 ps of the simulation time, the respective donor-acceptor oxygen atom distances were very low, decreasing the instantaneous barrier. The protons after the transfer underwent motions of large amplitude in the potential well at the acceptor side, and one of them briefly returned to the donor side at 4.1 ps. When the amplitudes of the donor-acceptor motions rose again, the back-transfer to the donor side took place at 4.34-4.39 ps, and this was a synchronous motion in one continuous slide along the potential energy path. The remaining correlated proton-sharing events for the crystalline phase of Compound 1 (at 7.8/8.1 ps, 9.1/9.2 ps) were again of the rare and rapid type. The data for Compound 2 differed slightly from those reported above: the proton-sharing events were more frequent, but less correlated (more independent). Indeed, only the events at 18.1/18.7 ps in the gas phase and at 10.4/10.6 ps in the crystal can be considered as correlated. No sustained periods of proton transfer events were recorded for this compound.

Conclusions
The current study was designed to characterize the intra-and inter-molecular interactions in two selected derivatives of naphtho-and anthra-quinone. The theoretical investigations were performed in three phases: in vacuo, in polar solvent reproduced by the application of the PCM, and in the crystalline phase. The molecular properties of the two compounds were studied in the electronic ground and excited states. Concerning the intramolecular interactions, a detailed description of the hydrogen bond dynamics was provided on the basis of static DFT and CPMD models. The intermolecular forces were explored in dimers and trimers formed by the discussed compounds. Our theoretical findings were compared with experimental structural data [1,2]. Summarizing, the static DFT models provided us with a detailed description of the bridged proton reaction path. The second energy minimum was detected in both compounds; moreover, the energy barriers for the proton transfer to the acceptor side were found not to be very high. The introduction of the polar solvent did not significantly change the metric and electronic structure parameters. We observed that the polar solvent moderately decreased the energy barrier for the proton movement. The application of the HOMA index showed that the bridged proton position had an impact on the aromaticity changes in the rings. The decrease of aromaticity was noticed for the hydroxyquinone moiety ring, while its increase was observed in the 1,4-benzoquinone ring. The electronic structure and topological analysis indicated a partially covalent character of the intramolecular hydrogen bond showing a strong property dependence on the bridged proton position. The SAPT energy partitioning showed the decisive role of dispersion for the intermolecular attraction. The bifurcated hydrogen bonds in the dimers of Compound 2, postulated experimentally [2], were much weaker than the intramolecular contact.
TD-DFT calculations confirmed the assumption of a negligible effect of the proton transfer on the electronic excitation energies in the investigated compounds. As the proton transfer did not involve the change in the electron density and the distance covered during the transfer was too small to induce a change in the dipole moment of the molecule, it cannot lead to a significant change in the electron configuration of the systems. This was observed as flat excitation energy curves. The experimental excitation energies did not agree well with the calculated ones; however, this was most likely related to the incorrect description of the electron configuration of the parent compound, naphthazarin. However, the qualitative description of small influence of the substituents on the excitation energies was in good agreement with the experimental results.
The CPMD simulations revealed that in both compounds, the bridged protons were very labile, but localized on the donor side. However, proton-sharing events and rare, short proton transfer phenomena were noticed. Many of them were correlated between the bridges of Compound 1, much less so in the case of Compound 2. The only longer proton transfer period for both bridges in Compound 1 began with an almost synchronous movement of the two bridge protons: the delay of the response of the second proton was only a few O-H oscillations. The photochemical structure-property relationship makes the transmission of the input "signal" (a photonic-driven proton transfer at one of the bridges) to the output "gate" (the proton transfer in the other bridge will shortly follow) possible. This interesting feature along with the excited state properties of the studied compounds make Compound 1 a good candidate for further studies on molecular electronics building blocks.