Computational study of aromaticity, 1H NMR spectra and intermolecular interactions of twisted thia-norhexaphyrin and its multiply annulated polypyrrolic derivatives

The recently synthesized twisted thia-norhexaphyrin and its multiply annulated polypyrrolic derivatives have been studied computationally. Gauge-including magnetically induced current calculations predict a global nonaromatic character of the initial thia-norhexaphyrin due to the highly-twisted conformation of the macrocycle. Upon the oxidation of the thia-norhexaphyrin four multiply annulated polypyrrolic aromatic macrocycles are formed for which the global aromatic character is confirmed in agreement with experimentally measured H NMR spectra. The calculation of the proton chemical shifts for the studied compounds by direct comparison with the tetramethylsilane standard leads to a significant mean absolute error. At the same time a linear regression procedure for the two selected groups of protons (CH and NH protons) provides much better values of calculated chemical shifts and tight correlation with experiment. The separate consideration of NH protons is motivated by the numerous intermolecular hydrogen bonds in which the protons are involved, which induce considerable upfield shifts, leading to a significant underestimation of the corresponding chemical shifts. Such a selected correlation can be used for accurate estimation of proton chemical shifts of the related porphyrinoids. Bader’s theory of Atoms in Molecules has been applied for the studied twisted thia-norhexaphyrin and its multiply annulated polypyrrolic derivatives to characterize intramolecular H-bonds and other non-covalent interactions.


Introduction
Porphyrinoids represent a large extended class of porphyrins that have found practical application as active materials for photovoltaic cells, field-effect transistors, and organic electronics devices. [1][2][3] Also, porphyrinoids can be used as sensitizers in photodynamic therapy and bioimaging. 4,5 It follows that the electronic structure and aromaticity of numerous porphyrinoids have been investigated intensively during the last decade. [6][7][8][9][10][11][12][13][14] Particularly, it has been shown that the isophlorines are strongly antiaromatic and are air-stable at the same time. 14 On the other hand, the chemical and tautomeric transformations of them can lead to purely aromatic species, 12 meaning that porphyrinoids can be used as molecular aromatic/antiaromatic switches. Generally, porphyrinoid systems can demonstrate aromatic, antiaromatic or nonaromatic character depending on the type of heteroatoms included in the macrocycle, linkage moieties and outer substitutes. 15,16 The (non/anti)aromatic properties of expanded porphyrins, carbaporphyrinoids and other porphyrinoids determine their specific spectroscopic and photophysical properties, [6][7][8] which is a good case for computational studies.
In the present work we focus on the electronic structure and 1 H NMR spectra of five recently synthesized 17 unique porphyrinoids ( Fig. 1) -core-modified norhexaphyrin (1) containing a thiophene ring and four products of its oxidation: thiopyrrolo-pentaphyrin with an embedded pyrrolo [1,2]isothiazole (2), a sulfur-free pentaphyrin incorporating an indolizine moiety (3), thiopyranyltriphyrinoid containing a 2H-thiopyran unit (4) and a fused pentaphyrin containing a pyrrolizine moiety (5) which is formed upon desulfurization reaction of compound (2). Due to the extended and complicated structure of the studied systems there are multiple pathways of magnetically induced ring currents which can exist and compete between each other. In such cases gauge-including magnetically induced current (GIMIC) calculations can provide a solid explanation of the real distribution of ring currents and make it possible to determine the global (non/anti)aromaticity. Accurate calculations of proton chemical shifts for extended porphyrins and porphyrinoids remain an important task because they help to assign the experimentally observed 1 H NMR signals, especially for close-positioned signals, and to support the (non/anti)aromatic behavior of these species. Such calculations are still a challenge for computational chemistry because of the large size of the porphyrinoids and the variety of proton types (aromatic, aliphatic, methine protons of CH bonds, NH protons of pyrrole rings, etc.). In this paper we have analyzed the applicability of a linear regression procedure for the accurate calculation of chemical shifts for two separate groups of protons (CH and NH protons). The proposed technique works well within both the gauge-independent atomic orbital (GIAO) and continuous set of gauge transformation (CSGT) computational approaches at the DFT level, demonstrating tight correlations with experiment (R 2 up to 0.995). Finally, we have analyzed the electronic density distribution function for the studied molecules (1)-(5) within Bader's quantum theory of ''Atoms in Molecules'' (QTAIM) in order to identify and characterize the intermolecular non-covalent interactions and H-bonds that significantly affect the 1 H NMR signal position for the corresponding protons involved in these contacts.

Computational details
The ground electronic singlet state (S 0 ) structures of the studied porphyrinoids 1-5 have been optimized using the B3LYP/6-31G(d) [18][19][20] method at the DFT level starting from the initial molecular geometries obtained by the single-crystal X-ray diffraction technique. 17 Based on the optimized geometries QTAIM analysis 21 of the electronic density distribution functions has been carried out. The non-covalent interactions have been identified by the presence of (3, À1) type critical points and bond paths between the interacting atoms. The energy of the found non-covalent contacts was estimated through the Espinosa-Molins-Lecomte (EML) relationship: 22,23 E bond = 0.5n(r), where n(r) is the potential energy density (a.u.) at the corresponding (3, À1) bond critical point. The calculations of nuclear magnetic shielding tensors for the hydrogen atoms have been performed within the GIAO 24 and CSGT 25 methods at the B3LYP/6-311++G(d,p) level 18,19,26 of DFT. The role of solvent in the shielding constants was considered using the polarizable continuum model (PCM) 27 taking CHCl 3 as a model solvent. The final proton chemical shift values (d i (theor), ppm) were estimated by: (1) the direct comparison of computed isotropic shielding constants (s i(iso) , ppm) with the TMS standard (s TMS , ppm) calculated at the same level of theory: d i (theor) = s TMS À s i(iso) ; and (2) the linear regression method according to the following expression: d i (theor) = (s TMS À s i(iso) À b)/a, where the slope (a) and intercept (b) coefficients are estimated from the linear correlation between d i (theor) = s TMS À s i(iso) and d i (exp.). 28,29 All the above discussed calculations at the DFT level were performed using Gaussian16 software. 30 The QTAIM calculations were carried out using the AIMAll program package. 31 Magnetically induced current densities and current strengths have been calculated using the GIMIC code. [32][33][34][35] The NMR shielding calculations required for the GIMIC computations have been performed at the B3LYP/def2-TZVP 18,19,36 level of theory using the Turbomole program package. 37 The strength of the global magnetically induced current for the case of planarized porphyrinoids 2, 3 and 5 was estimated by applying an external magnetic field along the z axis perpendicular to the molecular xy plane. In the case of compound 4, the magnetic field was oriented perpendicular to the plane of the cyclopenta[a]pyrrolo[2,1,5-de]quinolizine 14p-electron system, while the rest of the molecule is almost coplanar with this fragment. In the case of compound 1 we have used the approach described in ref. 38. The current strengths were computed by the rotation of the applied magnetic field in the xz plane (assuming that the macrocycle is oriented in the imaginary xy plane) in the range 0-1801 with a step of 301. The maximal values of total current strength and corresponding current pathways for compound 1 are only discussed in the text.
Ring-current strengths (I, nA T À1 ) and current strength pathways for selected chemical bonds have been obtained by integrating the current density that flows through planes intersecting the chemical bonds.

GIMIC calculations of global aromaticity
As follows from the optimized and X-ray structure, thianorhexaphyrin (1) 17 possesses a highly-twisted conformation Fig. 1 The structure of the porphyrinoids studied in this work including the principal scheme of their synthesis starting from norhexaphyrin (1). Details about the synthetic procedure are published in ref. 17. of the macrocycle that prevents efficient p-conjugation and electron delocalization. GIMIC calculations indeed provide a very weak net diatropic ring current (0.4 nA T À1 ) for the twisted macrocycle of thia-norhexaphyrin (1) (Fig. 2). At the same time each pyrrole and thiophene ring sustains significant local diatropic ring currents (the current strength varies in the range 5.2-7.2 nA T À1 ). It means that thia-norhexaphyrin (1) can be considered as a globally nonaromatic twisted porphyrinoid containing six locally aromatic coupled heterocycles (five pyrroles and one thiophene). Moving to the thiopyranyltriphyrinoid (4) containing a 2H-thiopyran unit, two different cyclic subsystems can be classified. The first one is a cyclopenta[a]pyrrolo[2,1,5-de]quinolizine 14p-electron system which possesses aromatic character (the net diatropic current strength is 14.7 nA T À1 ). The calculated currents along the three N-C bonds inside this fragment are very small (1.2, 0.3 and 0.9 nA T À1 ) meaning that the inner N atom is not involved in the aromatic system. The remainder of the molecule (second subsystem) containing the 2H-thiopyran cycle sustains a net non-aromatic character because of the broken p-conjugation within the saturated -CO-CH 2 -CHz fragment (Fig. 2). As a result only a weak local aromaticity can be postulated for the two pyrrole rings.
In contrast to compounds 1 and 4, porphyrinoids 2, 3 and 5 demonstrate clear global aromatic character and significant net diatropic ring currents (21.8, 22.1 and 28 nA T À1 , respectively) in agreement with the previously reported ACID and NICS calculations. 17 However, our GIMIC calculations make it possible to determine the multiple global aromatic pathways within each porphyrinoid molecule (2, 3 and 5) based on the absolute value of the current strength along the selected bond and quantitative splitting of this current at the bifurcation points (marked by turquoise spots in Fig. 3). For instance, for thiopyrrolo-pentaphyrin (2) there are two most important bifurcation points (I and II, Fig. 3). The first one (I) represents a splitting of the incoming current (21.8 nA T À1 ) into two flows that circulate through the outer perimeter (15.9 nA T À1 ) and through the inner subperimeter (5.9 nA T À1 ) (Fig. 3). Surprisingly, at the bifurcation point (II) the main part of the diatropic current circulates along the C-N bond (15.2 nA T À1 ) and only a small amount (6.6 nA T À1 ) goes to the C-C bond (Fig. 3).
This fact actually means that the nitrogen atom in the pyrrolizine fragment is involved in the global aromatic electronic structure of thiopyrrolo-pentaphyrin (2). The same situation (i.e. involvement of the nitrogen atom in the pyrrolizine fragment in the global aromatic pathway) is also observed for the fused pentaphyrin containing a pyrrolizine moiety (5), which is similar to thiopyrrolo-pentaphyrin (2) but does not contain the -S-linker ( Fig. 1). Despite the fact that the initial strong diatropic current (28.0 nA T À1 ) is almost equally split at the bifurcation point (V), at the next two points (VI) and (VII) this current (15.0 nA T À1 ) continues to split onto equivalent branches, so it is hard to distinguish the global aromatic pathway in this case. At the same time, the diatropic current in the outer perimeter (15.0 nA T À1 , green lines and arrows in Fig. 3) does not change significantly and becomes even higher because of the combination with the smaller currents from the inner p-delocalized system. Interestingly, compound 3, being very similar to compounds 2 and 4, demonstrates a short global aromatic pathway avoiding the multiply annulated fragment (Fig. 3). In the bifurcation point (III) the initial diatropic current (18.1 nA T À1 ) splits almost equally -the first pathway (8.7 nA T À1 ) gradually dissipates over the multiply annulated part, while the second pathway (9.4 nA T À1 ) circulates through the inner perimeter forming a short global aromatic pathway (Fig. 3, green perimeter in compound 3).
It is notable that the global aromatic pathways calculated by the GIMIC method (marked by green in Fig. 2 and 3) cannot be identified with the conjugation circuits usually presented for the porphyrinoid systems. 17,[39][40][41] The conjugation circuits are based on general chemical principles and single-crystal X-ray information about single/double bond alternation. Rather, the current pathways predicted by GIMIC demonstrate the real distribution of the magnetic currents due to the perturbation of the delocalized p-electronic molecular shell by the external magnetic field (usually applied perpendicular to the molecular plane). Actually, the number of electrons in the conjugation circuit and its satisfaction or unsatisfaction to the Hückel (4n + 2) rule need not be fulfilled strictly for complicated porphyrinoids because all of the delocalized p-electrons interact with the external magnetic field (not only those marked by the conjugation circuit). One should stress that our GIMIC  17 The GIMIC results are also in excellent agreement with the previously published NICS(0) local aromaticity indices for all studied porphyrinoids 1-5. 17

1 H chemical shift calculations
Generally, accurate calculations of the proton chemical shifts of porphyrinoid molecules still pose a challenge for computational chemistry. Actually, only a limited number of reports have been published 10,12,14,[42][43][44][45] in which the 1 H chemical shifts of porphyrin and some porphyrinoids are calculated and discussed. In most previous studies the direct comparison of the calculated isotropic shielding constants with the TMS standard was used to predict the 1 H chemical shifts. Generally, this approach provides significant mismatch between the experimental and theoretically estimated chemical shifts (up to 1 ppm for the CH protons and up to 4 ppm for NH protons). For our systems we observe the same trend (  (Table 1). 28 Especially, it relates to the low-and high-field regions (around (À5) to (À4) ppm and 10-14 ppm). The CSGT/B3LYP/6-311++G(d,p) approximation demonstrates a slightly smaller slope parameter (1.155) and thus a smaller contribution of systematic error, but the correlation coefficient is even smaller relative to the GIAO approach (R 2 = 0.975), meaning a larger random error contribution (Table 1). Additionally, the intercept parameter is almost two times higher in absolute value for the CSGT/B3LYP/6-311++G(d,p) approximation relative to the GIAO one (À2.2078 vs. À1.2072, Table 1), which also indicates a stronger systematic underestimation of 1 H chemical shifts computed by the CSGT scheme. The estimated MAEs for the combined number of 63 protons (53 CH and 10 NH atoms) equal 0.330 and 0.373 ppm, respectively, by the GIAO and CSGT schemes (Table 1 and Table S2, ESI †). However, the MAEs for the particular groups of CH and NH protons within these combined correlations are very different (0.208 vs. 0.980 ppm and 0.238 vs. 1.090 ppm, Table 1) meaning that only the linear relation between d i (exp) and d i (theor) cannot explain accurately the chemical shifts for the NH protons.
The correlation between d i (exp) and d i (theor) can be improved by separating NH and CH protons into two different groups (53 and 10 atoms, respectively, Table S3, ESI †) for which each particular correlation is more tight (Fig. 5) relative to those for the overall 63 hydrogen atoms. Indeed, for the CH protons the R 2 parameter within the GIAO and CSGT schemes is around 0.995, which is an acceptable value for a well-performing computational method. 28,29 The slope parameter is still around 1.15-1.16, demonstrating a considerable contribution of systematic error, but the intercept parameter becomes significantly smaller for CH protons, indicating a decrease of systematic underestimation (Table 1). At the same time, the R 2 , slope and intercept values for the linear correlation between d i (exp) and d i (theor) of NH protons are a bit worse compared to those for CH protons but much better compared to the generalized correlation for all 63 protons (Fig. 4). The resulting MAE values for the computed chemical shifts of NH protons are also much better relative to those obtained by the combined linear regression method as well as with a direct comparison with the TMS standard scheme (Table 1). It means that separating the different types of protons provides much better agreement between d i (exp) and d i (theor) values for the CH and NH protons in terms of a linear regression procedure, while the direct comparison with the TMS reference is much worse by the MAE criterion (Table 1, especially for NH protons). Accounting for solvent effects by the polarizable continuum model (PCM) almost does not affect the R 2 and slope parameters of the proposed correlations at all (Table 1), but definitely improves the intercept parameter. More precise calculations considering the influence of zero-point and thermal vibrations, and dynamic and relativistic effects 46-49 on the 1 H chemical shifts of porphyrinoid molecules should improve the correlation parameters in the future, however, such calculations are very time and resource consuming for large porphyrinoid systems.
Finally, we conclude that gas phase B3LYP/6-311++G(d,p) calculations of shielding constants by both the GIAO and CSGT schemes using the linear regression procedure for the different groups of protons (NH and CH in our case) demonstrate a sufficiently good accuracy of predicted chemical shifts as well as an optimal balance between accuracy and time/resource expenses. As an additional remark we can suppose that intramolecular non-covalent interactions affect significantly the chemical shifts of NH protons and generally provide a significant underestimation trend (the intercept parameter is large) in contrast to the CH protons, which are less sensitive and less active with respect to the inter-and intramolecular noncovalent interactions. Due to this reason, we will focus in the next section on the particular non-covalent interactions and  their electronic parameters calculated at the QTAIM level for the studied porphyrinoid molecules 1-5.

Validation of 1 H chemical shift calculations
Several previous studies have demonstrated that functionals with a low percentage of Hartree-Fock exchange (HFE) should not be employed for conjugated (and, particularly, aromatic) systems. [50][51][52] More recently, in the context of porphyrins, it has been also commented by Casademont-Reig et al. 53 and Torrent-Sucarrat et al. 54 That is why, in order to validate the applicability of the B3LYP/6-311++G(d,p) technique (20% HFE) with respect to the role of the HFE amount, we have recalculated 1 H chemical shifts by the hybrid meta generalized gradient M06-2X(54% HFE) approach using the same 6-311++G(d,p) basis set. The M06-2X/6-31G(d) optimized geometries were used for the 1 H chemical shift calculations discussed next in this section. All the results are summarized in Fig. S1 and S2 and in Table S4-S7 (ESI †). One should stress that one can see from Table S4 that the R 2 parameter for all correlations based on the M06-2X functional does not exceed 0.988, meaning the higher contribution of random error within the applied computational scheme. It should be noted that by the R 2 criterion none of the M06-2X based correlations satisfy the minimal acceptable value (R 2 = 0.995) for a well-performing computational method. 28,29 However, the slope and intercept values are much closer to the ideal 1 and 0 values, respectively, meaning the much smaller contribution of systematic error for the M06-2X based correlations. Comparing the MAE values in Table 1 and Table S4 (ESI †), one should stress that the CSGT approach demonstrates much better performance for the direct and generalized (CH + NH) schemes when combined with the M06-2X/6-311++G(d,p) method rather than with the B3LYP/6-311++G(d,p) approach. However, comparing the particular correlations for the 1 H chemical shifts within the CH and NH manifolds, B3LYP/6-311++G(d,p) provides the better performance by the MAE criterion (no more than 0.149 ppm for CH protons and no more than 0.573 ppm for NH protons) relative to the M06-2X/6-311++G(d,p) approach (the minimal MAE equals 0.211 and 0.585 ppm for CH and NH protons, respectively, Table S4, ESI †). It means that the proposed differentiation scheme for the evaluation of 1 H chemical shifts within separate sorts of CH and NH protons works adequately at the B3LYP/6-311++G(d,p) level and the considerable systematic error within this approach can be effectively reduced by the linear regression procedure.

Intramolecular non-covalent interactions and tautomerism
Before we start the analysis of non-covalent interactions in the studied porphyrinoids we have investigated the possible tautomerism of these molecules by changing the position of the NH proton within the porphyrinoid macrocycle. The results are collected in Fig. S3 (ESI †). We have found that for all cases the structure of porphyrinoids 1-5 established by single-crystal X-ray diffraction corresponds to the lowest energy tautomer Fig. 5 Comparison between the theoretically calculated and experimentally found 1 H chemical shifts for the two groups of CH (52 atoms, green spots) and NH (11 atoms, orange spots) types, calculated by both the GIAO and CSGT schemes with and without accounting for solvent effects at the PCM level. The R 2 parameter corresponds to the linear correlation (black line) between d i (exp) and d i (theor).
presented in Fig. 1. All other tautomers for compounds 1-3 and 5 are higher in energy by no less than 3.9 kcal mol À1 , which implies that the lowest ratio between the main tautomer and another one is only 1000 : 1 (if one considers the Boltzmann distribution at room temperature). Only in the case of compound 4 are the two possible tautomers (4 and 4a in Fig. S3, ESI †) quite close in energy, giving a distribution ratio 12 : 1, but accounting for the energy barrier for direct proton transfer along the NHÁ Á ÁN coordinate this ratio should be much higher. Indeed, in all cases for compounds 1-5 only the main tautomer was isolated in the single crystal state and no additional signals were identified in the 1 H NMR spectra of the studied porphyrinoids. 17 Therefore in this section we focus our attention only on the intramolecular interactions occurring in the main tautomer form of porphyrinoids 1-5 (Fig. 1). All the non-covalent interactions identified in the investigated porphyrinoids (Fig. 6) can be generally classified into two main groups: H-bonds and van der Waals (vdW) interactions. A few types of H-bonds can be identified in the studied porphyrinoids (NHÁ Á ÁN, NHÁ Á ÁS, CHÁ Á ÁN, CHÁ Á ÁO, CHÁ Á ÁF, and CHÁ Á ÁC) and many sorts of vdW contacts (FÁ Á ÁC, NÁ Á ÁC, FÁ Á ÁF, NÁ Á ÁS, OÁ Á ÁC, and HÁ Á ÁH) exist in these molecules. The electronic characteristics of all these interactions are collected in Table S4 (ESI †). We have performed QTAIM analysis both for the B3LYP/6-31G(d) optimized geometries and for the experimental geometries taken from X-ray analysis. Generally, intramolecular H-bonds of only the NHÁ Á ÁN and CHÁ Á ÁN types occur in molecules 2, 3 and 5, while S-containing porphyrinoids 1 and 4 also possess CHÁ Á ÁO and NHÁ Á ÁS interactions. The main parameters of these H-bonds only slightly depend on the initial geometry (optimized or X-ray, Table S4, ESI †), which is why only the data for B3LYP/6-31G(d) optimized geometries are presented in Table 2. Accounting for the positive values of the Laplacian (r 2 r(r)) and electron energy densities (h e (r)), one can formally classify all the found non-covalent interactions (Table 2 and Table S4, ESI †) as closed-shell interactions (in the framework of the QTAIM formalism), which are characterized by the minor sharing of electron density in the interatomic space.
However, the h e (r) values are very close to zero for most contacts (up to order 10 À4 ), meaning a close balance between the potential n(r) and kinetic g(r) energy densities, which are negative and positive, respectively, by definition. 55 In such cases, the absolute values of electron density r(r) and potential energy density n(r) at the corresponding bond critical points (BCPs) should be considered and used for the estimation of the bond strengths. The EML bond energies based on the n(r) values are predicted to be in the range 3.8-8.9 kcal mol À1 for NHÁ Á ÁN bonds, 3.7-5.2 kcal mol À1 for CHÁ Á ÁN bonds, and 1.9-4.9 kcal mol À1 for NHÁ Á ÁS and CHÁ Á ÁO contacts. It means that all the NH protons in the studied porphyrinoids are involved in a network of strong intramolecular bonds, which causes significant electron cloud depletion around these H atoms. This could be a reason why the conventional DFT approach strongly underestimates the corresponding shielding constants used for the 1 H chemical shift prediction. In this case, the linear regression scheme discussed above seems to be a reasonable way to exclude such systematic effects.
The fact itself that accounting for solvent effects within the PCM model improves the intercept parameters of the d i (exp) vs. d i (theor) correlations (Table 1) also supports the idea that the correct polarization of NH protons (but not only) is a proper way to accurately estimate d i (theor) 1 H chemical shifts. The vdW interactions presented in Table 2 and Table S4 (ESI †) are expectedly less strong than the H-bonds, but they are very important for the stabilization of the conformational structure of the studied porphyrinoids. Particularly, C 6 F 5 -groups form numerous weak vdW interactions (mainly of FÁ Á ÁC nature) with the central macrocycle, stabilizing the positions of these side substituents.
In order to validate the EML calculated energies of hydrogen bonds presented in Table 2 we have applied the technique proposed by Kaila et al., 56 which is based on the linear correlation between the diamagnetic part of the magnetically induced current density ( J dia ) along a selected hydrogen bond and the energy of this bond (E curr ). The results are summarized in Table 2. One can observe excellent agreement between the EML calculated bond energies and those estimated by the E(curr) p J dia dependence. It is interesting to note that the agreement is better for strong hydrogen bonds (generally, the deviation is no more than 1 kcal mol À1 ) but for weak interactions the deviation can reach 2.4 kcal mol À1 . It can be explained by the fact that the E(curr) p J dia relation by Kaila et al. 56 was calibrated using reference H-bond energies no less than 3 kcal mol À1 . Table 2 and Table S4 (ESI †) also consist of ellipticity (e) values, which characterize the curvature of the bond path between the interacting atoms. Actually, a high ellipticity (more than 0.5) represents a dynamic instability of these interactions, 21,57 i.e. these bonds are stable only upon the fixed position of atoms. It explains why in moving from B3LYP/6-31G(d) optimized to X-ray taken geometries the main electronic parameters (that are distance-dependent) are changed only slightly except for the ellipticity, which is interatomic disposition dependent. Concluding shortly, we should note that QTAIM analysis cannot be considered as a chemically accurate method for electronic structure calculations but with respect to the characterization of non-covalent interactions it seems to be the most powerful tool developed until now. Our QTAIM calculations also shed light on the very complicated electronic surroundings of the NH protons in the studied porphyrinoids, which become even more sophisticated in the presence of external magnetic field induced magnetic dia-or paratropic currents.

Conclusions
In the present paper the electronic structure, aromaticity and 1 H chemical shifts have been computationally studied for five recently synthesized porphyrinoids -multiply annulated polypyrrolic derivatives 2-5 of twisted thia-norhexaphyrin (1). The magnetically a d -interatomic distance, r(r) -electron density at the bond critical point (BCP), n(r) -potential energy density at the BCP, g(r) -kinetic energy density at the BCP, h e (r) -electronic energy density at the BCP (n(r) + g(r)), r 2 r(r) -Laplacian of the electron density at the BCP, E -EML bond energy, e -bond ellipticity; a 0 -Bohr radius. b E(curr) was estimated from the linear relation J dia = aE(curr) + b, where a and b are dimension coefficients equal to À0.078 (nA T À1 )/(kcal mol À1 ) and 0.345 nA T À1 , respectively, at the B3LYP/def2-TZVP level of theory. 56 induced current density calculations predict a non-aromatic character of thia-norhexaphyrin (1), partial aromatic character of thiopyranyltriphyrinoid (4), which contains a 2H-thiopyran unit, and strong global aromaticity of pentaphyrin with an embedded pyrrolo [1,2]isothiazole (2), a sulfur-free pentaphyrin (3) incorporating an indolizine moiety and fused pentaphyrin (5) containing a pyrrolizine moiety. These results are in excellent agreement with previously published experimental observations. Rather surprisingly, for compounds 2 and 5 we found that the pyrrolizine nitrogen atom is involved in the global aromatic pathway, implying significant electronic delocalization over the pyrrolizine fragment. Such a phenomenon has never been observed before despite the fact that a pyrrolizine fragment frequently appears in porphyrinoid systems. 39,40 Considering the groups of CH and NH protons as distinct, calculations using the GIAO/B3LYP/6-311++G(d,p) method, and analyzed using a linear regression scheme, made it possible to estimate the 1 H chemical shift values with high accuracy (the MAE is 0.138 and 0.391 ppm, respectively) and with an acceptable correlation criterion (R 2 = 0.995). Computations within the continuous set of gauge transformations (CSGT) scheme provide comparable results relative to the GIAO method but gave slightly worse results. Accounting for solvent effects on the isotropic shielding constants did not improve the correlation tightness, but reduced the systematic underestimation of the computed chemical shifts. We assign this effect to the fact that all the NH protons of the studied molecules are involved in moderate or strong H-bonds and in weak van der Waals interactions. This provides a stronger depletion of electron density around the NH protons compared to the CH protons, which usually are inactive with respect to inter-and intramolecular interactions.

Conflicts of interest
There are no conflicts to declare.