Partitioning Hückel–London Currents into Cycle Contributions

Abstract

Ring-current maps give a direct pictorial representation of molecular aromaticity. They can be computed at levels ranging from empirical to full ab initio and DFT. For benzenoid hydrocarbons, Hückel–London (HL) theory gives a remarkably good qualitative picture of overall current patterns, and a useful basis for their interpretation. This paper describes an implemention of Aihara’s algorithm for computing HL currents for a benzenoid (for example) by partitioning total current into its constituent cycle currents. The Aihara approach can be used as an alternative way of calculating Hückel–London current maps, but more significantly as a tool for analysing other empirical models of induced current based on conjugated circuits. We outline an application where examination of cycle contributions to HL total current led to a simple graph-theoretical approach for cycle currents, which gives a better approximation to the HL currents for Kekulean benzenoids than any of the existing conjugated-circuit models, and unlike these models it also gives predictions of the HL currents in non-Kekulean benzenoids that are of similar quality.

1. Introduction

Benzene was first isolated almost 200 years ago [1] and the term ‘aromatic’ came into use as a description for this and similar compounds soon afterwards [2]. Since Kekulé’s famous identification of the special structure of benzene [3], the importance, meaning and even existence of ‘aromaticity’ have been hotly debated, and these discussions show no sign of reaching a universally accepted conclusion [4,5,6,7,8,9,10,11,12]. However, one widely accepted working criterion for aromaticity is the manifestation in a cyclic system of global currents (ring currents) induced by application of an external magnetic field [13,14,15,16,17,18,19,20]. This definition of aromaticity appeals to the community of theoretical chemists who calculate molecular electric and magnetic response properties, and it has featured extensively in the scientific career of Riccardo Zanasi, from their early work with Paolo Lazzeretti in Modena, to their work over several decades with colleagues in Salerno. As a definition, it also has the desirable feature that the criterion is, at least in principle, clearcut: either there is a global current or not, and if there is one, it has a sense of circulation with respect to the axis of the external field, which leads to a natural division of (monocyclic) ring systems into disjoint aromatic, non-aromatic and anti-aromatic classes. This criterion is ideally suited to probing by theoretical methods that calculate induced current either directly, or via other response magnetic properties as proxies. The ring-current picture has a close connection to experiment, through ring-current effects on ${}^1$H NMR chemical shifts [16,17] and ‘exaltation of diamagnetism’ [13,14,15,21].

Over the last quarter of a century, the field has gained impetus from new possibilities for plotting physically realistic ab initio maps of the current density induced by an external magnetic field [22,23,24,25], and for interpreting these maps in terms of chemical concepts such as orbital energy, symmetry and nodal character [20,25]. Riccardo Zanasi has participated in all of these developments [26]. One paper from the Salerno group of particular relevance to the present topic is [27], where quantities from the Aihara model, to be discussed below, are used to aid interpretation of ab initio current maps.

In this paper, we concentrate on the oldest model for mapping induced currents in benzenoids and similar systems: Hückel–London (HL) theory [14,28], which can be formulated in several equivalent ways: as a finite-field method [29], a perturbation method based on bond-bond polarisabilities [30,31,32,33], or a treatment of current as the formal superposition of cycle contributions [34,35]. The purpose of the present paper is to draw attention to this third version of HL theory, which is associated with the name of the late Professor Jun-Ichi Aihara. His innovative reformulation of the HL problem has not always received the attention from other chemists that it deserves. Although the concepts that it generated, such as Topological Resonance Energy, Bond Resonance Energy and Magnetic Resonance Energy (TRE, BRE and MRE), are influential, it is rare to find examples of direct use by other chemists of the specifics of the method itself. This may be because the Aihara formalism employs a number of concepts from graph theory that are unfamiliar to most chemists, or because the defining equations are scattered over a long series of interlocking papers, so that their conversion to a workable algorithm has not always appeared straightforward.

Our aim here is to remedy this situation, by giving an explicit implementation. Our main motivation was not to calculate HL current maps (for which several easily implemented algorithms already exist), but to exploit the defining feature of Aihara’s approach: the emphasis on cycle contributions to current, where every cycle within the molecular graph, be it a chemical ring or larger, is taken into account. This feature has assumed new relevance over the last decade with the revival of interest in conjugated-circuit (CC) models [36,37,38,39,40,41]. A cycle C in a graph G is a conjugated circuit if both G and $G–C$ (the graph where all vertices of C and their associated edges have been deleted) have a perfect matching. In a CC model, each conjugated circuit contributes currents along its edges, with weights specific to the model [42].

Conjugated-circuit models have an attractive simplicity, but have crucial drawbacks for non-Kekulean systems, where they predict zero current, and for Kekulean systems with fixed bonds, where they predict ‘dead zones’ of vanishing current [43,44,45]. The current maps from conjugated-circuit models can be seen as approximate versions of HL current maps in which only certain ‘important’ cycles have been selected and given model-dependent weightings. The Aihara approach can be used as a toolkit to test these approximations, and to design better models.

Comparison of HL and CC currents in benzenoids by cycle size has allowed us to evaluate these selection and weighting schemes, and to propose a new model, also based on matchings, that gives an approximation to HL currents for both Kekulean and non-Kekulean benzenoids that is better than any of the published CC models [43]. The dual nature of HL theory as a graph theoretical method based on molecular-orbital theory, makes it interesting to compare HL results with conjugated-circuit models on the one hand, and with more sophisticated wavefunction and density functional approaches to electronic structure on the other. The relevance of the present graph-theoretical investigation to ab initio calculation is that HL currents already typically mimic pseudo-$\pi$ currents [43], which in turn are usually excellent mimics for current maps derived from full ab initio and density functional calculations. Some systematic exceptions to this statement are discussed in [43].

The symmetries and energies of HL molecular orbitals provide a useful basis for rationalising the frontier-orbital analysis of current maps obtained from ipsocentric calculations at these higher levels [20,25], even though HL and ipsocentric definitions of molecular-orbital contributions are markedly different. In delocalised $\pi$ systems, current maps calculated within the ipsocentric approach are dominated by the frontier orbitals. In contrast, as usually formulated, HL currents in these systems have significant contributions from lower-lying molecular orbitals.

Graph Theoretical Background

An undirected graph G consists of a set V of vertices and a set E of edges where each edge corresponds to an unordered pair of vertices from V. We use n to denote the number of vertices of a graph and m to denote the number of edges. A graph is planar if it can be drawn in the plane with no crossing edges. When traversing the faces of a graph, each edge $(u,v)$ is treated as the two arcs $(u,v)$ and $(v,u)$. A traversal of each face of the graph uses each arc exactly once.

The graphs considered in this paper are benzenoids. Benzenoids may be defined as simply connected subgraphs of the hexagonal lattice composed of edge-fused hexagons. Hence, they correspond to connected planar graphs having all internal faces of size 6. The vertices on the interior have degree 3. The vertices on the perimeter (external face) either have degree 2 or degree 3. As is well known, the $\pi$ systems of benzenoids support circulations of electrons induced by an external magnetic field with consequences for magnetic susceptibilities and ${}^1$H NMR chemical shifts [13,14,15,16,17,21]. The calculation of this magnetic response in HL theory requires an embedding of the molecular graph, with explicit coordinates for the atomic positions. The embedding used here for benzenoids idealises each carbon framework as planar and composed of regular hexagons of side $1.4$ Å, embedded without overlap in the hexagonal tessellation of the plane.

When representing current, the graph is converted to a directed graph. If there is a current of magnitude k on arc $(u,v)$ and a current of magnitude r on arc $(v,u)$ then the net current on arc $(u,v)$ is equal to $k-r$. A current of magnitude k on arc $(u,v)$ is equivalent to a current of magnitude $-k$ on arc $(v,u)$. In our depictions of currents, the current contributions and arc directions are shown so that all magnitudes are greater than or equal to zero. In our maps, diatropic currents, representing aromatic currents, are those in a counter-clockwise direction, and conversely paratropic currents, representing anti-aromatic currents, are those in a clockwise direction.

By convention, the ‘absolute’ currents obtained from HL theory are often reported on a scale where unit current is equal to the HL current along an edge of an isolated, neutral benzene ring with side length $1.4$ Å [46]. When comparing different models, it is more useful to consider scaled current, as empirical methods for approximating currents give relative and not absolute results. A scaled current is obtained from the current picture by dividing the current value of each edge by the maximum current value. Scaled currents have a current of one on each arc that bears maximum current.

2. The Hückel–London Model as a Superposition of Cycle Contributions

The Aihara formulation of Hückel–London theory was refined over a series of papers, and here we give the working equations needed for its implementation. As a practical check, our implementation was run on all the small benzenoids (both Kekulean and non-Kekulean) having up to ten hexagons and the computed results matched against HL currents from the standard finite-perturbation approach, giving computational verification that our interpretation of the equations is correct. Aihara’s basic formalism was presented in two papers from 1979 [34,35] in which the relationship to London’s approximations [14] was established. In London theory, the effect of an external magnetic field is to perturb the original Hückel secular matrix of the molecule, effectively converting the $+1$ entries in the adjacency matrix into exponentials that reduce to $+1$ in the limit of vanishing applied magnetic field. This gives an easily implemented finite-field version of HL theory, e.g., [29]. In contrast, the Aihara formalism is an analytic perturbation theory and hence the calculated current densities are simple functions of field-free characteristic polynomials [47].

The first step is to find the eigenvalues ${\lambda }_1,{\lambda }_2,\dots ,{\lambda }_n$ of the adjacency matrix $\mathbf{A}\left(G\right)$ of the graph G. The number of times that a value ${\lambda }_k$ appears as an eigenvalue is the multiplicity of ${\lambda }_k$, denoted by $m_k$. The multiplicity of the zero eigenvalue is the nullity of the graph, $\eta$. The characteristic polynomial, $P_G\left(x\right)$, for a graph G is equal to

$$P_G\left(x\right)=|x\mathbf{1}-\mathbf{A}\left(G\right)|=\prod _{k=1}^n(x-{\lambda }_k),$$

(1)

where $\mathbf{1}$ is the $n\times n$ identity matrix. If a graph has no vertices, then the characteristic polynomial is 1.

In the Hückel model, eigenvectors of the adjacency matrix correspond to molecular orbitals, and eigenvalues correspond to orbital energies. It is usual to choose $\alpha$ for the origin of the energy scale and $\left|\beta \right|$ for the energy unit, where $\alpha$ and $\beta$ are the (negative) Coulomb and Resonance integrals from Hückel theory. The energy of an electron occupying one of the shell of $m_k$ degenerate orbitals that have eigenvalue ${\lambda }_k$ in the field-free $\pi$-system is then $\alpha +{\lambda }_k\beta$, giving the correspondence between values ${\lambda }_k>0$, ${\lambda }_k=0$, and ${\lambda }_k<0$ and the bonding, non-bonding or antibonding character of the shell, respectively. Electrons are assigned to orbitals using the Aufbau and Pauli Principles and Hund’s Rule of Maximum Multiplicity. In brief: orbitals are filled sequentially in non-increasing order of eigenvalues (Aufbau), each orbital has a maximum occupancy of two electrons (Pauli), and where $m_k>1$, no orbital receives a second electron until each of the orbitals in the given shell contains at least one (Hund).

In practice, for calculations of HL currents the fractional occupation approximation is used, where the average number of electrons per orbital over a given shell is taken as the occupation number for every orbital in that shell. For the shell with eigenvalue ${\lambda }_k$, this occupation number is denoted by $n_k$. In cases where all $n_k$ are equal to 2 or 0, we have a closed shell; otherwise there is partial occupation of some orbital or orbitals and we have an open shell.

As benzenoids have bipartite molecular graphs, their Hückel eigenvalues and eigenvectors are subject to the Coulson-Rushbrooke Pairing Theorem [48]. If a bipartite graph has an eigenvalue satisfying $|{\lambda }_k|>0$, then the spectrum of the graph includes shells with eigenvalues ${\lambda }_k$ and $-{\lambda }_k$, with the same multiplicity. Moreover, to each of the $m_k$ orthonormal eigenvectors for the ${\lambda }_k$ shell there is an eigenvector in the $-{\lambda }_k$ shell derived by swapping the signs of all entries for vertices belonging to one partite set. Hence, a bipartite molecular graph has equal numbers of bonding and anti-bonding molecular orbitals, and the number of non-bonding orbitals ($\eta$) has the same parity as n, the number of conjugated carbon centres. Application of Aufbau, Pauli and Hund rules for the neutral molecule leads to an electron configuration with all bonding orbitals full ($n_k=2$), all anti-bonding orbitals empty ($n_k=0$), and any non-bonding orbitals half full ($n_k=1$).

We can simplify the calculation of HL current for a benzenoid by noting that an electron in a non-bonding orbital of a bipartite system makes no contribution to the total current (see Section 3), so we may deal with the equivalent closed-shell system where any non-bonding occupancies have been set to zero. This closed shell corresponds to the neutral molecule for a Kekulean benzenoid, and otherwise to the cation of charge $+\eta$. This simplification is specific to bipartite systems treated within the HL model. In ipsocentric ab initio approaches [20,25,49], non-bonding orbitals may contribute current, whether or not the molecular graph is bipartite. Consideration of the nearest closed-shell for a non-Kekulean benzenoid has a direct chemical application in the study of benzenoid cations; a modified CC approach has been used to give an account of the current maps for these species [50].

Characteristic polynomials are manipulated in order to calculate a resonance energy for each cycle C. The fundamental quantities in the Aihara analysis are Circuit Resonance Energies (CRE), denoted by the symbol $A_C$ (originally RE${}^{\left(C\right)}$ [51]) where C is any cycle in G, conjugated or not. The CRE represents a partition of the Topological Resonance Energy into cycle contributions, interpretable as a gain or loss in energy arising from conjugation along a cycle C [52,53]. It is defined in dimensionless form as [54]:

$$A_C=2\sum _k^{}{n}_k{f}_k\left({\lambda }_k\right),$$

(2)

where $n_k$ is the (shell average) number of electrons in the orbital with eigenvalue ${\lambda }_k$, $f_k$ is the function described below, and the sum in (2) runs over distinct values ${\lambda }_k$ (each eigenvalue counted once, irrespective of multiplicity). The value $A_C$ can be compared to the standard case of neutral benzene, for which $A_C=2/9$ (see Section 5.1).

The function $f_k\left({\lambda }_k\right)$ is defined in terms of the characteristic polynomials of graphs G and $G–C$, and is described in different ways for cases with $m_k=1$ and $m_k\ne 1$. In the simple case where $m_k$ = 1, the function $f_k$ is

$$f_k\left(x\right)=\frac{P_{G–C}\left(x\right)}{U_k\left(x\right)},$$

(3)

where the polynomial $U_k\left(x\right)$ is defined as

$$U_k\left(x\right)=\frac{P_G\left(x\right)}{{(x-{\lambda }_k)}^{m_k}}.$$

(4)

For $m_k=1$, $U_k\left({\lambda }_k\right)$ is equal to the first derivative of $P_G\left(x\right)$ evaluated at ${\lambda }_k$, i.e.,

$$U_k\left({\lambda }_k\right)={P_G}^{\prime }\left({\lambda }_k\right).$$

(5)

In systems with degenerate orbitals, the contribution from these orbitals to the CRE must account for the splitting induced by the external magnetic field [55]. To do so, the function $f_k$ is defined as

$$f_k\left(x\right)=\frac{1}{(m_k-1)!}{\left\{\frac{d^{m_k-1}}{d{x}^{m_k-1}}\frac{P_{G–C}\left(x\right)}{U_k\left(x\right)}\right\}}_{x={\lambda }_k},$$

(6)

If $d^0/d{x}^0$ is taken to be the identity, (3) is the formal limit of (6) for $m_k=1$.

The $A_C$ value for a given cycle can be converted into a Hückel–London cycle current, $J_C$, by accounting for the area of the cycle [56]. The cycle contribution to the total current-density map is defined, again in dimensionless form, as

$$J_C=\left(\frac{9}{2}\right)A_C{S}_C,$$

(7)

Written in this way, the equation gives the cycle contribution as a multiple of the unit HL current for neutral benzene [47]. The quantity $S_C$ is the area of cycle C in terms of benzene rings. In benzenoids, $S_C$ is therefore simply the number of hexagons enclosed by the cycle. In non-benzenoids, $S_C$ is the area of the cycle normalised to that of a hexagon, i.e., the faces inside the cycle are considered to be regular polygons and their areas are summed and divided by the area of a regular hexagon with the same side length. Hence, each polygonal ring of size p that is enclosed by cycle C contributes $p\sqrt{3}cot(\pi /p)/18$ to $S_C$.

The HL current-density map for a benzenoid is obtained edge by edge by summing contributions from all cycles that pass through the given edge to assign the bond current. A more compact representation uses ring currents assigned to the faces; current on a perimeter edge equates to the ring current for the face containing it; current on an interior edge is the vector sum of the ring currents flowing in the two faces that meet along that edge.

We denote the ring current of a face F by ${\widehat{J}}_F$. Note that the ring current on a face in a polycyclic system is not in general equal to the current contribution for that same cycle as given by the Aihara formula (7). The two are of course equal for benzene, and unscaled ring currents in dimensionless form are therefore also specified as ratios to the benzene value.

The sum of $A_C$ values over all cycles is used to define as a proposed aromaticity index, the magnetic resonance energy (MRE) of G [57]:

$$MRE=\sum _C{A}_C.$$

(8)

Aihara has argued that this index has a physical advantage over raw ring current as it is independent of cycle area, whereas ring currents are not. One of their most recent papers [58] is an encyclopaedic survey of the magnetic criteria of aromaticity, in which he concludes that MRE is for many purposes an ideal aromaticity index. This paper also gives a good working summary of all the basic equations of the Aihara approach.

A third cycle property related to aromaticity on the magnetic definition is the magnetic susceptibility of a cycle C, ${\chi }_C$, which has an even stronger dependence on cycle area and is defined, again in dimensionless form referred to the susceptibility of benzene (which is diamagnetic and therefore negative) as [35]:

$${\chi }_C=-\left(\frac{9}{2}\right)A_C{\left(S_C\right)}^2.$$

(9)

The $\pi$-electron contribution to the molecular magnetic susceptibility, $\chi$, is obtained by summing ${\chi }_C$ over all cycles.

Thus, the three quantities of circuit resonance energy ($A_C$), cycle current, ($J_C$), and cycle magnetic susceptibility (${\chi }_C$) all contain the same information, weighted differently. Aihara’s objection to the use of ring currents as a measure of aromaticity also applies to the magnetic susceptibility. A related point was made by Estrada [59], who argued that correlations between magnetic and energetic criteria of aromaticity for some molecules could simply be a result of underlying separate correlations of susceptibility and resonance energy with molecular weight.

3. A Pairing Theorem for HL Currents

As noted above, bipartite graphs obey the Pairing Theorem [48]. The theorem implies that when the eigenvalues of a bipartite graph are arranged in non-increasing order from ${\lambda }_1$ to ${\lambda }_n$, positive and negative eigenvalues are paired, with

$${\lambda }_k=-{\lambda }_{\overline{k}},$$

(10)

where $\overline{k}$ is shorthand for $n-k+1$. If $\eta$ is the number of zero eigenvalues of the graph, $n-\eta$ is even. Zero eigenvalues occur at positions ranging from $k=(n-\eta )/2+1$ to $k=(n+\eta )/2$. HL currents for benzenoids and other bipartite molecular graphs also obey a pairing theorem, as is easily proved using the Aihara Formulas (2)–(7),

We consider arbitrary elctron counts and occupations of the shells. Each electron in an occupied orbital with eigenvalue ${\lambda }_k$ makes a contribution $2f_k\left({\lambda }_k\right)$ to the Circuit Resonance Energy $A_C$ of cycle C (Equation (2)). The function $f_k\left({\lambda }_k\right)$ depends on the multiplicity $m_k$: it is given by Equation (3) for non-degenerate ${\lambda }_k$ and Equation (6) for degenerate ${\lambda }_k$.

Theorem 1.

For a benzenoid graph, the contributions per electron of paired occupied shells to the Circuit Resonance Energy of cycle C, $A_C$, are equal and opposite, i.e.,

$$f_k\left({\lambda }_k\right)=-f_k\left({\lambda }_{\overline{k}}\right).$$

(11)

Proof. 

The result follows from parity of the polynomials used to construct $f_k\left({\lambda }_k\right)$. The characteristic polynomial for a bipartite graph has well defined parity, as

$$P_G\left(x\right)={(-1)}^n{P}_G(-x).$$

(12)

On differentiation the parity reverses:

$$P_G{}^{\prime }\left(x\right)={(-1)}^{n-1}{P}_G{}^{\prime }(-x).$$

(13)

A benzenoid graph is bipartite, so all cycles C are of even size and $P_{G–C}\left(x\right)$ has the same parity as $P_G\left(x\right)$:

$$P_{G–C}\left(x\right)={(-1)}^n{P}_{G–C}(-x).$$

(14)

Therefore, for $m_k=1$,

$$f_k\left({\lambda }_k\right)=f_k\left(x\right){|}_{x={\lambda }_k}=\frac{P_{G–C}\left(x\right)}{P_G{}^{\prime }\left(x\right)}{|}_{x={\lambda }_k}=-\frac{P_{G–C}(-x)}{P_G{}^{\prime }(-x)}{|}_{-x={\lambda }_{\overline{k}}}=-f_k\left({\lambda }_{\overline{k}}\right).$$

(15)

The argument for the case for $m_k>1$ is similar. For a bipartite graph, the parity of $P_G\left(x\right)$ can equally be stated in terms of order or nullity:

$$P_G\left(x\right)=x^{\eta }\prod _{{\lambda }_k^{\prime }>0}(x^2-{{\lambda }_k^{\prime }}^2)={(-1)}^{\eta }{P}_G(-x).$$

(16)

Functions $U_k\left(x\right)$ and $U_{\overline{k}}(-x)$ are therefore related by

$$U_k\left(x\right)={(-1)}^{m_k+\eta }{U}_{\overline{k}}(-x),$$

(17)

as $P_G\left(x\right)={(-1)}^{\eta }{P}_G(-x)$ and $U_k\left(x\right)$ and $U_{\overline{k}}(-x)$ are formed by cancelling $m_k$ factors $(x-{\lambda }_k)$ and $(-x-{\lambda }_{\overline{k}})=(-1)(x+{\lambda }_k)$, respectively, from $P_G\left(x\right)$. Hence, the quotient function $P_{}\left(x\right)/U_k\left(x\right)$ behaves as

$$\frac{P_{G–C}\left(x\right)}{U_k\left(x\right)}={(-1)}^{n-m_k-\eta }\frac{P_{G–C}(-x)}{U_{\overline{k}}(-x)}.$$

(18)

Each differentiation flips the parity, and the pairing result for $m_k>1$ is therefore

$$f_k\left({\lambda }_k\right)={(-1)}^{n-m_k-\eta +m_k-1}{f}_{\overline{k}}\left({\lambda }_{\overline{k}}\right)=-f_{\overline{k}}\left({\lambda }_{\overline{k}}\right).$$

(19)

   □

Some straightforward corollaries are:

Corollary 1.

In the fractional occupation model, where all orbitals of a shell are assigned equal occupation, paired shells of a bipartite graph that contain the same number of electrons make cancelling contributions of current for every cycle C, and hence no net contribution to the HL current map.

Corollary 2.

In the fractional occupation model, all electrons in a non-bonding shell of a bipartite graph (${\lambda }_k={\lambda }_{\overline{k}}=0$) make no contribution to any cycle current $J_C$ and hence make no net contribution to the HL current map.

It should be noted that if a graph is non-bipartite, the non-bonding shell may contribute a significant current in the HL model. Furthermore, if G is bipartite but subject to first-order Jahn-Teller distortion, current may arise from the occupied part of an originally non-bonding shell; this can be treated by using the form of the Aihara model appropriate to edge-weighted graphs [58]. Corollary (2) also highlights a significant difference between HL and ipsocentric ab initio methods. In the latter, an occupied non-bonding molecular orbital of an alternant hydrocarbon can make a significant contribution to total current through low-energy virtual excitations to nearby ${\pi }^{*}$ shells, and can be a source of differential $\alpha$ and $\beta$ currents.

Corollary 3.

In the fractional occupation model, the HL current maps for the $q+$ cation and $q-$ anion of a π system that has a bipartite molecular graph are identical.

We can also note that in the extreme case of the cation/anion pair where the neutral system has gained or lost a total of n$\pi$ electrons, the HL current map has zero current everywhere. For bipartite graphs, this follows from Corollary (3), but it is true for all graphs, as a consequence of the perturbational nature of the HL model, where currents arise from field-induced mixing of unoccupied into occupied orbitals: when either set is empty, there is no mixing.

4. Implementation of the Aihara Method

4.1. Generating All Cycles of a Planar Graph

By definition, conjugated-circuit models consider only the conjugated circuits of the graph. In contrast, the Aihara formalism considers all cycles of the graph. A catafused benzenoid (or catafusene) has no vertex belonging to more than two hexagons. Catafusenes are Kekulean. For catafusenes, all cycles are conjugated circuits. All other benzenoids have at least one vertex in three hexagons, and have some cycles that are not conjugated circuits.

The size of a cycle is the number of vertices in the cycle. The area of a cycle C of a benzenoid is the number of hexagons enclosed by the cycle. One way to represent a cycle of the graph is with a vector $[e_1,e_2,\dots e_m]$ which has one entry for each edge of the graph where $e_i$ is set to one if edge i is in the cycle, and is set to 0 otherwise. When we add these vectors together, the addition is done modulo two. The addition of two cycles of the graph can either result in another cycle, or a disconnected graph whose components are all cycles.

A cycle basis B of a graph G is a set of linearly independent cycles (none of the cycles in B is equal to a linear combination of the other cycles in B) such that every cycle of the graph G is a linear combination of the cycles in B. It is well known that the set of faces of a planar graph G is a cycle basis for G [60]. The approach that we use for generating all the cycles starts with this cycle basis and finds the remaining cycles by taking linear combinations.

The cycles of a benzenoid that have unit area are the faces. The cycles that have area $r+1$ are generated from those of area r by considering the cycles that result from adding each cycle of area one to each of the cycles of area r. If the result is connected and is a cycle that is not yet on the list, then this new cycle is added to the list.

For the Aihara approach, a counterclockwise representation of each cycle of the graph is computed. It is easy to compute these as the cycles are generated. A face traversal algorithm [61] first provides the internal faces as traversed in counterclockwise order. If a new cycle $C_3$ is a linear combination of $C_1$ and $C_2$ then arcs that are in both $C_1$ and $C_2$ disappear and the remaining arcs should be oriented in the same way as they are in the cycle from which they came.

4.2. Efficient Computation of Necessary Derivatives

The derivative of a function f with respect to x is denoted here as $f^{\prime }\left(x\right)$. We first recall some elementary properties of the derivative. For a polynomial $p\left(x\right)$ of degree n that is equal to ${\sum }_{i=0}^n{c}_i{x}^i$, the derivative $p^{\prime }\left(x\right)$ is equal to ${\sum }_{i=1}^n{c}_{i}i{x}^{i-1}$. The product rule for a function $f\left(x\right)=p\left(x\right)q\left(x\right)$ states that $f^{\prime }\left(x\right)=p\left(x\right)q^{\prime }\left(x\right)+p^{\prime }\left(x\right)q\left(x\right)$. The quotient rule for a function $f\left(x\right)=p\left(x\right)/q\left(x\right)$ states that $f^{\prime }\left(x\right)=[q\left(x\right)p^{\prime }\left(x\right)-p\left(x\right)q^{\prime }\left(x\right)]/\left(q{\left(x\right)}^2\right)$.

In the set of small benzenoids we used for initial testing (Kekuléan benzenoids with at most seven hexagons) the maximum multiplicity of an eigenvalue is four (implying that the differentiation in the formula for $f_k\left(x\right)$ (Equation (6)) has to be applied three times). If the quotient rule is applied directly without additional simplification, then the degree of the denominator polynomial doubles. For example, starting with a polynomial of degree 30, results of one of degree 60. Differentiating a second time gives degree 120, and the third differentiation gives degree 240. Polynomials of such large degree resulted in numerical instability in the computations. In order to correct this problem, we changed the way that the differentiation was implemented. The new approach is as follows.

In the formula for $f_k\left(x\right)$ the two polynomials can each be expressed in the form ${\prod }_{i=1}^d(x-{\lambda }_i)$. For the numerator, $P_{G–C}\left(x\right)$, the ${\lambda }_i$ values are just the eigenvalues of $G–C$. For the denominator, $U_k\left(x\right)$, they correspond to the eigenvalues of G with each of the $m_k$ occurrences of an eigenvalue equal to ${\lambda }_k$ excluded. For a polynomial $p\left(x\right)=(x-{\alpha }_1)(x-{\alpha }_2)\dots (x-{\alpha }_{d_p})$ we use the notation $p\left(x\right)[-i]$ to denote the polynomial $(x-{\alpha }_1)(x-{\alpha }_2)\dots (x-{\alpha }_{d_p})/{(x-{\alpha }_i)}^{m_i}$ or in equivalent product form, $p\left(x\right)$ with the $m_i$ terms of the form $(x-{\alpha }_i)$ crossed out. (Eigenvalues ${\alpha }_i,{\beta }_i$ are not to be confused with the Hückel integrals $\alpha$, $\beta$.)

Suppose that the function that we want to differentiate is $g\left(x\right)=p\left(x\right)/q\left(x\right)$ for polynomials p and q with degrees $d_p$ and $d_q$, $p\left(x\right)=(x-{\alpha }_1)(x-{\alpha }_2)\dots (x-{\alpha }_{d_p})$ and $q\left(x\right)=(x-{\beta }_1)(x-{\beta }_2)\dots (x-{\beta }_{d_q})$. Applying quotient and product rules and cancelling out common terms in numerator and denominator gives this formula for $g^{\prime }\left(x\right)$:

$$g^{\prime }\left(x\right)=\frac{p\left(x\right)}{q\left(x\right)}\left[\sum _{i=1}^{d_p}\frac{1}{(x-{\alpha }_i)}-\sum _{j=1}^{d_q}\frac{1}{(x-{\beta }_j)}\right].$$

(20)

Note that, with this approach, the maximum degree increases by one each time instead of doubling. This results in better numerical stability.

For computing $f_k\left({\lambda }_k\right)$, it is not necessary to use a data structure that represents polynomials. Instead, vectors can be used. The recursive algorithm given below evaluates $f_k$ at $x={\lambda }_k$. The vectors (indexed starting from 0) are p[i]$=x-{\alpha }_{i+1}$ and q[i]$=x-{\beta }_{i+1}$. These are used to compute derivatives instead of computing characteristic polynomials explicitly.

The function eval_deriv differentiates $p\left(x\right)/q\left(x\right)$ power times, where the argument x at which to evaluate the derivative has already been chosen and the vectors have been pre-computed. If power is 0 then the answer is the product of the values p[0] to p[$d_p$-1] divided by the product of the values q[0] to q[$d_q$-1]. Otherwise the answer is computed from: $f^{\prime }\left(x\right)={\sum }_{i=1}^{d_p}\{p\left(x\right)[-i]/q\left(x\right)\}-{\sum }_{j=1}^{d_q}\{p\left(x\right)/\left[q\left(x\right)(x-i{\beta }_j)\right]\}$. MAX_DEG is the maximum degree of any polynomial.

double eval_deriv(int power, int dp, double p[MAX_DEG], int dq, double q[MAX_DEG])

{

   double r[MAX_DEG];

   double ans, top, bottom;

   int limit, pos, i, j;

// When power is 0, stop taking derivatives and evaluate.

   if (power == 0)

   {

      if (dp < dq) limit = dq;

      else limit = dp;

      ans = 1;

// The answer is the product of the p values divided by the product of the q values.

      for (i = 0; i < limit; i++)

      { 

         if (i < dp) top = p[i];

         else top = 1;

         if (i < dq) bottom = q[i];

         else bottom = 1;

         ans* = (top/bottom);

      }

      return(ans);

   }

   ans = 0;

// Compute qp’ / q^2 = p’/q.

// Ignore if dp = 0 since a polynomial of degree 0 has a derivative of 0. 

   if (dp > 0)

   {

//    If dp = 1 then the polynomial is x-a0 and the derivative of this is 1.

      if (dp == 1)

      {

         r[0] = 1;

         ans+ = eval_deriv(power-1, dp-1, r, dq, q);

      }

      else // dp > 1.

      {

         for (i = 0; i < dp; i++)

         {

//          Compute p(x)[-i]:

            pos = 0;

            for (j=0; j < dp; j++)

            {

               if (i != j)

               {

                  r[pos] = p[j];

                  pos++;

               }

            }

            ans+ = eval_deriv(power-1, dp-1, r, dq, q);

         }

      }

   }

// Now subtract off p q’ / q^2

   for (i = 0; i < dq; i++)

      r[i] = q[i];

   for (i = 0; i < dq; i++)

   {

      r[dq] = q[i];

      ans -= eval_deriv(power-1, dp, p, dq+1, r);

   }

   return(ans);

}

5. Some Examples of the Aihara Model

5.1. The Basic Case: Benzene

Benzene is the standard against which aromaticity of other molecules is judged, and is invoked in the dimensionless formulation of the Aihara Equations (2)–(9). For benzene, the characteristic polynomial and its derivative are

$$P_G\left(x\right)=(x^2-4){(x^2-1)}^2,$$

(21)

$$P_G{\left(x\right)}^{\prime }=6x(x^2-3)(x^2-1).$$

(22)

As benzene is a monocycle, $P_{G–C}\left(x\right)=1$. The eigenvalues are $\{+2,+1,+1,-1,-1,-2\}$, with occupation numbers in the neutral $6\pi$ system of $\{2,2,2,0,0,0\}$. Hence, the first shell has ${\lambda }_1=2$ and $n_1=2$ and, by (3),

$$f_1\left(2\right)={\left\{\frac{1}{P_G^{\prime }\left(x\right)}\right\}}_{x=+2}=\frac{1}{36}$$

(23)

and the second shell ${\lambda }_2=1$ and $n_2=2$ and, by (6),

$$f_2\left(1\right)={\left\{\frac{d}{dx}\frac{1}{{(x^2-4)}^{}{(x+1)}^2}\right\}}_{x=+1}=\frac{1}{36}.$$

(24)

Therefore, by (2), $A_C=2/9$. As $S_C=1$, the cycle contribution to current, which in this case is also the ring current, is 1 (by (7), and the (diamagnetic) susceptibility is $-1$. The value of $A_C$ for benzene is the reason for the factors of $9/2$ in the other Aihara equations.

Notice that in the HL model half of the ring current arises from the $2\pi$ LOMO and half from the $4\pi$ HOMO, in contrast to the ipsocentric picture where essentially the whole of the $\pi$ current arises from the HOMO [20].

5.2. An Analytical Example: The HL Current in Anthracene

Our approach is computational, but it is also interesting for interpretation purposes to see how the various quantities in the Aihara cycle decomposition of HL current can be worked out fully analytically in a simple case. The characteristic polynomial for anthracene is

$$\begin{array}{cc}\hfill P_G\left(x\right)& =x^{14}-16x^{12}+98x^{10}-296x^8+473x^6-392x^4+148x^2-16\hfill \\ \hfill & =(x-2)(x+2)\left(x^2+2x-1\right)\left(x^2-2x-1\right){(x-1)}^2{(x+1)}^2{\left(x^2-2\right)}^2,\hfill \end{array}$$

(25)

the roots of which are the eigenvalues of the adjacency matrix of the graph, split equally between bonding and anti-bonding shells. As anthracene is a catafusene, the graph is Kekulean and there are no non-bonding orbitals. The occupied orbitals of neutral anthracene correspond to eigenvalues $\left\{(1+\sqrt{2}),2,\sqrt{2},\sqrt{2},1,1,(-1+\sqrt{2})\right\}$. The unoccupied orbitals correspond to eigenvalues $\left\{(1-\sqrt{2}),-1,-1,-\sqrt{2},-\sqrt{2},-(1+\sqrt{2})\right\}$. Thus, anthracene has doubly degenerate pairs of orbitals at $\pm \sqrt{2}$ and $\pm 1$. In the Aihara formalism, every cycle within the graph is considered. For anthracene there are six possible cycles. Three are the individual hexagonal faces, two result from the naphthalene-like fusion of two hexagonal faces, and the final cycle is the result of the fusion of all three hexagonal faces. The cycles and corresponding polynomials $P_{G–C}\left(x\right)$ are displayed in

Table 1. Cycles and corresponding polynomials $P_{G–C}\left(x\right)$ in anthracene. Bold lines represent edges in C; removal of bold and dashed lines yields the graph $G–C$.

Cycle	Cycle Diagram	${\mathit{P}}_{\mathit{G}-\mathit{C}}\left(\mathit{x}\right)$
$C_1$		$\left(x^4-2x^3-2x^2+3x+1\right)\left(x^4+2x^3-2x^2-3x+1\right)$
$C_2$		$\left(x^4-2x^3-2x^2+3x+1\right)\left(x^4+2x^3-2x^2-3x+1\right)$
$C_3$		${\left(x^2+x-1\right)}^2{\left(x^2-x-1\right)}^2$
$C_4$		$\left(x^2+x-1\right)\left(x^2-x-1\right)$
$C_5$		$\left(x^2+x-1\right)\left(x^2-x-1\right)$
$C_6$		1

.

Individual circuit resonance energies, $A_C$, can now be calculated using Equation (2). For all occupied orbitals, $n_k=2$. Calculations can be reduced by accounting for symmetry-equivalent cycles. For anthracene, six calculations of $A_C$ reduce to four as $A_1$ = $A_2$ and $A_4$ = $A_5$.

First, the functions $f_k$ must be calculated for each cycle. For those eigenvalues with $m_k=1$, $f_k$ is calculated using Equation (3), where the appropriate form of $U_k\left(x\right)$ can be deduced from the factorised characteristic polynomial in Equation (25). For those occupied eigenvalues with $m_k=2$, $f_k$ is calculated using a single differentiation in Equation (6). This procedure yields the $A_C$ values in

Table 2. Circuit resonance energy (CRE) values, $A_C$, calculated using Equation (2) for cycles of anthracene. Cycles are labelled as shown in Table 1.

CRE	Formula	Value
$A_1=A_2$	$\left(\frac{53+38\sqrt{2}}{2128+1512\sqrt{2}}+\frac{19}{252}+\frac{-83\sqrt{2}}{392}+\frac{13}{36}+\frac{53-38\sqrt{2}}{1512\sqrt{2}-2128}\right)=\left(\frac{55}{126}-\frac{12\sqrt{2}}{49}\right)$	$\approx 0.0902$
$A_3$	$\left(\frac{153+108\sqrt{2}}{2128+1512\sqrt{2}}+\frac{-25}{252}+\frac{-113\sqrt{2}}{392}+\frac{17}{36}+\frac{153-108\sqrt{2}}{1512\sqrt{2}-2128}\right)=\left(\frac{47}{126}-\frac{43\sqrt{2}}{196}\right)$	$\approx 0.0628$
$A_4=A_5$	$\left(\frac{9+6\sqrt{2}}{2128+1512\sqrt{2}}+\frac{-5}{252}+\frac{85\sqrt{2}}{392}+\frac{-11}{36}+\frac{9-6\sqrt{2}}{1512\sqrt{2}-2128}\right)=\left(\frac{25\sqrt{2}}{98}-\frac{41}{126}\right)$	$\approx 0.0354$
$A_6$	$\left(\frac{1}{2128+1512\sqrt{2}}+\frac{-1}{252}+\frac{-57\sqrt{2}}{392}+\frac{5}{36}+\frac{1}{1512\sqrt{2}-2128}\right)=\left(\frac{17}{126}-\frac{15\sqrt{2}}{196}\right)$	$\approx 0.0267$

.

Circuit resonance energies, $A_C$, are converted to cycle current contributions, $J_C$, by Equation (7). These results are summarised in

Table 3. Cycle currents, $J_C$, in anthracene calculated using Equation (7) with areas $S_C$, and values $A_C$ from Table 2. Currents are given in units of the ring current in benzene. Cycles are labelled as shown in Table 1.

Cycle Current	Area, ${\mathit{S}}_{\mathit{C}}$	Formula	Value
$J_1=J_2$	1	$\left(\frac{55}{28}-\frac{54\sqrt{2}}{49}\right)$	≈0.4058
$J_3$	1	$\left(\frac{47}{28}-\frac{387\sqrt{2}}{392}\right)$	≈0.2824
$J_4=J_5$	2	$\left(\frac{225\sqrt{2}}{98}-\frac{41}{14}\right)$	≈0.3183
$J_6$	3	$\left(\frac{51}{28}-\frac{405\sqrt{2}}{392}\right)$	≈0.3603

.

The significance of these quantities for interpretation is that they allow us to rank the contributions to the total HL current, and see that even in this simple case there are different factors in play. Notice that the contributions $J_1$ and $J_3$ are not equal. The two cycles have the same area, and correspond to graphs $G–C$ with the same number of perfect matchings, so would contribute equally in a CC model. In the Aihara partition of the HL current, the largest contribution from a cycle is from a face ($J_1$ for the terminal hexagon), but so is the smallest ($J_3$ for the central hexagon). The contributions of the cycles that enclose two and three faces are boosted by the area factors $S_C$, in accord with Aihara’s ideas on the difference in weighting between energetic and magnetic criteria of aromaticity [57].

Finally, the ring currents in the terminal and central hexagonal faces of anthracene are calculated. They are listed in

Table 4. Ring currents, ${\widehat{J}}_F$, for the terminal and central rings of anthracene, calculated using the cycle currents from Table 3. Currents are given in units of the ring current in benzene. Cycles are labelled as shown in Table 1.

Face	Contribution	${\widehat{\mathit{J}}}_{\mathit{F}}$
Terminal hexagon	$J_1+J_4+J_6=J_2+J_5+J_6$	$\left(\frac{6}{7}+\frac{9\sqrt{2}}{56}\right)$$\approx 1.0844$
Central hexagon	$J_3+J_4+J_5+J_6$	$\left(\frac{18\sqrt{2}}{7}-\frac{33}{14}\right)$$\approx 1.2794$

and displayed in maps of ring and bond currents in Figure 1. As they must, the currents correspond exactly to the results of the finite-field numerical Hückel–London approach. Note that now the largest bond and ring currents appear in the central hexagon, not in the terminal hexagons. Although the local cycle contribution $J_1$ is larger than $J_2$, the ring current in the central hexagon has contributions from more of the large cycles. The same effect is seen in CC models. The profile of increasing ring current from the ends to the middle of a linear polyacene chain is also seen in ab initio calculations. It has given rise to the so-called ‘anthracene problem’ [42,62], which is seen as a difficulty for theories of local aromaticity, in itself a contentious concept.

5.3. A Numerical Example: An Non-Kekulean Case

As an illustration of how the Aihara version of the HL model deals with non-Kekulean benzenoids, we take the 5-ring dibenzo-derivative of phenalenyl that is shown as (I) in Figure 2a.

The graph (though not necessarily the molecule) has $C_{2v}$ symmetry, and three symmetry-distinct hexagons, $F_1$, $F_2$, and $F_3$, where the last two are related by symmetry to their images $F_2^{\prime }$ and $F_3^{\prime }$. The five hexagonal faces generate 19 cycles, which give 12 distinct cases, up to isomorphism, as listed in

Table 5. Cycle contributions to HL current in the non-Kekulean benzenoid I. D and P stand for diatropic and paratropic contributions, respectively.

Cycle	Size	${\mathit{S}}_{\mathit{c}}$	Composition	${\mathit{J}}_{\mathit{C}}$	Tropicity
$C_1$	6	1	F${}_1$	$+0.0795$	D
$C_2$	6	1	F${}_2$ ≅ F${}_2^{\prime }$	$+0.0852$	D
$C_3$	6	1	F${}_3$ ≅ F${}_3^{\prime }$	$+0.2386$	D
$C_4$	10	2	F${}_1$ + F${}_2$ ≅ F${}_1$ + F${}_2^{\prime }$	$+0.0795$	D
$C_5$	10	2	F${}_2$ + F${}_2^{\prime }$	$+0.0227$	D
$C_6$	10	2	F${}_2$ + F${}_3$ ≅ F${}_2^{\prime }$ + F${}_3^{\prime }$	$+0.1705$	D
$C_7$	12	3	F${}_1$ + F${}_2$ + F${}_2^{\prime }$	$-0.0170$	P
$C_8$	14	3	F${}_1$ + F${}_2$ + F${}_3$ ≅ F${}_1$ + F${}_2^{\prime }$ + F${}_3^{\prime }$	$+0.1193$	D
$C_9$	14	3	F${}_2$ + F${}_2^{\prime }$ + F${}_3$ ≅ F${}_2$ + F${}_2^{\prime }$ + F${}_3^{\prime }$	$+0.0341$	D
$C_{10}$	16	4	F${}_1$ + F${}_2$ + F${}_2^{\prime }$ + F${}_3$ ≅ F${}_1$ + F${}_2$ + F${}_2^{\prime }$ + F${}_3^{\prime }$	$-0.0227$	P
$C_{11}$	18	4	F${}_2$ + F${}_2^{\prime }$ + F${}_3$ + F${}_3^{\prime }$	$+0.0455$	D
$C_{12}$	20	5	F${}_1$ + F${}_2$ + F${}_2^{\prime }$ + F${}_3$ + F${}_3^{\prime }$	$-0.0284$	P

along with their respective contributions to current.

Collecting contributions, the ring currents in the unscaled map are ${\widehat{J}}_{F_1}=0.3864$, ${\widehat{J}}_{F_2}=0.5000$ and ${\widehat{J}}_{F_3}=0.5568$. Scaled to the maximum bond current, the ring currents are ${\widehat{J}}_{F_1}=0.6939$, ${\widehat{J}}_{F_2}=0.8980$ and ${\widehat{J}}_{F_3}=1.0000$. All are positive and hence diatropic, but arise from different balances of three terms: (i) the local contribution from the face itself (strongest for $F_3$), (ii) the diatropic contribution from the other cycles of size 2 mod 4 (strongest for face $F_2$) (iii) the summed paratropic contribution from the cycles of size 0 mod 4 (weakest for $F_3$). As Figure 2b shows, the terminal faces $F_3$ and $F_3^{\prime }$, which support the largest ring current, have the smallest contributions to local spin density in the neutral radical from the single electron in the non-bonding Hückel molecular orbital.

6. A New Cycle Current Model

The implementation of the Aihara version of the HL model described in Section 4 was used in our investigation of the possibility of improving existing conjugated-circuit models of induced current. This stage of the work has been described in detail in [43], and here we give only a brief sketch of the strategy and and an even shorter summary of the results.

6.1. Conjugated-Circuit Models of Current

As noted earlier, the conjugated-circuit (CC) model of aromaticity is based on the idea that the important cycles in a $\pi$-system are those cycles C for which both G and $G–C$ have perfect matchings. By analogy with aromatic benzene and anti-aromatic cyclooctatetraene, cycles of size 2 mod 4 or 0 mod 4 are taken to have positive or negative effects on aromaticity. Initially, this idea was used as a way of correlating $\pi$ resonance energy with counts of cyclic substructures [36,63,64]. The idea was extended to the magnetic criterion of aromaticity by associating contributions from conjugated circuits with paratropicity (antiaromaticity)/diatropicity (aromaticity), according to the divisibility of the cycle size [37,39,40,41]. This is in effect an importation of monocycle results from HL theory, in that Kekulé structures in themselves do not have an obvious connection to the direction of current induced in ring by a perpendicular external magnetic field.

In the chemical literature, a popular approach to calculation of energetics and currents was to construct all conjugated circuits by making pairwise unions of Kekulé structures (perfect matchings) of the molecular graph G [39]. The union consists of a set of disjoint $K_2$ units where these edges are in both matchings, and a collection of even cycles where the edges in each cycle come alternately from the first and second matching. Some chemical models distinguish between cases where the union of two Kekulé structures contains only one cycle and where it contains a set of two or more disjoint cycles; with opinion divided over whether the latter case should be included in CC models of current.

A simpler way to look at the array of different CC models, which leads to a framework for classifying the various models is given in [42]. All the published CC models of current are covered by a formula for the contribution of a cycle C to bond currents. Each CC model attributes a bond current of [43]

$$w_C(a,b)=\mp 2{\left|S_C\right|}^{a}m{(G–C)}^b{F}_C$$

(26)

to each edge of the cycle, where the ∓ sign allows for the diatropic/paratropic sense of the current and the factor of 2 arises from the two perfect matchings of an even cycle. The parameter a (equal to $-1,0$, or 1) defines the weighting by cycle area $S_C$ in the particular model. The parameter b (equal to 1 or 2) defines how the weights in the model depend on $m(G–C)$. The factor $F_C$ is usually unity but allows for the extra parameters that are introduced in some CC models. This factor has also been used to normalise CC currents in various ways, though this is not relevant for our comparisons of scaled current maps. The published CC models and their parameters $(a,b)$ are [43]:

Model R, $a=0,b=2$ [39], Model CKCDA, $a=1,b=2$ [40], Model M, $a=-1,b=1$ [41], and Model GM, originally with $a=1,b=1$, later used with $b=2$ [37]. The models that include non-trivial parameters $F_C$ are M and GM, as noted in [43]:

From the viewpoint of HL theory and the possibility of partitioning currents into cycle contributions, all CC models are approximations in which the true HL cycle current $J_C$ from Equation (7) has been replaced either by $w_C(a,b)$ when C is conjugated, and zero otherwise. In fact, despite the different physical rationales for the choices of a, b and $F_C$, all CC models give qualitatively similar accounts of the current maps of most benzenoids.

However, as noted in the Introduction, the CC model maps diverge widely from HL maps for benzenoids in two cases. The first is for the benzenoids with fixed bonds. These are either perylenoids (Kekulean with fixed single bonds) or zethrenoids (Kekulean with both fixed single and fixed double bonds). All fixed bonds carry zero current in the CC current maps, as a fixed bond cannot appear in a cycle in the union of two perfect matchings. In HL maps, however, such edges may carry significant current.

In the second case, the CC maps diverge from their HL counterparts in a more dramatic way. This is for the non-Kekulean benzenoids. All CC models necessarily predict zero current on all bonds of a non-Kekulean benzenoid, whereas the HL model gives a non-zero current map for all benzenoid species with $\pi$ electron counts between 1 and $2n-1$.

6.2. A Cycle-Current Model for All Benzenoids

Figure 3 gives an example of how CC and HL current maps can show major qualitative differences. The example is zethrene, the smallest of the zethrenoid family. In this case, the scaled map is identical for all CC models. This graph has five fixed single bonds and two fixed double bonds in the bridge region, and hence vanishing CC current in all bonds outside the terminal naphthalenoid units. However, the HL map shows currents on the bridge with $40\%$ of the maximum strength, and zero current in only one bond, where it is forced by symmetry. Our goal was to find a cycle current model that would be simple, give good current approximations overall, and avoids this ‘dead-zone’ limitation of conjugated-circuit models.

We use the notation $c_i\left(G\right)$ for the coefficient of $x^i$ in the characteristic polynomial of G. Note that for a benzenoid G, $c_0\left(G\right)$ is 0 for odd n, and otherwise ${(-1)}^{n/2}m{\left(G\right)}^2$. The value of $c_0\left(G\right)$ is also equal to the determinant of the adjacency matrix of G.

The new model was inspired by the Aihara partition of current, and the best form for the first term turned out to be identical with that of the CKCDA model, which itself was designed to have the same dependence on area as the Aihara expression (7) for cycle current. The CKCDA model [40] can be thought of as assigning a current of weight $2S_C{c}_0(G–C)/c_0\left(G\right)$ to each conjugated cycle C, with clockwise/anticlockwise sense for cycles of length $4n$ or $4n+2$, respectively. We took this as the first term of the new approximation formula.

In order to ensure that some current is predicted in the non-conjugated cycles, a second term was included. Experimentation suggested adding a term proportional to $2S_C{c}_2(G–C)/c_2\left(G\right)$. This term is analogous to the first term, but involves the next ($x^2$) coefficients in the characteristic polynomials for the bipartite graphs G and $G–C$.

For a nonsingular matrix A, the classical adjoint matrix for A, denoted Adj$\left(A\right)$ is equal to $A^{-1}/$Det$\left(A\right)$. Application of Jacobi’s theorem allows the value of $c_2\left(G\right)$ to be computed in $O\left(n^2\right)$ time given the classical adjoint matrix.

The relative weight of first and second terms gives a disposable parameter in the new model, and after further experimentation it was decided to settle for a weighting of 4. This gives a formula for cycle contributions in Kekulean benzenoids that depends on the tail coefficients (i.e., those of $x^0$ and $x^2$) in $P_G\left(x\right)$ and $P_{G–C}\left(x\right)$:

$$w_C\left(W\right)=2S_C\left\{\frac{c_0(G–C)}{c_0\left(G\right)}+4\frac{c_2(G–C)}{c_2\left(G\right)}\right\}.$$

(27)

As a last step, to ensure that some current is predicted for non-Kekulean benzenoids, we rewrite this formula in a more general way, replacing $c_0$ by c${}_{\eta }$ and $c_2$ by $c_{\eta +2}$, where $\eta$ is the nullity of the benzenoid graph G, so that in this case too we are using the tail coefficients of the characteristic polynomials. Thus, in final form the new model (Model W in [43]), has cycle contributions to current given by

$$w_C\left(W\right)=2S_C\left\{\frac{c_{\eta }(G–C)}{c_{\eta }\left(G\right)}+4\frac{c_{\eta +2}(G–C)}{c_{\eta +2}\left(G\right)}\right\}.$$

(28)

Whatever the number of non-bonding orbitals in the benzenoid, the cycle contribution is specified in terms of the lowest and next-to-lowest powers of x that occur in $P_G\left(x\right)$. The result of this change is that the formula now gives currents for both Kekulean and non-Kekulean benzenoids, offering a unified solution to the two problems of fixed bonds and open shells that beset CC models.

6.3. Testing the Model

An evaluation of Model W is reported in [43], where its ability to track HL current maps was compared to that of the four published CC models and four hypothetical variants. For this comparison, the test set of benzenoids on up to 10 hexagonal rings was used: it comprises 18,360 Kekulean benzenoids (of which 2388 are perylenoids and 2184 are zethrenoids) and 20,112 non-Kekulean benzenoids.

Two types of comparison were made. Overall statistical measures of model quality were based on the bond-current error function for an edge $uv$ of G, ${\Delta }_{uv}$. This function is calculated for two sets of scaled currents, $\left\{j_{uv}^A\right\}$ from the model under test and $\left\{j_{uv}^B\right\}$ from the HL reference, using the formula ${\Delta }_{uv}=|j_{uv}^A-j_{uv}^B|$, where each current is taken in the sense of the arc from u to v. Qualititative incorrectness of some maps is detected by counting misdirected graphs. A graph G is misdirected if at least one edge of G carries currents in $\left\{j_{uv}^A\right\}$ and $\left\{j_{uv}^B\right\}$ that are both non-negligible (magnitude $>{10}^{-7}$), run in opposite directions and give rise to ${\Delta }_{uv}>0.1$. Error norms $L_1$, $L_2$ and $L_{\infty }$ are computed for the set of bond-current errors $\left\{{\Delta }_{uv}\right\}$ for each model. ($L_1$ is the mean absolute error, $L_2$ is the root mean square error, and $L_{\infty }$ is the maximum absolute error, all averaged over the molecules in the given test set). For misdirected graphs, a simple count is made.

Extensive tabulations of the relative performances of eight CC models and Model W for the test set and various subsets are given in [43]. The main conclusions are as follows.

First, Model W performs better than the best CC model for the set of Kekulean benzenoids. The errors calculated with $L_1$, $L_2$ and $L_{\infty }$ norms are all reduced by factors of two or more compared to the best CC model. Model W has $L_1=4\%$, $L_2=5\%$ and $L_{\infty }=9\%$, expressed as percentages of the maximum scaled current in each molecule. This good performance is maintained when the test set is restricted to zethrenoids. Every CC model gives at least 2247 misdirected Kekulean benzenoid graphs, including at least 952 zethrenoids, whereas the new model gives only 110 in total, all of which are zethrenoids.

Secondly, the new model performs even better for non-Kekulean benzenoids. For the non-Kekulean benzenoids, Model W gives errors of $L_1=3\%$, $L_2=4\%$ and $L_{\infty }=7\%$, and no misdirected graphs at all.

7. Conclusions

The design of this useful effective new model of current maps in benzenoids benefitted considerably from the availability of an implementation of the Aihara cycle partition of HL theory. Other applications are planned.