On the magnetic criteria of aromaticity in topologically spherical molecules: theoretical tools for the magnetic response of giant fullerenes C₂₄₀ and C₅₄₀

Abstract

Magnetically induced three-dimensional current densities, net bond current strengths and divergence of the isotropically averaged Lorentz force density have been used to characterize the magnetic response of giant I_h-fullerenes containing 240 and 540 carbon atoms. The weakly diamagnetic response of C₆₀ is also confirmed for these larger fullerenes, which do not fulfil the 2(n + 1)² rule for spherical aromaticity. The presence of peculiar paratropic ring currents local to 6-membered rings is demonstrated. Their appearance is due to the superposition of migrating diatropic ring currents around each paratropic loop.

---

1 Introduction

Several efforts have been made by chemists to synthesize highly symmetrical molecules with the ideal shape of the five Platonic solids,^1–5 previously described by Pythagoreans and Theaetetus (according to Euclid¹ and the Byzantine Encyclopedia Suidas²), who also gave the definition of “perfect numbers” as those obtained by summing their divisors, e.g. 6 = 1 + 2 + 3.⁶ These polyhedra provide tangible models to guide group theorists in the study of molecular symmetry and crystal structure.⁷ The Archimedean solids,⁸ valued by mathematicians for their high symmetry, can be derived from Plato's using a series of geometric operations, e.g., truncation, snubbing, and other symmetry-preserving transformations like expanding or rectifying.

The most celebrated of the Archimedean solids is possibly the truncated icosahedron of buckminsterfullerene, C₆₀, first obtained by Kroto et al.⁹ It corresponds to one of the “medial” Goldberg polyhedra; see Fig. 1 of ref. 10. The chemical properties of Archimedean solids have been widely studied, especially within the context of “spherical aromaticity”,^11–16 with the rule 2(n + 1)² for n π-electrons. The extension to open-shell systems has been considered by Poater and Solà:¹⁷ according to these authors, a compound is aromatic if it has 2n² + 2n + 1 π-electrons with a spin of S = (n + 1/2), indicating a half-filled last energy level. A further extension to inorganic compounds has also been reported.¹⁸

In 2023, an important article appeared, possibly a milestone that represents a break from some ideas commonly accepted within the chemical community, because it marks a turning point between earlier and later literature on aromaticity,¹⁹ a “suspicious concept of low reputation”.²⁰ In light of this paper, reporting the attempts by a group of thirteen chemists to build resilience against the conceptual spread of aromaticity and foretelling a change of paradigm¹⁹ within the spirit of a Kuhn revolution,²¹ the notion of spherical aromaticity, frequently considered to explain the magnetic properties of fullerenes, probably deserves to be reconsidered,¹⁹ together with a plethora of other aromaticity types. In fact, many “geometrical shapes” of aromaticity have been proposed, e.g., “cubic”²² and “octahedral”.²³ Also inorganic clusters may exhibit cubic aromaticity,²⁴ with a 6n + 2 electron counting law. The present authors are convinced that the notion of “aromaticity shape” may also be re-examined within the spirit of this recent paper.¹⁹ Magnetic criteria for spherical aromaticity, particularly effective in the case of giant fullerenes, have been developed in this study: they are discussed in Sections 2 and 4. Two large fullerenes of I_h symmetry containing 240 and 540 carbon atoms will be examined, together with fullerene-60 as a reference compound.

1.1 Some details on Platonic and Archimedean solids

In Plato's dialogue Timaeus,^25–27 five perfect bodies are discussed: the tetrahedron, octahedron, cube, icosahedron and dodecahedron. The corresponding figures by Leonardo da Vinci, appearing in the treatise by the mathematician Luca Pacioli, De Divina Proportione,²⁸ are the tetracedron, octocedron, exacedron, ycocedron and duodecedron. The first four correspond, respectively, to the elements fire, air, earth and water. The dodecahedron is not identified by any element: it is customarily associated with the heavenly matter, i.e., ether. The book by Pacioli, composed in the years around 1498 in Milan, and published in Venice in 1509 by Paganino Paganini, had been preceded by that of the Italian painter and mathematician Piero della Francesca, De quinque corporibus regularibus on the geometry of polyhedra,²⁹ written in the 1480s or early 1490s. Piero gives an account of the Platonic solids, together with descriptions of five among thirteen Archimedean convex polyhedra, and of several other irregular polyhedra.^29,30

Plato explains how the Demiurge forms these ideal bodies by employing elementary triangles of two types: “one being the isosceles, and the other that which always has the square on its greater side three times the square on the lesser side” (Timaeus 54 B 5).^25–27 Such a triangle, “the fairest – ϰλλιστον – of the triangles, that triangle out of which, when two are conjoined, the equilateral triangle – σπλευρον τργωνον – is constructed as a third” (Timaeus 54 A 5),^25–27 is obtained by taking one half of the equilateral triangle. Thus, it is the right scalene triangle in which the square of the longer leg† is three times the square of the shorter leg.

A simple example of the second type is given by that having the shorter leg with length 1, the other with length and hypotenuse of length 2. According to Plato, this is the most beautiful right scalene triangle, the other angles being 30° and 60°; see Timaeus 54 A 7. The triangular faces of three regular solids, tetrahedron, octahedron and icosahedron, are obtained by connecting six elementary triangles of the first shape to form the equilateral triangles of their faces: by joining them in pairs along the hypotenuse, the Demiurge obtains three quadrilaterals which, when further joined, form the equilateral triangle.

Four right isosceles triangles, connected through their right angles for each face, are used by the Demiurge to shape the cube, which produces the element earth (Timaeus 55 B 3–55 D 8). Eventually, Plato affirms that the Demiurge employed the dodecahedron to embellish the structure of the Cosmos, Timaeus 55 C 4–6, i.e., “used it up for the Universe in his decoration thereof”.^25–27

2 Magnetic criteria of aromaticity

Despite the ongoing debate about the meaning of the aromaticity concept, as described in the previous section, the majority of chemists agree on the fact that, in terms of electron delocalization, aromaticity shows up very clearly in NMR because it changes the magnetic environment of the nuclei. Delocalized electrons create their own current circuit that causes a strengthening or weakening of the external field at different points in space (anisotropic effect). Aromatic molecules are characterized by a diatropic ring current that reinforces the external magnetic field outside the ring and reduces it inside the ring; see [18]annulene for an extraordinary example.³¹ In this regard, the value of NMR doubles, because it also accounts for the definition and recognition of anti-aromatic compounds, such as those having a paratropic current that reduces/reinforces the external field outside/inside the ring; see bisdehydro[12]annulene as one of the first examples of this kind.³² Therefore, the main criterion of magnetic aromaticity/antiaromaticity consists of checking for the presence of a diatropic/paratropic ring current. While this can be deduced indirectly from experimental data, it can be determined directly through a quantum mechanical calculation.³³

For planar systems, the magnetic response is well addressed by maps of induced current densities^34–37 and bond current strengths (BCSs)^38,39 sustained by a magnetic field at right angles to the molecular plane. In particular, BCSs give a measure of delocalization of the induced current, which can be compared in different molecules. In non-planar cases, such as carbon nanobelts and carbon nanotubes, very useful information can be obtained by studying the current density induced by a magnetic field parallel to the cylindrical axis.^40–44 Magnetically induced current densities in non-planar polycyclic aromatic hydrocarbons have also recently been presented for the infinitene molecule,⁴⁵ some planar and non-planar clarenes⁴⁶ and carborane^47,48 in the context of 3D aromaticity.^49,50

In general, for non-planar systems, the investigation is less straightforward, because, apart from the systems cited above, there is no single direction of the inducing field that is optimal to compute all the currents. To overcome this difficulty, we consider here the isotropically averaged Lorentz force density (IALFD),⁵¹ together with current-density and BCS maps for some selected orientations of the external field. The IALFD is a vector field that accounts for the energy change of a randomly tumbling closed-shell molecule placed in a uniform magnetic field. Its divergence (DIAL) provides a local definition of diamagnetic and paramagnetic responses, e.g., diamagnetic when DIAL < 0 and paramagnetic when DIAL > 0. Both IALFD and DIAL are independent of the orientation of the magnetic field. Furthermore, they are gauge-origin invariant and independent of the point-of-view; see ref. 51 for details.

As regards very popular, field-independent quantities defined previously, such as NICS⁵² and ACID,^53,54 it should be noted that: NICS requires integration over the whole space producing a value that is biased by a geometric factor,⁵⁵ and, as it has been shown,⁵⁶ it is poorly correlated with the current strength; ACID does not distinguish between aromatic and antiaromatic contributions and neglects the antisymmetric component of the current-density tensor, which should be the most relevant term in the case of a highly diamagnetic aromatic molecule.⁵⁷ On the other hand, it should also be noted that more sophisticated versions of NICS do lead to good correlations with the current strength.⁵⁶

3 The giant fullerenes C₂₄₀ and C₅₄₀

As is well documented,^58,59 in a fullerene, each carbon atom is bonded to three others to form a closed cage that is topologically equivalent to a sphere. In particular, fullerene structures are trivalent pseudospherical polyhedra that contain only pentagonal and hexagonal faces and satisfy the Euler's polyhedron formula, V − E + F = 2. For a fullerene C_n, the number of vertices is V = n, the number of edges is E = 3n/2 and the number of faces is F = n/2 + 2. Since E must be an integer, odd numbers of vertices are precluded, as well as n = 22.⁵⁹ Working a little with these relations, it can be shown that fullerenes contain P = 12 pentagonal faces and H = n/2 − 10 hexagonal faces.

The number of distinct isomers of a fullerene with a given number of carbon atoms increases rapidly as n increases. For n = 60, there are 1812 nonisomorphic structures, of which only one fulfils the isolated pentagon rule (IPR),⁶⁰i.e., the well-known buckminsterfullerene of full icosahedral symmetry and stable closed-shell electronic structure. In general, considering the nuclear framework, icosahedral fullerenes have a number of vertices given by n = 20(i² + ij + j²), where i > 0, j ≥ 0 and i ≥ j are Coxeter's integer parameters.^58,59 When i = j or j = 0, the fullerene belongs to the full icosahedral point group I_h; otherwise, it belongs to the rotational subgroup I, which contains no reflections. The two series n = 60i² and n = 20i² for i = j and j = 0, respectively, are connected by the so-called “leapfrog” rule. Unlike 60i² fullerenes, several neutral I_h-fullerenes from the n = 20i² series do not have enough electrons to fill the degenerate HOMOs and undergo first-order Jahn–Teller distortion to lower symmetry.

It is interesting to note that I_h-fullerenes from the n = 60i² series with 240, 540,…, vertices have, between any two pentagonal faces, 1, 2,…, hexagons inserted, as illustrated in Fig. 1. These patterns are interesting, because they entail the insertion of 1, 2,…, crowns of benzene rings around each pentagon, which, for a given fullerene, cannot host all together a Clar sextet in any of the resonance structures; otherwise, an unpaired electron would result in an inner pentagon; see Fig. 2. This raises at least two questions: (i) could the greater stability of structures containing the largest number of Clar sextets lead to an open-shell ground state? (ii) How is the overall aromaticity of these giant fullerenes affected by the presence of such “less/more aromatic” rings in the closed/open-shell case? Regarding the first question, the stability of the wavefunction has been checked for both C₂₄₀ and C₅₄₀; see the following section for a discussion of the stability results. The second question will be addressed in Section 5 using methods we have developed so far based on the magnetically induced current density. Fullerene-60 is also included in the study as a reference to make useful comparisons. None of these neutral I_h-fullerenes fulfil the 2(n + 1)² Hirsch rule for spherical aromaticity.¹¹

Fig. 1 Repeating motifs in fully icosahedral C₆₀ (top), C₂₄₀ (middle), and C₅₄₀ (bottom). Two symmetry planes perpendicular to each figure intersect in the middle point (red cross) of each of them. Symmetry-unique carbon atoms are enumerated.

Fig. 2 Six-membered rings around pentagonal faces in C₂₄₀ on the left and C₅₄₀ on the right, which cannot be occupied all together by a Clar sextet in any of the resonance structures of the molecules. Only those around blue pentagons are shown.

4 Methods

In our approach, the magnetically induced current-density vector is expanded into a power serieswhereis the α, β component of the second-rank asymmetric current-density tensor (CDT). With , we denote the α component of first-order current density (MAGIC), which is a continuous three-dimensional function of the vector position r. For our purposes, eqn (1) is truncated to first order.

To ensure origin invariance, the CDT is calculated according to the continuous translation of the origin of the current density (CTOCD) method,⁶¹ adopting the DZ2 variant.^62,63 The isotropically averaged Lorentz force density (IALFD) induced by a unitary magnetic field is calculated using⁵¹

Thus, its divergence (DIAL) is
Bond current strengths (BCSs), or net bond current susceptibilities,³⁸ are calculated by integration of the J^B(r) cross section for points over planes that bisect a selection of bonds, i.e., only symmetry-unique bonds; see ref. 39 for details. Here, in particular, the cross section of the current density is restricted within circles of radius equals to , where D is the length of the bond under examination. Then, integration is carried out using a slightly modified Simpson method. Of course, other methods of integration can be used; see, for example, ref. 64. BCS values are conveniently reported as a percentage of the benzene ring-current strength (BRCS) of 12 nA T⁻¹, taken as the yardstick.³⁹

Current-density calculations have been performed using the SYSMOIC program package,^65,66 which requires first-order perturbed molecular orbitals with respect to angular and linear momenta. These are obtained using the G-16 package of computer programs,⁶⁷ by means of the NMR = CSGT option.⁶⁸

Molecular geometries were optimized within the Born–Oppenheimer (BO) approximation for the (closed-shell) singlet electronic state at the B3LYP/6-31G(1d) level of theory for C₆₀ and C₂₄₀ and B3LYP/6-31G for C₅₄₀, to reduce the computer workload, assuming the I_h symmetry point group. For the optimized structures, the absence of imaginary vibrational frequencies was verified. The stability of the electronic configurations was also verified by means of the standard technique provided by the Gaussian package. It turns out that the wavefunction of both fullerenes is stable under the perturbations considered in G-16.⁶⁷

Coupled perturbed Kohn–Sham calculations for first-order perturbed molecular orbitals were performed using G-16 and adopting the BHandHLYP/6-311+G(2d) combination of functional and basis set for all molecules. This choice was based on our own experience⁶⁹ and on a recent benchmark,⁷⁰ in which the accuracy of calculated magnetizabilities, over a set of 27 small molecules, was compared among 51 different density functionals. Indicators of local aromaticity, HOMA and MCBO (see below) have been computed using Multiwfn.^71,72

4.1 σ/π orbital separation

Within LCAO-MO methods,⁷³ the CDT is computed using a sum of orbital contributions. In the CTOCD approach, the CDT is also given by a sum of orbital contributions.⁶¹ Therefore, whenever the orbital separation can be performed according to a specific criterion, for example by symmetry, the dissection of the current density into orbital contributions provides an extremely useful means for the analysis of the results. For a discussion on the partition of the current density into orbital contributions, see Steiner and Fowler.^37,74

However, there are at least a couple of points that require comment. The first point involves σ/π separation in non-planar molecules, for which a mixing contamination is always present. Such a contamination can be minimised, as we did, by selecting as π orbitals only those that provide the least contribution to the electron density around the nuclei and at bond centres. A second point of concern is the σ/π mixing of the perturbed molecular orbitals due to the applied magnetic field. It mainly regards the IALFD, which takes into account all possible orientations of the external magnetic field by definition. Whereas the all-electron calculation of the DIAL, in the limit of a complete basis set, is independent of the method used to distribute the gauge, the π-electron contribution to the DIAL does change. At any rate, the possibility of obtaining useful information from method-dependent dissections is not precluded, which partly justifies the many papers reporting this kind of calculation. Here, we will consider the contribution to the DIAL provided by π-electrons only, using the continuous transformation of the origin of the current density (CTOCD-DZ2) method.⁶³ The results reported in previous papers, see, for example ref. 75, show that this dissection, to be called π-DIAL, is effective in giving information on the onset of local or semiglobal ring currents,^76,77 with no need to adjust the orientation of the external magnetic field.

In the case of a pseudo-spherical cage, such as an I_h fullerene, there is no strict symmetry-enforced distinction between α/π orbitals, but an approximate separation persists between those orbitals locally tangential and those perpendicular to the curved spherical surface passing through the nuclei. The doubly occupied molecular orbitals of fullerenes C₆₀, C₂₄₀, and C₅₄₀ can be initially divided into 60, 240, and 540 core orbitals, 90, 360, and 810 σ-orbitals, and 30, 120, and 270 π-orbitals. The symmetries spanned by the core orbitals and σ-orbitals can be easily determined. For complete lists, see the Supplementary Information. Then, taking the orbital symmetries resulting from the CPKS calculations at the BHandHLYP/6-311+G(2d) level, the symmetries of the π-orbitals can be obtained by difference as:

Γ_π(C₆₀) = 1A_g ⊕ 1G_g ⊕ 2H_g ⊕ 1T_1u ⊕ 1T_2u ⊕ 1G_u ⊕ 1H_u (3)

Γ_π(C₂₄₀) = 2A_g ⊕ 2T_1g ⊕ 2T_2g ⊕ 4G_g ⊕ 6H_g ⊕ 4T_1u ⊕ 4T_2u ⊕ 4G_u ⊕ 4H_u (4)

Γ_π(C₅₄₀) = 4A_g ⊕ 5T_1g ⊕ 5T_2g ⊕ 9G_g ⊕ 13H_g ⊕ 1A_u ⊕ 8T_1u ⊕ 8T_2u ⊕ 9G_u ⊕ 10H_u (5)

For each occupied orbital ϕ_i, the partial electron density 2ϕ_iϕ_i at the center of each pair of bonded atoms was calculated using the TIPOMO program provided by the SYSMOIC package,^65,66 storing in two separate lists: (I) the largest value across all atom pairs, and (II) the total sum over all pairs. Then, the two lists are sorted from lowest to highest value to get π-orbital candidates in the first positions, as long as their symmetries satisfies (3), (4) and (5) for C₆₀, C₂₄₀, and C₅₄₀, respectively. According to criterion (I), orbitals that show a large local σ/π mixing are downgraded, while criterion (II) rewards orbitals that show a low σ/π mixing over all the molecular space. For C₆₀ and C₂₄₀, we found that the first 30 and 120 orbitals, respectively, collected according to criterion (II) match exactly the symmetries in (3) and (4). For C₅₄₀, the first 270 orbitals collected according to criterion (II) match the symmetries in (5), except that H_g and H_u, which are 14 (one more) and 9 (one less), respectively. The lists prepared according to criterion (I) only match exactly for C₆₀. The disagreement increases with the size of the fullerene. Evidently, the orbitals that should be excluded by criterion (I) display a locally large σ/π mixing. On this basis, we have considered the π-orbitals provided by criterion (II), excluding one set of H_g in the case of C₅₄₀.

5 Results and discussion

Results are provided as a set of figures, one for each fullerene, showing the current density induced in the π-electron cloud by a perpendicular magnetic field, alongside all-electron BCS and π-DIAL plots.

As regards the current-density maps, due to the pseudo-spherical shape of the molecules, J^B(r) was calculated at points 0.4a₀ apart, chosen to form two surfaces parallel to the faces of the polyhedra, one 1a₀ outside the cage, another 1a₀ inside the cage. For each point, an arrow parallel to the vector is plotted whose area is proportional to |J^B|. The inducing magnetic field is perpendicular to the plotting plane and points toward the reader. Diatropic circulations are clockwise, while paratropic circulations are anticlockwise. In this way, a semi-quantitative picture of the induced vector field is obtained.

A quantitative visualization of the induced currents is provided by BCS maps, which consist of a set of oriented arrows drawn in magenta along the bond directions and flanked by a number representing the BCS value as a percentage of the BRCS, i.e., 100 stands for 12 nA T⁻¹. Each arrow area is proportional to the corresponding BCS value. Also, in this case, the inducing magnetic field is perpendicular to the plotting plane and points toward the reader.

π-DIAL scalar fields are quite useful when there is no optimal single direction of the inducing magnetic fields. Fullerenes represent exemplary good cases of application. Here, in the following, π-DIAL is represented by means of surfaces of different colours, which are displayed for two isovalues: −0.003 (blue) diamagnetic regions, and +0.003 (red) paramagnetic regions.

In order to give some credible support to the results presented here in the following, let us first compare the experimental data of magnetizability (ξ) and ¹³C chemical shift (δ_I) available for C₆₀ with the predictions that are obtained integrating the corresponding density functions, written in terms of the CDT as

Chemical shifts with respect to TMS are in general calculated using the isotropic component of the shielding tensor σ_I = 1/3Tr(σ^I) by means of

δ_I = σ_ref − σ_I + δ_ref (8)

where “ref” stands for a proper reference compound, for example, benzene. In our cases σ^I has only one distinct component by symmetry.

Calculated magnetizabilities and ¹³C chemical shifts are reported in Tables 1 and 2, respectively, together with the results also calculated for C₂₄₀ and C₅₄₀.

Table 1 Calculated magnetizabilities at the BHandHLYP/6-311+G(2d,1p) level of theory (in cgs ppm), obtained viaeqn (6); the conversion factor from ppm cm³ mol⁻¹ to 10⁻³⁰ J T⁻² is ⁷⁸

ξ	C60	C240	C540
Calc.	−276	−2583	−7946
Expt^79	−260 ± 20
Expt^80	−252

---

Table 2 Calculated ¹³C chemical shifts at the BHandHLYP/6-311+G(2d,1p) level of theory with respect to TMS in ppm, obtained viaeqn (7) and (8). Symmetry-unique carbon atoms are reported with labels as in Fig. 1. In parenthesis is the number of equivalent nuclei

δ	C60	#	C240	#	C540	#
C1	142.2	(60)	135.6	(60)	139.0	(60)
C2			127.5	(60)	129.3	(60)
C3			126.3	(120)	123.2	(120)
C4					120.3	(120)
C5					119.8	(60)
C6					118.7	(120)
Expt^81	142.68
Expt^82	143.2

---

As can be observed, the agreement between theoretical and experimental data for C₆₀ is excellent for both ξ and δ. The improvement with respect to the HF estimates of magnetizability (−357 cgs ppm) reported much earlier^83,84 is particularly evident. Comparing the magnetizability of C₆₀ with that of benzene, which is measured to be −54.8 cgs ppm,⁸⁵ it is clear that the diamagnetism of C₆₀ is quite small, thus leading to a conclusion that “C₆₀ seems to be an aromatic molecule with a vanishingly small π-electron ring-current magnetic susceptibility”.⁷⁹

Now, looking at the current-density maps in the top row of Fig. 3, it can actually be seen that the π-electron contribution to the magnetically induced current density (MICD) is diffusely diamagnetic outside the cage and predominantly paratropic inside, where it is mainly located in the vicinity of the pentagonal rings. This is in agreement with previous visualization; see, for example, ref. 86. Even the all-electron BCS map reveals only net paratropic ring currents, whose strength is nearly equal to 70% of the BRCS for this orientation of the magnetic field.

Fig. 3 Top row: π-Electron contribution to the current density induced in C₆₀ by a magnetic field perpendicular to the plotting plane and pointing towards the reader, on surfaces that are outside (top-left) and inside (top-right) the nuclear framework at 1a₀ from it. Diatropic/paratropic circulations are clockwise/anticlockwise. Bottom row: Left, all-electron contribution to bond current strength as a percentage of the benzene ring-current strength (BRCS) of 12 nA T⁻¹, taken as a yardstick;³⁹ right, the π-DIAL; surface values are ±0.003 a.u., red is positive (antiaromatic) and blue is negative (aromatic).

However, it should be observed that π-electron diatropic ring currents in the vicinity of hexagonal faces on both sides of the cage cannot be excluded at all. This is clearly summarized by the π-DIAL map in the bottom right of Fig. 3, where disconnected paramagnetic regions local to the 5-membered rings (red surface) are surrounded by a continuous diamagnetic region (blue surface). In particular, we underline that the sum of diatropic ring currents of equal strength on each hexagonal face provides the calculated zero BCS value along the edges shared by any two hexagons. A sum of pyracylene fragments⁶³ seems quite appropriate to describe the calculated pattern of induced current densities.

According to the few-electron model, only a small subset of the high-lying π electrons dominate the more complex patterns of current in polycyclic π systems.³⁷ Moreover, orbital contributions to the current density obey a symmetry-based selection rule.³⁶ In this regard, it is interesting to observe that the symmetries of the HOMO, HOMO−1, and LUMO are H_u, H_g, and T_1u, respectively. Since T_1g ∈ H_u ⊗ T_1u and T_1u ∈ H_g ⊗ T_1u, a pure paratropic contribution is predicted for the HOMO → LUMO virtual transition and a pure diatropic contribution is predicted for the HOMO−1 → LUMO one. These two should prevail over all other virtual transition contributions. This is nicely confirmed by the plots in Fig. 4, as the sum of the HOMO and HOMO−1 contributions to the MICD outside and inside the C₆₀ cage is practically indistinguishable from that in the upper row of Fig. 3. The few-electron model is another way to show the prevalence of the diatropic/paratropic circulations outside/inside the cage, which is compatible with the small diamagnetism of C₆₀.

Fig. 4 MICD contributions from frontier orbital virtual transitions outside (top row) and inside (bottom row) the C₆₀ cage: left, HOMO paratropic; centre, HOMO−1 diatropic; right, their sum.

In the larger fullerenes discussed in the following, the frontier orbital symmetries are exactly the same as in C₆₀. However, it is not easy to provide a simple illustration like that given in Fig. 4, due to the large number of virtual orbitals contributing to the MICD for both the HOMO and HOMO−1.

5.1 Fullerene-240

The current-density maps, BCS map and π-DIAL surfaces of C₂₄₀ are shown in Fig. 5. The calculated current-density map, induced in the π-electron cloud by a magnetic field perpendicular to the figure plane and pointing toward the reader, is shown on the top for two sets of points outside and inside the cage at 1a₀ from the nuclear framework. All-electron BCSs are reported on the bottom right. The contribution to the BCS from all but π-electrons is negligible. Therefore, the two kinds of maps provide complementary information on a qualitative and quantitative basis, respectively.

Fig. 5 Top row: π-Electron contribution to the current density induced in C₂₄₀ by a magnetic field perpendicular to the plotting plane and pointing towards the reader, on surfaces that are outside (top-left) and inside (top-right) the nuclear framework at 1a₀ from it. Diatropic/paratropic circulations are clockwise/anticlockwise. Bottom row: Left, all-electron contribution to bond current strength as a percentage of the benzene ring-current strength (BRCS) of 12 nA T⁻¹, taken as a yardstick;³⁹ right, the π-DIAL; surface values are ±0.003 a.u., red is positive (antiaromatic) and blue is negative (aromatic).

The MICD maps outside and inside the C₂₄₀ cage show more similarities than those in C₆₀. However, some remarkable differences can be observed. On the outside of the cage, at least two concentric semiglobal diatropic circulations can be seen, all surrounding a paratropic ring current in the central 6-membered ring and surrounded by a faint global diatropic loop on the equator. Inside the cage, nine well-formed paratropic ring currents can be observed, four in the vicinity of 5-membered rings on the vertices of a rhombus and the other five over 6-membered rings, four on the edges of the rhombus and one in its centre. The presence of local paratropic benzene ring currents is of particular interest and for certain represents the most striking feature observed in the MICD maps. For convenience of the discussion, in the following we will refer to these paratropic 6-membered rings as P6R.

Examining the BCS map at the bottom left of Fig. 5, we observe that when a bond is shared by two different circuits, for example, a semiglobal diatropic outside and a local paratropic inside the cage, the strength of the current increases, as in the case of the two lateral pentagonal faces in the center, where the strength of the paratropic current increases from 97–98% up to 140–144% of the BRCS. The same happens for the bonds shared by the innermost semiglobal diatropic circulation and the P6Rs on the edges of the rhombus, where the current strength reaches 105–106% of the BRCS. The outer semiglobal diatropic current reaches 118% of the BRCS when it share the bonds forming the pentagonal faces at the rhombus vertices. The strength of the global diatropic circulation is estimated to oscillate between 14% and 58% of the BRCS.

All this is wonderfully represented by the π-DIAL map in the bottom right of Fig. 5: red paratropic islands emerging from a blue diatropic ocean. The shape of the paratropic surfaces is different for the 5- and 6-membered rings, thus allowing easy distinction of one type from the other. A diatopic surface runs over the entire molecule following the perimeters of the rings. Does this picture represent a diamagnetic or paramagnetic molecule? How does it compare with C₆₀? Considering the calculated values shown in Table 1, it can be observed that the magnetizability of C₂₄₀ is predicted to be about 9 times higher than that of C₆₀. Since C₂₄₀ has a radius that is almost twice the radius of C₆₀ (6.95 vs. 3.55 Å) and it contains four times as many carbon atoms, it can be concluded that this fullerene is also slightly diamagnetic.

This is confirmed by the numerical value of the global current passing a plane containing a half fullerene, from the centre of the molecule up to the outer region. For an inducing magnetic field parallel to the C₂ symmetry axis, i.e., perpendicular to the fragments in Fig. 1, the integrated global current sustained by the π-electrons of C₂₄₀ is calculated to be −11.3 nA T⁻¹, to be compared with the value of −8.4 nA T⁻¹ determined for C₆₀ (the minus sign indicates that both currents are diamagnetic).

The observation of local antiaromatic paratropic currents at the benzene rings in a neutral closed-shell molecule represents something truly unusual that deserves further attention. Usually, the induced current in benzene rings should be diatropic as in an ideal Clar system, namely, a resonance structure with the largest number of disjoint aromatic π-sextets, i.e., benzene-like moieties.⁸⁷ However, as already mentioned before, placing a Clar sextet in each of the P6Rs would result in an unpaired electron inside the pentagonal rings. Quite likely, P6Rs are the result of many adjacent circulations, which add together to provide the enhanced paratropic loops.

This last consideration opens an interesting question that leads to an equally interesting conclusion, i.e., how many Clar structures can we devise and how many of them present a Clar's sextet in the P6R? To simplify the analysis, we consider only a piece of C₂₄₀ consisting of a carboncone with one pentagon on the cone cap and three intact circles of fused hexagonal rings (see Fig. 6), which corresponds to carboncone[1,3].⁸⁸ Starting from the corannulene fragment, two separated sextets can be placed at once, and can migrate to the adjacent rings in 5 different ways. For each of them, it is possible to position at most 7 additional Clar sextets in 5 different ways, resulting in a total of 5 × 5 = 25 Clar structures with 9 sextets each; see Fig. 6 for a scheme of one of them. Now, one can count how many times a sextet is found in the P6R. Such a count turns out to be zero, a result that is not easy to guess a priori. With the exception of the rings on the corners, for which the number of sextets is also 0, all other rings do have Clar sextets for some of the 25 Clar structures. This means that P6Rs experience diatropic currents induced in the rings around them. In other words, one can think of a migrating diatropic ring current around each P6R, which manifests itself through the fixed paratropic loops.

Fig. 6 A fragment of C₂₄₀, delimited by green bonds, corresponding to carboncone[1,3]. One of its 15 Clar structures is represented.

The calculated ¹³C NMR chemical shifts of C₂₄₀, reported in Table 2, lie within the aromatic region. Carbon nuclei forming the pentagonal faces (C1) are shifted downfield by ∼9 ppm with respect to the carbon nuclei of the 6-membered rings (C2 and C3). The latter display a very similar chemical environment judging from their calculated δs.

5.2 Fullerene-540

The current-density maps, BCS map and π-DIAL surfaces of C₅₄₀ are shown in Fig. 7. As in the previous case, the calculated current-density map, induced in the π-electron cloud by a magnetic field perpendicular to the figure plane and pointing toward the reader, is shown at the top for two sets of points outside and inside the cage at 1a₀ from the nuclear framework. All-electron BCSs are reported on the bottom right.

Fig. 7 Top row: π-Electron contribution to the current density induced in C₅₄₀ by a magnetic field perpendicular to the plotting plane and pointing towards the reader, on surfaces that are outside (top-left) and inside (top-right) the nuclear framework at 1a₀ from it. Diatropic/paratropic circulations are clockwise/anticlockwise. Bottom row: Left, all-electron contribution to bond current strength as a percentage of the benzene ring-current strength (BRCS) of 12 nA T⁻¹, taken as a yardstick;³⁹ right, the π-DIAL; surface values are ±0.003 a.u., red is positive (antiaromatic) and blue is negative (aromatic).

The MICD maps outside and inside the C₅₄₀ cage are much more similar than those in C₂₄₀. However, on the outside of the cage, at least three concentric semiglobal diatropic circulations can be seen, all surrounding a central small rhombus with four P6Rs in its vertices. A faint global diatropic loop on the equator can be seen. Inside the cage, some concentric semiglobal diatropic circulations are also visible, which include sixteen well-formed paratropic ring currents; four are in the vicinity of 5-membered rings on the vertices of a large rhombus, eight are P6Rs along the edges of the large rhombus, and four are P6Rs on the vertices of a small rhombus in the map centre. Then, the presence of numerous P6Rs is also confirmed in C₅₄₀.

The BCS map at the bottom left of Fig. 7 shows the larger currents located at 5-membered rings, where strengths reach 150% of the BRCS. For 6-membered rings, BCSs overcoming the BRCS can also be seen, particularly in the vicinity of the bonds shared by P6R and the semiglobal diatropic circulations; see, for example, the values of 127% and 115% of the BRCS along the larger rhombus edges. The strength of the global diatropic circulation is estimated to oscillate between 27% and 67% of the BRCS.

The virus-like shape of the π-DIAL map of C₅₄₀, at the bottom right of Fig. 7, is really impressive: paratropic regions emerge prominently from the underlying diatropic framework. The relatively small magnetizability of the molecule, i.e., −7946 cgs ppm, in comparison to those of C₆₀ and C₂₄₀, (see Table 1) also document that this large fullerene is very weakly diamagnetic. Nonetheless, the half-fullerene global current sustained by the π-electrons is calculated to be −47 nA T⁻¹, which is large compared to that of the intermediate fullerene-240, but it accounts for the bigger size of the cage and the greater number of π-electrons.

The calculated ¹³C NMR chemical shifts (see Table 2) are in the aromatic range between 139–119 ppm. Also, in this case the most deshielded nuclei (C1) are those of the pentagonal rings.

5.3 Local aromaticity

The hypothesis that a local paramagnetic current can result from a summation of nearby stronger local diamagnetic currents rests on the possibility of defining local patterns of circulation, and can be placed in the framework of discussions on local aromaticity.⁸⁹ Of the several tools available in the literature to address local aromaticity, we chose one of the oldest, the Polansky index, which is a measure of the similarity of the Hückel density matrix of a given six-membered ring and that of the isolated benzene. This index is limited in the [0,1] interval, and has been proved to be higher for rings endowed with Clar sextets in a larger fraction of resonance structures.⁹⁰ As can be seen from Fig. 8, the rings endowed with a paratropic ring current turn out to have a Polansky index smaller than all the six abutting rings, indicating that they have a smaller local aromaticity: the six abutting rings are more similar to benzene. The same result is obtained with two other indices: the HOMA (Harmonic Oscillator Model of Aromaticity), which is based on the similarity of bond lengths with those of benzene,⁹¹ and the Multi Center Bond Index,⁹² later called the Kekulé Six Center Index (KSCI)⁹³ or Multi Center Bond Order (MCBO),⁷² which is a generalization of the Mayer two-center bond order⁹⁴ to the six-center case. As can be seen from Table 3, both indicators give a lower local aromaticity than all abutting rings for the P6Rs, consistently with the Polansky indices.

Fig. 8 Polansky indices of the benzenoid rings of C₂₄₀ and C₅₄₀ and labels for symmetry-unique rings.

Table 3 Indicators of local aromaticity for the symmetry-unique rings of C₂₄₀ and C₅₄₀

Molecule	Ring	Polansky Index	HOMA	MCBO
C240	A	0.7710	0.5353	0.0120
	B	0.7754	0.5375	0.0183
	C	0.7364	0.2756	0.0099
C540	A	0.7676	0.6034	0.0132
	B	0.7685	0.7011	0.0188
	C	0.7479	0.2561	0.0098
	D	0.7657	0.5315	0.0151
	E	0.7519	0.3248	0.0098

---

6 Conclusions

The magnetic response of the giant I_h-fullerenes C₂₄₀ and C₅₄₀, for the (closed-shell) singlet electronic state, is not dissimilar to that of the iconic C₆₀. The study of the magnetically induced current density, net bond current strengths and divergence of the isotropically averaged Lorentz force density led us to the conclusion that both fullerenes, which do not fulfil the 2(n + 1)² rule for spherical aromaticity, are only weakly diamagnetic.

However, we show the presence, in each large fullerene, of an abundant number of peculiar strong paratropic ring currents local to 6-membered rings, which can be rationalized as a result of the superposition of migrating diatropic ring currents around each paratropic ring. Therefore, the local magnetic aromaticity of the rings can be different from what is observed in the total current-density pattern.

Consistently, indicators of local aromaticity based on the unperturbed density matrix or the geometry assign a lower local aromaticity to these rings.

Author contributions

All authors contributed equally to this work.

Conflicts of interest

There are no conflicts to declare.

Data availability

The code for calculating the magnetically induced current density, as well as bond current strengths and isotropically averaged Lorentz force densities, can be found at https://sysmoic.chem.unisa.it/DISTRIB/. Part of the data supporting this article has been included as part of the supplementary information (SI). See DOI: https://doi.org/10.1039/d6nj00094k.

Acknowledgements

Financial support from FARB 2022 and FARB 2023 is gratefully acknowledged.

Notes and references

1. Selections illustrating the history of Greek mathematics. 1: From Thales to Euclid, The Loeb Classical Library, ed. Heiberg, 2006. At p. 379, there is scholion i. of Eucl. Elements xiii: “Eν τοτ τ βιβλ, τουτστι τ τ, γρφεται τ λεγμενα Пλτωνος πντε σχματα, ατο[small upsilon, Greek, tilde] μν οϰ στιν, τρα δ τ[small omega, Greek, tilde]ν προειρημνων πντε σχημτων τ[small omega, Greek, tilde]ν Пυαγορεων στν τεϰβος ϰα πυραμ ϰα τ δωδεϰεδρον. Θεαιττου δ ττε ϰτεδρον ϰα τ εϰοσεδρον. τν δ προσωνυμαν λαβεν Пλτων δι τ μμνησαι ατν ν τ Tιμα περ ατ[small omega, Greek, tilde]ν”, with English translations by Ivor Thomas: “In this book, that is, the thirteenth, are described the five Platonic figures, which are however not his, three of the aforesaid five figues being due to the Pythagorean, namely, the cube, the pyramid and the dodecahedreon, while the octahedron and the icosahedreon are due to Theaetetus. They received the name Platonic because he discourses in the Timaeus about them.”.
2. Suidas, 2025, https://www.physics.ntua.gr/mourmouras/ilias_odysseia/Lexiko/souida.html, [Byzantine literary lexicon and encyclopedia dating back to the 10th century]. The Greek voice Theaetetus translates to: “Theaetetus, an Athenian, astronomer, philosopher, a pupil of Socrates, taught in Heraclea. He was the first to describe the five solids so-called. He lived after the Peloponnesian wars.”.
3. H. P. Schultz, J. Org. Chem., 1965, 30, 1361–1364.
4. H. Hopf, Classics in hydrocarbon chemistry: syntheses, concepts, perspectives, Wiley-VCH, Weinheim, Im Kolophon: Milton Keynes: Lightning Source, 2011.
5. N. Sherden and K. Allan, Platonic Solids, 2005, https://stoltz2.caltech.edu/seminars/2005_Allan_Sherden.pdf, [Online; accessed 23-August-2025].
6. F. Acerbi, Arch. History Exact Sci., 2005, 59, 319–348.
7. M. Hamermesh, Group theory and its application to physical problems, Dover publ, New York, Facsim, 1989.
8. Archimedean Solid, 2025, https://mathworld.wolfram.com/ArchimedeanSolid.html, [Online; accessed 23-August-2025].
9. H. W. Kroto, J. R. Heath, S. C. O’Brien, R. F. Curl and R. E. Smalley, Nature, 1985, 318, 162–163.
10. M. Goldberg, Tôhoku Math. J., 1934, 40, 226–236.
11. A. Hirsch, Z. Chen and H. Jiao, Angew. Chem., 2000, 39, 3915–3917.
12. M. Bühl and A. Hirsch, Chem. Rev., 2001, 101, 1153–1184.
13. Z. Chen and R. B. King, Chem. Rev., 2005, 105, 3613–3642.
14. A. Muñoz Castro, Int. J. Quantum Chem., 2018, 118, e25645.
15. C.-C. Shu, D. W. Szczepanik, A. Muñoz-Castro, M. Solà and Z.-M. Sun, J. Am. Chem. Soc., 2024, 146, 14166–14173.
16. A. Muñoz-Castro, Chem. – Asian J., 2024, 19, e202400670.
17. J. Poater and M. Solà, Chem. Commun., 2011, 47, 11647.
18. M. Kulichenko, N. Fedik, N. V. Tkachenko, A. Muñoz Castro, Z.-M. Sun and A. I. Boldyrev, Aromaticity, Elsevier, 2021, pp. 447–489.
19. G. Merino, M. Solà, I. Fernández, C. Foroutan-Nejad, P. Lazzeretti, G. Frenking, H. L. Anderson, D. Sundholm, F. P. Cossío, M. A. Petrukhina, J. Wu, J. I. Wu and A. Restrepo, Chem. Sci., 2023, 14, 5569–5576.
20. M. Solà, Front. Chem., 2017, 5, 22.
21. T. S. Kuhn and I. Hacking, The structure of scientific revolutions, The University of Chicago Press, Chicago, London, 4th edn, 2012.
22. O. López-Estrada, E. Selenius, B. Zuniga-Gutierrez, S. Malola and H. Häkkinen, J. Chem. Phys., 2021, 154, 204303.
23. O. El Bakouri, M. Duran, J. Poater, F. Feixas and M. Solà, Phys. Chem. Chem. Phys., 2016, 18, 11700–11706.
24. P. Cui, H.-S. Hu, B. Zhao, J. T. Miller, P. Cheng and J. Li, Nat. Commun., 2015, 6, 6331.
25. Plato, Timaeus, Critias, Cleitophon, Menexenus, Epistles, Translated by R. G. Bury, Loeb Classical Library 234, Cambridge, MA, Harvard University Press, 1929. Internet Archive, https://archive.org/details/bwb_S0-EHE-311/page/n9/mode/2up.
26. Plato, Timaeus, Translated by Donald J. Zeyl, Indianapolis, Hackett Publishing Company, 2000.
27. Plato, Timaeus, Translated by Peter Kalkavage, Newburyport, MA, Focus Publishing, 2nd edn, 2016, https://dn710003.ca.archive.org/0/items/plato-keith-whitaker-parmenides-focus-1996/Plato%20-%20Timaeus%20%5Bkalkavage%5D.pdf.
28. L. Pacioli and D. D. Proportione, https://https://archiviostorico.mediobanca.com/wp-content/uploads/2021/01/De-Divina-Proportione.pdf, 1956, [Online; accessed 18-August-2025].
29. Wikipedia contributors, De quinque corporibus regularibus – Wikipedia, The Free Encyclopedia, 2024, https://en.wikipedia.org/w/index.php?title=De_quinque_corporibus_regularibus&oldid=1217884125.
30. J. V. Field, Arch. Hist. Exact Sci., 1997, 50, 241–289.
31. F. Sondheimer, Pure Appl. Chem., 1963, 7, 363–388.
32. J. A. Pople and K. G. Untch, J. Am. Chem. Soc., 1966, 88, 4811–4815.
33. P. Lazzeretti, Prog. Nucl. Magn. Reson. Spectrosc., 2000, 36, 1–88.
34. P. W. Fowler, E. Steiner, B. Cadioli and R. Zanasi, J. Phys. Chem. A, 1998, 102, 7297–7302.
35. A. Ligabue, U. Pincelli, P. Lazzeretti and R. Zanasi, J. Am. Chem. Soc., 1999, 121, 5513–5518.
36. E. Steiner and P. W. Fowler, Chem. Commun., 2001, 2220–2221.
37. E. Steiner and P. W. Fowler, J. Phys. Chem. A, 2001, 105, 9553–9562.
38. J. Jusélius, D. Sundholm and J. Gauss, J. Chem. Phys., 2004, 121, 3952.
39. G. Monaco, R. Zanasi, S. Pelloni and P. Lazzeretti, J. Chem. Theory Comput., 2010, 6, 3343–3351.
40. G. Monaco and R. Zanasi, J. Phys. Chem. Lett., 2017, 8, 4673–4678.
41. G. Monaco, E. L. T. Scott and R. Zanasi, ChemistryOpen, 2020, 9, 616–622.
42. F. F. Summa, G. Monaco, L. T. Scott and R. Zanasi, J. Phys. Chem. Lett., 2020, 11, 7489–7494.
43. H. Kono, Y. Li, R. Zanasi, G. Monaco, F. F. Summa, L. T. Scott, A. Yagi and K. Itami, J. Am. Chem. Soc., 2023, 145, 8939–8946.
44. N. Kai, H. Kono, T. Stünkel, D. Imoto, R. Zanasi, G. Monaco, F. F. Summa, L. T. Scott, A. Yagi and K. Itami, Chem. Sci., 2025, 16, 9195–9202.
45. G. Monaco, R. Zanasi and F. F. Summa, J. Phys. Chem. A, 2022, 126, 3717–3723.
46. G. Monaco, F. F. Summa, R. Zanasi and P. Lazzeretti, Chem. – Eur. J., 2024, 30, e202401167.
47. Z. Noori, M. Solà, C. Viñas, F. Teixidor and J. Poater, Beilstein J. Org. Chem., 2025, 21, 412–420.
48. Z. Noori and J. Poater, Org. Chem. Front., 2025, 12, 4321–4331.
49. J. Aihara, J. Am. Chem. Soc., 1978, 100, 3339–3342.
50. R. B. King, Chem. Rev., 2001, 101, 1119–1152.
51. G. Monaco and R. Zanasi, J. Chem. Phys., 2020, 153, 104114.
52. P. V. R. Schleyer, C. Maerker, A. Dransfeld, H. Jiao and N. J. R. V. E. Hommes, J. Am. Chem. Soc., 1996, 118, 6317–6318.
53. R. Herges and D. Geuenich, J. Phys. Chem. A, 2001, 105, 3214–3220.
54. D. Geuenich, K. Hess, F. Köhler and R. Herges, Chem. Rev., 2005, 105, 3758–3772.
55. S. Pelloni, G. Monaco, R. Zanasi and P. Lazzeretti, AIP Conf. Proc., 2012, 1456, 114–118.
56. S. Radenković and S. Đorđević, Phys. Chem. Chem. Phys., 2021, 23, 11240–11250.
57. G. Monaco and R. Zanasi, J. Phys. Chem. A, 2018, 122, 4681–4686.
58. P. W. Fowler, Philos. Trans. R. Soc., A, 1993, 343, 39–52.
59. P. W. Fowler and D. E. Manolopoulos, An atlas of fullerenes, Clarendon Press, Oxford [u.a.], 1995.
60. R. Sure, A. Hansen, P. Schwerdtfeger and S. Grimme, Phys. Chem. Chem. Phys., 2017, 19, 14296–14305.
61. P. Lazzeretti, M. Malagoli and R. Zanasi, Chem. Phys. Lett., 1994, 220, 299–304.
62. R. Zanasi, J. Chem. Phys., 1996, 105, 1460.
63. P. Fowler, R. Zanasi, B. Cadioli and E. Steiner, Chem. Phys. Lett., 1996, 251, 132–140.
64. T. J. P. Irons, L. Spence, G. David, B. T. Speake, T. Helgaker and A. M. Teale, J. Phys. Chem. A, 2020, 124, 1321–1333.
65. G. Monaco, F. F. Summa and R. Zanasi, J. Chem. Inf. Model., 2020, 61, 270–283.
66. G. Monaco, F. F. Summa and R. Zanasi, SYSMOIC Version 3.3, 2022, on-line manual: https://SYSMOIC.chem.unisa.it/MANUAL/.
67. M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, G. Scalmani, V. Barone, G. A. Petersson, H. Nakatsuji, X. Li, M. Caricato, A. V. Marenich, J. Bloino, B. G. Janesko, R. Gomperts, B. Mennucci, H. P. Hratchian, J. V. Ortiz, A. F. Izmaylov, J. L. Sonnenberg, D. Williams-Young, F. Ding, F. Lipparini, F. Egidi, J. Goings, B. Peng, A. Petrone, T. Henderson, D. Ranasinghe, V. G. Zakrzewski, J. Gao, N. Rega, G. Zheng, W. Liang, M. Hada, M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, T. Vreven, K. Throssell, J. A. Montgomery, Jr., J. E. Peralta, F. Ogliaro, M. J. Bearpark, J. J. Heyd, E. N. Brothers, K. N. Kudin, V. N. Staroverov, T. A. Keith, R. Kobayashi, J. Normand, K. Raghavachari, A. P. Rendell, J. C. Burant, S. S. Iyengar, J. Tomasi, M. Cossi, J. M. Millam, M. Klene, C. Adamo, R. Cammi, J. W. Ochterski, R. L. Martin, K. Morokuma, O. Farkas, J. B. Foresman and D. J. Fox, Gaussian16 Revision C.01, Gaussian Inc., Wallingford CT, 2016.
68. J. R. Cheeseman, G. W. Trucks, T. A. Keith and M. J. Frisch, J. Chem. Phys., 1996, 104, 5497–5509.
69. F. F. Summa, G. Monaco, P. Lazzeretti and R. Zanasi, Phys. Chem. Chem. Phys., 2021, 23, 15268–15274.
70. S. Lehtola, M. Dimitrova, H. Fliegl and D. Sundholm, J. Chem. Theory Comput., 2021, 17, 1457–1468.
71. T. Lu and F. Chen, J. Comput. Chem., 2012, 33, 580–592.
72. T. Lu, J. Chem. Phys., 2024, 161, 082503.
73. R. M. Stevens, R. M. Pitzer and W. N. Lipscomb, J. Chem. Phys., 1963, 38, 550–560.
74. E. Steiner and P. W. Fowler, Phys. Chem. Chem. Phys., 2004, 6, 261–272.
75. G. Monaco, F. F. Summa, R. Zanasi and P. Lazzeretti, Chem. – Eur. J., 2024, 30, e202401167.
76. D. Sundholm, H. Fliegl and R. J. Berger, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2016, 6, 639–678.
77. A. Stanger, G. Monaco and R. Zanasi, ChemPhysChem, 2019, 21, 65–82.
78. P. J. Mohr, D. B. Newell, B. N. Taylor and E. Tiesinga, Rev. Mod. Phys., 2025, 97, 025002.
79. R. C. Haddon, L. F. Schneemeyer, J. V. Waszczak, S. H. Glarum, R. Tycko, G. Dabbagh, A. R. Kortan, A. J. Muller, A. M. Mujsce, M. J. Rosseinsky, S. M. Zahurak, A. V. Makhija, F. A. Thiel, K. Raghavachari, E. Cockayne and V. Elser, Nature, 1991, 350, 46–47.
80. R. S. Ruoff, D. Beach, J. Cuomo, T. McGuire, R. L. Whetten and F. Diederich, J. Phys. Chem., 1991, 95, 3457–3459.
81. R. Taylor, J. P. Hare, A. K. Abdul-Sada and H. W. Kroto, J. Chem. Soc., Chem. Commun., 1990, 1423.
82. H. Ajie, M. M. Alvarez, S. J. Anz, R. D. Beck, F. Diederich, K. Fostiropoulos, D. R. Huffman, W. Kraetschmer and Y. Rubin, et al. , J. Phys. Chem., 1990, 94, 8630–8633.
83. P. Fowler, P. Lazzeretti and R. Zanasi, Chem. Phys. Lett., 1990, 165, 79–86.
84. R. Zanasi and P. Fowler, Chem. Phys. Lett., 1995, 238, 270–280.
85. J. Hoarau, N. Lumbroso and A. Pacault, C. R. Acad. Sci., 1956, 242, 1702.
86. M. P. Johansson, J. Jusélius and D. Sundholm, Angew. Chem., Int. Ed., 2005, 44, 1843–1846.
87. M. Solà, Front. Chem., 2013, 1, 22.
88. Z.-Z. Zhu, Z.-C. Chen, Y.-R. Yao, C.-H. Cui, S.-H. Li, X.-J. Zhao, Q. Zhang, H.-R. Tian, P.-Y. Xu, F.-F. Xie, X.-M. Xie, Y.-Z. Tan, S.-L. Deng, J. M. Quimby, L. T. Scott, S.-Y. Xie, R.-B. Huang and L.-S. Zheng, Sci. Adv., 2019, 5, eaaw0982.
89. M. Bühl, Chem. – Eur. J., 1998, 4, 734–739.
90. O. E. Polansky and G. Derflinger, Int. J. Quantum Chem., 1967, 1, 379–401.
91. T. M. Krygowski, J. Chem. Inf. Comput. Sci., 1993, 33, 70–78.
92. M. Giambiagi, M. S. De Giambiagi and K. C. Mundim, Struct. Chem., 1990, 1, 423–427.
93. P. Bultinck, R. Ponec and S. Van Damme, J. Phys. Org. Chem., 2005, 18, 706–718.
94. I. Mayer, Chem. Phys. Lett., 1983, 97, 270–274.

---

Footnote

† The Greek word πλευρν, usually translated with “side” when speaking of triangles, should be “leg” in this case, as evident from the context.

---