Ring-Current Maps for Benzenoids: Comparisons, Contradictions, and a Versatile Combinatorial Model

As a key diagnostic property of benzenoids and other polycyclic hydrocarbons, induced ring current has inspired diverse approaches for calculation, modeling, and interpretation. Grid-based methods include the ipsocentric ab initio calculation of current maps, and its surrogate, the pseudo-π model. Graph-based models include a family of conjugated-circuit (CC) models and the molecular-orbital Hückel-London (HL) model. To assess competing claims for physical relevance of derived current maps for benzenoids, a protocol for graph-reduction and comparison was devised. Graph reduction of pseudo-π grid maps highlights their overall similarity to HL maps, but also reveals systematic differences. These are ascribed to unavoidable pseudo-π proximity limitations for benzenoids with short nonbonded distances, and to poor continuity of pseudo-π current for classes of benzenoids with fixed bonds, where single-reference methods can be unreliable. Comparison between graph-based approaches shows that the published CC models all shadow HL maps reasonably well for most benzenoids (as judged by L1-, L2-, and L∞-error norms on scaled bond currents), though all exhibit physically implausible currents for systems with fixed bonds. These comparisons inspire a new combinatorial model (Model W) based on cycle decomposition of current, taking into account the two terms of lowest order that occur in the characteristic polynomial. This improves on all pure-CC models within their range of applicability, giving excellent adherence to HL maps for all Kekulean benzenoids, including those with fixed bonds (halving the rms discrepancy against scaled HL bond currents, from 11% in the best CC model, to 5% for the set of 18 360 Kekulean benzenoids on up to 10 hexagonal rings). Model W also has excellent performance for open-shell systems, where currents cannot be described at all by pure CC models (4% rms discrepancy against scaled HL bond currents for the 20112 non-Kekulean benzenoids on up to 10 hexagonal rings). Consideration of largest and next-to-largest matchings is a useful strategy for modeling and interpretation of currents in Kekulean and non-Kekulean benzenoids (nanographenes).


INTRODUCTION
This article reports an approach to comparing and improving models for mapping induced ring currents in benzenoids. The global currents induced in conjugated π systems by a perpendicular external magnetic field have long been regarded as indicators of aromaticity. 1−6 Techniques for estimating and calculating ring-current strength and direction also have a long history; 1,4,6−19 they vary from purely graph-theoretical and empirical models to ab initio quantum-mechanical perturbation theories. As usual in theoretical chemistry, there is a balance between the requirement for quantitative information on specific systems and the search for patterns in classes of molecules sharing chemically significant substructures, hence the emergence of a variety of methods.
Ab initio work on currents using the ipsocentric approach 12,17−19 gives physically realistic maps for many individual molecules 20,21 and provides a transferable conceptual analysis of ring currents in terms of symmetry, nodal character, and energy of molecular orbitals. 18,19,22 These advantages are retained in the pseudo-π model (also ipsocentric), where, in the spirit of London's pioneering calculation on benzene, 1 a cluster of H atoms acts as a surrogate for a conjugated carbon framework. The in-plane current is taken as a mimic of the out-of-plane π ring current of the real system. The pseudo-π model reduces computational cost by orders of magnitude, while reproducing the overall pattern of current flow and retaining the chemical interpretation of the ipsocentric approach. 23 Other ab initio approaches that can produce maps for visualization of magnetic response include GIMIC, 15 which concentrates on total induced current density rather than orbital contributions, ACID, 16 which plots anisotropy of the induced current density tensor, and other methods which are based on induced local magnetic fields. 24 This last quantity is related to the spatial distribution of local shieldings, and hence to the widely used NICS (nucleus-independent chemical shift) method, as developed by Schleyer and his school. 25 All these ways of producing maps are essentially grid based, and associate a current−density vector (or some proxy) with points defined on a mesh in molecular space, consistent with the through-space character of current flow in the electron cloud. Another approach is possible. In graph-theoretical methods, the current−density map is replaced by a directed graph decorated with induced currents flowing only along arcs. Of these methods, the oldest and arguably the best founded in terms of quantum mechanics is the Huckel-London (HL) theory, 1,7−11 in which magnetic response is calculated using Huckel theory with a perturbed adjacency matrix. The perturbation is carried by position-dependent gauge factors attached to notional atom-centered basis functions, hence indirectly introducing a dependence of current on the molecular geometry, although the use of idealized equilateral hexagons for all rings restores the graph-theoretical nature of the results. 10 A vast amount of literature on the HL method has shown it to be informative, often giving maps of current in striking overall agreement with those from more sophisticated calculations. Ipsocentric ab initio, pseudo-π, and Huckel-London methods are all based on molecular-orbital (MO) theory.
More recently, attention has turned to graph-theoretical models of current. A cycle C in a graph G is a conjugated circuit (CC) if both C and G−C have at least one perfect matching (Kekuléstructure). Conjugated circuits, identified by comparing pairs of Kekuléstructures, were first used in calculations of resonance energies of conjugated systems. 26 Conjugatedcircuit models of current appeared in the literature as early as 1979, 27 and then seemed to be forgotten. They were revived by Randić2 8,29 and others, and now exist in at least four related 30 variants. 27,29,31,32 In these models, bond currents arise by superposition of ring currents contributed from conjugated circuits. Signs for the ring currents are assigned ad hoc according to divisibility of the cycle length (aromatic/diatropic for (4N + 2) cycles, and anti-aromatic/paratropic for (4N) cycles), and are fixed to agree with MO results on single cycles. Claims for plausibility of the resulting maps have been made for each particular CC model, with our pseudo-π and ab initio maps often taken as the "gold standard" (e.g., in refs 29,31,32). The comparisons are at the level of visual inspection and typically made for only a few molecules. In contrast, the current models studied in this paper are evaluated analytically on a much larger set of molecules. Figure 1 shows the current− density maps obtained for the coronene molecule with a variety of MO and CC methods; all are in qualitative agreement, to the extent that they all predict concentric diatropic and paratropic rim and hub currents for this molecule, although with different estimates of relative strengths.
Hence, we arrive at the first aim of the present work: to devise a quantitative method of comparison between calculated current maps. This gives firmer ground for deciding whether one or another CC method is "best", absolutely or for some class of applications. It turns out that the automated comparison procedure detects a hitherto unnoticed problem with grid methods, leading to a reassessment of the pseudo-π method for specific classes of benzenoids. Finally, the analysis also suggests an adaptation of the CC method to give improved correlation with HL for both Kekulean and non-Kekulean benzenoids. Non-Kekulean benzenoids were previously not amenable to pure conjugated-circuit strategies, because in graphs lacking perfect matchings, such models necessarily give no current. The new model gives excellent agreement with (and a combinatorial interpretation of) HL current maps for all types of benzenoid. In turn, the new maps match well with those computed for closed-shell ions of non-Kekulean benzenoids with the pseudo-π method and rationalized with a hybrid model combining modified CC schemes and MO delocalization indices. 33,34 The structure of the paper is as follows. Section 2 (Methods) presents strategies for reduction of grid currents to currents on a graph, and discusses appropriate data sets for analysis of the pseudo-π method. Section 4 (Results) presents results on conservation of current for the graphs derived from pseudo-π maps. Benzenoids with fixed bonds are identified as the class that gives rise to poor conservation of current and systematic errors in specific pseudo-π maps. Graph-reduced pseudo-π and HL maps are compared and the essential Figure 1. Current−density maps for the π system of the coronene molecule, as obtained with different published methods. (a) Ab initio ipsocentric calculation. 20 (b) Pseudo-π model. 23 (c) Graph-based methods. When scaled so that the largest bond current (on the perimeter) is 1, all graphbased methods give the same picture, differentiated only by the value j, the size of the inner paratropic current. For the respective models discussed in this paper, j is 0.289 (Huckel-London, HL); 0.524 (RandićCC model 29 ); 0.333 (CC model due to Ciesielski et al. 31 ); 0.286 (CC model due to Mandado 32 ); 0.251 (CC model due to Gomes and Mallion 27 ). For the models to be discussed later in the text, j is 0.182 (R*), 0.122 (CKCDA*), 0.690 (M*), 0.713 (GM*), and 0.305 (Model W). similarity of maps from these two MO methods is demonstrated. An extensive comparison is then made of maps from HL theory with all published CC models and variants. Section 4 (Discussion) then sets out a new combinatorial model that gives excellent agreement with HL maps for both Kekulean and non-Kekulean benzenoids (about which pure CC methods have nothing to say). Finally, some conclusions and brief comments on future directions are presented.
All graph-based comparisons in this paper are made for normalized currents, where the maximum edge current for each molecule is set to 1. CC methods are not expected to reproduce accurate magnitudes of bond currents.

METHODS
2.1. Strategies for Conversion of Grid to Graph Currents. Ab initio calculations are usually performed at the equilibrium molecular geometry, which typically has a distribution of bond lengths and angles. However, pseudo-π calculations for benzenoids normally assume a molecular graph drawn on a regular hexagonal lattice (with fixed bond length 1.4 Å). We use these pseudo-π maps as proxies for full ab initio computation in what follows, since we need to survey hundreds of systems. Huckel-London calculations are typically carried out with this idealized geometry, and CC methods that use weightings based on area 10,11,27,31,32 also normally assume regular hexagonal rings. All maps in the rest of this paper are calculated under these constraints.
The first hurdle is to devise tactics for comparing maps from grid-and graph-based methods. In the grid methods, the primary output is a vector for each grid point (this is j, the first derivative of induced current density with respect to the applied field, here called "the current" for short). Our grid maps 18,22 show arrows representing the projection of j in the plotting plane, distributed on a sparse grid for clarity, and superimposed on a denser contour plot of the magnitude |j|, (see, for example, Figure 1a). The task is to condense this information onto the molecular graph, for which the nongrid (HL and CC) methods deliver a single arc per edge representing magnitude and sense of current along the bond.
Some basic preliminary considerations are connected with the physical nature of current. It is clearly desirable that a map should show continuity of current. The HL and CC maps obey conservation of current automatically, as current flow at each node of the graph follows Kirchhoff's First Law in that the sums of scalar in-and out-currents are equal. An MO method uses explicit basis functions, however, and maintenance of current conservation is dependent on operator sum rules that are fulfilled only approximately with finite basis sets. 35 In the pseudo-π method, this can become a problem. The basis is minimal, consisting of a single STO-3G H 1s function per site, and there are departures from continuity, although they are not always evident in maps with a sparse grid of arrows. When unusually large, these deviations may indicate problems with the representation of the ground-state wave function. However, the errors for this simplified version of the ipsocentric method are still typically much smaller than the glaring discrepancies that occur in conventional single-origin calculations with medium-sized basis sets. 21 A second requirement is that the map should show the full symmetry of the molecule in the magnetic field. For a planar molecule, this implies that the vector field of induced current should transform as an in-plane rotation in the point group of the molecular framework. This is automatic for graph-based The Journal of Physical Chemistry A pubs.acs.org/JPCA Article methods. However, symmetry is easily lost in grid-based methods, simply because of the mismatch between the predefined square grid and the hexagonal array of atoms and bonds. Figure 2 shows magnitudes of pseudo-π currents (as integers 1 to 99 that express the current relative to the maximum for the molecule) for benzene plotted on a square grid with origin at the center of the hexagon. Currents are seen to be only roughly symmetric. Clearly, for proper comparison with the graph-theoretical picture, smoothing or a denser grid is needed. To achieve a clean comparison, we use here a symmetry-adapted triangulated grid with an odd number of points per graph edge, which treats every hexagon equally and maintains symmetry for all benzenoids. Figure 3 confirms that this benzenoid-adapted grid for benzene (with nine grid points along each carbon−carbon bond) now displays the full natural symmetry. Next, to make a grid-to-graph conversion of the pseudo-π map on the hexagonal grid, we investigated three options for calculating bond currents along graph edges: (1) Midpoint: To choose the sense, from vertex u to vertex v, or vice versa, and calculate the magnitude of the equivalent bond current, take the j vector for the grid point nearest the midpoint of the line between the nuclei. (2) Catchment: Take all grid points and make a crude sum of |j| by assigning every grid-point value to the nearest edge center. (3) Cross-stream: Take a perpendicular bisector of each edge and sum |j| over 2k + 1 points symmetrically distributed with respect to the bond midpoint, attributing the total to the edge midpoint. This gives a local sampling of current in the region of the bond. The value of 2k + 1 = 7 was chosen here.
Versions of each option using vector summation of j with resolution onto the bond direction were also investigated, giving marginal improvements in the statistics. All options were first evaluated according to how well the currents that they produced obeyed conservation of flow. The next step was to compare the grid currents to HL currents. As comparisons are made for scaled maps, currents were normalized in the converted pseudo-π map so that the largest bond current in each graph was a strength one.
2.2. Note on the Computational Approach. In the course of this research, many different current models were considered and, except for those described below, rejected. The initial phase started with Kekulean benzenoids having up to seven hexagons. 36,37 To facilitate examination of subsets, the comparison program used an auxiliary file that listed for each molecule the graph number and a code 0 (exclude) or 1 (include). A program was created to make auxiliary files according to various rules. As discussed below, the first restriction was exclusion of cases with proximity problems. Later, classes of benzenoids with fixed bonds arose as problematic cases. These are detected by the program that makes the auxiliary files by computing the number of matchings per edge. Use of an auxiliary file was much easier than attempting to generate each class of interest (especially as in the initial stages it was not known which molecules would turn out to be problematic). For the second phase (consideration of non-Kekulean benzenoids), the f usgen module of the CaGe program 37 was used to construct all the fusenes, and a second program was used to screen out the cases with at least one perfect matching. We could alternatively have continued using a larger set (all fusenes) with an auxiliary file, but this way saved space.
The first step in all computing the currents in the various graph-based models is to generate all the cycles of the benzenoid. We did this by starting with a list of the internal faces. It is well-known that the set of faces is a cycle basis for the set of all cycles of a planar graph (meaning that every other cycle is a linear combination of these: when two cycles are added, if an edge appears an even number of times it does not appear in the sum). Then for h going from one to the number of faces, each of the faces is added to each of the cycles containing h hexagons. Cases having one new cycle that has not already been generated are added to the list of cycles.

Generation and Choice of Data Sets.
The initial test set of benzenoids (simply connected subgraphs of the hexagonal lattice composed of edge-fused hexagons) consisted of the 266 Kekulean benzenoids composed of seven or fewer fused hexagons (Test Set 1). CC models make no prediction at all for non-Kekulean systems, and all our pseudo-π calculations are performed for closed shells, so only Kekulean benzenoids were included in the initial comparison. However, even this natural set of benzenoids presents a subtlety. The pseudo-π method has a particular limitation. 38 In this method, there is no concept of a bond. There is just a set of centers bearing basis functions and giving a σ simulation of the carbon π system; 23 bonds are present only implicitly through intervertex distances. Thus, the method can never distinguish between bonded and nonbonded pairs of carbon atoms with the same separation, as may occur in the regular hexagonal geometry. (In real molecular systems, nonbonded distances are larger, and interactions are damped by steric and electronic effects of intervening hydrogen atoms, which may also lead to out-ofplane distortion.) In fjord regions on the perimeter of a benzenoid inscribed on the regular hexagonal lattice, extra "fake edges" will lead to a pseudo-π map appropriate to a different benzenoid (see Figure 4). Benzenoids that give rise to this problem were removed from the set, to leave 232 that are more fairly representative of the correct use of the pseudo-π approach (Test Set 2).

RESULTS
3.1. Conservation of Flow in Pseudo-π Maps. Table 1 summarizes the results on conservation of flow in pseudo-π maps for Test Set 2. In these comparisons, the continuity error at a vertex v, Δ v , was defined as the absolute value of the difference between total current in and total current out. Three measures of error were then applied. The L 1 norm (mean absolute error) takes the total error as the average sum of Δ v taken over all vertices. The more usual L 2 norm (root-meansquare error) takes the total error as the square root of the sum of squares of Δ v averaged over all vertices. The L ∞ norm (maximum absolute error) takes the error as the maximum value of Δ v over all vertices.
The main conclusion from the error measures is that the resolved cross-stream method gives the best continuity overall. All further results will be reported for graph-reduced pseudo-π maps obtained by this method. The comparisons also give more information about the source of continuity errors. A useful parameter for highlighting cases of poor continuity is D v , the percentage current drop across vertex v, defined by expressing Δ v as a percentage of the larger of incoming and outgoing currents. The value D for the graph is the largest value of D v . In the tables, D̅ is the average of D over all graphs in the data set and D max is the largest value for any graph in the set. The large differences between average and maximum values of D (final columns of Table 1) indicate that the errors are not evenly distributed across the whole set.
Much of the remaining continuity error in maps obtained with the cross-stream method is attributable to violations of current conservation in the pseudo-π model for specific types of benzenoid. This is made clear by the partitioning of the data set shown in Table 2. To anticipate the discussion below, we note that the very high values of current drop occur for the two classes that we call zethrenoids and perylenoids. Zethrenoids are those benzenoids that have a fixed double bond (an edge present in all perfect matchings); perylenoids are those benzenoids that have no fixed double bond but at least one fixed single bond (an edge present in no perfect matchings).
The smallest examples of each of these types are shown in Figure 5. All the examples with D max ≥ 15% in Test Set 2 belong to one of these types. If we exclude all cases with D max ≥ 30%, we remove exactly the 7 zethrenoids; if we exclude all cases with D max ≥ 15%, we also remove all but 2 of the 17 perylenoids in the Test Set (and no other graphs). . Two non-isomorphic benzenoid frameworks that give necessarily identical pseudo-π maps for idealized geometries because they have the same set of 1.4 Å pairwise distances on the ideal hexagonal-grid.  a The results are calculated using graph reductions from the three sampling methods described in the text, each in a vector (res) and scalar version. The dataset is Test Set 2, consisting of the 232 Kekulean benzenoids in the range that do not have the short nonbonded contacts that would give rise to proximity errors in pseudo-π maps. Δ v is the difference between incoming and outgoing scaled currents at vertex v. The error norms are described in the text.
D is the percentage current drop at the least well conserved vertex of the graph. D̅ is the average of D over the dataset and D max is the largest value for any graph in the set.  a The continuity errors noted in Table 1 are recomputed for a partition of the Kekulean benzenoids of Test Set 2 into zethrenoids (z), perylenoids (p), and remaining (r) benzenoids (with no fixed bonds). Notation as in Table 1.
The Journal of Physical Chemistry A pubs.acs.org/JPCA Article 3.2. Comparison of Normalized Pseudo-π and HL Maps. Table 3 reports a full comparison using Test Set 2 of graph-reduced pseudo-π maps from the cross-stream sampling method with the HL maps. Note that when comparing maps that consist of two sets of scaled bond currents {j a } and {j b }, the error on edge uv is calculated as Δ uv = |j uv a − j uv b |, where each current is taken in the sense of the arc from u to v. Error norms L 1 , L 2 , and L ∞ are then applied to this edge function. Similar comparisons were made for all the sampling methods, and it turns out that the resolved cross-stream sampling method gives the closest match of pseudo-π to Huckel-London maps, possibly because it best reflects the conservation property implicit in the HL currents.
The main conclusion is that there is evidence of a reasonable degree of correlation between maps from pseudo-π and HL models (rms error <7%), as might have been hoped for methods both derived from minimal-basis molecular-orbital theory. However, one striking observation from the initial survey reported in both Tables 2 and 3 is that the errors are not distributed uniformly. Difficulties with nonconservation of current in the reduced pseudo-π maps and poor fit with HL maps are both concentrated in certain subsets of the tested benzenoids: the zethrenoids and perylenoids defined above. This observation provoked further inspection of the original unreduced maps, where an anomalously severe lack of conservation was immediately apparent in these "difficult" cases, as discussed below.
3.3. Benzenoids with Fixed Bonds. The common factor for the aberrant benzenoids is that all contain some fixed bond(s). In chemical terms, some bonds (edges of the graph) have either unit or zero Pauling bond order (the fraction of all perfect matchings of the graph that contain that edge), and so cannot be part of any conjugated circuit. A fixed double bond implies fixed single bonds, but not vice versa. Our definition of the zethrenoids and perylenoids captures all benzenoids that have fixed bonds: those benzenoids that have one or more edges with Pauling bond order 1 are the zethrenoids, and those that have one or more edges with Pauling bond order 0 that are not zethrenoids are the perylenoids.
As an example of the nature of the problem with zethrenoids, Figure 6 shows pseudo-π maps for zethrenoids in Test Set 2. These maps all show what we call "eyebrow" features, which are plainly unphysical. Current is launched onto the central bridge from each naphthalene moiety, but peters out almost immediately. The zethrene molecule belongs to a class of alternant hydrocarbons in which Kekulean fragments (here naphthalenoid) are connected to the rest of the graph via two vertices at even distance in the isolated fragment (and hence in the same partite set) through formal single bonds to the rest of the graph. The zethrene molecule and its analogues have inspired a rich chemistry. 39 The bridging region of zethrene necessarily contains only fixed bonds. It is clear that the pseudo-π method does not cope with the description of currents in this intermediate region.
The pseudo-π maps for zethrene were not noticeably improved by repeating the calculation in tests at the ab initio CHF or DFT levels. Problems with zethrene and similar molecules have been noted before, and are often attributed to incipient radical character, 40,41 implying the need for a multireference treatment. Similar difficulties have been noted in treatments of ballistic conduction, where currents in some conjugated molecules can be badly misrepresented in singlereference treatments. 42 Whatever the source of the problem, there is a clear subset of bad benzenoid graphs that give major continuity discrepancies for the pseudo-π method in bridge regions of fixed bonds.
Three further data sets were constructed to explore this problem. Test Set 2z consists of the seven zethrenoids in Test Set 2 (one further zethrenoid from Test Set 1 is excluded for proximity reasons); Test Set 2p consists of the 17 perylenoids in Test Set 2 (two further perylenoids from Test Set 1 are excluded for proximity reasons); Test Set 2r, with 208 members, consists of the remaining benzenoids in Test Set 2, in other words, the set of Kekulean benzenoids on up to seven hexagons that suffer neither from the proximity limitation nor from fixed bonds. Table 3 shows the effect on the quality of the match between pseudo-π and HL maps. The overall conclusion is that HL currents and pseudo-π currents show a clearly improved match when zethrenoids in particular, and less significantly, perylenoids are excluded from the comparison.
Looking ahead, we note that the different methods in our comparative study treat the zethrene case in different ways. In the Huckel-London (HL) method, bonds have intermediate character and carry current except where precluded by symmetry. For zethrenoids, HL gives a map with flow around the perimeter of the bridging region (Figure 6b). This is at least physically plausible, though it is not obvious that the computed currents are "correct", given the absence of a clean ab initio reference. In contrast, any pure CC method must allocate zero flow through all fixed bonds, and hence all CC methods give maps for zethrene itself that consist of disconnected naphthalenoid islands, as noted in the earliest paper in this area. 27 In CC methods, perylenoids also have islands of current separated by one or more empty bridge regions. HL and pseudo-π maps for small perylenoids start off with this qualitative general pattern, but with a minor component of current allowed on the bridge. In the larger systems, the bridge current increases and eventually becomes dominant. For   Datasets and notation as in Table 2.
The Journal of Physical Chemistry A pubs.acs.org/JPCA Article example, pseudo-π and HL maps for the symmetrically extended perylenes with up to four hexagons in the central belt are illustrated in Figure 7.The two methods agree, at least qualitatitively, in their progressive departure from the very simple CC island pattern: in both sets of maps, as the systems become larger, significant current begins to circulate in the central region. The increasing complexity of ipsocentric current maps for perylenes has been noted before, and led the authors of one paper 41 to make the sweeping recommendation that CC models should be abandoned altogether. We would argue that such models are often useful zero-order pictures, and that efforts should instead concentrate on improving them and extending their range of applicability (see Discussion of an improved model below). For zethrenes and perylenes, we should, however, be aware of the significant shortcomings of single-reference treatments, which can have drastic effects on predictions of molecular conduction, for example. 42 Other authors have suggested that some wholesale "renormalization" of the MO results may be needed. 43 Here, we can simply note that there is a well-defined class of benzenoids (those with fixed bonds) for which the results of pseudo-π and HL methods are likely to deviate from those of pure CC models.  . This set is partitioned as before into three subsets: 3z (zethrenoids), 3p (perylenoids), and 3r (the rest). The choice of HL for comparison with CC methods is motivated by several considerations. These include our acceptance of HL as a reliable proxy for all the MO-based methods, the shared graph-theoretical nature of HL and CC models, the implicit neglect of full electron interaction common to HL and CC, and the fact that HL currents themselves can be written exactly as sums of cycle contributions. On this last point, Aihara 44−50 has championed the cycle partition of Huckel energies and magnetic properties. He developed a formulation of Huckel theory in which current is calculated as a sum over cycles (not just the conjugated circuits). This description of HL theory informed the construction of a CC model by Ciesielski, Krygowski, Cyranśki, Dobrowolski, and Aihara (CKCDA). 31 In some sense, we can regard all the CC models as approximations to the Aihara formulation of HL theory in which model weights (including zero for the cycles that are not conjugated circuits) have replaced the exact HL cycle weights. A paper in preparation will give a detailed description of a new implementation of Aihara's method that clarifies this comparison with CC models and uses it for further analysis. Meanwhile, Aihara's approach 44−50 is used here as motivation for construction of the improved current model that we propose in the following section.
As mentioned earlier, there are at least four CC models in the chemical literature. In each of them, currents are assigned to the cycles that are conjugated circuits, with a direction that depends on divisibility of the cycle size (diatropic for (4N + 2)-CC, paratropic for (4N)-CC). The strength of the CC current is computed with a different weighting scheme in each model. Bond currents are calculated by superposition of CC contributions. As shown elsewhere, 30 the different weighting schemes can be represented by a unified formula where the ± sign allows for the diatropic/paratropic orientation, and the factor of 2 comes from the two perfect matchings of the cycle. The parameter a (equal to −1, 0, 1) defines the weighting by cycle area (defined with reference to a standard hexagon) in the particular model. The parameter b (equal to 1 or 2) shows how the weights in the model depend on m(G − C), the number of perfect matchings of the subgraph induced by deletion of the cycle from the molecular graph G. The factor F C is usually just unity, but allows for the extra empirical multiplicative factors introduced in two of the models (see below). The weights may also be normalized in different ways, but as we use scaled maps, this is not relevant. The literature models are due to Randić2 9 (Model R, a = 0, b = 2), the group of Cieselski, Krygowski, Cyranśki, Dobrowolski, and Aihara 31 (Model CKCDA, a = 1, b = 2), Mandado 32 (Model M, a = −1, b = 1), and Gomes and Mallion 27 (Model GM, originally 30 used with a = 1, b = 1, though more recently 11 with a = 1, b = 2). Of the two models that include an empirical weighting, Mandado's model M includes a factor of (0.78) −2 ≈ 1.644 that is applied to boost the cycle current for those CC formed of linear chains of 3 or more hexagons; for simplicity, we neglect this in our implementation. The other is the GM model, which has sizedependent multipliers derived from earlier semiempirical calculations (see Table 1  The different values of a reflect different physical intuitions of the modelers. The b parameter arises because, in the literature, conjugated circuits were originally constructed and counted by considering the superpositions of ordered pairs of Kekuléstructures of the molecular graph. The union of two perfect matchings is a set of disjoint cycles of size 4 or more (each of which is a CC) and a set of matching edges (which can be regarded formally as trivial 2-cycles). Only the CC contribute to ring current in conjugated-circuit models. In some models, 29,31 contributions to current of ∓2|A| a are included for every nontrivial cycle for each ordered pair of perfect matchings; this corresponds 30 to weighting each CC by m(G − C) 2 . In other models, 27,32 a contribution to the current arises only for ordered pairs of perfect matchings with exactly one CC; this corresponds 30 to weighting each CC by m(G − C). Calculation by eq 1 is more efficient than considering each ordered pair of perfect matchings, especially for highly branched benzenoids with very large numbers of Kekuleś tructures. 51,52 For Kekulean benzenoids, only diatropic ring currents appear in CC models, as all conjugated circuits are of type 4N+2. Of course, for non-Kekulean benzenoids all currents are zero for all CC models, because non-Kekulean benzenoids do not have conjugated circuits in the strict sense.
We note that for every published CC current model there is a hypothetical twin with the same value of a but the opposite value of b (recall that b = 1, or 2). We will call these twin models R*, CKCDA*, M*, and GM*, respectively.
A word about electron counting in HL and CC calculations is needed here. HL theory includes electrons explicitly and the predicted currents can be written as sums of contributions from occupied orbitals. As nonsingular bipartite graphs, Kekulean benzenoids in Huckel theory have no nonbonding orbitals, and the neutral molecules have closed-shell π electronic configurations in which all n 2 bonding π orbitals are filled and all n 2 antibonding π orbitals are empty, as implied by the Pairing Theorem. 53 Each Kekuléstructure implies the same π electron count, based on localized two-electron π bonds on edges in the perfect matching. Therefore, in comparing the maps from nπ-electron HL calculations on benzenoids with those from CC models we are indeed comparing like with like. Table 4 reports a statistical comparison of HL and CC maps for Kekulean benzenoids using the four published CC models and the variants derived by filling out the array of possible a and b values. Since the HL and CC methods do not have proximity limitations and do not suffer from discontinuity errors, we take the full set of Kekulean benzenoids (Test Set 3) and its partition into zethrenoids, perylenoids, and the rest. The various error norms defined earlier are used to compare models according to their edge error functions, Δ uv . We also introduce a simple criterion for how far the description of a single molecule differs qualitatively between models: a misdirected graph for a given model is one where, on at least one edge, after scaling so that the maximum current in both maps is 1, the model and HL currents are each non-negligible (exceed a nominal threshold of 10 −7 ), run in opposite directions, and give rise to a value of Δ uv ≥ 0.1.
The table shows that all the published CC models have some claim as rough approximations to the HL maps for  43 Although extant CC models all achieve some success in shadowing HL theory, it is natural to wonder if their accuracy could be improved. The range of systems to which they apply is also limited. By definition, CC models make dramatically oversimplified predictions about currents in molecular regions where some bonds are fixed, and no prediction at all about the currents in non-Kekulean systems. As we will show in the next section, both problems can be addressed with a relatively simple modification of the CC approach. The aim is to augment the CC model in a way that preserves its combinatorial spirit but improves the fit to HL currents for Kekulean benzenoids, and removes the blanket prohibition on current in fixed bonds and non-Kekulean structures. A hint comes from Aihara's partition of HL current into cycle contributions. The Aihara formulation replaces the standard HL sum over orbital contributions with a sum over orbitals and all cycles present in the graph. The orbital contribution to the current flowing in cycle C arising from a shell of m k degenerate orbitals with common eigenvalue λ k and n k electrons per orbital is 48 for m k > 1. Here, P G-C (x) is the characteristic polynomial of G − C and U k (x) is the reduced polynomial P x x ( )/( ) , with P G (x) the characteristic polynomial of G, and A(C) is the area of the cycle (i.e., in a benzenoid, the number of hexagons that it encloses). For a single nondegenerate orbital, the orbital contribution is The characteristic polynomial for a Kekulean benzenoid G is A link can be made between characteristic polynomials of benzenoids and CC models via the Leibniz formula for the determinant. 54 The characteristic polynomial, P G(x) , of a bipartite graph G (the molecular graph of an alternant conjugated π system) is Hence, g n = +1 and g 0 = (−1) n det A. The Leibniz formula for a determinant implies that the coefficient g i in the characteristic polynomial of a bipartite graph (a graph that has no odd cycles) is where the first summation runs over all subsets of size i of the vertices of G, and for each such subset the second summation runs over ordered pairs of perfect matchings of the graph G−S (including the s of each perfect matching with itself). The sign to be attributed to P is where the product runs over all cycles of the graph induced by M 1 ∪ M 2 . For this computation, all even cycles contribute a factor −1 to the product. A similar formula can be derived for graphs that are not bipartite, but then the sum no longer ranges over ordered pairs of perfect matchings, and there can also be odd cycles that contribute +1 to the sign. If an edge is in both perfect matchings then it is considered to be a 2-cycle, and hence makes a contribution of −1 to the sign. Equation 7 gives a straightforward derivation of a useful standard result for the coefficient g 0 in the characteristic polynomial of a benzenoid. The result is that g 0 = (−1) n/2 m(G) 2 , where m(G) is the number of Kekuléstructures of benzenoid G. This follows from the special case of eq 7 with i = 0, which gives The value of sign(P) is (−1) c , where c is the number of cycles (including 2-cycles) in the graph induced by the pair of perfect matchings. Equivalently, sign(P) = (−1) s , where s is the number of connected components of the graph induced by M 1 ∪M 2 .The graph induced by M 1 ∪M 2 consists of 2-cycles from common edges and larger even cycles from CC. As all CC of a benzenoid are of size 2 mod 4, the number of even cycles (including 2-cycles) in G is odd for n ≡ 2 mod 4 and even for n ≡ 0 mod 4, yielding the desired sign(P) = (−1) n/2 for all P.
We can also find c 0 , the tail coefficient of the characteristic polynomial of G − C where G is a Kekulean benzenoid and C is a CC of G. The derivation in this case is as follows. Let H = G − C, where C is a CC of Kekulean benzenoid G. If the number of vertices of G is n ≡ r mod 4, then the number of vertices of H is (r + 2) mod 4. The constant term in the characteristic polynomial P H (x) is then Note that if C is a cycle in G that is not a conjugated circuit, c 0 in P H (x) is 0. This follows from the fact that a cycle C is a CC if and only if G − C has at least one perfect matching.
As an aside, a straightforward general observation can be made about the relationship between perfect matchings of G − C and G for any bipartite graph G and any even cycle C. Let C be an (even) cycle of G. Any perfect matching of H = G − C can be extended to a perfect matching of G by adding in alternating edges from C.
The relevance of the above expressions for g 0 and c 0 to CC current models is that for a Kekulean benzenoid they give The simplest parameter-free CC model with the same explicit dependence on cycle area as the Aihara HL formulation is Model CKCDA. Equation 11 allows the weights in this model to be expressed as products of cycle area and the tail coefficient of P G−C (x). The proposal for an improved current model is to use this as a starting point and to take into account further coefficients from the characteristic polynomials. This can be done in a general way that applies to both Kekulean and non-Kekulean benzenoids.
The adjacency matrix of a general benzenoid has η zero eigenvalues, where η is the nullity of the graph. The general form of the characteristic polynomial for benzenoid G is then As benzenoids are bipartite, n and η have the same parity (by the pairing theorem 53 ), so the polynomial in eq 12 contains only odd powers of x when n and η are both odd, and only even powers of x when n and η are both even. Likewise, deletion of an even cycle C of size r from a general benzenoid leads to a polynomial with coefficients c i = 0 either for all i ≡ 0 mod 2 or all i ≡ 1 mod 2. With all this taken into account, the new current model (Model W) has cycle contributions based on the tail coefficients of eq 12 and the corresponding coefficients of each P G−C (x): where α is an adjustable parameter. Rough optimization suggests α = 4 as a good value for benzenoids. The versions tested had α = 0, 1, 4. When the weight of a cycle is negative, a current of magnitude |w(α, C)| is assigned to flow anticlockwise on the cycle; when the weight is positive, the current is routed in the clockwise (paratropic) direction. For benzene, formula 13 evaluates to −1/2. From the deduction in eq 7 from the Leibniz formula, it follows that the coefficients appearing in the model all have combinatorial interpretations in terms of pairs of matchings. For Kekulean benzenoids, g 0 and c 0 are related to pairs of perfect matchings in G and G − C, and the second term of the model involves g 2 and c 2 . The coefficient g 2 is where the first summation runs over pairs of vertices in G and the second over ordered pairs of perfect matchings of the graph induced by deletion of these vertices. The graph G − v 1 − v 2 is in general no longer a benzenoid, and pairs of matchings may have even cycles of types 4N and 4N + 2, leading to differences The Journal of Physical Chemistry A pubs.acs.org/JPCA Article in sign(P). Figure 8 shows an example of a Kekulean benzenoid that has some terms contributing +1 to the second sum and others that contribute −1. Similar considerations apply to the graphs involved in the expression for the coefficient c 2 . 4.1.2. Model W for Kekulean Benzenoids. Table 5 shows the results for model W (eq 13) with parameters α = 0, 1, 4 compared to the best results from pure CC models. Model W with α = 0 is equivalent to model CKCDA. The Test Set 3 and the subsets 3z, 3p, and 3r are as used in Table 4. For Kekulean benzenoids, the addition of the extra term in Model W improves all error measures. With α = 4, the computed maps of scaled currents on all 18 360 Kekulean benzenoids on up to 10 hexagons have very few "misdirected" cases, and these disappear altogether when zethrenoids and perylenoids are removed from the test set. All published CC models have much larger percentages of misdirected graphs, and also higher error norms. Model W reproduces HL currents to within 5% rms error, showing that inclusion of next-to-largest matchings is a key element in building a realistic current map for these benzenoids.

Model W for Non-Kekulean Benzenoids.
Here the test set is Test Set 4, consisting of the 20 112 non-Kekulean benzenoids on 10 hexagons or fewer. The comparison for non-Kekulean benzenoids is extremely satisfactory. Taking into account tail coefficients based on nullity gives current maps that are superior in two ways: they not only have nonvanishing current, but they are in good agreement with HL predictions. In fact, the error measures in the last section of Table 5 show that the model is performing even better for non-Kekulean cases.
Again, it is important to be sure that we are comparing like with like in terms of electron count. The Hund's Rule πelectron configuration for a neutral molecule with a bipartite molecular graph has double occupation of all (n − η)/2 bonding orbitals. The remaining η electrons are distributed over the η orbitals of the nonbonding shell. The (n − η)/2 antibonding orbitals are empty. The HL calculation is carried out in the fractional-occupation approximation, where each orbital in a shell of p orbitals with a total occupation of q electrons is assigned occupation number q/p. Hence for neutral benzenoids, the nonbonding shell (if η ≠ 0) has occupancy of one electron per orbital. However, it is easy to show that the total contribution to the HL current from a nonbonding shell with uniform occupation is exactly 0. Hence, the HL maps that we compute for benzenoid B would apply to all species from B η+ to B η− . A similar logic dictates the use of nearest cationic species in calculations of pseudo-π maps for non-Kekulean benzenoids. 33 This would be significant for experimental investigation of ring currents, as conventional NMR measurements require closed-shell, diamagnetic species. The ability to treat larger high-spin species such as the triangulenes, 55−57 which have been proposed as prototypes of spintronic systems, is an important advantage of a post-CC model.
As a concrete example, Figure 9 sets out the steps in the analytical calculation of the current according to Model W for the smallest non-Kekulean benzenoid, the phenalenyl radical. Here η = 1, and the three cycle types are C 6 , C 10 , and C 12 , enclosing 1, 2, and 3 hexagons. From eq 13, using the characteristic polynomials listed in the figure caption, the cycle weights are 10 207 w( , C ) 6 ( 1) 12 The net current in a perimeter bond has unscaled value w w w ( , C ) 2 ( , C ) ( , C ) 1 3 28 207 leading to a scaled map with a diatropic perimeter current of 1, and no current in the interior region. This is in agreement with HL maps, and fully consistent with the results of ab initio calculations for the 12π and 14π ions of this molecule. 58 Contributions to the currents in eqs 15 to 17 are also readily interpreted in terms of perfect matchings of vertex-deleted subgraphs of G and G − C. The graphs G − C are paths (Figure 9b), and their analysis is particularly straightforward: all give negative contributions to c 1 (G − C) and the first two give positive contributions to c 3 (G − C)). Figure 9c−h gives the slightly more involved analysis needed for the graphs G − v and G − u − v − w, as simplified by the use of symmetry. One interesting case is the graph G − v when v is the central black vertex (left-hand graph in Figure 9c). This graph is a 12-cycle with two matchings, which give complete cancellation of the contributions to g 1 from their four pairings. In the other graph of type G − v (the right-hand graph in Figure 9c), the contributions to g 1 from pairs of matchings are all positive. Likewise, all contributions to g 3 arising from pairs of matchings of graphs G − u − v − w are negative. The two terms of the model give reinforcing diatropic currents. The first is that, as with all cyclebased models, Model W automatically obeys conservation of current, and as with all graph-based models, it does not suffer from proximity limitations caused by indirect definitions of bonded contacts. As presently formulated, the model is intended to apply to bipartite graphs/alternant molecules: for non-bipartite graphs, odd cycles may contribute current, and although the coefficients c 0 , c 2 , ..., g 0 , g 2 , ... can still be calculated, their interpretation purely in terms of perfect matchings no longer applies.
The second observation is that the key to the success of Model W in dealing with current through regions of the molecule that have fixed bonds is the second term in eq 13. To get flow of the required type, it is necessary to assign current to some cycles that are not CC. As Figure 6c shows, for example, this allows the model to mimic the HL flow of current in the problematic fixed-bond region of a zethrenoid. For non-Kekulean systems, the success of the model stems from the change from considering perfect matchings, as in Kekulean systems, to maximum matchings. For these systems, the second term in eq 13 typically acts to improve agreement with the HL model.
Third, our comparisons between MO and combinatorial models (including Model W) are at the level of the scaled maps, and this is done for two reasons: the models often have no settled normalization scheme, 43 and they have no secure claim to reproduce absolute currents.
A fourth point relates to the mathematical nature of the model. The characteristic polynomial coefficients needed for the calculation of current for benzenoids in this model are all familiar objects, which can be determined by standard methods in a number of ways. They are related to pairs of matchings in a way that is expressed precisely by the Leibniz-related formula 7. They can also be calculated directly from eigenvalues of the adjacency matrix A(G). Given this set of values, the tail coefficient g 0 is found from their product, g 2 is a sum of (n − 2) subsets, and so on. As a third possibility, Jacobi's equality for complementary minors 59 can be used to find the (n − i) × (n − i) subdeterminants of A needed for a Leibniz-type calculation, by expressing them in terms of smaller i × i determinants. This simplification avoids the need to calculate eigenvalues, but requires matrix inversion, as the formula requires the inverse of A for a Kekulean benzenoid, or inversion of adjacency matrices for vertex-deleted subgraphs for non-Kekulean cases.
Finally, although the model is a significant extension of the CC approach in terms of applicability and accuracy, there is also continuity with the CC picture and the underlying 4N/4N + 2 classic rules for magnetic response of individual cycles. This is illustrated by a simple inspection of the signs of η and η + 2 terms of eq 13. As particular coefficients of the Calculations are for benzenoids with at most 10 hexagonal rings. Test Set 3 is the set of all Kekulean benzenoids in this range. Test Set 4 is the set of all non-Kekulean benzenoids in this range. The row labeled "Improvement" shows the equivalent result from the best-performing CC model for the given test set, and (in parentheses) the factor by which model W with α = 4 improves on that result. Entries (−) indicate that model W has eliminated a particular error. Other definitions as in Table 4. In other words, if the two-term model predicts current in a given cycle, it will be in the standard diatropic/paratropic sense predicted by the MO method.
Consistently with the conjectured rule, the combinations (+,−) and (−,+) are not seen at all in the survey, though all others are. The combination (0, 0) is rare, seen only in one Kekulean case of a cycle with size 2 mod 4 ( Figure 11), but of course this may be an effect of small sample size. Two combinations that are impossible for Kekulean benzenoids are (+,+) and (+,0) for cycles of size |C| ≡ 0 mod 4 (because for a 4N-cycle C in a Kekulean benzenoids, G − C has no perfect matchings); (+,+) is possible for cycles with |C| ≡ 0 mod 4 in non-Kekulean benzenoids, and indeed appears to be the most common case. For cycle size |C| ≡ 2 mod 4, the case where the second term reinforces the first, (−,−) predominates for both Kekulean and non-Kekulean benzenoids.

CONCLUSIONS
5.1. Beyond Benzenoids? This paper has concentrated on models for mapping induced currents in benzenoids, where CC, HL and the new combinatorial Model W give qualitatively similar maps for Kekulean systems. Unanimity is not to be expected to persist in general, even for Kekulean systems containing nonhexagonal rings. The example of a molecule as simple as butalene shows how disparate predictions can arise ( Figure 10). Butalene has a bipartite molecular graph G consisting of two edge-fused 4-cycles. This graph has three perfect matchings, and three conjugated circuits consisting of the two 4-cycles and the perimeter 6-cycle, having cycle currents running in paratropic, paratropic, and diatropic senses, and relative areas 1/2, 1/2, and 1, respectively. By symmetry, all perimeter currents are equal and the central edge carries no current. As the matching counts m(G − C) are equal to 1 for all three conjugated circuits, the parameter b is not relevant. However, the weighting of cycle areas is critical to the qualitative prediction for the global current from each CC model: models with a = −1 (M, M*) predict a paratropic Figure 9. Example of the use of Model W to predict the current map of a non-Kekulenoid benzenoids (in this case, the phenalenyl radical). (a) Numbering scheme for the molecular graph G and the list of orbits of equivalent vertices. The characteristic polynomial of G is P G (x) = x(x 12 − 15x 11 + 84x 9 − 226x 7 + 309x 5 − 207x 3 + 54). (b) The three non-isomorphic graphs G − C obtained by deletion of an even cycle C from G. They are the paths on 7, 3, and 1 vertices, with respective characteristic polynomials x(x 6 − 6x 4 + 10x 2 − 4), x(x 2 − 2), and x. (c) Three non-isomorphic graphs G − v obtained by deletion of a vertex belonging to the larger partite set in G, showing total contributions to the coefficient of x in P G (x). (d−h) Graphs G − u − v − w obtained from the five non-isomorphic choices of two black vertices (one per row) and any compatible white vertex, showing total contributions to the coefficient of x 3 in P G (x). Conventions for graph sets (c) to (h): contributions to coefficients of the characteristic polynomial P G (x) are given in the form s × t where s is the number of symmetry copies of the graph, and t is the summed contribution of all pairs of perfect matchings of a single copy. One perfect matching is illustrated for each graph. Deleted portions of graphs are indicated in gray.  . Butalene, a molecule where CC models give contradictory predictions for global ring current. This π system is variously predicted to be (a) antiaromatic (with a paratropic ring current), (b) localized, or (c) aromatic (with a diatropic ring current), respectively, depending on whether cycle contributions are taken to be inversely proportional to, independent of, or directly proportional to cycle area in the particular CC model employed to predict the current map.
The Journal of Physical Chemistry A pubs.acs.org/JPCA Article perimeter, models with a = 0 (R, R*) a vanishing perimeter current, and models with a = 1 (GM, GM*, CKCDA, CKCDA*) a diatropic perimeter. The HL and W models both predict a diatropic perimeter, as does a coupled Hartree−Fock ab initio calculation, 60 all in agreement with the CC models with linear dependence of current on cycle area. CC models have also been used for systems containing even numbers of rings of odd size. 11,27,61 In this context, GM, for example, has so far shown a good correlation with HL results, though admittedly for a small sample of graphs. 5.2. General Comments. A protocol has been defined for reduction of the grid-based pseudo-π maps of current density in benzenoids to graphs with flow only along edges, allowing direct comparisons with graph based models for current maps. These comparisons revealed an overall qualitative agreement between pseudo-π and Huckel-London current−density maps for most benzenoids. They also exposed disagreement in two specific circumstances. The first is when the benzenoid graph has short nonbonded contacts, which give rise to false currentbearing edges in the pseudo-π map. This "proximity limitation" is an inevitable consequence of the idealized hexagonal lattice geometry. The second was where pseudo-π maps show poor conservation of current in the neighborhood of formal fixed bonds, as in zethrene, perylene, and analogues, which may be a limitation caused by the single-reference nature of the wave function in the standard pseudo-π method. In all other cases, Huckel-London maps appear to be qualitatively similar to pseudo-π maps, which in turn are excellent proxies for ab initio mapping of currents.
The same ideas have been used to compare Huckel-London current−density maps with those derived from conjugatedcircuit (CC) models, using the sample set of all Kekulean benzenoids with at most ten hexagonal rings (18 360 molecules). The comparison is motivated by the realization that both HL and CC are effectively cycle decompositions of current, albeit with different weighting schemes. All published CC models give a reasonable simulation of HL results on average for most benzenoids, as judged by rms error norms, for example. Limiting ourselves to pure CC approaches, the two models, called GM and R here, would seem to be the best: the GM model because it has the closest match to HL current maps, and R for its simplicity. For GM, the parametrized model in the original form 27 gives an rms error of 8%. Its nearest competitor is the parameter-free version of the model proposed by Mandado,32 which achieves an rms error of 10% in the same test. On the other hand, if simplicity and a transparent relationship to well-studied combinatorial properties of graphs are priorities, then the model proposed by Randic (R) 29 is the best of the CC set; these features allow, for example, an analytical model of perimeter currents in benzenoids based on Pauling Bond Order. 52 Of those discussed here, Model W with α = 0 is the simplest that can treat both Kekulean and non-Kekulean benzenoids.
Finally, there are distinct advantages to going beyond the pure CC model. A new model has been proposed here for combinatorial calculation of current maps of benzenoids. Model W uses tail and next-to-tail coefficients of the characteristic polynomial of the benzenoid and corresponding coefficients for the graphs induced by cycle deletions, all of which can be expressed in terms of perfect matchings (Kekuleś tructures) of vertex-deleted subgraphs. For Kekulean benzenoids, the new model outperforms all CC models by a factor of 2 in reproducing HL maps. It also gives excellent statistics and realistic current maps for non-Kekulean benzenoids, which cannot be modeled at all by conventional CC methods.