Gleb V.
Baryshnikov
*ab,
Rashid R.
Valiev
cd,
Qizhao
Li
e,
Chengjie
Li
e,
Yongshu
Xie
e and
Hans
Ågren
af
aDivision of Theoretical Chemistry and Biology, School of Engineering Sciences in Chemistry, Biotechnology and Health, KTH Royal Institute of Technology, 10691, Stockholm, Sweden. E-mail: glibar@kth.se
bDepartment of Chemistry and Nanomaterials Science, Bohdan Khmelnytsky National University, 18031, Cherkasy, Ukraine
cResearch School of Chemistry & Applied Biomedical Sciences, National Research Tomsk Polytechnic University, Lenin Avenue 30, Tomsk 634050, Russia
dDepartment of Chemistry, University of Helsinki, FIN-00014, Helsinki, Finland
eKey Laboratory for Advanced Materials and Joint International Research Laboratory of Precision Chemistry and Molecular Engineering, Feringa Nobel Prize Scientist Joint Research Center, School of Chemistry and Molecular Engineering, East China University of Science & Technology, 130 Meilong Road, Shanghai 200237, P. R. China
fCollege of Chemistry and Chemical Engineering, Henan University, Kaifeng, Henan 475004, P. R. China
First published on 24th October 2019
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 1H 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.
In the present work we focus on the electronic structure and 1H NMR spectra of five recently synthesized17 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).
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. |
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 1H 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 (R2 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 1H NMR signal position for the corresponding protons involved in these contacts.
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.
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 π-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–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 π-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 π-electrons interact with the external magnetic field (not only those marked by the conjugation circuit). One should stress that our GIMIC results do not contradict but excellently confirm the rise of the aromaticity degree in the series of multiply annulated molecules 2 < 3 < 5 based on a comparison of the total diatropic current strength (21.8 < 22.1 < 28 nA T−1) and the difference in the chemical shifts of the inner and outer protons Δδout–in (7.19 < 10.70 < 12.51 ppm).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
Fig. 4 Comparison between the theoretically calculated and experimentally determined 1H chemical shifts for the single group of protons containing 63 atoms of NH and CH types. The R2 parameter corresponds to the linear correlation (green line) between δi(exp) and δi(theor). The idealized correlation δi(theor) = δi(exp) (black line) is presented for comparison. Experimental data are taken from ref. 17. |
R 2 | Slope (a) | Intercept (b) | MAE, ppm | MAE(CH), ppm | MAE(NH), ppm | |
---|---|---|---|---|---|---|
CH + NH (GIAO gas phase) | 0.981 | 1.164 | −1.207 | 0.330 | 0.208 | 0.980 |
CH + NH (CSGT gas phase) | 0.975 | 1.155 | −2.208 | 0.373 | 0.238 | 1.090 |
CH + NH (GIAO PCM) | 0.980 | 1.170 | −1.148 | 0.333 | 0.206 | 1.003 |
CH + NH (CSGT PCM) | 0.974 | 1.161 | −2.156 | 0.386 | 0.249 | 1.111 |
CH (GIAO gas phase) | 0.995 | 1.152 | −0.905 | — | 0.138 | — |
NH (GIAO gas phase) | 0.995 | 1.102 | −2.086 | — | — | 0.391 |
CH (CSGT gas phase) | 0.994 | 1.157 | −1.985 | — | 0.149 | — |
NH (CSGT gas phase) | 0.991 | 1.073 | −3.122 | — | — | 0.512 |
CH (GIAO PCM) | 0.995 | 1.162 | −0.864 | — | 0.126 | — |
NH (GIAO PCM) | 0.993 | 1.103 | −2.038 | — | — | 0.437 |
CH (CSGT PCM) | 0.994 | 1.166 | −1.952 | — | 0.143 | — |
NH (CSGT PCM) | 0.989 | 1.075 | −3.072 | — | — | 0.573 |
CH + NH (GIAO gas phase) direct | — | — | — | 0.553 | 0.345 | 1.654 |
CH + NH (CSGT gas phase) direct | — | — | — | 1.194 | 0.889 | 2.812 |
CH + NH (GIAO PCM) direct | — | — | — | 0.612 | 0.426 | 1.602 |
CH + NH (CSGT PCM) direct | — | — | — | 1.104 | 0.838 | 2.766 |
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 δi(exp) and δi(theor) cannot explain accurately the chemical shifts for the NH protons.
The correlation between δi(exp) and δ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 R2 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 R2, slope and intercept values for the linear correlation between δi(exp) and δ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 δi(exp) and δ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 R2 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 effects46–49 on the 1H 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 non-covalent 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.
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 (∇2ρ(r)) and electron energy densities (he(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.
Fig. 6 The structure of porphyrinoids 1–5 with the identified non-covalent interactions (dashed lines) identified by QTAIM analysis. |
Bond | d, Å | ρ(r), e a0−3a | ν(r), a.u. | g(r), a.u. | h e(r), a.u. | ∇2ρ(r), e a0−5 | ε | E, kcal mol−1 | E(curr),b kcal mol−1 [Jdia, nA T−1] |
---|---|---|---|---|---|---|---|---|---|
a d – interatomic distance, ρ(r) – electron density at the bond critical point (BCP), ν(r) – potential energy density at the BCP, g(r) – kinetic energy density at the BCP, he(r) – electronic energy density at the BCP (ν(r) + g(r)), ∇2ρ(r) – Laplacian of the electron density at the BCP, E – EML bond energy, ε – bond ellipticity; a0 – Bohr radius. b E(curr) was estimated from the linear relation Jdia = 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 | |||||||||
1 | |||||||||
N1H⋯N2 | 2.083 | 0.0241 | −0.0185 | 0.0188 | 0.0003 | 0.0763 | 0.029 | −5.8 | −6.1 [0.81] |
N3H⋯N2 | 2.058 | 0.0254 | −0.0196 | 0.0198 | 0.0002 | 0.0804 | 0.034 | −6.1 | −6.2 [0.82] |
2 | |||||||||
N2H⋯S | 2.360 | 0.0199 | −0.0157 | 0.0165 | 0.0008 | 0.0697 | 0.193 | −4.9 | −5.9 [0.79] |
N2H⋯N3 | 2.068 | 0.0246 | −0.0193 | 0.0197 | 0.0004 | 0.0807 | 0.063 | −6.1 | −6 [0.8] |
N3⋯S | 3.027 | 0.0166 | −0.0105 | 0.0114 | 0.0009 | 0.0495 | 0.063 | −3.3 | −2.1 [0.5] |
3 | |||||||||
N3⋯HN2 | 1.909 | 0.0348 | −0.0283 | 0.0278 | −0.0005 | 0.1094 | 0.056 | −8.9 | −9.5 [1.06] |
N1⋯HC13 | 2.288 | 0.0167 | −0.0118 | 0.0145 | 0.0027 | 0.0689 | 0.457 | −3.7 | −2 [0.5] |
4 | |||||||||
N1H⋯N2 | 1.946 | 0.0315 | −0.0253 | 0.0251 | −0.0002 | 0.0996 | 0.025 | −7.9 | −8.1 [0.96] |
N1H⋯S | 2.751 | 0.0106 | −0.0062 | 0.0080 | 0.0018 | 0.0391 | 0.430 | −1.9 | −4 [0.65] |
O⋯HC13 | 2.225 | 0.0191 | −0.0155 | 0.0172 | 0.0017 | 0.0754 | 0.709 | −4.9 | −6.1 [0.81] |
N5⋯HC6 | 2.126 | 0.0227 | −0.0166 | 0.0175 | 0.0009 | 0.0735 | 0.051 | −5.2 | −7 [0.88] |
5 | |||||||||
N1⋯HN2 | 2.352 | 0.0149 | −0.0109 | 0.0139 | 0.0030 | 0.0676 | 0.769 | −3.4 | −1 [0.43] |
N2H⋯N3 | 2.268 | 0.0166 | −0.0120 | 0.0130 | 0.0010 | 0.0561 | 0.143 | −3.8 | −1.9 [0.49] |
N3⋯HC16 | 2.234 | 0.0185 | −0.0129 | 0.0135 | 0.0006 | 0.0571 | 0.024 | −4.0 | −4 [0.65] |
However, the he(r) values are very close to zero for most contacts (up to order 10−4), meaning a close balance between the potential ν(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) and potential energy density ν(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 ν(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 1H 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 δi(exp) vs. δ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 δi(theor) 1H 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, C6F5-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 (Jdia) along a selected hydrogen bond and the energy of this bond (Ecurr). 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) ∝ Jdia 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) ∝ Jdia 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 (ε) 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,57i.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.
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 1H chemical shift values with high accuracy (the MAE is 0.138 and 0.391 ppm, respectively) and with an acceptable correlation criterion (R2 = 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.
Footnote |
† Electronic supplementary information (ESI) available. See DOI: 10.1039/c9cp04819g |
This journal is © the Owner Societies 2019 |