On the aromaticity of uracil and its 5-halogeno derivatives as revealed by theoretically derived geometric and magnetic indexes

The problem of aromaticity in heterocyclic rings of uracil and its 5-halogenoderivatives (5XU) was analyzed theoretically by calculating modified harmonic oscillator model of aromaticity (HOMA) for Heterocycle Electron Delocalization (HOMHED), nucleus-independent chemical shift parameters (NICS) and the so-called scan experiments, using helium-3 atom as a magnetic probe. The impact of halogen electronegativity on C5 atom’s NBO charges was also investigated. Water, as a polar environment, has a negligible impact on 5XU aromaticity. The most stable diketo tautomer shows a very low aromaticity while the “rare” dihydroxy form (tautomer No 6) is aromatic and resembles benzene. This is in agreement with traditional drawing of chemical formula of uracil’s six-membered ring, directly showing three alternating single and double bonds in its tautomer No 6. No good correlation between magnetic and geometric indexes of aromaticity for the studied 5XU tautomers was found. Linear correlation between the magnitude of NICS minimum, as well as the distance of the minimum above uracil ring plane center from 3He NMR chemical shift scan plot with respect to halogen electronegativity were observed. A strong linear dependence of magnetic index of aromaticity and the electronegativity of 5X substituent was observed.


Introduction
Uracil is an important component of living systems and belongs to small molecules which are responsible for information transfer [1]. Its structure could be derived from pyrimidine (1,3-diazine). Detection of the occurrence of uracil and pyrimidine in space [2] indirectly supports the fact that such molecules, containing two nitrogen atoms in a six-membered ring are stable and present in both extraterrestrial space matter and in living systems. For over 50 years, its 5-halogeno derivatives have been studied [3,4]. Among them, 5fluorouracil has been widely used as an anticancer drug [3,4].
The stability of uracil and its derivatives reflect their specific structure. Recently, on the basis of theoretical modeling using density functional theory (DFT at B3LYP-D3/aug-cc-pVQZ level of theory), it was postulated [5] that mutual arrangement of local C=O and N-H functional groups in uracil, as well as in 1,3-diazines, influences the electron density and local dipole moments distribution in such a way that it is energetically more favorable than in both 1,2-and 1,4diazines.
Stability of many organic molecules is often related to aromaticity [6][7][8][9][10]. Among them are cyclic organic compounds, with double bonds and heteroatoms with lone electron pairs included in their structure, which could be aromatic. The concept of aromaticity is based on equilibration of planar ring bond lengths due to the presence of shared and delocalized electrons. The number of delocalized electrons could be determined according to Hückel's formula 4n + 2, where n = 1 for benzene. There are two commonly used theoretical indexes of aromaticity, the geometric one (HOMA, and closely related parameters, derived from structural dimensions [11][12][13][14]) and magnetic [7,9,15,16] (NICS). For benzene, HOMA is close to 1. HOMHED is a reparametrized version of HOMA and is dedicated to heterocyclic molecules. Magnetic indexes of aromaticity are usually calculated in the middle of ring as NICS(0), 1 Å above that point, NICS(1), or its zzcomponent, NICS (1) zz . Negative values of NICS are calculated for aromatic compounds while positive for anti-aromatic molecular systems. Another magnetic probe of aromaticity is nuclear shielding and chemical shift of helium-3 atom [17][18][19][20][21]. Studies of aromaticity using 3 He are not so popular but they possess one advantage: in principle the predicted 3 He NMR chemical shifts could be verified by experiment [17,18,22].
For example, a good agreement was observed between a position of a narrow resonance signal in a 3 He NMR spectra of He@C60, as well as for other fullerenes and theoretically predicted isotropic value of helium-3 nuclear magnetic shieldings, converted to chemical shift [20,21,[23][24][25]. In many experimental studies, the fullerene cage size, symmetry, and aromaticity in the presence of the encapsulated 3 He atom [20,21,23,24] was studied by noble gas NMR technique. In such works, the 3 He atom was treated as another, independent magnetic probe of aromaticity.
In Fig. 1 are presented six optimized tautomeric structures of 5-fluorouracil, which is the most important 5-halogeno derivative of uracil. The structures are labeled with the commonly used abbreviations [5]. From the six potential uracil and its 5XU tautomers [26][27][28][29][30][31][32], the diketo form (tautomer No 1 in Fig. 1) is the most stable, as verified by experiment and theoretical calculations. On the other hand, it is possible to guess that tautomer No 6 should be the most aromatic one [5]. This tautomer is formed by three alternating single and double bonds, respectively. However, it is difficult "to guess" quantitatively the amount of aromaticity for tautomers 1 and 6 from the formulas, presented in Fig. 1, only. Besides, the nature of 5X substituent should also alter the aromaticity of the corresponding uracil derivatives. Moreover, the magnitude of magnetic shielding in the vicinity of uracil plane is related to the degree of electron delocalization and it should change perpendicularly to the plane with the distance R from the ring center. These phenomena could be modeled by so-called scanning experiment based on calculation of NICS(R) value as function of varied distance R [16,33,34] or by the change of 3 He isotropic nuclear magnetic shielding σ R , often converted to chemical shift δ R , when moving the helium probe by a distance R above the ring center [17,18].
Recently, the aromaticity of acenes was correlated with the HOMO-LUMO energy gap E g [22]. In case of 5XU, by changing the electron-donating properties of X substituent, the magnitude of E g and electron delocalization within the uracil plane should also vary [5]. Thus, by changing the halogen substituent one could expect a direct impact on E g , as well as on the aromaticity of 5-halogenouracils. However, the above mentioned problems were not discussed explicitly in the literature [5] or used to explain some structural and spectroscopic properties of 5XU molecules.
The aim of the current study is to characterize the aromaticity of uracil and its 5X-derivatives using geometric index of aromaticity HOMHED [35]. Additionally, a profile of magnetic field above the uracil ring will be modeled in separate scan experiments for NICS indexes and helium-3 chemical shift. The electron-donating or electron-accepting nature of X substituent could also influence the magnitude of charge at C5 atom, as well as ring aromaticity. Thus, the changes of charge at C5 and X atoms will be studied using natural bond orbital [36] concept (NBO). Finally, we will look for same correlations between uracil HOMO/LUMO gap, electronegativity of halogen substituent and aromaticity. Since solvent affects aromaticity indexes very little, we will model our compounds in vacuum, and only in some cases compare the results to those, obtained with the presence of water, incorporated by polarized continuum model (PCM) [37].

Computational part
All calculations were performed with Gaussian 16 program [38]. In this study, we selected a popular and fairly reliable B3LYP [39][40][41] hybrid density functional improved by inclusion of dispersion interactions via added Grimme's GD3 term [42], for modeling of uracil and its four 5-halogeno derivatives in vacuum and in the condense phase. The latter calculations were conducted using a polar surrounding (water) approximately described by PCM model [37]. The presence of all positive harmonic frequencies confirmed that the calculated structures were the true energy minima. Large Dunning-type correlation-consistent valence basis set of quadruple-zeta quality, augmented with a diffuse function (aug-cc-pVQZ) [43] was selected for better characterization of multiple bonds and lone electron pairs. Due to the problems with SCF convergence and optimization of 5-iodouracil tautomers, calculated with larger basis sets, we tested several combinations of available basis sets. Finally, all tautomers of 5 IU were optimized using 6-31 + G(d) for C, N, O, H, and 6-311G for I. We also managed to optimize the 5 IU1 and 5 IU6 tautomers with two larger basis sets, aug-cc-pVQZ for C, H, N, and O, and aug-cc-pVDZ-PP for I. The latter basis set was downloaded from EMSL basis set exchange [44,45]. Magnetic index of aromaticity NICS [7,9] was used to determine the aromaticity of the studied tautomers and monohalogenoderivatives of benzene as reference. In addition, geometric indexes of aromaticity, HOMA [12][13][14] and HOMHED [35], were calculated. NBO charges [36] were calculated at optimized geometries. Gauge including atomic orbital [46,47] (GIAO) NMR calculations of NICS parameters and 3 He isotropic nuclear magnetic shieldings at the ring center and above, at distance of 1 Å and R, from tautomers 1 and 6 of uracil and 5XU were also performed at the B3LYP/aug-cc-pVQZ level of theory (only for 5 IU, a previously mentioned 6-31 + G(d) for C, N, O, H, and 6-311G for I basis sets were used). The distance R from the middle of rings was varied from 0 to 2 Å with steps of 0.05 Å, from 2 to 5 Å with step of 0.5 Å and from 5 to 10 Å with steps of 1 Å. In particular, in agreement with our earlier studies [17,18], the theoretical chemical shifts of 3 He in the gas phase (δ He ) were obtained using the isotropic value of free probe atom nuclear magnetic shielding (σ o ) as reference and calculating its shielding σ R at a distance R from the middle of ring, according to the following formula At that point, we want to stress that the geometry of uracil and 5XU derivatives was optimized at B3LYPD3/aug-cc-pVQZ level of theory, but the subsequent GIAO NMR calculations were conducted without Grimme correction term since the results obtained with, and without dispersion corrections, were the same.
For graphical correlations, both linear and nonlinear least square fits were performed and their quality assessed from the corresponding R 2 parameter. The latter one used second order polynomial formula, abbreviated as polynomial_2.

HOMA, HOMHED, and NICS parameters of uracil and 5XU
The calculated HOMHED parameters for the studied structures are gathered in Table 1. It is apparent from Table 1 that all tautomers of uracil and its 5-halogeno derivatives are to some extent aromatic, as reflected by their fairly high HOMHED index (between 0.7 and 1). Following the suggestion of an anonymous reviewer, we also checked the aromaticity of 5FU6 tautomer with OH group at C4 position rotated and directed toward the X atom. This molecule was only negligibly less aromatic than the one with opposite direction of hydroxyl substituent (HOMHED of 0.984 vs. 0.987 for 5FU6). However, comparing the magnitudes of calculated aromaticity indexes, there is a clear difference between the most stable tautomer 1 and 6. Thus, lower values of  6 and are close to the value calculated for benzene (HOMEHED~1). The same order of aromaticity for the studied tautomers was earlier reflected by the calculated NICS and HOMA values in ref. [5] In order to test a potential correlation between geometric parameters of aromaticity and the magnetic ones, in Fig. 2 we plotted the HOMA and HOMHED values for all tautomers of uracil and 5FU vs. the calculated NICS(1) values. The HOMHED parameters were calculated in the current study and both NICS and HOMA ones were taken from our earlier work [5].
Interestingly, it is apparent from Fig. 2 that no good correlation between geometric HOMA, HOMHED, and magnetic index of aromaticity NICS(1) exists. However, the corresponding R 2 parameter for a dependence between NICS(1) calculated for all studied 5XU tautomers and geometric indexes is higher for HOMHED than for HOMA (0.82 vs. 0.67). Probably, this also illustrates the multidimensional character of aromaticity derived by using various parameters and points out to some advantage of using HOMHED combined with NICS(1) for characterization of heterocyclic compounds instead of traditional HOMA descriptor.

NICS scan above molecular planes of uracil and 5XU
A more clear picture of aromaticity profile above 5FU ring is shown in Fig. 3A and it somehow represents a "magnitude" of ring current at different distances from the molecular plane, e.g., above the center of the ring. As result of "scan experiment [16-18, 33, 34]" it is possible to show a characteristic profile of NICS(R) dependence on the distance R above the ring center of tautomers 1 and 6 of uracil, as well as benzene, taken as a model compound with a NICS showing a minimum at R min distance, markedly below 1 Å. The "scan" curve resembles a "Morse curve" and is related to the aromaticity profile above (and below) the benzene ring [10,[16][17][18]34]. However, in case of the diketo tautomer of uracil (U1), the minimum is very shallow what indicates a very small aromaticity. Interestingly, on the contrary to tautomer 1, tautomer 6 shows a scan curve closely resembling benzene. Thus, it is apparent that the significant depth of scan minimum, NICS(R min ) of about − 7.5 ppm for tautomer 6 indicates its high aromaticity while tautomer 1 is significantly less aromatic (NICS(R min ) is about − 1 ppm). Evidently, the curve depth of about − 11 ppm and pattern shape depends on the πelectron current density, being larger in benzene (NICS(R min ) than in tautomer 6 of uracil. In Fig. 3B, the results of similar "scan" experiments for tautomer 1 of all the studied 5XU1 compounds are shown. Analogous NICS(R) curve shapes with minimum in the range of R = 0.5 and 1 Å are observed for all 5XU except the tautomer 1 of 5-   (Fig. 3B). In the latter case the minimum is observed at R = 0 Å. The same situation was observed for pyrrole [18]. Besides, in that study we demonstrated that the position (R min ) and magnitude of NICS(R min ) could be roughly correlated with the aromaticity degree of ring compounds [18]. Analyzing curve patterns from Fig. 3B, one could observe some systematic changes of R min and NICS(R min ) parameters with the type of X substituent. Thus, the dashed red line, placed in the figure for the reader's convenience, starts at R min = 0 Å and goes nicely through all the curves minima. At this point, we also tested the impact of water, modeled by PCM, on the aromaticity of uracil's tautomer 1 and 6. The scan experiment for R changing from 0 to 2 Å was performed in vacuum and in water (see Fig. S1A and B in the supplementary material). The NICS(R) scan curves calculated in the presence of water nearly coincide with those, obtained in the gas phase. Thus, the subsequent calculations were performed in the gas phase only since the used model indicates a lack of effect of polar solvent on the aromaticity of uracil.
One could also expect a direct correlation between electronegativity of halogen substituent and the density of π-electron ring current, as reflected in "NICS scan experiments". In fact, nice linear (R 2 > 0.9) dependences of R min and NICS(R min ) vs. Pauling's electronegativity of X substituent are observed in Fig. 4A and B. The increase of X substituent electronegativity results in decreasing R min and a small increase of ring aromaticity (see Fig. 4B), as reflected by more negative value of NICS(R min ).
To support our data for 5XU molecular systems we also checked the impact of halogen electronegativity on monohalogenobenzenes C 6 H 5 X, where X = H, I, Br, Cl, and F. The NICS(R) scan experiment produced similar curve profiles to those, presented in Fig. 3B (see Fig. S2 in the Supplementary Material). However, the minima for halogenobenzenes were significantly lower (more negative) than for 5XU1 indicating a higher degree of aromaticity for these compounds (NICS(R min ) from about − 9 to − 11 ppm). In addition, for C 6 H 5 F, the calculated R min~0 .7 Å. Besides, strong linear correlations between R min , NICS(R min ) and electronegativity were observed for monohalogenobenzenes (Fig.  S3 in the supplementary material).

Electronegativity of X substituent and NBO charge at C5 atom
It is also important to check the impact of X substituent electronegativity on NBO charge at C5 atom in the most stable tautomer of 5XU1 molecules. It is apparent from Fig. 5 that changes of NBO charge of X substituent is reflected in formation of almost equal but opposite sign charge at carbon No 5. The observed changes of charge are very regular and point to the electronegativity as the main factor controlling the electronic impact on the halogeno-substituted uracil ring. The difference of charges (X -C5) also decreases in a regular, nonlinear way. Here, we want to stress, that in this case, no Nonlinear dependence between NBO charge of X substituent, C5 atom, their difference (X -C5), and electronegativity of X substituent in tautomer 1 of 5XU linear dependence exist. Instead, a second-order polynomial nicely fits the calculated points (R 2 > 0.99).
The observed patterns of NBO charge variations at C5 are similar to those for benzene ring carbon atom, directly bonded to halogen (see Fig. S4 in the supplementary material).
Our results are also in agreement with earlier studies on electronic effect of X substituent, in particular its donating and accepting properties, changing only to a little extend the aromaticity of benzene ring [48][49][50][51]. On the other hand, olefins are significantly more sensitive to such effects, including "push-pull cases") [52]. Therefore, the above data strongly point out to halogen electronegativity (and fluorine as an extreme case), reflected as an inductive effect (I), being responsible for small changes of aromaticity in 5-halogenouracils. On the other hand, halogen atoms also produce a resonance (mezomeric) effect (R) which is absent for X = H. The order of this effect is F > Cl > Br > I. Thus, each substituent shows a competing I and R effects and fluorine is an extreme example. It is also worth mentioning that in case of borazine derivatives, the fluorine substituent also resulted in NICS(R) scan curve with a minimum at R = 0 Å [53] (the reported curve patterns were similar to that in Fig. 3B). The behavior of F-substituent somehow differs from the other halogens probably because of its extremely low electronegativity. However, the fluorine anomaly is well known in the literature (see its anomalous reactivity in ref. [54]). 3 He atom as a magnetic probe of uracil and 5XU aromaticity Despite general acceptance of aromaticity as a very important description of structural, energetic and spectroscopic properties of chemical compounds and their reactivity, there exists many controversies in this field [10,51]. Aromaticity is often regarded as a multidimensional property. Thus, sometimes no direct correlation between the individual indexes of aromaticity could be seen. In this regard, in the current study we were interested to see, if there is some correlation between the geometric indexes of aromaticity for uracil and its 5-halogeno derivatives and their two, apparently not related magnetic indexes-NICS and helium-3 nuclear shieldings.
At this stage of our study, the helium-3 atom has been used as another magnetic probe of aromaticity. Isotropic values of 3 He nuclear magnetic shielding constants for the helium atom located 1 Å above the ring center for uracil and its tautomers, as well as for the remaining 5XU molecules were calculated. For brevity, in Table S1 in the supporting material are shown 3 He isotropic nuclear magnetic shieldings and chemical shifts of helium-3 atom 1 Å above the ring center (calculated with respect to isolated probe in the gas phase). In order to test a correlation between two magnetic indexes of aromaticity, δ( 3 He) and the NICS(1) for all tautomers of uracil and 5XU, in Fig. 6 were plotted the calculated parameters and fitted using a linear least square (LSQ) procedure. As result, a very good linear correlation between the two magnetic parameters of aromaticity, calculated in vacuum was produced (R 2 = 0.98888). This result suggests that in the future carefully designed experiments, 3 He NMR spectroscopy could be used to determine experimentally the "magnitude of aromaticity" of the studied molecular systems. In other words, we postulate to use helium-3 not only as a computational but also an experimental magnetic probe of aromaticity.
On the contrary, a significantly worse linear correlation between the chemical shift of the 3 He atom 1 Å above the middle of ring and the geometric index of aromaticity, HOMHED is apparent from Fig. 7A (R 2 = 0.84). This correlation is even worse for HOMA (not shown, R 2 = 0.71).

HOMO-LUMO energy gap and aromaticity of uracil and 5XU1
Recently, we reported on exponential correlation between E g and NICS index calculated for a series of acenes [22]. In the current studies, we also expected a similar correlation between E g and aromaticity indexes of uracil tautomers and their derivatives. In Table S2 are gathered B3LYP-D3/aug-cc-pVQZ calculated HOMO, LUMO, and E g energies for uracil tautomers, as well as their derivatives in the gas phase, and in Fig. 7B, a dependence between NICS(1) and E g is plotted.
The results presented in Fig. 7B are very scattered and once more demonstrate the fact that aromaticity is a very complicated and multidimensional phenomenon [10]. Contrary to the earlier work on regularly growing size of conjugated acenes [22], no correlation could be observed between E g and the calculated indexes of aromaticity for all tautomers of uracil Fig. 6 Linear dependence between B3LYP/aug-cc-pVQZ calculated δ( 3 He) for helium atom located at a distance of 1 Å above the middle of 5XU heterocyclic ring and NICS (1) and 5-halogenuracils (an attempt of the linear fit produced RMS of only 0.14304). However, by analyzing E g and NICS(1) of tautomer 1 only a nice linear correlation is observed with R 2 of 0.92 (See Fig. S5 in the supplementary material).

Conclusions
DFT modeling of molecular structure of six tautomers of uracil and its 5-halogeno derivatives was performed at DFT (B3LYPD3/aug-cc-pVQZ) level of theory. The aromaticity of the studied heterocyclic compounds was determined using geometric indexes of aromaticity, HOMA, and HOMHED. The most stable diketo tautomer 5XU1 showed very small aromaticity (HOMA and HOMHED about 0.5 and 0.76, respectively) while the "rare" U6 one was fairly aromatic (HOMA and HOMHED close to 1). GIAO NMR calculations of magnetic indexes of aromaticity, NICS(0), NICS(1), NICS(1)zz, and NICS(R), for R = 0 to 10 Å with step of 0.05 Å were also conducted on the previously optimized structures at B3LYP/aug-cc-pVQZ level of theory, e.g., without dispersion effects included. Electronegativity of X substituent and its induction effect was the main factor controlling the charge at C5 atom in 5XU tautomers. Regular changes of NBO charges at C5 (nicely fitted with second order polynomial formula), with respect to electronegativity of X substituent were observed for tautomer 1 of 5XU structures.
The last magnetic probe of aromaticity used was isotope of helium-3. Its nuclear shielding at R = 1 Å, as well as scans up to 10 Å were also calculated and converted to helium chemical shifts.
Both the geometric indexes of aromaticity and the magnetic ones reflected low aromaticity of tautomer 1 and high, comparable to benzene for tautomer 5XU6. High aromaticity of tautomer 6 confirms the suitability of its traditional formula, presented with three alternated single and double bonds. A linear increase of 5XU aromaticity, expressed by R min and NICS(R min ), obtained from scan experiments, with the electronegativity of halogen substituent was observed (R 2 about 0.92 and 0.93).
Linear correlations were observed for δ( 3 He) and NICS(1) with R 2 close to 1 (R 2 = 0.99). However, the linear correlation between δ( 3 He) and HOMHED with R 2 of 0.83 was less efficient (correlation for HOMA was even worse, with R 2 = 0.71). On the contrary, no direct correlation between NICS(1) for all tautomers and E g was observed (R 2 about 0.14). On the contrary, a strong linear correlation between these parameters was observed for tautomer 1.
It was also observed that the presence of a polar environment (water modeled by PCM) increased only negligible the aromaticity of 5XU.
Acknowledgments Calculations have been carried out using resources provided by the Wrocław Centre for Networking and Supercomputing (http://wcss.pl).
Code availability Gaussian 16 at WCSS Wroclaw.  Conflict of interest The authors declare that they have no competing interests.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.