Orbital contributions to magnetically induced current densities using gauge-including atomic orbitals

Abstract

We have developed a method to calculate orbital contributions to magnetically induced current density (MICD) susceptibilities in molecules using gauge-including atomic orbitals (GIAO). The methods implemented in the GIMIC program have been used for analyzing orbital contributions to magnetically induced ring-current (MIRC) strengths. We have studied five aromatic, one nonaromatic, and four antiaromatic molecules. We show here that the contributions to the MIRC strength of all orbitals belonging to a given irreducible representation of the molecular point group in the presence of an external magnetic field are divergence free, whereas the MICD susceptibility of the individual orbitals are generally not divergence free. The largest contribution to the MIRC strength of antiaromatic molecules originates from the transition between the highest occupied molecular orbital (HOMO) and the lowest unoccupied molecular orbital (LUMO), whereas aromatic molecules have significant contributions involving many occupied orbitals. The MIRC contributions of σ orbitals are significant for planar molecules with a strained molecular structure.

---

1 Introduction

The aromaticity concept was originally introduced to characterize molecules whose electron-delocalization properties resemble those of benzene.^1,2 The electron delocalization of aromatic molecules leads to enhanced thermodynamic stability, reduced bond-length alternation, particular reactivity, and typical magnetic and spectroscopic properties.³ The aromatic nature according to the magnetic criterion can be determined by measuring or calculating proton nuclear magnetic resonance (¹H NMR) chemical shifts or by calculating various magnetic response properties such as nucleus-independent chemical shifts (NICS),⁴ magnetically induced current density (MICD) susceptibilities,^5–11 and magnetically induced ring currents (MIRC).^11–15

These computational approaches have to some extent been used for estimating contributions to the aromatic nature of individual molecular orbitals.^16–23 A general and reliable method for calculating orbital contributions to MICD susceptibilities and molecular properties would be useful, since decomposition into orbital contributions provides insights into the electronic structure and can be used for predicting changes due to substitution, structural variations, or charge redistribution. Common approaches to qualitatively understand the role of the frontier molecular orbitals (MO) in aromatic bonds comprise analyses of orbital energies and spatial distributions of orbitals as well as using symmetry selection rules.^24,25

Symmetry-based selection rules of the transitions between occupied and virtual orbitals have been used for assessing which occupied orbitals contribute to the MICD susceptibilities and magnetic properties. When the electronic transition between the highest-occupied molecular orbital (HOMO) and the lowest-unoccupied molecular orbital (LUMO) is electric-dipole (translationally) allowed, ring-shaped molecules may sustain a diatropic MIRC in the presence of an external magnetic field. Aromatic molecules can be identified using group theory. Molecules can be aromatic when the product of the irreducible representations of the HOMO, the LUMO and the translation operator contains the total symmetric irreducible representation. Analogously, molecules can be antiaromatic when the irreducible representation of the rotational operator multiplied by the irreducible representations of the HOMO and LUMO contains the total symmetric irreducible representation. The transition is then magnetically dipole (rotationally) allowed, which eventually leads to a paratropic MIRC.^26–28 Nonaromatic molecules are in this context conjugated molecular ring that are neither aromatic nor antiaromatic. Even though a transition is formally translation or rotationally allowed, a molecule may still be nonaromatic due to a large spatial separation of the occupied and virtual orbitals or because the energy separation between them is large.

The symmetry analysis provides a qualitative interpretation of molecular aromatic properties, whereas the MIRC strengths provide quantitative information about the aromatic nature because molecules sustaining a net diatropic MIRC are aromatic, whereas antiaromatic molecules have a net paratropic MIRC.²⁹ The MIRC strength can be used for determining how strongly aromatic or antiaromatic the molecule is.¹¹

Orbital contributions to the MICD susceptibility have previously been calculated using the continuous transformation of the origin of the current density (ipsocentric) approach.^7,10,30 In the ipsocentric approach, the origin of the vector potential of the external magnetic field is the position where the MICD susceptibility is calculated. Steiner and Fowler performed the first calculations of orbital contributions to the MICD susceptibility using the few-orbital model.^24,31 They studied the MICD susceptibility of benzene with the ipsocentric approach using an early version of the SYSMOIC package.¹⁴ The ipsocentric approach employs atomic basis functions that do not depend on the magnetic perturbation, whereas gauge-including atomic orbitals (GIAO) also known as London atomic orbitals incorporate the gauge origin in the basis functions^32,33 leading to a fast basis-set convergence of the MICD susceptibility and NMR shielding constants.^34,35

The gauge-including magnetically induced currents (GIMIC) method employs GIAOs.^{11–13,15,28} GIMIC calculations need the one-electron density matrix and the three magnetically perturbed density matrices as input. The density matrices can be calculated at various levels of theory using common quantum chemistry program packages. Cartesian coordinates of the molecular structure and basis-set information are also needed as input data. The GIMIC method has been successfully applied in MICD susceptibility studies of aromatic and antiaromatic molecules including planar molecular rings,²⁹ Möbius-twisted molecular rings,^36,37 toroidal molecules,³⁸ organometallic molecules,^39–41 and inorganic molecular clusters.^42,43

MICD susceptibilities calculated with GIMIC suffer from a tiny charge leakage, i.e., the divergence of the MICD susceptibility does not vanish when using finite basis sets. Divergence-free MIRC strengths can be calculated by performing numerical line integration of the induced magnetic field along the stagnation line of the MIRC vortex using the Ampère-Maxwell law.⁴⁴ However, it has been shown that the charge-leakage error of GIMIC calculations is negligibly small already when using commonly used moderate-size basis sets and that it practically vanishes when using large basis sets.¹¹ Recently, an expression for correcting the divergence was derived enabling calculations of divergence-free MICD susceptibilities using any basis set.⁴⁵

NMR shielding tensors have been partitioned into contributions of various orbitals by performing natural chemical shielding (NCS) analyses.⁴⁶ However, unitary transformations of the orbitals lead to spurious magnetic shielding contributions,²⁵ which can be identified as unphysical core and valence-σ contributions to the MICD susceptibility and to the NMR shielding tensors due to the mixing of occupied orbitals belonging to the same irreducible representation.^43,47 Orbital contributions to the MIRC strengths can also be estimated by partitioning the MICD susceptibility.^25,48 Orbital contributions and atomic orbital contributions to NMR shielding tensors can then be calculated by numerical integration.^49,50

In this work, we report an extension of the GIMIC method that enables calculations of orbital contributions to the MICD susceptibility. This is achieved by transforming the density matrices in the basis-function (atomic) basis to the canonical orbital basis where molecular orbital contributions can be identified. The selected orbital contributions to the density matrices are transformed back to the basis-function basis. The orbital contributions to the MICD susceptibility can be analyzed and visualized in the same way as the total MICD susceptibility.

We have studied eleven planar molecules and one nearly planar monocyclic molecule with different types of aromatic character. Six of them are aromatic (benzene (D_6h), borazine (D_3h), 1,5-dibora-2,4-diazabenzene (C_2v), porphin (D_2h), cyclopropane (D_3h) and the cyclopropenium cation (D_3h)), one is nonaromatic (1,4-cyclohexadiene, (D_2h)), and four are antiaromatic (cyclooctatetraene (D_4h and D_2d), hexadehydro[12]annulene (D_3h), tetraoxa-isophlorin (D_2h) and cyclobutadiene (D_2h)), where the molecular point group is given within parenthesis.

Cyclobutadiene and the cyclopropenium cation were studied because they have bond angles that significantly deviate from the standard bond angles of sp² hybridized hydrocarbons leading to ring strains. These strained molecules, along with molecules like cyclopropane, cyclobutane, and cubane, sustain an MIRC in the valence σ orbitals leading to what is called σ aromaticity or σ antiaromaticity. The direction (tropicity) of the MIRC is determined by the number of electrons in the σ orbitals of the C–C bonds.^29,51–59

The results for some of the molecules are compared to available data calculated using the ipsocentric approach. Orbital contributions to the MICD susceptibility of most of the studied molecules are reported here for the first time to the best of our knowledge.

The method we use for calculating orbital contributions to the MICD susceptibility is described in Section 2. We describe the employed computational levels in the Section 3. The obtained results are presented in Section 4. The main results are summarized in Section 5, which also contains the most important conclusions of our study.

2 Methodology

2.1 Decomposition of density matrices

The one-electron density matrix in the atomic-orbital basis (D^AO) consist of contributions of the individual occupied orbitalswhere N_occ is the number of occupied orbitals and n_i is the occupation number of MO i, and C_i contains the coefficients of MO i. The orbital contributions to D^AO (D^AO_i) can be identified by transforming it to the canonical-orbital basis D^MO. The basis transformation is done by multiplying from right with SC, where S is the overlap matrix and C is matrix containing the MO coefficients, and with its transpose C^TS from left

D^MO = C^TSD^AOSC. (2)

The first-order magnetically perturbed density matrices in the AO basis for each Cartesian direction of the magnetic field (B_λ, λ ∈ {x, y, z}) can be obtained from the MO coefficients and the first-order change of the MO coefficients (C_λ) due to the external perturbation B_λ

C_λ can be expressed as a unitary (U_λ) transformation of the unperturbed occupied MO coefficientswhich is the orbital response due to the external magnetic field B_λ in the limit of vanishing perturbation.^60–62 The occupied–virtual block of U_λ (U_λ,iv) is obtained by solving the orbital response equations. The virtual–virtual block of U_λ does not contribute to the perturbed density matrices.⁶³ The occupied–occupied block of U_λ (U_λ,ij) does not vanish when perturbation-dependent GIAOs are used.⁶⁴ However, it can be determined i.e., by differentiating the orthonormalization condition^60,61

Using the relation in eqn (4), the final expression for the unitary transformation matrix can be written as⁶⁵

where S^B_λ is the first derivative of the overlap matrix with respect to the external magnetic field (B_λ) in the limit of vanishing perturbation. The contributions of orbital i to the perturbed density matrices can also be obtained by performing the same unitary transformation and orbital selection procedure as applied on the one-electron density matrix. The orbital contribution to the perturbed density matrices in the MO basis consists of one row and one column of the matrices. The common elements of two orbitals in the occupied–occupied block are shared equally between the two orbitals.

The orbital contributions to the density matrices are transformed to the AO basis. The dimension of the AO density matrices is N_basis², where N_basis is the number of basis functions. The orbital contributions to the MICD susceptibility can then be calculated with GIMIC in the same way as the total MICD susceptibility by using the corresponding orbital contributions to the density matrices D^AO_λ,i and D^AO_i as input. Contributions to the MIRC strength of the i^th MO are obtained by modifying the occupation matrix n_i and performing the linear transformations and thereby avoiding time-consuming calculations.

The computational costs for calculating orbital contributions to the MICD susceptibility are comparable to those of all-electron GIMIC calculations because the computational time for preparing the density matrices in the AO basis for selected orbitals is much smaller than for the GIMIC calculation, whose computational time is independent of whether it is the total MICD susceptibility or orbital contributions to it.

The perturbed density matrices contain information about the response of all orbitals through the mixing of occupied and virtual orbitals, whereas in calculations of orbital contributions using the few-orbital model, only selected virtual orbitals are considered.^24,25,31 The selection might lead to wrong conclusions because many small contributions from pairs of occupied and virtual orbitals may affect the interpretation. We can perform calculations using the few-orbital model by identifying the orbital response involving selected occupied and virtual orbitals, whereas the rest of the response terms are set to zero. We have not implemented the approximate few-orbital model because the computational costs are the same as when considering the response of all orbitals.

2.2 Group theoretical analysis

The external magnetic field reduces the symmetry of the molecular point group when symmetry operations of the molecular point group do not map the magnetic field onto itself.^66,67 The resulting molecular point group in the presence of the external magnetic field is a subgroup of the molecular point group in the absence of the field and of C_∞h,⁶⁸ which is the point group of the external magnetic field. The order of the subgroup is usually larger than one, when the magnetic field coincides with a symmetry axis of the molecular point group.⁶⁶ The point group of the molecular structure in the presence of an external magnetic field can conveniently be determined by using the flow chart reported in ref. 66.

Orbital contributions to the MICD susceptibility are not in general divergence free, whereas the orbital contributions of all orbitals in a given irreducible representation of the combined point group are divergence free. Thus, the sum of the orbital contributions to the MIRC of all orbitals in each irreducible representation is then charge conserving, which leads to an MIRC whose strength is independent of the orientation of the integration plane.

The point groups of the studied molecules in the presence of an external magnetic field along the main symmetry axis, which is perpendicular to the molecular ring, are C_6h (benzene), C_4h (planar cyclooctatetraene), S₄ (bent cyclooctatetraene), C_3h (borazine, cyclopropenium cation, and hexadehydro[12]annulene) , C_2h (porphin, 1,4-cyclohexadiene, tetraoxa-isophlorin and cyclobutadiene), and C_s (1,5-dibora-2,4-diazabenzene).

2.3 Spatial averaging of the orbital contributions

The MICD (J(r)) is obtained by contracting the MICD susceptibility tensor with the direction vector of the external magnetic field. For a given direction of the external magnetic field, the circulation direction (tropicity) of the MIRC vortex can be in the classical diatropic (aromatic) direction, which is here the positive direction or in the negative paratropic (antiaromatic) direction. The tropicity can be assessed by visualizing the vector field. The profile of the MIRC and its strength (I) can be obtained by numerically integrating J(r) that passes through a surface (S) cutting the molecular ring.where the dot represents the scalar product of J(r) with the normal of the infinitesimal cross-section area dS. The two-dimensional integration domain begins in the center of the MIRC vortex and extends to a sufficiently large distance from the molecule in the other spatial directions. Some inaccuracy is introduced when the stagnation line of the MIRC vortex and the edge of the integration plane do not coincide. However, the studied molecules have high symmetry implying that the magnetic field, the stagnation line of the MIRC vortex and the main symmetry axis coincide. The integration plane is usually positioned so that it cuts the ring in the middle of a chemical bond. However, since charge conservation is almost fulfilled when using GIAOs and standard basis-set sizes, I is only weakly dependent on the position of S.⁶⁹ The tropicity of J(r) can be identified and separated into diatropic and paratropic contributions by following the vector field around the vortices of J(r) using e.g., the Runge–Kutta algorithm.⁷⁰

Since contributions to the MIRC strength of individual orbitals (I_i) may depend on the position of S, unique orbital contributions to the MIRC strength can be obtained by calculating the average of the MIRC strength for various orientations of S around the ring. This can be done by rotating S around the center of the ring with one of its edges along the MIRC vortex as shown in Fig. 1A. The average MIRC strength 〈I_i〉 is then obtained as

which corresponds to numerical integration around a circle using the trapezoidal rule. We use an integration step of Δα = 1° to ensure high accuracy.

Fig. 1 The orientations of the integration plane for calculating the angular dependence of the MIRC strength of benzene (A). The angular dependence of the MIRC strength of benzene for all orbitals in each irreducible representation of the C_6h point group (B). The angular dependence of the MIRC strength of the core, valence σ, and π orbitals of benzene (C).

3 Computational details

The molecular structures were optimized with Turbomole^71–73 at the density functional theory (DFT) level using the B3LYP functional,^74,75 the empirical D4 dispersion correction⁷⁶ and the def2-TZVP⁷⁷ basis sets. The one-electron density matrix and the three perturbed density matrices were obtained by calculating the nuclear magnetic resonance (NMR) shielding constants at the DFT level^78,79 using the B3LYP^74,75 functional, which is a popular functional that has been used in many applications. We have also calculated NMR shielding constants with the strong-correlation local hybrid functional (scLH22t).^80,81 The scLH22t functional was chosen because an accurate magnetizability of the ozone molecule is obtained with it. Further geometry optimizations and calculations of density matrices were performed for C₃H₃⁺ to investigate the convergence of the MIRC strength and orbital contributions to the MIRC strengths with respect to the basis-set size. The B3LYP functional and the def2-QZVP, cc-pVTZ, cc-pVQZ, cc-pV5Z*, and cc-pV6Z* basis sets were used.^77,82–84 Basis functions with a larger angular momentum number than l = 4 were removed from the cc-pV5Z and cc-pV6Z basis set because the present implementation cannot handle h and higher angular momentum functions.

The MICD susceptibility tensor was calculated with the GIMIC program.⁸⁵ The integration plane begins 8 bohr above the molecular plane and ends 8 bohr below it. In the other direction it begins at the symmetry axis and ends 15 bohr from the axis. The plane is divided into elements whose dimension is 0.05 × 0.05 bohr. We use 9 × 9 Gauss integration points in each element to ensure a high accuracy. The J(r) vector field was obtained by contracting the tensor with the direction of the external magnetic field. The external magnetic field was aligned along the main symmetry axis. The MIRC strengths were obtained by integrating the J(r) passing through a plane that cuts through the molecule.^11–13 The molecular structures and orbitals are visualized with Chemcraft,⁸⁶ and the J(r) pictures were made with ParaView.⁸⁷ The plots were generated using Matplotlib, a Python library for data visualization.⁸⁸

4 Results and discussion

4.1 Benzene

The orbital contributions to the MIRC strength calculated for the core, valence σ and π orbitals of benzene are given in Table 1. The total orbital contributions of all orbitals in each irreducible representation of the C_6h point are divergence free in the limit of complete basis set when using the GIMIC method. The MICD susceptibility calculated with the GIMIC program is practically divergence free when using standard basis sets (def2-TZVP) as seen in Fig. 1B. The MIRC contributions of all orbitals in the irreducible representations are also given in Table 1.

Table 1 The average orbital contributions to the MIRC strength (〈I_i〉 in nA/T) and the total MIRC strengths of the core, σ, and π orbitals of benzene calculated at the B3LYP and scLH22t levels. The orbital contributions and the contributions of all orbitals in each irreducible representation of the C_6h point group are compared to the orbital contributions of the isotropic shielding constant (σ(0) in ppm) previously calculated at the Hartree–Fock level in the center of the ring using the ipsocentric method.²⁵ We also compare to orbital contributions to the σ_zz(0) value calculated at the DFT level using GIAOs⁸⁹

Contribution
D6h	C6h	Type	B3LYP	scLH22t	σ(0)^25	σ zz(0)^89
a The core contributions are not included.
1a1g	1ag	Core	0.06	0.03	0.10	—
1e1u	1e1u	Core	0.09	0.15	−0.64	—
1e2g	1e2g	Core	1.19	1.22	0.73	—
1b2u	1bu	Core	0.77	0.84	−0.06	—
2a1g	2ag	σ	3.03	3.15	−0.40	14.4
2e1u	2e1u	σ	4.88	5.26	0.09	11.4
2e2g	2e2g	σ	4.35	4.56	8.94	8.8
3a1g	3ag	σ	2.81	2.73	0.34	12.3
2b2u	2bu	σ	−0.51	−0.36	−2.25	−2.9
3e1u	3e1u	σ	−3.01	−3.39	−7.47	−14.0
1b1u	3bu	σ	1.79	1.62	0.01	4.2
1a2u	1au	π	4.12	4.14	2.92	12.8
3e2g	3e2g	σ	−15.16	−15.46	−11.34	−67.6
1e1g	1e1g	π	7.53	7.59	18.90	23.4
Core			2.12	2.24	0.13	10.8
σ			−1.82	−1.90	−12.06	−33.3
Core + valence σ			0.30	0.34	−11.94	−22.5
π			11.65	11.73	21.83	36.2
	Ag		5.90	5.91	0.04	26.7^a
	Bu		2.05	2.10	−2.30	1.3^a
	E1u		1.96	2.02	−8.02	−2.6^a
	Au		4.12	4.14	2.92	12.8
	E2g		−9.62	−9.68	−1.67	−58.8^a
	E1g		7.53	7.59	18.90	23.4
Total			11.95	12.07	9.89	13.7

---

The MIRC strengths (in nA/T) calculated with GIMIC and the isotropic shielding constant in the center of the ring (σ(0) in ppm) calculated with the ipsocentric approach²⁵ correlate suggesting that similar orbital contributions to the MICD susceptibility are obtained using GIMIC and the ipsocentric method. Orbital contributions to the MICD susceptibility and other magnetic properties of benzene have also been studied using the ipsocentric method,^24,25 where they found that the HOMO contribution to J(r) of benzene is diatropic and mainly determines the MIRC strength.²⁵ The GIMIC calculations yield a total π contribution of 11.65 nA/T, whereas the core + valence σ contribution is only 0.30 nA/T. These values are close to the ones of 11.7 nA/T and 1.1 nA/T, respectively, which were previously calculated at the Hartree–Fock level with the SYSMOIC program using the ipsocentric approach.⁴⁸

The diatropic π contribution to the MIRC of benzene originates from the π orbitals in the A_u and E_1g irreducible representations of the C_6h point group. The orbitals in the E_2g irreducible representation sustain a strong paratropic MIRC of −9.62 nA/T, which is almost canceled by contributions of the valence σ and core orbitals in the other irreducible representations of the C_6h point group. Since the separation into the core and valence σ contributions to the MIRC strength depends on the orientation and position of the integration plane, we compare average contributions that are obtained as described in Section 2.3. The average contribution of the core and valence σ orbitals are 2.12 nA/T and −1.82 nA/T, respectively. The contributions of the valence σ and core orbitals are tiny when the integration plane cuts the ring in the middle of a C–C bond as seen in Fig. 1C.

In Table 1, one can see that all occupied valence orbitals contribute significantly to the MIRC strength. The contribution of the core orbitals is diatropic, whereas the contributions of the valence σ orbitals are both positive and negative. The 〈I_i〉 values of the valence orbitals have in general the same sign as the corresponding orbital contributions to σ(0), since diatropic orbital contributions to the MIRC result in diamagnetic contribution to σ(0) and paratropic orbital contributions to the MIRC result in paramagnetic shielding contributions.^49,50 Even though the orbital contributions to the MIRC strength and to σ(0) have the same trend, the relative size of the contributions significantly differs. The discrepancy is more severe when comparing the formally divergence-free contributions of all orbitals in each irreducible representation. The orbitals in the A_g irreducible representation sustain an MIRC of 5.90 nA/T at the B3LYP level, whereas the corresponding σ(0) value of 0.04 ppm is tiny.

For the B_u and E_1u irreducible representations, the MIRC strength and the σ(0) contributions have opposite signs. The MIRC contribution of the orbitals in the E_2g irreducible representation is strong and paratropic, whereas the σ(0) contribution is weakly paramagnetic. The π contribution to the MIRC strength of the only orbital in the A_u irreducible representation is about 35% of the MIRC, which is of the same relative size as the A_u contribution to σ(0). The MIRC contribution of the two π orbitals in the E_1g irreducible representation is about 65% of the MIRC strength, whereas the corresponding contribution to σ(0) is a factor of two larger than the total σ(0) value. The three π orbitals (HOMO−4, HOMO−1 and HOMO) and their MIRC contributions as a function of the orientation angle of the integration plane are shown in Fig. 2. The MIRC strength of HOMO−4 with respect to the orientation of the integration plane is almost independent of the angle because HOMO−4 is the only orbital in the irreducible representation A_u of the group of points C_6h. The oscillating behavior of the contributions of the two π orbitals in E_1g and that the sum of them is constant are also seen in Fig. 2. The small waves of the horizontal lines near the nuclei are due to a tiny J(r) leakage because finite basis sets have been used.¹¹

Fig. 2 The π orbitals of benzene (top), their MICD (middle) and the dependence of their contributions to the MIRC strength (in nA/T) as a function of the angle of the integration plane (bottom). The orbital pictures are made with ParaView.⁸⁷

There is in general no clear correlation between the orbital contribution to the MIRC strength and the orbital contributions to the σ_zz(0) and σ(0) values of benzene. However, the MIRC strength of the π orbitals and the π-orbital contributions to the σ_zz(0) values correlate. In a study on 43 molecular rings, the best linear correlation was obtained for NICS_π,zz(z) (−σ_π,zz(z)) values and MIRC strengths of the π orbitals when z is 1.25 Å.⁹⁰

Even though the σ(0) contributions are calculated at the Hartree–Fock level using the ipsocentric approach,⁴⁸ the comparison of the MIRC strengths and σ(0) values casts a shadow over the ability of NICS calculations to estimate orbital contributions to the MIRC strength. We discuss in the text only MIRC strengths calculated at the B3LYP level, since calculations on aromatic molecules using the scLH22t functional yield MIRC values that agree well with the ones obtained at the B3LYP level. See Table 1.

4.2 Borazine

The contributions to the MIRC strength of borazine of all orbitals of each irreducible representation of the C_3h point group are given in Table 2, whereas the contributions of each occupied orbital are given in the ESI.† The average MIRC contributions of the core, valence σ and π orbitals are also reported in Table 2. The orbital contribution to the MIRC strength of the HOMO is −0.47 nA/T. The HOMO is a doubly degenerate π orbital and the only orbital belonging to the E′′ irreducible representation. For benzene, the HOMO–LUMO transition is translationally allowed, resulting in a strong diatropic contribution to the MIRC of the HOMO, whereas the HOMO–LUMO transition of borazine is translationally and rotationally allowed,²⁵ resulting in a weak MIRC contribution. The J(r) of the HOMO of borazine consists of local vortices circulating around the nitrogen atoms.⁹¹ The MIRC strength of the HOMO is independent of the orientation of the integration plane because it is the only occupied orbital belonging to the E′′ irreducible representation.

Table 2 The average contributions to the MIRC strength of borazine (〈I_i〉 in nA/T) of all orbitals in each irreducible representation of the C_3h point group. The calculations were performed at the B3LYP and scLH22t levels. The contributions of the core, valence σ and π orbitals are also given. Orbital contributions are reported in the ESI

Contribution	B3LYP	scLH22t
A′	7.14	7.06
E′	−6.94	−6.87
A′′	3.52	3.50
E′′	−0.46	−0.63
Core	3.04	2.80
Valence σ	−2.85	−2.61
Core + valence σ	0.20	0.19
π	3.06	2.87
Total	3.25	3.06

---

The MIRC contribution of HOMO−2 is 3.52 nA/T, which is the third π orbital. Borazine is weakly aromatic⁹¹ sustaining an MIRC of 3.25 nA/T, which is at the borderline between the MIRC strength of aromatic and nonaromatic molecules. Molecules with an absolute MIRC strength of <3 nA/T i.e., 25% of the benzene value can be considered nonaromatic. The average contribution of the core orbitals of 3.04 nA/T almost cancels the contribution of the valence σ orbitals of −2.85 nA/T. The contributions of the core + valence σ orbitals of 0.2 nA/T and of the π orbitals of 3.06 nA/T are close to the ones of 0.3 nA/T and 1.8 nA/T, respectively, which were obtained in SYSMOIC calculations at the Hartree–Fock level using the ipsocentric approach.⁴⁸

4.3 1,5-Dibora−2,4-diazabenzene

1,5-Dibora-2,4-diazabenzene (C₂B₂N₂H₆) is isoelectronic with benzene and aromatic sustaining an MIRC strength of 7.67 nA/T, which is more than 60% of the one for benzene. Since it belongs to the C_2v point group, C_s is its highest possible point group in the presence of an external magnetic field. The contributions to the MIRC strength of C₂N₂B₂H₆ of all orbitals in each irreducible representation of the C_s point group are given in Table 3, whereas the contributions of each occupied orbital are given in the ESI.† The MICD susceptibility in the A′ (σ) and A′′ (π) irreducible representations is divergence free. The MIRC strength of the π orbitals is 7.59 nA/T. Thus, almost no MIRC is sustained in the σ orbitals. The average MIRC strength of the core orbitals of 2.67 nA/T is canceled by the contributions of the valence σ orbitals. The MIRC strengths of the core and valence σ orbitals depend on the orientation of the integration plane because the orbitals belong to the same irreducible representation of the C_s point group.

Table 3 The average contributions to the MIRC strength of C₂B₂N₂H₆ (〈I_i〉 in nA/T) of the orbitals in the two irreducible representations of the C_s point group. The calculations were performed at the B3LYP and scLH22t levels. The contributions of the core, valence σ and π orbitals are also given. Orbital contributions are reported in the ESI

Contribution	B3LYP	scLH22t
Core	2.67	2.52
Valence σ	−2.59	−2.42
A′ (core + valence σ)	0.08	0.10
A′′ (π)	7.59	7.55
Total	7.67	7.65

---

4.4 Porphin

For porphin, the average contribution of the core orbitals of 11.22 nA/T is almost canceled by the contribution of the valence σ orbitals of −10.75 nA/T. The resulting contribution of the core + valence σ orbitals of 0.47 nA/T is much smaller than the average contribution of 27.2 nA/T of the π orbitals. The total MICD strength is 27.68 nA/T. The contributions to the MIRC strength of porphin of all orbitals in each irreducible representation of the C_2h point group are given in Table 4. The contributions of the core, valence σ, and π orbitals are also reported.

Table 4 The average contributions to the MIRC strength (in nA/T) of porphin, tetraoxa-isophlorin, 1,4-cyclohexadiene and cyclobutadiene of the orbitals in each irreducible representation of the C_2h point group. The calculations were performed at the B3LYP and scLH22t levels. The contributions of the core, valence σ and π orbitals are also given. Orbital contributions are reported in the ESI

Contribution	Porphin		Isophlorin		1,4-Cyclohexadiene		Cyclobutadiene
	B3LYP	scLH22t	B3LYP	scLH22t	B3LYP	scLH22t	B3LYP	scLH22t
Ag	0.06	0.14	−0.26	−0.34	−4.08	−4.16	8.89	8.85
Bg	12.50	12.65	−76.80	−60.45	3.62	3.46	−19.88	−20.56
Au	14.70	14.90	13.62	12.95	−4.41	−4.13	3.84	3.85
Bu	0.41	0.48	−0.56	−0.59	4.23	4.26	−12.79	−12.91
Core	11.22	10.56	10.86	10.31	2.03	2.14	1.32	1.41
σ	−10.75	−9.94	−11.69	−11.24	−1.88	−2.04	−5.22	−5.47
Core + valence σ	0.62	−0.83	−0.93	0.15	0.11	−3.90	−4.06
π	27.20	27.55	−63.17	−47.50	−0.79	−0.67	−16.05	−16.71
Total	27.68	28.17	−64.00	−48.44	−0.64	−0.56	−19.95	−20.76

---

About half of the π contribution to the MIRC strength (14.70 nA/T) originates from the orbitals in the A_u irreducible representations of the C_2h point group and the contribution of the B_g orbitals is 12.50 nA/T. The averaged contributions of the individual orbitals in the ESI† show that nine π orbitals contribute more than 1 nA/T to the total MIRC. The individual MIRC contribution of these π orbitals is almost the same for seven of them (∼3.5 nA/T). The only exceptions are the contributions of HOMO of 2.20 nA/T and of HOMO−1 of 1.28 nA/T. The sum of their contributions is also ∼3.5 nA/T. The orbital contributions of four of the π orbitals are tiny, that is, in the range of ± 1nA/T. The contributions of HOMO−2 and HOMO−3 are −1.01 nA/T and −0.01 nA/T, respectively. A previous study of orbital contributions to the J(r) of porphin using the ipsocentric method suggested that four frontier π orbitals contribute to the MIRC strength,³¹ which does not agree with the results obtained in our study. The total MIRC strength is almost entirely determined by contributions of the π orbitals. However, many energetically low-lying π orbitals contribute to the total MIRC.

The diatropic contribution to the MIRC of the HOMO can be explained by the translationally allowed transition from HOMO (b_1u) to LUMO (b_3g), whereas the diatropic contributions of the energetically low-lying π (b_1u) orbitals using selection rules may involve transitions to many virtual orbitals belonging to the same irreducible representation as the LUMO.

4.5 Nonaromatic molecule

The absence of continuous electron delocalization is generally the reason why molecular rings are nonaromatic. Here, we study 1,4-cyclohexadiene as an example of a nonaromatic molecular ring. The electron delocalization of 1,4-cyclohexadiene ring is disrupted by the out-of-plane hydrogen atoms at the two sp³-hybridized carbon atoms. The MIRC strengths of the σ and π orbitals of 1,4-cyclohexadiene can be separated because its molecular structure belongs to the D_2h point group. In the presence of a magnetic field along the main symmetry axis, the core and valence σ orbitals belong to the A_g and B_u irreducible representations of the C_2h point group. The sum of their contributions to the MIRC strength is 0.15 nA/T. An MIRC strength of the π orbitals of −0.79 nA/T is obtained as the sum of the contributions of the orbitals in the A_u and B_g irreducible representations. Since 1,4-cyclohexadiene does not sustain any strong MIRC in the σ nor in the π orbitals it is nonaromatic, even though the transitions between occupied and unoccupied frontier orbitals are translationally allowed. The MIRC strengths of the orbitals in each irreducible representation are given in Table 4 where we also report the average contributions of the core, valence σ and π orbitals. Orbital contributions are reported in the ESI.†

4.6 Antiaromatic molecules

Tetraoxa-isophlorin is a strongly antiaromatic porphyrinoid. Previous studies showed that it has a large magnetic transition moment from the electronic ground state to the lowest excited state resulting in a significant paramagnetic contribution to the magnetizability.^92,93 Hexadehydro[12]annulene is also an antiaromatic molecule sustaining a strong paratropic MIRC.^94–96 The average contributions to the MIRC strength of tetraoxa-isophlorin for all orbitals in the irreducible representations of the C_2h point group are reported in Table 4. The MIRC strength of hexadehydro[12]annulene for all orbitals in the irreducible representations of the C_3h point group are given in Table 5. Orbital contributions are reported in the ESI.†

Table 5 The average contributions to the MIRC strength (in nA/T) of hexadehydro[12]annulene (HDH[12]A), C₃H₃⁺ and C₃H₆ of the orbitals in each irreducible representation of the C_3h point group. The calculations were performed at the B3LYP and scLH22t levels. The contributions of the core, valence σ, and π orbitals are also given. Orbital contributions are reported in the ESI

	HDH[12]A		C3H3^+		C3H6
	B3LYP	scLH22t	B3LYP	scLH22t	B3LYP	scLH22t
A′	−0.01	0.14	3.33	3.35	3.02	3.07
E′	−0.02	−0.26	3.14	3.04	6.42	6.42
A′′	−36.86	−31.05	3.90	3.92	2.84	2.81
E′′	12.77	12.20	—	—	−2.36	−2.30
Core	11.18	9.26	1.07	1.30	0.85	0.79
σ	−11.22	−9.40	5.40	5.09	8.59	8.70
Core + σ	−0.04	−0.14	6.47	6.40	9.44	9.49
π	−24.10	−18.86	3.90	3.92	0.49	0.50
Total	−24.13	−18.99	10.37	10.32	9.92	9.99

---

The average contributions to the MIRC of the core + valence σ orbitals are nearly zero. The total MIRC strength originates mainly from the π orbitals, which belong to the B_g and A_u irreducible representations for tetraoxa-isophlorin, and to the A′′ and E′′ irreducible representations for hexadehydro[12]annulene,.

The contributions of the individual orbitals show that the MIRC strength of nine π orbitals of tetraoxa-isophlorin is diatropic. Each of them contributes more than 1.5 nA/T to the total MIRC strength. Their total diatropic MIRC strength is 23.81 nA/T. The diatropic contribution of four other π orbitals is 2.22 nA/T. The total contribution of all π orbitals except one is 26.03 nA/T, which is almost as large as for porphin. The large paratropic contribution of −89.97 nA/T of the HOMO calculated at the B3LYP level makes tetraoxa-isophlorin antiaromatic. Calculations using the scLH22t functional yield a HOMO contribution of −73.54 nA/T showing the importance of significant long-ranged Hartree–Fock exchange in the functional to properly describe magnetic response of strongly antiaromatic molecules, whereas for aromatic molecules the MIRC strengths obtained with the two functionals are almost equal.

Five of the six π orbitals of hexadehydro[12]annulene sustain a total diatropic MIRC of 16.67 nA/T, whereas the MIRC contribution of the HOMO is −40.76 nA/T making the molecule antiaromatic. The strong paratropic MIRC contribution of the HOMO can be explained by using symmetry selection rules. The transition of hexadehydro[12]annulene is purely rotationally allowed as also for the HOMO (b_2g)–LUMO (b_3g) transition of tetraoxa-isophlorin.

The planar structure of cyclooctatetraene (COT), which belongs to the D_4h point group, is reduced to the C_4h point group in the presence of an external magnetic field along the main symmetry axis. The average contributions to the MIRC strength of COT of all orbitals in each irreducible representation are given in Table 6. The average contributions of the core, valence σ, and π orbitals are also given. Orbital contributions are reported in the ESI.†

Table 6 The average contributions to the MIRC strength (〈I_i〉 in nA/T) of planar cyclooctatraene for all orbitals in the irreducible representations of the C_4h point group and of the bent cyclooctatraene for all orbitals in the irreducible representations of the S₄ point group are given. The calculations were performed at the B3LYP and scLH22t levels. Average contributions of the core, valence σ and π orbitals are also given. Orbital contributions are reported in the ESI

Contribution (C4h)	B3LYP	scLH22t
Ag	6.45	6.46
Bg	−5.65	−5.83
Eg	7.48	7.30
Au	4.10	4.11
Bu	−51.68	−46.47
Eu	−1.11	−1.17
Core	4.06	3.74
σ	−4.37	−4.28
Core + valence σ	−0.30	−0.54
π	−40.10	−35.06
Total	−40.41	−35.60

---

Contribution (S4)	B3LYP	scLH22t
A	−7.47	−6.37
B	−2.02	−2.04
E	6.67	6.46
Total	−2.82	−1.96

---

Planar COT is antiaromatic and characterized by a strong paratropic HOMO contribution of −51.68 nA/T to the MIRC strength, which can be explained by the rotationally allowed HOMO (b_2u)–LUMO (b_1u) transition.⁹⁷ The three other π orbitals sustain diatropic contributions to the MIRC strength. The sum of their contributions is 11.58 nA/T, which is close to the MIRC strength of benzene. The core + valence σ contribution of only −0.30 nA/T can be ignored. Planar COT is strongly antiaromatic sustaining a paratropic MIRC of −40.41 nA/T at the B3LYP level. The MIRC strength of −35.60 nA/T calculated with the scLH22t functional is slightly smaller.

The MIRC contributions of the core + valence σ and π orbitals of planar cyclooctatetraene have been previously calculated at the Hartree–Fock level with the SYSMOIC program. The obtained core + valence σ contribution of −1.1 and the π contribution of −18.0 nA/T also show that the π contribution dominates making the molecule antiaromatic.⁴⁸

Bent cyclooctatetraene belonging to the D_2d point group is energetically more stable than the planar structure.⁹⁸ The average MIRC contribution of the HOMO of −12.77 nA/T for the bent cyclooctatetraene structure shows that bending the molecule significantly reduces the HOMO contribution to the MIRC strength. The total MIRC strength of the bent COT is only −2.82 nA/T suggesting that it is nonaromatic. The paratropic MIRC contribution of the HOMO is the largest one. Many other orbitals contribute significantly. The average contribution of the core orbitals is 2.77 nA/T. Orbital MIRC strengths are reported in the ESI.† The divergence-free contributions of all orbitals in each irreducible representation of the S₄ point group are given in Table 6. The spatial distributions of the J(r) of the HOMO of the planar and bent cyclooctatetraene are shown in the ESI.†

4.7 Strain

Molecular properties of planar ring-shaped molecules, whose ∠CCC angles significantly smaller than the standard ∠CCC angle of sp² hybridized hydrocarbons, differ from those of ordinary molecular rings due to the ring strain.⁹⁹ The strain leads to delocalization of the electrons in the σ orbitals and often to σ aromaticity. However, all strained molecular rings do not exhibit σ aromaticity.

We study orbital contributions to the J(r) of the planar cyclopropenium cation (C₃H₃⁺) and cyclobutadiene (C₄H₄) because their magnetic properties have significant contributions of the σ orbitals. Previous studies of the delocalization energy of the C₃H₃⁺ and C₄H₄ suggested that C₃H₃⁺ is σ + π aromatic and that C₄H₄ is σ + π antiaromatic.^100,101 Calculations of the orbital contributions to the MIRC strength of C₄H₄ show that 20% of the MIRC strength is sustained by the core + valence σ orbitals. Since the core contribution is diatropic, the paratropic MIRC contribution of valence σ orbitals is even 25% of the total MIRC strength. The MIRC strength of all orbitals of C₄H₄ belonging to the irreducible representations of the C_2h point group are given in Table 4. The total MIRC of the antiaromatic C₄H₄ is −19.95 nA/T at the B3LYP level and −20.76 nA/T at the scLH22t level. The HOMO (b_2g) → LUMO (b_3g) and the HOMO−1 (b_2u) → LUMO+1 (b_3u) transitions of C₄H₄ are rotationally allowed leading to σ + π antiaromaticity. Thus, a significant paratropic MIRC contribution of −12.25 nA/T arises from the HOMO−1 σ orbital due to the strain.

The core + valence σ contribution to the MIRC strength of C₃H₃⁺ is 6.47 nA/T, which is even larger than the π contribution of 3.90 nA/T. Since the average core contribution is 1.07 nA/T, the MIRC is dominated by the valence σ contribution. The MIRC contributions of C₃H₃⁺ for all orbitals of the irreducible representations of the C_3h point group are given in Table 5. The total MIRC strength of 10.37 nA/T is almost as larger as the MIRC strength of benzene. All orbitals except the HOMO sustain a diatropic MIRC. The MIRC of the HOMO is unexpectedly paratropic with a MIRC strength of −3.49 nA/T because the HOMO (e′) to LUMO+1 (e′) transition is rotationally allowed.

Previous computational studies on cyclopropane^54,56 (C₃H₆) suggested that it is not σ aromatic, whereas other MICD susceptibility studies on C₃H₆ yielded a strong MIRC in the σ orbitals and a weak MIRC contribution in the π orbitals.^48,55 The present calculations of the orbital contributions to the MIRC of C₃H₆ confirm that the MIRC is mainly sustained in the σ orbitals as shown in Table 5. Even though the total MIRC of C₃H₃⁺ (10.37 nA/T) and of C₃H₆ (9.92 nA/T) are almost equal, the contribution of the core + valence σ orbitals to the MIRC strength of C₃H₆ of 9.44 nA/T dominates. It is an order of magnitude larger than the π-orbital contribution of 0.49 nA/T. The contribution of the HOMO to the MIRC strength of C₃H₆ is 1.57 nA/T. The HOMO (e′) → LUMO+2 (e′) transition is translationally allowed, whereas the HOMO (e′) → LUMO transition is both translationally and rotationally forbidden. The contributions to the MIRC strength of core + valence σ and π orbitals obtained in this study are also in good agreement with the previously reported ones of 9.7 nA/T and 0.5 nA/T, respectively, which were obtained using the SYSMOIC program.⁴⁸

The extended Hückel rule for σ orbitals⁵¹ suggests that C₄H₄ is σ antiaromatic and C₃H₃⁺ is σ aromatic, since C₄H₄ has 4n core + valence σ electrons and C₃H₃⁺ has 4n + 2 core + valence σ electrons.

The angular dependence of the MIRC contributions of the core, valence σ, and π orbitals of C₃H₃⁺ and C₄H₄ is shown in Fig. 3. The core contribution is diatropic for all directions of the integration plane and vanishes when the plane passes through the middle of the C–C bond. The core contribution increases in the vicinity of the nuclei, where the valence σ contribution to the MIRC strength decreases. The sum of the two contributions is independent of the orientation angle of the integration plane. The core + valence σ contribution is tiny for most of the studied molecules, whereas for C₃H₃⁺ and C₄H₄ the sum of the two contributions is independent of the angle of the integration plane and it is of the same size as the π contribution to the MIRC strength.

Fig. 3 The angular dependence of the contributions of the core, valence σ, and π orbitals of C₃H₃⁺ (top) and of cyclobutadiene (bottom). The MIRC strengths of C₃H₃⁺ were calculated using the def2-QZVP basis sets to show that the divergence is smaller when using larger basis sets. The positions of the nuclei are indicated by the vertical dashed lines.

4.8 Basis-set effects

The MICD susceptibility is formally divergence free. The charge conservation condition of the orbital contributions to J(r) is fulfilled when considering all orbitals in a given irreducible representation of the point group of the molecular structure in the presence of an external magnetic field. However, the condition is not fulfilled in every point in space when finite basis sets are used in the GIMIC calculation. Even though one can correct the divergence of the MICD susceptibility,⁴⁵ it is generally not necessary because the divergence of the MICD susceptibilities obtained in GIMIC calculations is small and does not affect most interpretations. The MIRC strengths were calculated as a function of the orientation angle of the integration plane as discussed in Section 2.3. The difference between the maximum and minimum values of the total MIRC strength of C₃H₃⁺ and the orbital contributions to the MIRC strength of all orbitals in each irreducible representation of the C_3h point group are reported in Table 7. The angular dependence of the MIRC strength decreases when using systematically larger basis sets from cc-pVTZ to cc-pV6Z, while the mean values are almost converged with the triple-ζ quality basis sets. The orbital contributions as a function of the orientation angle of the integration plane are not constant even with the largest basis set suggesting that larger basis sets than cc-pV6Z must be employed to obtain a MICD susceptibility that is divergence free in all points of space. The cc-pV6Z basis sets are large but they are probably not flexible enough near the nuclei.

Table 7 The mean (nA/T) and the range (Δ, nA/T) of the MIRC strength, calculated as the difference between the minimum and maximum values, are obtained from the angular dependence of the MIRC contributions of C₃H₃⁺ for all orbitals in each irreducible representation of the C_3h point group. The MIRC strengths are calculated at the B3LYP level using the cc-pVTZ, cc-pVQZ, cc-pV5Z and cc-pV6Z basis sets

	cc-pVTZ		cc-pVQZ		cc-pV5Z		cc-pV6Z
	Mean	Δ	Mean	Δ	Mean	Δ	Mean	Δ
A′	3.33	0.29	3.32	0.29	3.32	0.21	3.32	0.16
E′	3.15	0.92	3.21	0.47	3.22	0.27	3.22	0.17
A′′	3.90	0.10	3.91	0.05	3.91	0.03	3.91	0.01
Total	10.37	0.81	10.44	0.27	10.44	0.11	10.45	0.04

---

5 Summary and conclusions

We have developed and implemented a novel approach in the GIMIC program for calculating contributions of occupied orbitals to the magnetically induced current density (MICD) susceptibility tensor. The method has been applied to various aromatic, nonaromatic and antiaromatic monocyclic molecules. The orbital contributions to the MIRC strength of benzene, borazine, 1,4-cyclohexadiene, 1,5-dibora-2,4-diazabenzene, the cyclopropenium cation, cyclobutadiene, planar cyclooctatetraene, bent cyclooctatetraene, hexadehydro[12]annulene, porphin and tetraoxa-isophlorin have been thoroughly investigated. The strength of the magnetically induced ring current (MIRC) has been obtained by numerically integrating the MICD susceptibility passing through an integration plane that crosses the molecular ring using accurate integration grids. We have used various orientations of the integration plane to assess how well the charge-conservation condition of the MICD susceptibility is fulfilled. The present results agree with previously reported data that have been calculated for some of the studied molecules using alternative methods. The present work is more comprehensive presenting for the first time detailed investigations of the orbital contributions to the MICD susceptibility and the MIRC of the eleven studied molecules.

The orbital contributions to the MICD susceptibility are not divergence free implying that the orbital contributions to the MIRC strength vary around the molecular ring. Unique average orbital contributions to the MIRC strength were obtained by calculating the MIRC strength for various orientations of the integration plane. Average MIRC strengths were obtained by numerically integrating over the orientation angle around the whole ring. The total MIRC and the MIRC contributions of all orbitals belonging to a given irreducible representation of the molecular point group in the presence of an external magnetic field are formally divergence free. The charge leakage of the MICD susceptibility calculated using the GIMIC method is tiny and can in most applications be ignored. The leakage decreases when increasing the size of the basis set and vanishes in the complete basis-set limit.

For planar antiaromatic molecules, the MIRC contribution of the HOMO is strongly paratropic and much larger than the contributions of the other orbitals, which is not the case for aromatic molecules. The total MIRC of aromatic molecules consists of orbital contributions of many orbitals. The orbital contributions to the MIRC strengths of C₄H₄, C₃H₃⁺ and C₃H₆ showed that the ring strain leads to a strong MIRC strength in an energetically high-lying σ orbital. The MIRC of the σ orbital of C₄H₄ is paratropic and for C₃H₃⁺ and C₃H₆ it is diatropic.

The presented algorithm works well for molecules whose ground-state wave function is dominated by a single Slater determinant. The size of the molecules that can be studied is limited only by the calculations of the NMR shielding constants. Since molecular magnetic properties can be obtained with the GIMIC program by numerical integration of the MICD susceptibility multiplied by the vector potential of the magnetic property, the present work also opens new possibilities for investigating orbital contributions to molecular magnetic properties.

Data availability

The data supporting this article have been included as part of the ESI.† Density functional theory calculations were performed with Turbomole version 7.8. The Turbomole webpage is https://www.turbomole.org/. GIMIC, version 2.0 can be freely downloaded from https://github.com/qmcurrents/gimic and https://zenodo.org/record/8180435. The new features of the GIMIC program are not yet available to everyone. ParaView can be downloaded from https://www.paraview.org/.

Author contributions

RTN implemented the methods, performed the calculations and prepared the pictures. DS proposed the project and developed the theory together with RTN. RTN and DS wrote the first version of the manuscript. All authors contributed to writing the article.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

This work has also been supported by the Academy of Finland through projects 340582 and 340583. The CSC – the Finnish IT Center for Science is acknowledged for computer time.

References

1. R. Hoffmann, Am. Sci., 2015, 103, 18–22.
2. 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.
3. V. I. Minkin, Pure Appl. Chem., 1999, 71, 1919–1981.
4. P. von Ragué Schleyer, C. Maerker, A. Dransfeld, H. Jiao and N. J. R. van Eikema Hommes, J. Am. Chem. Soc., 1996, 118, 6317–6318.
5. J. A. Pople, Mol. Phys., 1958, 1, 175–180.
6. R. McWeeny, Mol. Phys., 1958, 1, 311–321.
7. T. A. Keith and R. F. W. Bader, Chem. Phys. Lett., 1993, 210, 223–231.
8. P. Lazzeretti, M. Malagoli and R. Zanasi, Chem. Phys. Lett., 1994, 220, 299–304.
9. R. Zanasi, P. Lazzeretti, M. Malagoli and F. Piccinini, J. Chem. Phys., 1995, 102, 7150–7157.
10. P. Lazzeretti, Prog. Nucl. Magn. Reson. Spectrosc., 2000, 36, 1–88.
11. J. Jusélius, D. Sundholm and J. Gauss, J. Chem. Phys., 2004, 121, 3952–3963.
12. H. Fliegl, S. Taubert, O. Lehtonen and D. Sundholm, Phys. Chem. Chem. Phys., 2011, 13, 20500–20518.
13. D. Sundholm, H. Fliegl and R. J. F. Berger, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2016, 6, 639–678.
14. G. Monaco, F. F. Summa and R. Zanasi, J. Chem. Inf. Model., 2021, 61, 270–283.
15. M. Dimitrova and D. Sundholm, in Aromaticity: Modern Computational Methods and Applications, ed. I. Fernández López, Elsevier, 2021, pp. 155–194.
16. P. von Ragué Schleyer, M. Manoharan, Z.-X. Wang, B. Kiran, H. Jiao, R. Puchta and N. J. R. van Eikema Hommes, Org. Lett., 2001, 3, 2465–2468.
17. T. Heine, P. v. R. Schleyer, C. Corminboeuf, G. Seifert, R. Reviakine and J. Weber, J. Phys. Chem. A, 2003, 107, 6470–6475.
18. T. Heine, R. Islas and G. Merino, J. Comput. Chem., 2007, 28, 302–309.
19. A. C. Tsipis, I. G. Depastas and C. A. Tsipis, Symmetry, 2010, 2, 284–319.
20. R. Islas, T. Heine and G. Merino, Acc. Chem. Res., 2012, 45, 215–228.
21. N. D. Charistos, A. G. Papadopoulos, T. A. Nikopoulos, A. Muñoz-Castro and M. P. Sigalas, J. Comput. Chem., 2017, 38, 2594–2604.
22. R. Gershoni-Poranne and A. Stanger, in Aromaticity: Modern Computational Methods and Applications, ed. I. Fernández López, Elsevier, 2021, ch. 4, pp. 99–154.
23. I. Pérez-Juste, M. Mandado and L. Carballeira, Chem. Phys. Lett., 2010, 491, 224–229.
24. E. Steiner and P. W. Fowler, J. Phys. Chem. A, 2001, 105, 9553–9562.
25. E. Steiner and P. W. Fowler, Phys. Chem. Chem. Phys., 2004, 6, 261–272.
26. E. Steiner and P. W. Fowler, Chem. Commun., 2001, 2220–2221.
27. R. J. F. Berger and A. Viel, Z. Naturforsch., B:J. Chem. Sci., 2020, 75, 327–339.
28. D. Sundholm, M. Dimitrova and R. J. Berger, Chem. Commun., 2021, 57, 12362–12378.
29. H. Fliegl, D. Sundholm, S. Taubert, J. Jusélius and W. Klopper, J. Phys. Chem. A, 2009, 113, 8668–8676.
30. A. Soncini and P. W. Fowler, Chem. Phys. Lett., 2004, 396, 174–181.
31. E. Steiner and P. W. Fowler, ChemPhysChem, 2002, 3, 114–116.
32. F. London, J. Phys. Radium, 1937, 8, 397–409.
33. R. Ditchfield, Mol. Phys., 1974, 27, 789–807.
34. J. R. Cheeseman, G. W. Trucks, T. A. Keith and M. J. Frisch, J. Chem. Phys., 1996, 104, 5497–5509.
35. T. Helgaker, M. Jaszuński and K. Ruud, Chem. Rev., 1999, 99, 293–352.
36. L. N. Wirz, M. Dimitrova, H. Fliegl and D. Sundholm, J. Phys. Chem. Lett., 2018, 9, 1627–1632.
37. H. Fliegl, D. Sundholm and F. Pichierri, Phys. Chem. Chem. Phys., 2011, 13, 20659–20665.
38. K. Reiter, F. Weigend, L. N. Wirz, M. Dimitrova and D. Sundholm, J. Phys. Chem. C, 2019, 123, 15354–15365.
39. R. R. Valiev, T. Kurten, L. I. Valiulina, S. Y. Ketkov, V. N. Cherepanov, M. Dimitrova and D. Sundholm, Phys. Chem. Chem. Phys., 2022, 24, 1666–1674.
40. D. Blasco and D. Sundholm, Dalton Trans., 2024, 53, 10150–10158.
41. A. Rabe, Q. Wang and D. Sundholm, Dalton Trans., 2024, 53, 10938–10946.
42. Y. C. Lin, J. Jusélius, D. Sundholm and J. Gauss, J. Chem. Phys., 2005, 122, 214308.
43. M. Orozco-Ic, L. Soriano-Agueda, D. Sundholm, E. Matito and G. Merino, Chem. Sci., 2024, 15, 12906–12921.
44. R. J. F. Berger, M. Dimitrova, R. T. Nasibullin, R. R. Valiev and D. Sundholm, Phys. Chem. Chem. Phys., 2022, 24, 624–628.
45. G. Monaco, F. F. Summa, R. Zanasi and R. J. F. Berger, J. Chem. Phys., 2024, 161, 194105.
46. J. A. Bohmann, F. Weinhold and T. C. Farrar, J. Chem. Phys., 1997, 107, 1173–1184.
47. M. Orozco-Ic, N. D. Charistos, A. Muñoz-Castro, R. Islas, D. Sundholm and G. Merino, Phys. Chem. Chem. Phys., 2022, 24, 12158–12166.
48. G. Monaco, R. Zanasi, S. Pelloni and P. Lazzeretti, J. Chem. Theory Comput., 2010, 6, 3343–3351.
49. R. K. Jinger, H. Fliegl, R. Bast, M. Dimitrova, S. Lehtola and D. Sundholm, J. Phys. Chem. A, 2021, 125, 1778–1786.
50. H. Fliegl, M. Dimitrova, R. J. F. Berger and D. Sundholm, Chemistry, 2021, 3, 1005–1021.
51. M. J. S. Dewar, Bull. Soc. Chim. Belg., 1979, 88, 957–967.
52. M. J. S. Dewar, J. Am. Chem. Soc., 1984, 106, 669–682.
53. D. Moran, M. Manoharan, T. Heine and P. V. R. Schleyer, Org. Lett., 2003, 5, 23–26.
54. S. Pelloni, P. Lazzeretti and R. Zanasi, J. Phys. Chem. A, 2007, 111, 8163–8169.
55. P. W. Fowler, J. Baker and M. Lillington, Theor. Chem. Acc., 2007, 118, 123–127.
56. W. Wu, B. Ma, J. I-Chia Wu, P. von Ragué Schleyer and Y. Mo, Chem. - Eur. J., 2009, 15, 9730–9736.
57. R. Carion, B. Champagne, G. Monaco, R. Zanasi, S. Pelloni and P. Lazzeretti, J. Chem. Theory Comput., 2010, 6, 2002–2018.
58. Y. Hao, J. Wu and J. Zhu, Chem. - Eur. J., 2015, 21, 18805–18810.
59. M. Orozco-Ic and D. Sundholm, Phys. Chem. Chem. Phys., 2023, 25, 12777–12782.
60. P. Pulay, in Advances in Chemical Physics: Ab Initio Methods in Quantum Chemistry Part 2, ed. P. K. Lawley, John Wiley & Sons Ltd, Singapore, 1987, pp. 241–286.
61. M. Häser, R. Ahlrichs, H. P. Baron, P. Weis and H. Horn, Theor. Chim. Acta, 1992, 83, 455–470.
62. K. Wolinski, J. F. Hinton and P. Pulay, J. Am. Chem. Soc., 1990, 112, 8251–8260.
63. J. Gauss, Chem. Phys. Lett., 1992, 191, 614–620.
64. J. Olsen, K. L. Bak, K. Ruud, T. Helgaker and P. Jørgensen, Theor. Chim. Acta, 1995, 90, 421–439.
65. J. Gerratt and I. M. Mills, J. Chem. Phys., 1968, 49, 1730–1739.
66. A. Pausch, M. Gebele and W. Klopper, J. Chem. Phys., 2021, 155, 201101.
67. S. Pelloni and P. Lazzeretti, Int. J. Quantum Chem., 2011, 111, 356–367.
68. A. J. Ceulemans, Group theory applied to chemistry, Springer, 2nd edn, 2024.
69. H. Fliegl, J. Jusélius and D. Sundholm, J. Phys. Chem. A, 2016, 120, 5658–5664.
70. Q. Wang, J. Pyykkö, M. Dimitrova, S. Taubert and D. Sundholm, Phys. Chem. Chem. Phys., 2023, 25, 12469–12478.
71. R. Ahlrichs, M. Bär, H. S. Plitt and H. Schnöckel, Chem. Phys. Lett., 1989, 161, 179.
72. S. G. Balasubramani, G. P. Chen, S. Coriani, M. Diedenhofen, M. S. Frank, Y. J. Franzke, F. Furche, R. Grotjahn, M. E. Harding, C. Hättig, A. Hellweg, B. Helmich-Paris, C. Holzer, U. Huniar, M. Kaupp, A. Marefat Khah, S. Karbalaei Khani, T. Müller, F. Mack, B. D. Nguyen, S. M. Parker, E. Perlt, D. Rappoport, K. Reiter, S. Roy, M. Rückert, G. Schmitz, M. Sierka, E. Tapavicza, D. P. Tew, C. van Wüllen, V. K. Voora, F. Weigend, A. Wodyński and J. M. Yu, J. Chem. Phys., 2020, 152, 184107.
73. Y. J. Franzke, C. Holzer, J. H. Andersen, T. Begušić, F. Bruder, S. Coriani, F. Della Sala, E. Fabiano, D. A. Fedotov, S. Fürst, S. Gillhuber, R. Grotjahn, M. Kaupp, M. Kehry, M. Krstić, F. Mack, S. Majumdar, B. D. Nguyen, S. M. Parker, F. Pauly, A. Pausch, E. Perlt, G. S. Phun, A. Rajabi, D. Rappoport, B. Samal, T. Schrader, M. Sharma, E. Tapavicza, R. S. Treß, V. Voora, A. Wodyński, J. M. Yu, B. Zerulla, F. Furche, C. Hättig, M. Sierka, D. P. Tew and F. Weigend, J. Chem. Theory Comput., 2023, 19, 6859–6890.
74. A. D. Becke, J. Chem. Phys., 1993, 98, 5648–5652.
75. C. Lee, W. Yang and R. G. Parr, Phys. Rev. B, 1988, 37, 785–789.
76. E. Caldeweyher, C. Bannwarth and S. Grimme, J. Chem. Phys., 2017, 147, 034112.
77. F. Weigend and R. Ahlrichs, Phys. Chem. Chem. Phys., 2005, 7, 3297–3305.
78. K. Reiter, F. Mack and F. Weigend, J. Chem. Theory Comput., 2018, 14, 191–197.
79. Y. J. Franzke, J. Chem. Theory Comput., 2023, 19, 2010–2028.
80. C. J. Schattenberg, A. Wodyński, H. Åström, D. Sundholm, M. Kaupp and S. Lehtola, J. Phys. Chem. A, 2023, 127, 10896–10907.
81. C. J. Schattenberg and M. Kaupp, J. Phys. Chem. A, 2024, 128, 2253–2271.
82. J. Dunning and H. Thom, J. Chem. Phys., 1989, 90, 1007–1023.
83. D. E. Woon, J. Dunning and H. Thom, J. Chem. Phys., 1994, 100, 2975–2988.
84. A. Wilson, T. van Mourik and T. H. Dunning Jr, J. Mol. Struct., 1997, 388, 339–349.
85. GIMIC, version 2.0, a current density program, can be freely downloaded from https://github.com/qmcurrents/gimic and https://zenodo.org/record/8180435 (assessed 1.1.2025), 2019.
86. Chemcraft, Chemcraft - graphical software for visualization of quantum chemistry computations, https://www.chemcraftprog.com, assessed: 23.01.2025.
87. J. Ahrens, B. Geveci and C. Law, ParaView: An End-User Tool for Large Data Visualization, Visualization Handbook, Elsevier, 2005, ISBN-13: 978-0123875822, see also: http://www.paraview.org (assessed 18.1.2025).
88. J. D. Hunter, Comput. Sci. Eng., 2007, 9, 90–95.
89. N. D. Charistos, A. G. Papadopoulos and M. P. Sigalas, J. Phys. Chem. A, 2014, 118, 1113–1122.
90. S. Radenković and S. Ðorđević, Phys. Chem. Chem. Phys., 2021, 23, 11240–11250.
91. R. Báez-Grez and R. Pino-Rios, RSC Adv., 2022, 12, 7906–7910.
92. R. R. Valiev, H. Fliegl and D. Sundholm, Phys. Chem. Chem. Phys., 2017, 19, 25979–25988.
93. R. R. Valiev, G. V. Baryshnikov, R. T. Nasibullin, D. Sundholm and H. Ågren, J. Phys. Chem. C, 2020, 124, 21027–21035.
94. J. A. Pople and K. G. Untch, J. Am. Chem. Soc., 1966, 88, 4811–4815.
95. J. Jusélius and D. Sundholm, Phys. Chem. Chem. Phys., 2001, 3, 2433–2437.
96. J. Jusélius and D. Sundholm, Phys. Chem. Chem. Phys., 2008, 10, 6630–6634.
97. P. Fowler and E. Steiner, Chem. Phys. Lett., 2002, 364, 259–266.
98. R. W. Havenith, P. W. Fowler and L. W. Jenneskens, Org. Lett., 2006, 8, 1255–1258.
99. R. Huisgen, Angew. Chem., Int. Ed., 1986, 25, 297–311.
100. D. Cremer and J. Gauss, J. Am. Chem. Soc., 1986, 108, 7476–7477.
101. G. Hohlneicher, L. Packschies and J. Weber, Phys. Chem. Chem. Phys., 2007, 9, 2517–2530.

---

Footnote

† Electronic supplementary information (ESI) available: Visualization of magnetically induced current densities, average orbital contributions to the magnetically induced ring (MIRC) currents, the angular dependence of the MIRC strength of the core, valence σ and π orbitals, the angular dependence of the sum of the MIRC strengths of all orbitals in each irreducible representation, pictures of the molecular structures of the studied molecules, and the Cartesian coordinates of the studied molecules (zipped). See DOI: https://doi.org/10.1039/d5sc00627a

---