Experimental and theoretical characterization of chelidonic acid structure

Chelidonic acid (4-oxo-4H-pyran-2,6-dicarboxylic acid) is present in plants of Papaveraceae family, especially in Chelidonium majus. Due to its anticancer, antibacterial, hepatoprotective, and antioxidant properties, it has been used in medical treatments. In this work, the X-ray structure of methanol solvate of chelidonic acid was determined. Layers of chelidonic acid are held by hydrogen bonds via COOH and C = O fragments and additionally bridged by methanol. The formed H-bond network between two acid units is different from typical –COOH dimers observed, e.g., in crystals of isophtalic acid. The molecular structure of 2,6-dimethyl-γ-pyrone (2Me4PN) and chelidonic acid, a 2,6-dicarboxylic derivate of γ-pyrone (4PN), was verified in silico using density functional theory (DFT-B3LYP) combined with large correlation-consistent basis sets. The impact of –CH3 and –COOH substituents on 4PN ring structure, dipole moments, geometric/magnetic indexes of aromaticity, and NBO charges was assessed following unconstrained geometry optimization in the gas phase, chloroform, methanol, DMSO, and water with solvent effect introduced using the polarized continuous model (PCM). H-bond network formed in chelidonic acid–methanol complex was analyzed and their interaction energy estimated. Theoretical modeling enabled prediction of accurate structural parameters, dipole moments, and geometric/magnetic indexes of aromaticity of the studied 4PN, 2Me4PN, and chelidonic acid molecules.


Introduction
Chelidonic acid is a biogenic dicarboxylic acid, which has a γ-pyrone core [1], and is present as a secondary metabolite in different Papaveraceae family plants (e.g., Cassia minosodia, Sorghum vulgare, Cassia spectabilis), but mainly in Chelidonium majus L. [2,3]. For brevity, in the Supplementary Information, we mentioned some general information on plants containing chelidonic acid and other bioactive metabolites.
Recently, metabolites from C. majus have been studied using new techniques, such as NMR [4][5][6][7][8]. The current study aims at extending the basic knowledge about chelidonic acid, which is an interesting bioactive substance from C. majus.
Despite numerous bioactive compounds, present in C. majus L., their structure and spectroscopic properties are generally not known in details. On the other hand, many excellent NMR studies on these compounds were reported by Marek Radek group [4][5][6][7][8].
Unfortunately, due to experimental obstacles, the crystal structure of chelidonic acid or its solvates has been unknown. As a result of many attempts, we succeeded to raise the unstable crystal of methanol solvate of chelidonic acid and determined its room-temperature structure.
In the current work, we begin with a full X-ray characterization of methanol solvate of chelidonic acid, determined in this study, and compare it with known structures of isophtalic acid [9] and methyl monoester of chelidonic acid [10]. The schematic formula of chelidonic acid and its atom numbering, as given by X-ray studies, is shown in Fig. 1A. The molecular structure of chelidonic acid shows C 2 symmetry and apparently looks very simple. To learn about population of different H-bonding motifs in crystals of carboxylic acid dimers, a search in the Cambridge Structure Database (CSD [11]) was performed.
Some old theoretical studies on chelidonic acid, without recovering electron correlation energy, reported [12] on its geometrical structure, as predicted by the Hartree-Fock method, combined with very small Pople type basis sets (RHF/6-31G*). Later, DFT calculations at B3LYP/6-31+ G* level of theory [13,14] provided more accurate results. Ab initio studies on 4PN (Fig. 1B), which is a parent molecule for chelidonic acid, were also reported [15,16]. Recently, the molecular structure and vibrational spectra of chelidonic acid monomer and dimer were studied theoretically at B3LYP/6-311+ + G** level of theory by Malaganvi et al. [17]. The authors also reported on experimental FT-IR and Raman spectra of solid chelidonic acid.
In this study, the title molecule has been studied together with two closely related molecules, 4PN and 2Me4PN, selected as references ( Fig. 1B and C).
We have paid a special attention to very accurate theoretical prediction of all structural parameters in various environments. The calculated structures include precise hydrogen atoms localization which hopefully could be compared with future experiment [18], as well as the assessment of dipole moments and distribution of charges in chelidonic acid and its two derivatives with H and CH 3 substituents at ring positions C2 and C6.
In general, theoretical results predicted with more complete and flexible basis sets provide better reproduction of experimental parameters [19][20][21]. Thus, we have compared the experimentally determined structure of methanol solvate of chelidonic acid with the structures of two closely related molecules, calculated with very large basis sets in the gas phase, chloroform, methanol, DMSO, and water solutions, modeled with the polarizable continuum model (PCM) using the integral equation formalism variant [22] (IEF-PCM). Comparison between theoretically predicted gas-phase geometries of molecules with data from solid samples is commonly used to study crystal structures of biologically active molecules [23][24][25]. The geometrical changes between the reported gas phase structure of 4PN, and available crystal structures of 2Me4PN and methanol solvate of chelidonic acid, determined in this work, have been analyzed. Hopefully, this will allow an observation of the substituent effect on the geometric properties of the parent molecule (4PN) which are not obscured by any intermolecular interactions, including solvent or crystal forces [19]. Special attention has been paid to C = O and C = C bonds variations. Thus, the heterocyclic rings of 4PN and 2Me4PN molecules served as a reference for geometry changes upon H/CH 3 and H/COOH replacements at C2 and C6 positions.
Moreover, in this study, an advanced complete basis set [26][27][28][29][30][31][32][33] (CBS) estimated structural parameters have been directly compared with the experimental geometry of the smallest investigated molecule in the gas phase, e.g., with 4PN and next with X-ray structures of its dimethyl-and dicarboxylate derivatives. Theoretical gauge-including atomic orbital [34] (GIAO)-calculated isotropic nuclear magnetic shieldings have been used to predict the magnetic indexes of aromaticity, NICS (nuclear independent chemical shift) [35]. In addition, geometric index of aromaticity (harmonic oscillator model of aromaticity, HOMA [36][37][38]) was derived and the degree of electron delocalization within the ring discussed. Finally, chelidonic acid associates with methanol were modeled in order to get some insight into the H-bonding patterns in the crystal structure and binding energetic.

Experimental
Chelidonic acid was purchased from Aldrich. For brevity, additional information on chelidonic acid purification and synthesis was included in the Supplementary Material (SI). Single crystal X-ray measurement of chelidonic acid mono crystal, obtained as a methanol solvate, was conducted at room temperature. Some important information on X-ray measurement and structure determination [39,40] has been gathered in Table S1 in SI. Crystal structures were visualized using Mercury program [41].

Density functional theory (DFT) structure prediction
Accurate DFT [42] unconstrained geometry optimization with efficient B3LYP hybrid density functional [13,14,16], combined with Dunning's type aug-cc-pVXZ basis set series [43,44], where X = D, T, Q, and 5, has been conducted with Gaussian 16 program package [45]. In this study, we have selected the B3LYP hybrid density functional because of its fairly good performance in predicting structure, energy, and spectroscopic parameters of organic systems. Very tight optimization criteria and large grid size (150,590) were used. All positive harmonic vibrational frequencies supported prediction of true potential energy minima structures. Subsequent estimation of selected structural parameters in the complete basis set limit [26][27][28][29][30][31][32][33] using the 2-parameter fit [26,30] to individual results was conducted. It also avoids acceptance of inaccurate results due to fortuitous error cancelation, often obscuring results of a single calculation with smaller basis sets. The use of systematically growing basis sets helps in noticing trends in changes of structural (and other) parameters upon increasing the cardinal number X. Moreover, to clearly notice changes in trends, we present our results using three and often four decimal points. One of the simple ways, leading to estimation of CBS parameters [26][27][28][29][30][31][32][33], for example the electronic energy (E), is to perform several calculations of the selected parameter using specially designed basis sets, for example, correlation-consistent ones of Dunning type [43,44], and next fit the obtained individual energies with 2-parameter function [26,30]: where E(CBS) and B are fitting parameters and X is the cardinal number of the basis set. Apart from energy, other atomic and molecular parameters, related to energy, could be derived in the CBS limit [26][27][28][29][30][31][32][33].
The selected molecules of 4PN, 2Me4PN, and chelidonic acid were initially optimized in vacuum, and using the PCM approach, in chloroform, methanol, DMSO, and water. In addition, for verification of calculation accuracy, an unconstrained optimization of 4PN in vacuum was performed at the coupled cluster with single and double excitations (CCSD/cc-pVTZ) level of theory [46] in Gaussian program [45]. Due to its huge computational cost only one coupled cluster optimization of structure, combined with relatively small size basis set, has been conducted.
Besides, a fragment of methanol solvate of chelidonic acid crystal structure, containing a dimer formed by chelidonic acid and methanol, was optimized in vacuum at B3LYPD3/ aug-cc-pVTZ level of theory. In this case, the Grimme GD3 dispersion correction [47] was used and interaction energy, corrected for basis set superposition error using the counterpoise method [47] (CP), was determined. This allowed an analysis of the formed H-bond network and reliable assessment of the interaction energy. To estimate the population of carboxylic acid dimers forming different motifs, a search in Cambridge Structure Database [11] was performed.
Individual charges on atoms have been predicted with the aid of natural population analysis (NBO) at the B3LYP/augcc-pV5Z level of theory [48]. Finally, magnetic and geometric indexes of aromaticity, NICS, and HOMA were calculated for the studied molecules at B3LYP/aug-cc-pVTZ level of theory. For an n-membered ring, HOMA was calculated as follows [36][37][38]: where R i and R opt are predicted and reference bond lengths, and α is normalization constant selected for a particular bond type (α is 257.7, 157.38, and 93.52 for CC, CO, and CN bonds). The reported reference CC, CO, and CN bond lengths are 1.388, 1.265, and 1.334 Å [36]. For these calculations, benzene, pyridine, pyrrole, and furan were additionally selected as reference molecules and their HOMA and NICS values calculated at the same level of theory and with geometry optimized using the same methodology.

Crystal structure of chelidonic acid
Due to the stability and crystallization problems, there have been no available X-ray structures of chelidonic acid reported in the literature (see the CSD basis [11]). Instead, some of its derivatives, including its mono methyl ester, were reported by other authors [10]. However, as a result of extensive attempts in the current study, we have managed to obtain a monocrystal of chelidonic acid with one solvent molecule (methanol) suitable for X-ray examination. Thus, the structure of methanol solvate of chelidonic acid will be compared with our theoretical predictions in the gas phase and solution in the subsequent parts of the work. Obviously, the gas phase measured interatomic distances (and angles) between carbon, oxygen, or nitrogen atoms should be close to those, determined by the solid-state X-ray technique [49,50]. In addition, the electron diffraction (ED) and microwave studies (MW) are also able to accurately measure C-H, N-H, and O-H bond lengths where traditional X-ray technique fails (see also ref. [18] for an enormous progress in localization of hydrogen atoms in low temperature X-ray experiments). Unfortunately, ED and MW techniques are only applicable to relatively small molecules in the gas phase [49]. For brevity, the details of X-ray measurements of methanol solvate of chelidonic acid in the current study are included in Table S1 in the Supplementary Material. Figure 2 shows the X-ray-determined structure of a single molecule of chelidonic acid together with the closest molecule of methanol, and with atom labeling.
The room-temperature crystal structure of chelidonic acid solvated with methanol, determined in the current work, is deposited in CCDC crystal basis under the number of 2,144,108 (see also some X-ray-determined interatomic distances in Table S2 in the Supplementary Material). The corresponding packing diagram of the studied acid and the accompanying methanol molecules, indicating potential interatomic interactions, is shown in Fig. 3.
One of -COOH groups in chelidonic acid is rotated by 180° to better accommodate the solvent molecule. As result, the molecular structure of chelidonic acid in the crystal state loses its C2 symmetry. Methanol molecules play an essential role in the crystal network of chelidonic acid and act as a kind of bridging units between two layers of acid (Figs. 3 and 4). It is also apparent that each two chains of acid molecules has no direct contacts via H-bonds. Instead, molecules of methanol play this indirect role. The molecules of acid in the chain (or layer) are kept by the H-bond formed between one -COOH group and the carbonyl group from the neighboring chelidonic acid molecule, which is located at C4 position in the ring. The other carboxyl group binds two solvent molecules. The methanol molecules also form H-bonds with a second layer of chelidonic acid what results in a two-layer structure (Fig. 3).
As shown in Fig. 4, hydrogen bonds between molecules of chelidonic acid and methanol control its two-layered arrangement in the crystal structure. Of particular interest are two contacts: a non-typical H-bond between H atom in methanol and acid O4 atom (1.983 Å, sum of vdW distance = 2.72 Å and a typical H-bond between acid H3 and methanol O atom being 1.691 Å). These bonds are very important in formation of two-layered crystal structure of chelidonic acid since they connect both layers via methanol molecules.

Crystal structure of monomethyl ester of chelidonic acid
Monomethyl ester of chelidonic acid [10] also forms a similar crystal pattern to monomethanol solvate of chelidonic acidthere are H-bonds between the -COOH group of one molecule and carbonyl group at C4 position from the neighboring molecule (Fig. 5). Packing of methyl ester of chelidonic acid [10] is shown in Fig. S1 in the Supplementary Material. Typical -COOH bonds in dimers of carboxylic acids, for example present in formic acid dimer [51,52], are also not visible in case of mono methyl ester of chelidonic acid.
It is important to notice that Malaganvi et al. [17] performed theoretical calculations and experimental IR/ Raman studies in KBr, concluding that chelidonic acid dimer forms a typical structure with two carboxylic groups held by strong (and short) intermolecular hydrogen bonds what resembles formic acid (or isophthalic acid) dimers. However, their proposed dimeric structure is not confirmed by our X-ray structure of crystalline methanol solvate of chelidonic acid or methyl ester of chelidonic acid [10].

Crystal structure of isophthalic acid
It has been interesting to compare the determined in this work H-bonding network in the crystal of methanol mono solvate of chelidonic acid, as well as, in its earlier reported monoester with another organic molecule, bearing two symmetrically arranged carboxylic groups at meta-positions of benzene ring. The monocrystal of isophtalic acid is a good example of such a system from the literature [9,53]. The crystal of isophtalic acid, measured at 90 K, where molecules form chains with typical two carboxylic groups held together, differs essentially from crystalline methanol solvate of chelidonic acid, measured in this study, and is similar to dimers of formic acid [51,52,54]. In this case, the H-bonds also allow a two-layered arrangement of acid molecules (Fig. 6) and leads to their corresponding crystal packing (Fig. S2).
Thus, different crystal structures of both acids with two COOH groups or with one blocked carboxylic group due to ester formation are clearly seen in Figs. 4, 5, and 6.

A short analysis of population of observed motifs of carboxylic acids H-bond network in CSD
The above discussed three clearly distinguished motifs of H-bonding in carboxylic acid crystals, similar to those, shown in Figs. 4, 5, and 6, appear in CSD [11] and are shown schematically in Fig. 7.
It seems that CSD could shed some light on the abundance of different H-bond motifs holding together crystalline carboxylic acids. However, one should keep in mind some reservation about subsequent accuracy of reliable hydrogen atoms localization [18], mainly based on hybridization considerations, and subsequent recalculation of O-H bond length using ab initio/DFT methodology. Nevertheless, it us useful to explore the available CSD structures according to three selected motifs, shown in Fig. 7. Obviously, the observed H…O distances give an important information about the H-bond strength. H-bonds between isophtalic acid molecules in a two-layer crystal pattern [9] It was interesting to check the population of crystal structure motifs, reported in CSD and forming typical carboxylic dimers presented in Fig. 7A and other kind of dimers between a molecule bearing the -COOH group and the neighboring one, with C = O substituent (Fig. 7B). The third set of H-bond network is formed by two acid molecules with their carboxylate groups separated by two methanol molecules (Fig. 7C).
It is apparent from CSD data [11] that typical dimeric structures (A) form a very large number of crystals (7497) with an average O…H distance of 1.775 Å. Additionally, it is expected that due to crystal packing and intermolecular interactions (and within the experimental error), this set could be divided into symmetric (5520 structures for d1 = d2) and asymmetric dimers (1977 crystals, or 26.37%, including 977 cases with d1 < d2 and 1000 crystals with d1 > d2). However, the differences between d1 and d2 distances are very small but still could reflect some features of a real picture.
A significantly smaller number of dimers (469), formed by interacting -COOH and O = C functional groups with O-H…O = C distances ranging between 1.5 and 2.7 Å, are also gathered in CSD. The corresponding distance in methanol solvate of chelidonic acid is 1.837 Å. Moreover, molecules with carboxylic groups connected (bridged) with each other by two methanol molecules are very rare in CSD. There are only 35 crystals with this type of very strong H-bond network (d1-d4 from 1.620 to 1.815 Å). Out of them, four (about 10.8% of this crystal set), with refcodes BANMUZ, CILLUE, GENMES, and QIQYEW, are held by very strong bonds (about 1.5 Å). However, this unique pattern is also observed in the methanol solvate of chelidonic acid, analyzed in the current study with two different H-bonds (1.685 and 2.142 Å), and indeed, this is very rare case.
Thus, from the total of 8001 dimeric structures in CSD [11], 93.70% are typical carboxylic dimers, 5.83% are B type structures, and only 0.44% are carboxylic dimers bridged by methanol molecules (C structures). Comparing the H-bond lengths, it is evident that the methanol solvate of chelidonic acid (Fig. 4) forms less stable acid -methanol cluster and the corresponding bonds are longer by about 0.07-0.33 Å.

Theoretical structures of 4PN, 2Me4PN, and chelidonic acid
As expected, the three calculated molecules in the gas phase and solution are planar. In the following parts of our study, the predicted interatomic distances and bond angles in these molecules are analyzed and compared with available experimental data in the gas phase or in the crystalline state. Thus, our approach allows an observation of the influence of substituents on structural variations between the studied molecules without any additional perturbations. Thus, it is important to stress that such studies are free from environmental effects. On the other hand, we also have included the implicit solvent effect on their structures. As expected, the C-C and C-O single bonds in crystals closely resemble those in the gas phase while C = O ones in the crystal are elongated due to H-bonding [49].

4PN
For brevity, all the B3LYP/aug-cc-pVXZ calculated bond lengths and interatomic angles for an isolated molecule of 4PN are gathered in Tables S3 and S4 in the SI. It is apparent that the bond distances only very slightly change with the improvement of basis set completeness and size. Thus, the bonds calculated with aug-cc-pVDZ and augcc-pV5Z basis sets differ from experiment by 0.005 to 0.010 Å, i.e., are within experimental error bars. A very good performance of B3LYP density functional, expressed by the complete basis set limit estimated bond lengths, in prediction of 4PN structure is apparent from a very close reproduction of experimental parameters in the gas phase, as revealed by microwave studies [50] (RMS = 0.003 Å, see Table 1).
Additionally, in order to verify the quality of our DFT predicted bond lengths of 4PN, we have performed optimization of the free molecule using the state-of-the-art coupled cluster with single and double excitations [46] level of theory (CCSD). The CCSD/cc-pVTZ-optimized structure was also very accurate; however, the RMS deviation of bond lengths from experiment (0.005 Å) was slightly worse than the B3LYP/CBS result (0.003 Å). Obviously, further optimization of 4PN using CCSD method combined with better quality basis set was impractical.
There is no significant impact of solvent effect inclusion within PCM model on the optimized bond lengths of 4PN (see Table 1 with RMS raising only from 0.003 in vacuum to 0.006 Å in water). As expected, only a slightly better Fig. 7 Three H-bond network motifs observed in carboxylic acid dimers [11] agreement between theory and experiment in the gas phase is observed. Thus, moving from an isolated 4PN molecule to water, no change of the C-H bond is observed and the most important effect is related to elongation of C4 = O2 bond by 0.013 Å. At the same time, the adjacent C3-C4 bond decreases by 0.009 Å and the C2 = C3 double bond elongates by 0.003 Å. In addition, the C2-O1 bond shortens by 0.004 Å. These small but not negligible changes indicate a redistribution of electron density between alternating double and single bonds. As a result, a more efficient delocalization between the two terminal oxygen atoms could exist.
The B3LYP/CBS bond lengths, estimated in the complete basis set limit, agree with experiment in the gas phase nearly ideally (RMS = 0.003 Å) but the angles are slightly less accurate (Table S5, RMS = 0.314°). The angles are also more sensitive to the quality of the basis set (the corresponding differences between angles obtained with cardinal numbers X = D and X = 5 are from 0.016 to 0.150°). In addition, the CCSD/cc-pVTZ-optimized structure of 4PN also predicts marginally worse bond angles than the B3LYP/CBS calculations (RMS of 0.318 vs. 0.314°). Solvent polarity also influences bond angles in 4PN more than the bond lengths (0.3 to 0.6° difference between water and the gas phase is observed [50]).
2Me4PN The B3LYP/aug-cc-pVXZ calculated bond lengths and interatomic angles for an isolated molecule of 2Me4PN and in the selected solvents are gathered in Tables S6 and S7 in the SI. In Table 2, the theoretical CBS estimated bond lengths in vacuum and the studied solvents with the reported experimental X-ray data are compared [55]. Besides, the RMS deviations of CBS(T-5) values from experimental bond lengths are also given in Table 3.
Similarly to 4PN, the bonds only slightly change with the improvement of the basis set completeness and size (Table S4). A very good performance of B3LYP density functional, expressed by the complete basis set limit estimated bond lengths, in prediction of 2Me4PN structure is also apparent from a close reproduction of the experimental parameters of the parent molecule, 4PN in the gas phase, as revealed by microwave studies [50] (compare Tables 2  and 3). However, it is important to notice a longer C4 = O2 bond length in the crystal due to formation of H-bonding (1.244 Å) than our distance for a free molecule (1.225 Å). In analogy to crystal effects, water modeled by PCM increases this bond to 1.239 Å.
The angles in 2Me4PN molecule also depend on the quality of the basis set (the corresponding differences between angles obtained with cardinal numbers X = D and  Table S7). As before, polar solvents have stronger impact on the CBS predicted bond angles and the changes of up to 0.7° are calculated (Table S8). Considering the fact that the 4PN structure in the gas phase is perfectly predicted by B3LYP/CBS calculations, we believe that the corresponding 2Me4PN and chelidonic acid parameters in the gas phase are also predicted with a very high accuracy. However, in the crystal state, their structural parameters (mainly C = O) could differ due to intermolecular interactions.

Chelidonic acid
B3LYP/aug-cc-pVXZ calculated structural parameters of free (isolated) molecule of chelidonic acid are gathered in Tables S9 and S10 in the SI. The corresponding B3LYP/ CBS(T-5) calculated bond lengths in free chelidonic acid are compared with parameters, extracted from our crystal structure of chelidonic acid with methanol (Table 3).
Obviously, the experimental carbonyl bond length of 1.240 Å differs from the value calculated in the gas phase (1.219 Å) due to the formation of H-bonding in the crystal state. The corresponding CBS estimated bond angles of chelidonic acid are compared with experiment and shown together with the corresponding RMS values in Table S11 in SI.

Comparison of theoretically predicted bond lengths of chelidonic acid with experimental data
In Table 3, the corresponding B3LYP/CBS(T-5) bond lengths with experimental X-ray results, obtained in this research, for crystal network of chelidonic acid H-bonded with one molecule of solvent (CH 3 OH) are compared. An elongation of C = O bond by about 0.04 Å as a result of H-bonding in crystal was observed between the calculated and experimental data of chelidonic acid. In Fig. 8, we compare the theoretically predicted C2 = C3, C3-C4, C2-O1, and C4 = O2 bond lengths in three studied molecules with available experimental data. In general, theoretical modeling fairly accurately reproduce the experimental data in the three studied molecules, in particular those, which are determined in the gas phase. There are also some small bond lengths variations due to the presence of different substituents at C2 and C6 positions. The agreement between theoretical and experimental bond lengths, measured by X-ray, is negligible lower but still accurate.

Chelidonic acid -methanol dimer optimized from the crystal structure
The entire structure of mono methanol solvate of chelidonic acid is held together in the crystal by H-bonds and significantly weaker, long-range intermolecular interactions. Thus, it was interesting to test a possibility of free standing 1:1 molecular dimer formation starting from its crystal geometry.
In Fig. 9, the result of B3LYP/aug-cc-pVTZ optimization of chelidonic acid -methanol -is shown. The obtained gas-phase structure of this dimer is stable and no significant deformation of geometry from the initial crystal structure is visible. Thus, the formed H-bond lengths are fairly short

Dipole moment
High dipole moment indicates a possibility of strong interaction of solute with polar solvent or in the crystalline state. 4PN is a very polar molecule with an experimental dipole moment of 3.79 ± 0.02 D [50]. In the current study, the dipole moments of 4PN were calculated for a free molecule and in four solvents with increasing polarity, as characterized by their dielectric constant ε. The calculated dipole moment for a free 4PN molecule overestimates experiment by about 0.3 D (deviation from experiment is about 8.4%, see Table S12). Interactions with solvents modeled within the PCM approach result in a large increase of 4PN dipole moments (from 1.5 to 2.1 D). Its 2,6-dimethyl derivative and chelidonic acid are also predicted as highly polar structures (dipole moments of both molecules exceed 5 D, see Table 4).

NBO charges
Molecular dipole moment, reactivity, and electronic structure of the studied compound [19] are related to atomic charge distribution. In order to get some information about charge distribution in the studied molecules and the impact of two substituents at positions 2 and 6 in their six-membered ring, we have calculated a distribution of natural bond orbital charges [19]. Obviously, one should be aware that only differences between atomic charges and not their absolute values are important and should be interpreted [19]. The corresponding NBO values for the selected atoms are graphically presented in Fig. 10. As expected from structural considerations [56][57][58], the C2 and C4 atoms in 4PN ring bear positive charges while C3 and both oxygen atoms are negative. The attachment of methyl groups at positions 2 and 6 causes an increase in the polarity of all bonds in the molecule, i.e., the charges on C2 and C4 atoms increase while on C3, O1, and O2 atoms decrease. It is very interesting that the carboxylic groups at the 2 and 6 positions have a significantly lower effect than the methyl groups on the charge distribution in the 4-pyrone ring, and for most atoms, it introduces changes in the opposite direction. This apparently is due to the fact that the carboxyl group is an electron-withdrawing group while the methyl group has a positive inductive effect.

NICS
Aromaticity is an important concept in chemistry. In particular, aromaticity is very popular in organic synthesis and reactivity concept and is related to electron delocalization, planarity of ring structures, and equilibration of bond lengths, as well as to substituent effect [59]. Among different indexes characterizing aromatic (or anti-aromatic) nature of compounds are magnetic and geometric ones. Calculated NICS and HOMA values of the studied molecules and benzene, pyridine, pyrrole, and furan, selected as reference molecules are gathered in Table 5. Benzene with negative NICS(0), NICS(1), and NICS(1)zz values and HOMA parameter close to unity is a typical example of an aromatic compound. According to negative NICS values, pyridine, pyrrole, and furan selected as additional examples are also aromatic. However, a significantly lower value of HOMA (0.33) indicates lower aromaticity of furan. The last molecule, furan, as an example of oxygen containing heterocyclic ring, seems to be less aromatic according to its lower geometric index of aromaticity. This is in agreement with earlier report [60] on furans and oxazoles and could indicate a multidimensional aspect of aromaticity.
As seen from Table 5, all NICS(0) values of the three studied γ-pyrone molecules are slightly positive what indicates a very small anti-aromaticity in the middle of their ring plane. However, 1 Å above the ring plane a reverse sign of magnetic index of aromaticity is calculated. In addition, small and negative HOMA parameters are predicted what also points out a weak (or lack of) aromaticity of these compounds. Our results agree with earlier reported moderate aromaticity, calculated for α-and γ-pyrones [15]. In conclusion, only some electron delocalization 1 Å above the molecular plane exists in the studied molecules. Thus, these molecules are planar but not really aromatic and significantly differ from benzene, pyridine, pyrrole, or furan, calculated as reference molecules.

Conclusions
For the first time, the X-ray structure of methanol solvate of chelidonic acid was determined in this study and critically compared with reported structures of monomethyl ester of chelidonic acid and isophthalic acid. An important role of methanol molecules connecting two layers of  chelidonic acid was observed. This produced a unique crystal structure of dicarboxylic acids, held by H-bond motifs in the presence of bridging methanol molecules. From the total of 8001 dimeric structures of carboxylic acids available in CSD, 93.70% were typical carboxylic dimers (R-COOH…HOOC-R), 5.83% formed R-COOH… O = C motifs, and only 0.44% were carboxylic dimers bridged by methanol molecules with general formula R-COOH… (MeOH) 2 …HOOC-R.
B3LYP/aug-cc-pVXZ calculations in the gas phase and solution modeled by the PCM approach were successfully used to predict very accurate structure of γ-pyrone as a parent molecule of chelidonic acid and its 2,6-dimethyl derivative. Structures of 2,6-dimethyl-and 2,6-dicarbooxylic derivatives were also fairly accurately predicted and closely reproduced their X-ray structures. Large dipole moments (over 4 D) of the studied compounds were calculated in vacuum. In addition, a significant impact of solvent polarity on the predicted dipole moment was observed. According to the calculated geometric and magnetic indexes of aromaticity, the studied molecules are non-aromatic despite their planar ring structures.
The obtained structural parameters and dipole moments could serve as first step supporting detailed spectroscopic IR, Raman, and NMR characterization of numerous biologically active molecules, present in Chelidonium majus and similar plants. In addition, these results could provide a deeper insight for future experimental and theoretical spectroscopic studies of chelidonic acid, its complexes with metal ions and similar natural products, and their bioactive properties in hydrophobic and hydrophilic environments of living organisms.