Magnetic Aromaticity of Cycloporphyrin Nanorings

: The ascertainment of magnetic aromaticity is not necessarily straightforward, especially for large and bent systems, such as the cycloporphyrin nanorings recently synthesized by the group of Anderson. Six of these cycloporphyrin nanorings were studied here computationally. Indirect methods, based on nuclear shielding and magnetizabilities, and direct methods, based on standard quantum mechanics, were both used effectively to determine their magnetically induced current strength, which mostly conﬁrmed Anderson’s classiﬁcation. However, in the case of hexanions, and in particular for cyclohexaporphyrin hexacations, a signiﬁcant cancellation of delocalized diatropic and paratropic ﬂow occurred, showing that the resultant faint aromatic character was a result of competing aromatic and antiaromatic contributions, as also evidenced by the ipsocentric method. A warning is renewed on the use of isotropic shielding to determine the tropicity of the magnetically induced current.


Introduction
The research on aromaticity is kept always active by novel theoretical endeavors to manage its fuzziness [1] and by the synthesis of new conjugated polycyclic molecules of increasing complexity [2][3][4], which force the extension of the concept of aromaticity far beyond pioneering, yet still useful landmarks as the aromatic sextet [5]. Of the many tools used to assess aromaticity, those based on the magnetic response are currently the most used. That response was first investigated experimentally, first through measurements of magnetizabilities and later of nuclear magnetic shieldings [6]. Then, models were developed to deduce the tropicity of the magnetically induced current from the experimental values. The terms diatropic and paratropic have been proposed in an experimental overview of chemical shifts of annulenes to describe molecules disclosing a diamagnetic or a paramagnetic ring current, respectively [7]. Eventually, theoretical and computational advances allowed reliable computations, not only of the experimental values, but also of the current density tensor field itself, which is generated by the external magnetic field in the molecular domain and determines by integration the experimentally observable magnetic properties [8][9][10][11][12][13]. Perhaps, the availability of reliable and handy software to compute the induced current can gradually limit the studies based on its indirect determination, but as a matter of fact, as a followup of the old indirect way of retrieving the tropicities, the computations of the nucleus-independent chemical shift (NICS [14]) continue to be widely used as a tool to grasp the main features of the magnetic aromaticity of a molecule. This practice is not straightforward, neither at the computational nor the qualitative level [15][16][17]; however, it must be considered that quantitative agreement has been obtained between integrated current strengths and NICS scans or values at appropriate heights in monocycles [18][19][20], and an overall qualitative agreement has been obtained for planar polycyclics [21].
Albeit not the only one, magnetic aromaticity assessed through indirect methods has certainly been a basic tool in the classification of the beautiful nanorings of cycloporphyrins, recently synthesized by Anderson and coworkers [22][23][24]. Large nanorings of tetracations and hexacations of cycloporphyrins have been classified as aromatic/antiaromatic, and the classification agrees with the Hückel rule when the electrons are counted along the main circuit of current flow. Those nanorings are likely the largest experimental rings to follow the rule so far. The classification as aromatic for one of the species discussed by Anderson has been recently challenged by Matito et al. [25], who stated that an inadequate computational method led to an erroneous conclusion.
Considering the technological relevance of small-gap organic systems, such as the cycloporphyrin nanorings, and the widespread use of the magnetic response in addressing their characterization, the understanding of the above divergent results seemed to us of primary importance. Here, we report our results on the magnetic aromaticity of the six experimentally synthesized [24] cycloporphyrin nanorings c-P6 Q , c-P7 Q , and c-P8 Q , where the charge Q is either +4 or +6.

Materials and Methods
Geometries from [24], obtained at the LC-ωhPBE [26] (ω = 0.1)/6-31G*), were taken as the starting points. After the computation of the magnetically perturbed wavefunction (the CSGT method [8,27]) on those geometries with the same method of [24], they were reoptimized with the same basis set, but using the BHandHLYP functional [28], ranked best in a recent benchmark of computations of magnetizabilities [29] and among the best in a benchmark based on hypervirial relationships [30]. Magnetizabilities and the magnetically perturbed wavefunction have been computed once again with the CSGT method. The newly obtained geometries are all true minima of symmetry D 2h , D 3h , C s , C s , D 2h , and C s for c-P6 4+ , c-P6 6+ , c-P7 4+ , c-P7 6+ , c-P8 4+ , and c-P8 6+ , respectively.
These computations were performed using Gaussian 16 [31]. NICS scans have also been computed with Gaussian16, but with GIAO, rather than CSGT, to avoid the spurious dependency of the results on the number of computed points in the scan, resulting from the present implementation in Gaussian 16 [32]. The computation of the magnetically induced current density and its integration and visualization were performed using SYSMOIC [13] with the CTOCD-DZ2 [33,34] method (CTOCD-DZ1 [34,35] for orbital decomposition).

From Magnetic Properties to Ring Current Models
Ring current models (RCMs) have been traditionally used to interpret the values of magnetic properties, such as magnetizabilities and chemical shifts [36]. The simplest of such RCMs, called infinitely thin circular loop of current (ICLOC) [18,19], can be considered an archetypal reference model. Even for this simplest model, the induced magnetic field cannot be generally expressed in terms of elementary functions: computations of the ICLOC model using the needed elliptic functions have been rarely reported in the literature [24,37,38]. The expression of the field simplifies at large distances, where it becomes the field of a magnetic dipole, and along the axis of the loop. Such a loop of radius s, placed in the xy plane, generates a contribution to the shielding at the center of the loop given by: where the vacuum permeability µ 0 = 4π × 10 −7 N A −2 = 4π × 10 −6 Å(nA) −1 T, is the contribution given by the loop to the parallel component of the magnetizability tensor and I B is the signed current strength, which is the amount of current induced by a unitary magnetic field perpendicular to the loop in the linear approximation and is positive or negative for a paratropic or a diatropic circulation, respectively [19]. While the magnetic dipole approximation was rather common in the older literature, the possibility to compute NICS and magnetizabilities has given the occasion to consider the ICLOC model in more recent times [18,19,[39][40][41][42]. The ICLOC model has two parameters: the radius of the loop s and the signed current strength I B . The determination of the current strength by a single observable (magnetizability or proton chemical shift, preferably the single parallel components ξ zz or σ zz ) assuming the geometric radius of the ring often leads to significant errors. On the other hand, the combined use of two observables ξ zz and σ zz (0) has given a considerable improvement for a set of diverse monocycles [43,44]. The coupling of these two observables can be used to define an effective average radius: which when plugged into Equation (1) or (2) gives the signed current strength. The cases where s av turns out to be slightly larger than the geometrical radius s can be interpreted in terms of a model that has a finite width, the toroidal circular loop of current (TCLOC). From a mathematical point of view, the equations developed for this model should not be considered for s av > (1 + 2/3)s = 1.29s (α < 2/3 in Equation (26) of [18]). However, whenever the difference between s av and the geometric radius exceeds the physical spread of the current (roughly smaller than 1 Å), there is a hint of the inadequacy of both the ICLOC and TCLOC models.
The definition of the more advanced RCMs developed to date, requiring three or four parameters, calls for more input data, typically taken from NICS zz scans. In an investigation of a family of diverse monocycles, particular merit has been credited to (i) the ICLOC2 model, which consists of two identical ICLOCs equally displaced off the molecular plane (three parameters: I B , s, and the off-plane displacement z) and (ii) the ICLOC2C model, which consists of two concentric coplanar ICLOCs (four parameters: I B 1 , I B 2 , s 1 , and s 2 ). The ICLOC2 and ICLOC2C models are the representative models for the π and σ currents in simple monocycles [19]. Remarkably, at odds with statements based on a small set of computational results [45], both of them allow the onset of a maximum or a minimum in the NICS zz scan, although for ICLOC2, this nonmonotonous plot only occurs for small cycles, such that the displacement z of the loops along the axis is such that s/z < 2 [19]. Tables 1 and 2 report relevant values for simple RCM modeling of the nanorings according to the computational protocols used in this paper. In the last column of the tables, we anticipate the DFT signed current strength, to be discussed below. Starting from Table 1, we notice that at least one of the two ICLOC-based estimates of the current strength reported in [24] (obtained through a comparison of the off-axis field to either experimental chemical shifts or computationally demanding three-dimensional NICS zz scans) matches reasonably well with the DFT current strength (Appendix A). We also note that the computations based on the s av method basically have a better accuracy, although they are far simpler to perform. Moreover, the s av computations has the added value that it indicates that an ICLOC modeling for c-P6 6+ is suspicious, which is a relevant anticipation, as we shall see.
The same approach performed with the BHandHLYP functional has several similarities. The signs of the signed currents strengths are all preserved, and thus, the tropicities continue to the follow Hückel rule. However, the magnitude of the currents are definitely smaller, even by an order of magnitude. The largest variations occur for c-P6 6+ and c-P8 8+ . Comparing the last columns of Tables 1 and 2, it is quite apparent that the choice of the functional is a matter of importance, as stated by Matito et al. [25].
With the change of the functional, the difference between the effective radius s av and the geometric radius s tends to increase for all hexacations, especially for c-P6 6+ , where s av is as large as 29.1 Å, more than twice the geometrical radius. This unphysical value definitely calls for the inadequacy of both the ICLOC and TCLOC models. Table 1. Estimation of the signed current strength of the nanorings from the magnetic properties computed at the LC-ωhPBE (ω = 0.1)/6-31G* level. Nuclear magnetic shielding and magnetizability components σ zz (0) and ξ zz in ppm and 10 −30 J T −2 , respectively. Belt radii s and s av in Å. Signed current strengths I B in nA T −1 . Among more advanced RCMs, we considered the reference ones for monocycles, ICLOC2 and ICLOC2C. Figure 1 compares the shielding computed for these models with that obtained with an ICLOC model and with that obtained from DFT. It is apparent that the best fit is obtained by the ICLOC2C model. This model also gives the signed current strength and magnetizability best matching the DFT values. Moreover, the best-fit parameters for the other two models are unphysical: the radius of the ICLOC is almost threetimes the geometrical radius 12.9 Å, and the off-plane displacement of the two identical loops of the ICLOC2 model by far exceeds the height of the porphyrin methines above the ring symmetry plane (Table 3). Eventually, also the individual parameters of the ICLOC2C model have a consistent deviation from the DFT values; however, these parameters are highly correlated, and different solutions are possible [19]. Fixing the values I 1 = 5.7 and I 1 = −7.7 nA T −1 , we still obtain a good fit of the scan, with a worsening of the computed magnetizability. The increased values of the radii s 1 and s 2 can be partly understood because what matters in a simplified RCM is an average radius, which as discussed before is typically greater than the geometric radius. Table 3. Parameters and derived quantities of the nonlinear RCMs used to fit the magnetic shielding scan of Figure 1. The DFT values are also reported: the radii correspond to the locations of the extrema in the plane used to compute the current strength ( Figure 3). The set of parameters of the ICLOC2C* model were obtained by fixing the current strengths at the values retrieved from DFT.

Current Strengths
To easily grasp the main features of the three-dimensional current density field, Sundholm and coworkers introduced net bond current strengths, or net bond current susceptibilities [10,12,46], which are obtained by the integration of the induced current densities crossing planes bisecting selected bonds. When current delocalization is present, a sizeable net bond current flow is detected, whose strength can be used to make a direct comparison among different molecular systems. To this end, it is customary to take the benzene C-C current strength as the reference and to compare the relative current strength. This procedure allows the comparison of the results coming from very different methods, even non-ab initio [47]. Absolute current strengths require extending the integration domain toward the onset of a circulation close to a different bond. This is best done using contour levels [13,48], but the procedure can become nonstraightforward for small cycles or congested molecules [49]. In these cases, the use of square [46] or circular [50] integration domains can be convenient. At any rate, the relative current strengths are less sensible to the tail of the currents. Figure 2 shows sketches of a cycloporphyrin unit of the rings. The center of the ring is below, and the magnetic field is directed upwards, so that arrows directed to the right/left indicate paratropic/diatropic currents. Units are percents of the benzene current (for this figure, squares of 4 au sides centered on the midpoints of the bonds were used; with the same method, the benzene signed current strength is −9.3 nA T −1 ). The continuity of the current is pretty well respected.
The higher values of current strength predicted for tetracations are here clearly appreciated. In four cases out of six, the current along the inner pathway passing through the nitrogen atoms is 2-3-times the current passing on the outer pathway, at odds with what happens in zinc-porphyrin for a magnetic field perpendicular to the porphyrin plane [51]. Exceptions are c-P8 6+ and c-P6 6+ . In the latter case, in particular, the current on the outer pathway even exceeds the one on the inner pathway, and both of them are very small, not reaching 10% of the reference benzene current. Tropicities along the inner pathway are consistent with the ansatz of a Hückel counting of 14×N − Q electrons for a c-PN Q species. Electron counting along the outer pathway would give the opposite prediction in the case of c-P6 6+ : 15×N − Q = 84. This Hückel counting argument turns out to be false as the small values of the total current strength along the outer perimeter are always diatropic. However, as anticipated by the ICLOC2C model, and as will be discussed below, the currents are actually the result of a competition of sizeable paratropic and diatropic currents.  The magnetic field is oriented upwards and the center of the nanoring is behind the paper, so that arrows directed towards the right denote a paratropic current and arrows directed towards the left denote a diatropic current.
Further insight into the different behavior of the tetracations and hexacations can be obtained considering in more detail the pattern of the current along a bond in the circuit. Figure 3 compares the contour levels of the current on the plane bisecting the middle bond of the butadienelink. The integration, carried up to 10 −3 au, is the source of the values reported in Tables 2 and 3. As can be seen in all systems (apart from two additional minor diatropic domains in c-P6 4+ ), there are two integration domains, indicating two concentric loops of current, which are corotating for tetracations and contrarotating for hexacations. There is considerable cancellation of the currents in the determination of the net current. Indeed, in the case of c-P6 6+ , the net current of −2.0 nA T −1 stems from the cancellation of two contributions as large as −7.7 and 5.7 nA T −1 . Less important cancellation happen in c-P7 6+ (+12.1 − 2.1 = +10.0) and c-P8 6+ (−11.8+1.7 = −10.1). These cancellations, especially that of c-P6 6+ , are not a consequence of a localized flow, but mainly stem from the presence of two contrarotating currents flowing inside and outside the ring (Figure 4), as was also found in ultrashort [5,5] carbon nanotubes [52] and much earlier in fullerene [53].   Tables 1 and 2. The current in the average molecular plane is also shown. A reference arrow with the size and direction of the maximum value of the delocalized benzene current is shown in blue. The magnetic field is oriented upwards and the center of the nanoring is behind the paper, so that arrows directed towards the right denote a paratropic current and arrows directed towards the left denote a diatropic current.
A B Figure 4. Current density plotted 1.2 au inside (A) and outside (B) the c-P6 6+ nanoring on a surface having the molecular shape. The magnetic field is directed upwards and the point of view is in the center of the nanoring, so that paratropic (diatropic) currents are directed towards the left (right).

Orbital Contributions to the Current
The ipsocentric method allows for a decomposition of the total magnetically induced current into a sum of rotationally allowed and translationally allowed orbital contributions [8,[54][55][56]. Remarkably, only a few contributions coming from orbitals close to the Fermi level are generally sufficient to recover the pattern of delocalized currents in an almost quantitative way. In paradigmatic simple cycles, the paratropic currents come from a Jahn-Teller split pair, while a full pair is responsible for the diatropic currents. These two-and four-electron rules often allow predicting current patterns at a low level of computation [57,58]. Table 4 reports this orbital decomposition for the cycloporphyrins. c-P7 6+ is the only system whose current is obtained semiquantitatively from the single HOMO; in all other cases, either four or six orbital contributions are needed for a semiquantitative model. Not only are c-P6 6+ and c-P8 6+ the only systems requiring up to six orbitals, they are also peculiar because their HOMO contribution has the opposite tropicity of the all-orbitals signed current strength. In effect, the selection rules lead to both rotationand translation-allowed virtual transition for orbitals close to the Fermi level ( Figure 5).   It is remarkable that, with the single exception of c-P7 6+ , the few-electron models are not built taking sequentially four or six orbitals down from the Fermi level, because of the intervention of orbitals that do not give any contribution to the current. These orbitals are linear combinations of porphyrin orbitals with a nodal plane in the average plane of the nanoring, and they cannot couple with the orbitals responsible for the current, which are all symmetric upon reflection in the molecular plane. In the case of c-P6 6+ , these orbitals are just the HOMO, HOMO−1, . . . HOMO−5, so that the first orbital giving a contribution to the current is the HOMO−6. Interestingly, several works [59,60] have shown that the orbital energies depend heavily on the amount of HF exchange incorporated into the functional, to the point that, for long-range corrected functionals, a fine-tuning of the ω parameter is required in order to obtain a reliable description of both ground and excited states [61], in particular with increasing system size [59,60]. Indeed, the cited paper by Matito et al. [25] showed that a slight variation of the ω parameter leads to completely different energy minima. In line with this analysis, we found that LC-ωhPHE (ω = 0.1) predicts a significantly different orbital distribution for the molecules under study, leading in turn to different contributions to the total current.

Off-Axis Shielding
One of the principal arguments to classify the c-PN Q as aromatic or antiaromatic is their NMR spectrum showing significant upfield and downfield shifts for protons attached to a pendant silyl group and facing either inside or outside the ring. The ICLOC model is a standard interpretative tool to cope with these shifts. Using the elliptic functions to compute the off-axis induced magnetic field [62], we computed the contribution to the isotropic NICS, stemming from a single ICLOC or from two concentric ICLOCs, taking the parameters obtained above for c-P6 6+ (Table 3, ICLOC2C* parameters). As can be seen from Figure 6, for a representative location of inner protons (e.g., z = R = 10 Å), the presence of two concentric contrarotating loops leads to shifts, which are 2-3-times those of a single ICLOC, even without changes of the total current strength. Notably, a two-dimensional BHandHLYP/6-31G* NICS scan gives results that compare better to the ICLOC2C* model than to the simple ICLOC. In the bottom-right panel of Figure 6, the DFT computation of the contribution to the isotropic NICS coming from a magnetic field perpendicular to a porphyrin unit is also reported. It can be seen that NICS zz and NICS xx have comparable magnitudes, so that they are both relevant to determine the isotropic chemical shift, observed in solution.

Discussion
The magnetic response used to define magnetic aromaticity is a nontrivial property due to its tensor character and to the fact that it can result from opposing contributions coming from different subsets of electrons. For planar molecules, it is customary to refer to the π electrons and to the "zz"-response, i.e., the induced field parallel to the inducing external field. This basic ansatz is no longer straightforward if the molecules are nonplanar, as the cycloporphyrins studied here. We showed that, as the porphyrin rings are roughly perpendicular to the average plane, a sizeable delocalized current can be expected for a non-zz perturbing magnetic field, and once more, this indicates that care must be used in interpreting the average shielding as a result of the currents delocalized in the average molecular plane. Although we did not attempt a full disentangling of the isotropic shielding, as this would require computations of the conformationally averaged alkyl pendants, we succeeded in modeling the DFT current strength in terms of the ICLOC2C model: two infinitely thin concentric circular loops of current. This modeling, needed for hexacations and not for tetracations requires an extension of the standard "trannulene" model [63] with two corotating concentric loops of currents [64]. The current is always generated by the "π in " orbitals, perpendicular to the nanoring perimeter, but due to the presence of virtual transitions of a different nature, the currents on the inside and on the outside of the molecule have opposite directions, as previously found in fullerene [53] and nanotubes [52]. This competition is not the result of the contributions of different subsets of active orbitals, such as the "π in " and "π out " in some carbomers [65,66] or the σ and π in periodoannulenes [67], although it must be recognized that the π in orbital set can be divided into those localized on the butadiene links and on the porphyrin units. The presence of only these two subsets of π in orbitals is not sufficient to require an extension of the basic trannulene model, as evidenced by the strong and corotating currents found in tetracations. As a matter of fact, mere electron counting to predict the current patterns according to the Hückel rule is rather effective in five times out of the six studied here, but it could lead to conjecture a highly overestimated magnetic response in the case of c-P6 6+ .
At any rate, proper credit should be given to the topology of Anderson's cycloporphyrins, which allows investigating the magnetic aromaticity of very large rings. This possibility is hindered in annulenes, due to steric interactions between hydrogen atoms pointing towards the middle of the cycle, as well as in polycyclics, where the dissection of a perimetral current in local and global contributions is far from obvious and has been rigorously performed only at the tight-binding level [68]. At that level, circuits give an energy contribution that decreases exponentially with their size [69], but their influence on magnetic properties decreases less rapidly, due to the increased area of the circuit [68,70]. Apart from the faint aromaticity of c-P6 6+ , the present investigation gave an independent confirmation of the onset of delocalized current over very large conjugated circuits, of up to 112 atoms.

Conclusions
We reported the computations of the current density and magnetic properties for the six nanorings c-PN Q , with N = 6, 7, 8 and Q = 4+ or 6+. The analysis of a few magnetic properties, such as the magnetizabilities and central shielding, was sufficient to determine semiquantitatively the magnetically induced current strengths, as demonstrated by the direct computation and visualization of the latter. We confirmed the finding of Matito's group, concerning the need for a proper functional to address delocalization in these systems [25]. Indeed, current strengths at the BHandHLYP level were very smaller to those computed at the LC-ωhPBE (ω = 0.1) level, even by a full order of magnitude. As for the debated case of c-P6 6+ , its current strength turned out to be diatropic, but very small, suggesting classifying it as very weakly aromatic. However, a proper RCM for c-P6 6+ was not the basic ICLOC, but the ICLOC2C, consisting of two concentric and coplanar infinitely thin circular loops of current, and the very small magnetic aromaticity of c-P6 6+ resulted from the cancellation of two non-negligible delocalized heterotrophic currents along the main conjugation pathway, which is also consistent with the interpretation derived from the ipsocentric method [55,71]. In this case, the fundamentally unsolved problem of retrieving the current strength from experimental isotropic nuclear shieldings not only has to face the dependence of these isotropic shieldings from sizeable perpendicular contributions (not the "zz" contribution considered for a basic ICLOC model) [72], but it should also take care of the unusual shape of the desired delocalized current (the one causing the parallel magnetic responses). GIAO Gauge including atomic orbital ICLOC Infinitely thin circular loop of current ICLOC2 2 off-plane displaced infinitely thin circular loops of current ICLOC2C 2 concentric coplanar infinitely thin circular loops of current NICS Nucleus-independent chemical shift RCM Ring current model TCLOC Toroidal circular loop of current