Magnetically induced current density: a semiempirical approach

Abstract

This study examines the calculation of magnetically induced current densities within the CTOCD-DZ (continuous transformation of origin of current density-dimagnetic zero) framework using various semiempirical methods based on the neglect of diatomic differential overlap (NDDO) approximation. The performance of these methods was assessed for a series of benzenoid hydrocarbons, by comparing their results with the B3LYP/def2-SVP reference. Additionally, it was proposed how the original NDDO-based semiempirical methods for calculating magnetically induced current densities can be modified to further reduce computational costs. The modified NDDO-based approach was shown to significantly decrease computational time while maintaining accuracy.

---

Introduction

Calculations of the magnetically induced current densities represent a widely recognized approach for evaluating aromaticity and electron delocalization pathways in polycyclic conjugated molecules.^1–6 However, standard visualizations of the current densities, typically calculated on a properly chosen 2D grid placed near the molecule, or numerical integration of the current density passing through selected planes, can be computationally demanding and impractical for large systems. Quantum-chemical semiempirical methods^7–10 are commonly used to reduce computational costs and to extend the method's applicability. The philosophy of semiempirical methods involves starting from the first-principles formalism, followed by the introduction of certain approximations that speed up calculations. To correct errors caused by these approximations, empirical parameters are incorporated and calibrated against reliable experimental or theoretical reference data.

Historically, the π-electron Hückel molecular orbital theory was one of the first semi-empirical methods developed.¹¹ These methods have since become widely used in modeling of ring currents and other molecular magnetic properties. The extension of the Hückel–London theory¹² by Pople¹³ and McWeeny¹⁴ was the first semiempirical method designed to overcome computational limitations in calculating ring currents. The so-called topological ring currents, which combine the graph-theoretical approach and the Hückel–London theory, continue to be recognized today for providing an accurate description of ring currents in polycyclic conjugated molecules.^2,15

It is also worth noting that there are other approximate methods used for modeling ring currents. One such method, known as the pseudo-π approach,¹⁶ can provide very accurate predictions of π-electron current densities induced in polycyclic aromatic compounds.^17–19 This model relies on very crude approximations, where only carbon atom skeleton is considered, with carbon atoms replaced by hydrogen. This simplified model should mimic the real π-electron system of the corresponding conjugated molecule. Our research group has previously shown that the pseudo-π approach can be applied to both planar and non-planar benzenoids, providing both qualitative and quantitative agreement with full DFT results.²⁰

The success of all these approximate approaches for current density calculations can be understood through the work of Fowler and Steiner, who demonstrated that within the ipsocentric framework,^21,22 current densities can be interpreted and predicted by analyzing a few frontier orbitals. Their analysis showed that the molecular magnetic response is primarily determined by the symmetry, nodal character and energies of these orbitals. Since the nodal properties and symmetry are basically independent of the computational level, this explains why ring currents and, more generally, molecular magnetic response properties can be predicted at low levels of computation.

In the present work, the performance of several semiempirical methods based on the neglect of diatomic differential overlap (NDDO) approximation¹⁰ was evaluated for calculating magnetically induced current densities in a series of benzenoid hydrocarbons. These current densities were calculated using the continuous transformation of origin of current density (CTOCD) method,^23,24 specifically its diamagnetic-zero (DZ)²⁵ or ipsocentric version.²² Further details on the CTOCD method can be found in reviews^1,26 and elsewhere.²⁷ It is worth noting that the NDDO-methods can be parameterized for reliable calculations of NMR chemical shifts.¹⁰ As demonstrated below, additional approximations can be incorporated into the original NDDO-based semiempirical methods for magnetically induced current density calculations. It will be shown that these modifications further reduce computational time without sacrificing accuracy.

Methodology

Consider a closed-shell molecule exposed to a static homogenous magnetic field B. This field is described by the magnetic vector potential A(r), such that B = ∇ × A(r). The vector potential A(r) adopts the following form:where r₀ denotes the explicit origin of the vector potential. The first-order current density is origin-independent in the complete basis set limit. However, this is not the case with a finite basis set. In 1993, Keith and Bader introduced the continuous set of gauge transformation (CSGT)^23,24 method for solving the origin-dependence issue in practical applications with finite basis sets. Later, it was suggested that the CTOCD method is a more appropriate name for the CSGT approach.^25,26 In this method, instead of using a fixed common origin r₀, a continuous function f(r) is used to set the origin of the vector potential. Within the CTOCD framework, the diamagnetic-zero (DZ)²⁵ version is the simplest approach, where the origin at the point where the current density is calculated is set to be that point itself, i.e., r₀ = f(r) = r. For this reason, the CTOCD-DZ method is also known as the ipsocentric method.²² This choice of f(r) formally cancels the diamagnetic contribution to the first-order current density:^22,24,25where the summation runs over all occupied molecular orbitals, N stands for the number of electrons, p̂ is the linear momentum operator for a single electron, and ϕ⁽⁰⁾_i and ϕ⁽¹⁾_i are molecular orbitals and their first-order corrections, respectively. The CTOCD-DZ method is a perturbative approach, and the diamagnetic term reappears in the expression for the current density.²² Within a standard perturbative formalism, the first-order corrections to the occupied molecular orbitals are given as an extension of the basis of the unoccupied, unperturbed molecular orbitals:where K is the number of basis functions. In both Hartree–Fock (HF) and density functional theory (DFT) methods, the coefficients C⁽¹⁾_pi are determined by solving a set of coupled equations:where ε⁽⁰⁾ denotes molecular orbital energy, 〈qi|pj〉 and 〈ji|pq〉 are two-electron integrals, and Ĥ⁽¹⁾ is the first-order correction to the one-electron part of the Fock operator. The factor ξ takes a value of 1 for the HF level, and for DFT, it represents the portion of the HF exchange in the given functional. This implies that for pure generalized gradient approximation (GGA) functionals, ξ is set to zero.²⁷

Now, let us consider the form of Ĥ⁽¹⁾ from eqn (4). For the vector potential given by eqn (1) the operator Ĥ⁽¹⁾ has the following form:

where l̂ is the angular momentum operator for a single electron with respect to the origin r₀. However, in the CTOCD approach, the origin is continuously distributed. It has been shown that the current density is invariant to gauge transformations.²⁸ In the simplest case if the origin shifts by a constant vector d, the vector potential changes from A(r) to A′(r):

The change in vector potential from A(r) to A′(r) leaves the magnetic field B unchanged, but this alters the first-order correction of the one-electron operator Ĥ⁽¹⁾:

Substituting eqn (7) into eqn (4) leads to six coupled equations to be solved. These equations involve six different operators, l̂_x, l̂_y, l̂_z, p̂_x, p̂_y and p̂_z, which results in six different orbital corrections: and which are needed to calculate the current density by eqn (2). In most practical applications, the magnetic field is applied along a single axis, simplifying the problem to solving only a subset of coupled equations. For instance, if the field is oriented along the z-axis, only the corrections and are required to determine the first-order current density.

Computational details

Geometry optimization and frequency calculations were performed at the B3LYP^29–31/def2-SVP³² level of theory. Starting from these optimized geometries, single-point energy calculations were performed using several NDDO-based semiempirical methods (MNDO, AM1, PM3, PM6, PM7),^9,10 which provided the corresponding wavefunctions. All electronic structure calculations were carried out using the Gaussian 09 program.³³

Magnetically induced current densities (MICDs) were calculated at the B3LYP/def2-SVP level of theory, as well as using various NDDO-based semiempirical methods (MNDO, AM1, PM3, PM6, PM7). To verify the accuracy of the referent data, the MICDs of selected molecules were calculated using the B3LYP/def2-TZVP method. To calculate the MICD according to eqn (2), both the unperturbed molecular orbitals and their first-order corrections are required. The unperturbed molecular orbitals were obtained from the Gaussian09 output files (formatted chk-files), while the needed first-order orbital corrections were calculated using eqn (4). In all calculations, the external magnetic field was set to be aligned along the z-axis, with the molecule lying in the xy-plane. MICD and bond current strength³⁴ calculations were performed using an in-house FORTRAN program. The current densities were visualized using the Paraview program.³⁵ In all current density maps, clockwise circulations represent diatropic current densities. Bond current strengths (I, in nA T⁻¹) were obtained by numerical integration of the current densities passing through a rectangle that perpendicularly bisects the given bond. For C–C bonds located on the perimeter of the studied systems, this rectangle starts from the ring center, passes through the bond center, and extends 5 bohr outside the molecule. In the case of shared C–C bonds, the integration rectangle connects the centers of the two rings that share the given bond. In all cases, the integration rectangle extends 5 bohr above and 5 bohr below the molecular plane.

Results and discussion

The performance of the NDDO-based methods in calculating MICDs was evaluated by comparing them to the B3LYP/def2-SVP method as a reference. Fig. 1 and Fig. S1 (ESI†) display the total and π-only current density maps of benzene, dibenzo[e,l]pyrene and coronene, as obtained with the B3LYP and MNDO methods. Fig. 1 illustrates that the MNDO method correctly predicts a strong diatropic current in benzene, as well as very pronounced local diatropic circulations in the four peripheral hexagonal rings in dibenzo[e,l]pyrene, in agreement with predictions based on Clar's theory.^36,37 In coronene, the MNDO method predicts a strong diatropic circulation along the perimeter and weaker paratropic currents within the central hexagon, consistent with previous results.³⁸ These counter-rotating ring currents in coronene contradict the predictions of the “annulene-within-an-annulene” model,^38–40 which suggests that both the rim and the hub, containing 18 and 6π electrons, respectively, should exhibit diatropic currents. The intensity of the currents can be compared using the j_max quantity, which represents the maximum values of the modulus of the induced current density in the plotting plane. As illustrated in Fig. 1, the MNDO method tends to overestimate the j_max values when compared to those obtained from the B3LYP/def2-SVP method. For the π-electron current density maps, the agreement between the B3LYP and MNDO results is even more consistent than for the total current densities (Fig. S1, ESI†). The current density maps presented in Fig. S2–S4 (ESI†) further show that all the other semiempirical methods (AM1, PM3, PM6 and PM7) provide very similar results, with reasonably good agreement to the reference B3LYP approach. This indicates that the NDDO-based semiempirical methods are able to capture the key features of the induced currents in the studied molecules.

Fig. 1 Total current density maps plotted 1 bohr above the molecular plane of benzene (a), dibenzo[e,l]pyrene (b) and coronene (c) obtained using the B3LYP/def2-SVP (left) and MNDO (right) methods. The maximum magnitudes of the current in the plotting plane (j_max) are expressed in a.u.

The current density maps shown in Fig. 1 and Fig. S1–S4 (ESI†) demonstrate that while the employed NDDO-based methods predict the same general patterns of current density distribution, they differ in predicting the intensity of the induced current density in these molecules. Notably, the PM3 method consistently underestimates the current density intensities in these compounds. To quantitatively assess the performance of the used semiempirical methods, current strengths were calculated for selected carbon–carbon bonds in a series of 11 benzenoid hydrocarbons, as depicted in Fig. 2. Table S1 (ESI†) summarizes the calculated current strengths obtained using various methods. First, to assess the reliability of the B3LYP/def2-SVP results, the current strengths were also computed using the more extended def2-TZVP basis set. A good linear correlation is observed between both the total and π-only current strengths obtained with the def2-TZVP and def2-SVP basis sets (Fig. S5, ESI†). Additionally, there is excellent numerical agreement with the results from the two basis sets, as the average absolute errors (AAE) of the B3LYP/def2-SVP current strengths relative to the B3LYP/def2-TZVP results are found to be 0.4 and 0.3 nA T⁻¹ for the total and π-only current strengths, respectively. These findings justify the use of the B3LYP/def2-SVP method as the reference in this study.

Fig. 2 Molecules under study with labeled bonds for which current density strengths were calculated.

With the exception of PM3, all other semiempirical approaches predict slightly stronger current strengths than the B3LYP/def2-SVP method (Table S1, ESI†). The AAEs of the semiempirical methods relative to the B3LYP/def2-SVP results were found to be 4.2, 2.0, 2.3, 6.4 and 4.4 nA T⁻¹ for MNDO, AM1, PM3, PM6 and PM7, respectively. These errors are relatively large, though the best numerical agreement with the DFT results was obtained with the AM1 method. However, cross-correlation analysis presented in Table 1 reveals that despite relatively large AAE, there exist good correlations between the DFT and NDDO-based current strengths. The strongest correlation is found between the current strengths calculated using B3LYP and MNDO as shown in Fig. 3, while analogous correlations for other methods are presented in Fig. S6 (ESI†). It is important to note that all plots in Fig. 3 and Fig. S6 (ESI†) contain a negative-valued outlier, corresponding to paratropic currents circulating along the inner hexagon in coronene. The correlation coefficients calculated after removing this outlier further confirm a relatively good correlation between the DFT and NDDO-based current strengths (Table 1). This analysis indicates that all the employed NDDO-based methods provide reasonably good agreement with the full B3LYP/def2-SVP results. Given that the best correlation was found for the total bond current strengths obtained with the MNDO method, this method was also employed to compute the π-electron current density strengths for the same set of bonds in the studied molecules (Table S1, ESI†). The correlation between the π-electron current density strengths obtained with B3LYP/def2-SVP and MNDO (Fig. 3) shows that the correlation is nearly identical to that of the total currents. Moreover, the AAE between the MNDO and B3LYP π-electron current strengths is 3.5 nA T⁻¹, which represents a reduced error compared to the corresponding AAE for the total currents. Additionally, the current strength profiles⁵ obtained by stepwise current density integration along the line connecting carbon atoms in the para position of benzene provide further insights into the performance of the NDDO-based methods (Fig. S7, ESI†). In general, it can be observed that the MNDO method underestimates the intensity of both the diatropic global circulation involving H-atoms and the paratropic circulation around the benzene ring center. These discrepancies highlight the limitations arising from the minimal basis set and other approximations incorporated into the MNDO method. On the other hand, the MNDO provides a more accurate description of π-electron currents.

Table 1 Correlation coefficient (R) between the total current density strengths for the selected C–C bonds in molecules 1–11 (as specified in Fig. 2), obtained with different methods. The values in parentheses correspond to those calculated by removing the negative outlier point (see Fig. 3 and Fig. S6, ESI)

	B3LYP	MNDO	AM1	PM3	PM6	PM7
B3LYP	1.000 (1.000)
MNDO	0.975 (0.930)	1.000 (1.000)
AM1	0.968 (0.900)	0.988 (0.968)	1.000 (1.000)
PM3	0.944 (0.819)	0.982 (0.956)	0.977 (0.928)	1.000 (1.000)
PM6	0.971 (0.917)	0.989 (0.971)	0.999 (0.999)	0.971 (0.918)	1.000 (1.000)
PM7	0.972 (0.916)	0.989 (0.970)	0.999 (0.999)	0.973 (0.920)	1.000 (1.000)	1.000 (1.000)

---

Fig. 3 Correlation between the total (a) and π-electron (b) current density strengths (I, in nA T⁻¹) for the selected C–C bonds in molecules 1–11 (as shown in Fig. 2), obtained using the B3LYP/def2-SVP and MNDO methods.

Overall, the semiempirical-method-based MICD calculations provide reasonably accurate predictions of current density intensities and distribution in benzenoid hydrocarbons. On the positive side, these methods significantly reduce the computational resources required. For instance, for coronene (11), the number of basis functions used in the B3LYP/def2-SVP and MNDO calculations is 420 and 108, while the corresponding CPU times for the MICD calculations needed to generate the MICD plots are 27.32 and 0.45 seconds (using a single Intel Xeon CPU E5-2630 v3 @ 2.40GHz processor).

In what follows, we explore whether further improvements can be made to the efficiency of the MICD calculations using a semiempirical approach. The most computationally demanding aspect of current density calculations for larger molecules lies in the evaluation of the two electron-integrals required for eqn (4). Additionally, these integrals must be transformed into the molecular orbital basis, further increasing computational costs. As mentioned above, for pure GGA functionals, the calculation of two-electron integrals is not required (ξ = 0 in eqn (4)). Following the philosophy of the semiempirical approach, here we aimed to investigate the possibility of bypassing the calculations of two-electron integrals. In particular, we repeated all calculations by setting ξ = 0 in eqn (4), and to maintain accuracy, an additional empirical parameter α was introduced: where C⁽¹⁾_pi is obtained from eqn (4) for ξ = 0, and the corrected is then used for current density calculation by eqn (2). The parameter α is aimed at compensating for the omission of the two electron integrals in eqn (4).

The MNDO wavefunction was used in the remainder of this study. The selection of MNDO, among other semiempirical methods showing similar performances, was based on the quality of the correlation between the bond current strengths calculated with B3LYP (Table 1 and Fig. 3). By setting ξ = 0 and α = 1, total and π-only current density maps were calculated for benzene, dibenzo[e,l]pyrene and coronene as shown in Fig. S8 (ESI†). It can be observed that the modified MNDO approach gives significantly weaker current densities compared to the original MNDO method. However, even with this simplified approach, the method is still capable of capturing the key features of current density distribution in these molecules.

To optimize the value of the parameter α, an extended set of molecules was employed. In addition to the molecules shown in Fig. 2, this set also includes those from Fig. 4. The molecules were divided into a training set (molecules 1–15) and a test set (molecules 16–18). There are 52 and 23 carbon–carbon bonds considered in the training and test sets, respectively. Current strengths for the selected bonds of the compounds from the training set were calculated by varying the α value. Fig. 5 shows how the AAE of the current strengths, relative to the referent B3LYP results, depends on α. The minimum AAE was observed at α_opt = 3.1. Remarkably, for α_opt = 3.1, the AAE for the training set was 1.7 nA T⁻¹, which is significantly lower than the AAE of 4.0 nA T⁻¹ obtained for the original MNDO approach using molecules 1–11. The AAE for bond current strengths in molecules 1–11, calculated using the α_opt value, was found to be 1.8 nA T⁻¹ (Table S1, ESI†). This suggests that the modified MNDO approach (α-MNDO), reduces computational time while, counterintuitively, enhancing accuracy compared to the original MNDO method.

Fig. 4 Continuation of the set of studied molecules from Fig. 2, with labeled bonds for which current density strengths were calculated.

Fig. 5 Dependence of AAE (in nA T⁻¹) on the parameter α for the selected C–C bonds in molecules 1–15.

The performance of the α-MNDO approach was further tested using the compounds from the test set. The calculated current strengths for the selected bonds are provided in Table S2 (ESI†). The α-MNDO bond current strengths were plotted against the B3LYP/def2-SVP results in Fig. 6. A good correlation between the B3LYP and α-MNDO results was observed, with an R² value of 0.901 (R = 0.949). Furthermore, the AAE for the current strengths of the test set, relative to the B3LYP/def2-SVP results, was found to be 1.6 nA T⁻¹.

Fig. 6 Correlation between the total current density strengths (I, in nA T⁻¹) for the selected C–C bonds in molecules from the test set (16–18), obtained using the B3LYP/def2-SVP and α-MNDO methods.

Current density maps for compounds 17 and 18 were calculated using both the B3LYP and α-MNDO approaches (Fig. 7 and Fig. S10, ESI†). It is evident that, in addition to good quantitative agreement with the DFT results, the α-MNDO approach can efficiently reproduce current density maps of large benzenoid molecules while using very modest computational resources. For example, for compound 18 (C₁₁₄H₃₀), the number of basis functions used in the B3LYP/def2-SVP and α-MNDO calculations is 1860 and 486, respectively, while the corresponding CPU times for the MICD calculations needed to generate the MICD plots are 43 327.8 and 274.3 seconds (using a single Intel Xeon CPU E5-2630 v3 @ 2.40 GHz processor). This data clearly demonstrates computational efficiency of the α-MNDO method, which on the other hand, can provide reasonably good accuracy.

Fig. 7 Total current density maps plotted 1 bohr above the molecular plane of 18 obtained using the B3LYP/def2-SVP (left) and α-MNDO (right) methods. The maximum magnitude of the current in the plotting plane (j_max) are expressed in a.u.

It is worth noting that compounds 12–18 (Fig. 4) belong to the class of fully (or totally-resonant) benzenoid hydrocarbons, which can be represented as a single Clar formula in which all π-electrons are grouped in aromatic sextets.^36,41 Previous studies have shown that full benzenoids exhibit strong local diatropic circulations that follow the Clar-structure model, along with a global diatropic current around the molecular perimeter.^20,37 Both the current density maps (Fig. 7 and Fig. S9, S10, ESI†) and the integrated current density strengths (Tables S1 and S2, ESI†) reveal a consistent trend: the α-MNDO method tends to overestimate the intensity of local circulations within hexagonal rings, while underestimating the significance of the global circulation along the perimeter.

To investigate the extent to which the applicability of the α-MNDO method can be extended, the additional compounds depicted in Fig. 8 were considered. Total and π-only current density maps for biphenylene (20) are shown in Fig. 9, while Fig. S11 (ESI†) displays the maps for the other compounds from Fig. 8. The bond current strengths are summarized in Table S3 (ESI†). The aromaticity of biphenylene has been much studied in the past, with previous research showing that it sustains a strong paratropic circulation within the four-membered ring and a somewhat weaker diatropic circulation along the six-membered rings.^2,42 The α-MNDO method correctly predicts the tropicity of the currents in six- and four-membered rings of 20. Similarly, the α-MNDO accurately depicts the counter-rotating ring currents along the outer rim and the inner hub of the examined circulenes (19 and 21) in agreement with previous findings.^43,44 The best agreement with the B3LYP and α-MNDO methods is found for Zn–porphyrin (22).⁵ These results suggest that the α-MNDO method has potential to be successfully applied to non-benzenoid molecules, including those with heteroatoms or with opposing aromaticity characteristics. Although based on a limited set of molecules (Fig. 8), it can be concluded that the present α-MNDO method systematically underestimates the intensity of paratropic currents. It is important to note that the α-MNDO method was parameterized for benzenoid molecules, and its applicability to non-benzenoid molecules will be addressed in future studies.

Fig. 8 Non-benzenoid molecules studied, with labeled bonds for which current density strengths were calculated.

Fig. 9 Total (a) and π-electron (b) current density maps plotted 1 bohr above the molecular plane of 20 obtained using the B3LYP/def2-SVP (left) and α-MNDO (right) methods. The maximum magnitude of the current in the plotting plane (j_max) is expressed in a.u.

Conclusions

This study evaluates the performance of NDDO-based semiempirical methods for calculating magnetically induced current densities using the CTOCD-DZ approach, focusing on a series of benzenoid hydrocarbons. It was shown that all the semiempirical methods employed (MNDO, AM1, PM3, PM6 and PM7) successfully capture the key features of induced currents, in good agreement with the B3LYP/def2-SVP results. In addition to the qualitative agreement observed in current density maps, a good correlation was found between the bond current strengths obtained from the NDDO-based and those from DFT calculations.

Furthermore, a modification to the original MNDO method for calculating magnetically induced current densities was proposed to further enhance computational efficiency. Specifically, the proposed α-MNDO approach eliminates the two-electron integral calculations, introducing an additional empirical parameter α to compensate this approximation. The optimal value of α was determined based on a training set of 15 benzenoid hydrocarbons, which includes 52 distinct carbon–carbon bonds. The α-MNDO approach yielded an average absolute error of 1.6 nA T⁻¹ in bond current strength calculations for the test set. As a result, the α-MNDO method proved to be significantly more computationally efficient than the original MNDO method, while also offering even better agreement with full DFT results.

Data availability

The data supporting this article have been included as part of the ESI.†

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

This work is dedicated to the students and educators in Serbia who, during the 2024–2025 academic year, courageously opposed corruption and advocated integrity in the educational system. This work was supported by the Serbian Ministry of Education, Science and Technological Development (Agreement No. 451-03-136/2025-03/200122).

Notes and references

1. P. Lazzeretti, Prog. Nucl. Magn. Reson. Spectrosc., 2000, 36, 1–88.
2. J. A. N. F. Gomes and R. B. Mallion, Chem. Rev., 2001, 101, 1349–1384.
3. P. Lazzeretti, Phys. Chem. Chem. Phys., 2004, 6, 217–223.
4. P. W. Fowler, C. M. Gibson and D. E. Bean, Proc. R. Soc. A, 2014, 470, 20130617.
5. D. Sundholm, H. Fliegl and R. J. F. Berger, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2016, 6, 639–678.
6. D. Sundholm, M. Dimitrova and R. J. F. Berger, Chem. Commun., 2021, 57, 12362–12378.
7. J. J. P. Stewart, Semiempirical Molecular Orbital Methods, in Reviews in Computational Chemistry, ed. K. B. Lipkowitz and D. B. Boyd, 1990, pp. 45–81.
8. T. Clark, J. Mol. Struct.: THEOCHEM, 2000, 530, 1–10.
9. T. Bredow and K. Jug, Theor. Chem. Acc., 2005, 113, 1–14.
10. W. Thiel, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2014, 4, 145–157.
11. C. A. Coulson, B. O’Leary, B. O’Leary and R. B. Mallion, Hückel Theory for Organic Chemists, Academic Press, London, 1978.
12. F. London, J. Phys. Radium, 1937, 8, 397–409.
13. J. A. Pople, Mol. Phys., 1958, 1, 175–180.
14. R. McWeeny, Mol. Phys., 1958, 1, 311–321.
15. T. K. Dickens and R. B. Mallion, MATCH Commun. Math. Comput. Chem., 2016, 76, 297–356.
16. P. W. Fowler and E. Steiner, Chem. Phys. Lett., 2002, 364, 259–266.
17. M. Orozco-Ic, M. Dimitrova, J. Barroso, D. Sundholm and G. Merino, J. Phys. Chem. A, 2021, 125, 5753–5764.
18. M. Orozco-Ic and D. Sundholm, Phys. Chem. Chem. Phys., 2022, 24, 22487–22496.
19. A. Mahmood, M. Dimitrova and D. Sundholm, J. Phys. Chem. A, 2023, 127, 7452–7459.
20. M. Antić, S. Đorđević, B. Furtula and S. Radenković, J. Phys. Chem. A, 2020, 124, 371–378.
21. E. Steiner and P. W. Fowler, Chem. Commun., 2001, 2220–2221.
22. E. Steiner and P. W. Fowler, J. Phys. Chem. A, 2001, 105, 9553–9562.
23. T. A. Keith and R. F. W. Bader, Chem. Phys. Lett., 1993, 210, 223–231.
24. T. A. Keith and R. F. W. Bader, J. Chem. Phys., 1993, 99, 3669–3682.
25. P. Lazzeretti, M. Malagoli and R. Zanasi, Chem. Phys. Lett., 1994, 220, 299–304.
26. P. Lazzeretti, Theor. Chem. Acc., 2012, 131, 1222.
27. A. Soncini, A. M. Teale, T. Helgaker, F. De Proft and D. J. Tozer, J. Chem. Phys., 2008, 129, 74101.
28. A. Soncini, P. Lazzeretti and R. Zanasi, Chem. Phys. Lett., 2006, 421, 21–26.
29. A. D. Becke, Phys. Rev. A: At., Mol., Opt. Phys., 1988, 38, 3098–3100.
30. C. Lee, W. Yang and R. G. Parr, Phys. Rev. B: Condens. Matter Mater. Phys., 1988, 37, 785–789.
31. B. Miehlich, A. Savin, H. Stoll and H. Preuss, Chem. Phys. Lett., 1989, 157, 200–206.
32. F. Weigend and R. Ahlrichs, Phys. Chem. Chem. Phys., 2005, 7, 3297–3305.
33. M. J. Frisch, et al., Gaussian 09, Revision B.01, Gaussian Inc., Wallingford CT, 2009.
34. J. Jusélius, D. Sundholm and J. Gauss, J. Chem. Phys., 2004, 121, 3952–3963.
35. U. Ayachit, The ParaView Guide: A Parallel Visualization Application, Kitware, 2015.
36. E. Clar, The Aromatic Sextet, John Wiley & Sons Ltd, London, 1972.
37. E. Steiner, P. W. Fowler, A. Soncini and L. W. Jenneskens, Faraday Discuss., 2007, 135, 309–323.
38. E. Steiner, P. W. Fowler and L. W. Jenneskens, Angew. Chem., Int. Ed., 2001, 40, 362–366.
39. G. Monaco and R. Zanasi, Phys. Chem. Chem. Phys., 2013, 15, 17654–17657.
40. T. K. Dickens and R. B. Mallion, Phys. Chem. Chem. Phys., 2013, 15, 8245–8253.
41. M. Solà, Front. Chem., 2013, 1, 22.
42. S. Radenković, J. Tošović, R. W. A. Havenith and P. Bultinck, ChemPhysChem, 2015, 16, 216–222.
43. S. Radenković, I. Gutman and P. Bultinck, J. Phys. Chem. A, 2012, 116, 9421–9430.
44. G. V. Baryshnikov, R. R. Valiev, N. N. Karaush, D. Sundholm and B. F. Minaev, Phys. Chem. Chem. Phys., 2016, 18, 8980–8992.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: https://doi.org/10.1039/d5cp00638d

---