Scaling Factors for Models of Currents in Benzenoids

(1)

by

Parminder Kaur

B.Tech. Lovely Professional University 2013

A Master’s Project Submitted in Partial Fullfilment of the Requirements for the Degree of

MASTER OF SCIENCE

in the Department of Computer Science

©Parminder Kaur, 2017 University of Victoria

(2)

Scaling Factors for Models of Currents in Benzenoids

by

Parminder Kaur

B.Tech. Lovely Professional University 2013

Supervisory Committee

Dr. Wendy Myrvold (Department of Computer Science)

Supervisor

Dr. Frank Ruskey (Department of Computer Science)

(3)

Abstract

Conjugated-circuit (CC) methods give purely graph theoretical em-pirical models for estimation of molecular ring currents. This project considers the question of how to scale CC models to give better repro-duction of the quantum mechanical H¨uckel–London model, in order to combine the transparency of the CC approach with the physical inter-operability of the H¨uckel–London method. Scaling is carried such that L1-errors and L∞-errors are minimised for families of linear polyacenes and zigzag fibonacenes.

(4)

List of Tables

1 Conjugated-circuits that result from each pair of matchings. . . . 11

2 Pairs of matchings that contribute to the Mandado current model 15 3 Conjugated circuit models and their currents . . . 18

4 Unscaled current for each cycle for the Randi´c, Ciesielski et.al. and Mandado models . . . 19

5 Two new models proposed by Myrvold and Fowler . . . 21

6 The gi values for the 3-hexagon linear polyacene. . . 22

7 The ci values for the G − C graphs for cycles C1 to C6. . . 23

8 The current contributions for each cycle for the Two-term model. 23 9 The current contributions for each cycle for the 3-hexagon linear polyacene. . . 24

10 The gi values for the 13-vertex benzenoid from Figure 20. . . 28

11 The ci values for the G − C graphs for cycles C1 to C7 for the 13-vertex benzenoid. . . 28

12 The current contributions for each cycle for the Two-term model for the 13 vertex benzenoid. . . 28

13 Scaling factors for linear polyacenes (L∞). . . 34

14 Scaling factors for linear polyacenes (L1). . . 34

15 Scaling factors for zigzag fibonacenes (L∞). . . 35

16 Scaling factors for zigzag fibonacenes (L1). . . 35

17 L∞-errors for linear polyacenes. . . 35

18 L1-errors for linear polyacenes. . . 36

19 L∞-errors for zigzag fibonacenes. . . 36

(7)

List of Figures

1 Linear polyacenes. . . 2

2 Zigzag fibonacenes. . . 2

3 Some equivalent currents. . . 5

4 Examples of currents. . . 5

5 Current c1 on G. . . 6

6 Current c2 on G. . . 6

7 The current difference c1− c2. . . 7

8 Matchings for the 3-hexagon zigzag fibonacene. . . 8

9 Cycles for the 3-hexagon zigzag fibonacene. . . 9

10 Randi´c current weights for each cycle in the 3-hexagon zigzag fibonacene found by vector addition of bond currents. . . 10

11 Sum of currents for all cycles on the 3-hexagon zigzag fibonacene, using the Randi´c current model. . . 11

12 Final currents on the 3-hexagon zigzag fibonacene, using the Randi´c current model. . . 11

13 Current weights in the Ciesielski et. al. model on each cycle in a 3-hexagon zigzag fibonacene. . . 13

14 Sum of currents on each cycle of 3-hexagon zigzag fibonacene using the Ciesielski et. al. current model. . . 14

15 Final current on a 3-hexagon zigzag fibonacene using the Ciesiel-ski et. al. current model. . . 14

16 Mandado current weights on each cycle in the 3-hexagon zigzag fibonacene. . . 16

17 Sum of currents on each cycle of the 3-hexagon zigzag fibonacene, using the Mandado current model. . . 17

18 Final currents on the 3-hexagon zigzag fibonacene, using the Mandado current model. . . 17

19 Perfect matchings for G − C graphs. . . 20

(8)

21 The Two-term model currents for the 3-hexagon linear polyacene. 23 22 The Adjusted Two-term current for the 3-hexagon linear polyacene. 24

23 No perfect matching for G − C7graph. . . 24

24 Cycles in a 13-vertex benzenoid. . . 26

25 Two-term current for Cycles in 13-vertex benzenoid. . . 27

26 The Final Two-term current for 13-vertex benzenoid. . . 29

27 Scaling Mandado currents to H¨uckel–London currents. . . 33

28 L∞-error for zigzag fibonacenes. . . 37

29 L1-error for zigzag fibonacenes. . . 38

30 L1-error for linear polyacenes. . . 38

31 L∞-error for linear polyacenes. . . 39

32 Currents for 3-hexagon Linear Polyacene. . . 40

33 The Two-term model gives a smaller L∞-error than the Adjusted Two-term model for 3-hexagon Linear Polyacene. . . 41

(9)

Acknowledgements

This work has been encouraged and become possible under the guidance of Dr. Wendy Myrvold. Thank you Wendy, for consistently encourage me to do better. I would also like to thank Patrick Fowler for his feedback on project, it was helpful.

I would like to express my gratitude to my best friend Navpreet who uncon-ditionally took care of me during this journey and help me in accomplishing the dream of my parents.

(10)

1 Introduction

Ring currents have been at the heart of quantum chemistry for many decades. Ring currents are circulations of electrons induced in unsaturated carbon frame-works on application of a perpendicular external magnetic field. The Huckel-London model uses an empirical model for the response of a molecular electron distribution to external magnetic fields to calculate the currents that flow in the bonds of unsaturated carbon frameworks. The calculations require specification of both the molecular graph and the positions in space of the carbon centres. Currents are usually reported as ratios of the derivative of current with respect to field to the same quantity calculated for a benzene molecule. Huckel-London currents are typically in qualitative agreement with the results of more sophis-ticated quantum-chemical calculations, but are computed much more cheaply. This paper considers the Hückel–London current model [13] and several models based on conjugated-circuits; the model due to Randić [12], the model due to Ciesielski et. al. [5], the model due to Mandado [9] and the Two-term and the Adjusted Two-term models due to Myrvold and Fowler. The goal of our re-search is to find an approach for scaling conjugated-circuit models so that they provide approximations to Hückel–London currents.

We are considering two families of benzenoids, linear polyacenes which are polycyclic aromatic hydrocarbons consisting of linearly fused benzene rings. The smaller molecules (benzene, naphthalene, anthracene, tetracene,....) are shown in Figure 1. Zigzag fibonacenes (Figure 2) consist of chains of hexagons with a zigzag shape such that no straight segment is of length greater than two. There can be a number of ways zigzags can be plotted, our research considers only the ones which follow the pattern such that the ladder goes up as shown in Figure 2.

(11)

(a) (b) (c)

(d)

Figure 1: Linear polyacenes.

(a) (b) (c) (d)

(e) (f)

Figure 2: Zigzag fibonacenes.

The objective is to find an optimal value of a constant s which could be used to scale the Randi´c, Ciesielski et. al., Mandado, Two-term and Adjusted Two-term to the H¨uckel–London current model for linear polyacenes and zigzag fibonacenes. To achieve this, we use the mathematical optimization technique known as linear programming. Dantzig’s Simplex algorithm [4] is used to extract optimal scaling factors for the different current models.

(12)

2 Background

2.1 Basic Graph Theory Definitions

This paper takes definitions from West’s book [14]. An undirected graph G with n vertices and m edges consists of a vertex set V(G) = { v1, v2...., vn }

and an edge set E(G) = {e1, e2, ..., em}, where each edge corresponds to

an unordered pair of vertices. The notation (u, v) is used for an edge between vertex u and vertex v. A simple directed graph or digraph G consists of a vertex set V(G) and an arc set E(G), where each arc corresponds to an ordered pair of vertices. We write (u, v) ∈ E(G), if there is an arc from u to v. A weighted directed graph is a directed graph with weights assigned to its edges.

A cycle C is a graph with n vertices connected in a closed chain such that the number of vertices in C equals the number of edges, and every vertex has degree two. If G is a graph, then the subgraph induced by S ⊆ V (G) is the graph H where V (H) = S and E(H) = {(u, v) : (u, v) ∈ E(G) and u, v ∈ S}.

A matching in an undirected graph G is a set of pairwise disjoint edges. The vertices belonging to the edges of a matching are saturated by the matching; the others are unsaturated. if a matching saturates every vertex of G, then it is a perfect matching or complete matching. In Chemistry, perfect matchings are called Kekul´e structures. The number of perfect matchings of a graph G is denoted by m(G). For all conjugated-circuit models for currents, m(G) is assumed to be equal to one for a graph with no vertices.

A graph is planar if it can be drawn in the plane in such a way that none of its edges cross each other. If G is a planar graph, then any plane drawing of G divides the plane into regions, called faces. The graphs that correspond to benzenoids that are considered in this paper are planar graphs with no cut vertices, the vertices have degree two or three and all the internal faces are of size six, and also vertices not on the external face have degree three.

Given a molecular graph G, a current in G is represented by a weighted directed graph H such that V (G) = V (H) and if G has edge (u, v) then H has

(13)

either arc (u, v) or (v, u). The weight on arc (u, v) represents the magnitude of the current from u to v.

Eigenvalues are a special set of scalars associated with a linear system of equations (i.e., a matrix equation) that are sometimes also known as character-istic roots. For a square matrix A of order n, the number λ is an eigenvalue if and only if there exists a non-zero vector x such that

Ax = λx.

Using the matrix multiplication properties, we obtain

A − λIn

x = 0.

The equation has a non trivial solution if

detA − λIn

= 0.

The characteristic polynomial of a graph is the characteristic polynomial of its adjacency matrix. If the eigenvalues of the adjacency matrix for G are λ1, λ2, λ3, ..., λn then the characteristic polynomial for G can be written as:

(λ − λ1)(λ − λ2)...(λ − λn).

The characteristic polynomial of G with n vertices is given by:

n

X

i=0

giλi.

Or equivalently, it can be written as:

g0+ g1λ1+ g2λ2+ ... + gnλn

where g0, g1, g2, g3...gn are coefficients of characteristic polynomial.

2.2 Conjugated-Circuits

In graph theory terms, a cycle C is a conjugated-circuit in a molecular graph G if both C and G − C have a perfect matching [6].

(14)

2.3 Error between two current models

In a graph with arc (u, v), a current of magnitude cu,v from u to v is

math-ematically equivalent to a current of magnitude -cv,u from v to u. Figure 3

shows some equivalent currents. The net current on (u, v) is cu,v - cv,u. The

net current can be negative (see Figure 4).

Figure 3: Some equivalent currents.

Figure 4: Examples of currents.

For two models A and B, let cAu,v be the current on arc (u, v) for model A

and let cB

u,vbe the current for model B. The difference A-B of these two currents

has current cA

u,v- cBu,von arc (u,v). The magnitude of the difference on arc (u,v)

is defined to be | cA

u,v− cBv,u | and the direction is from u to v if cu,v > cv,u

and the direction is from v to u if cv,u > cu,v. For a current of zero value,

an arbitrary direction is assigned. Current values c1 and c2 on a graph G are

(15)

- c2.

In conjugated-circuit models, the computed current is the sum of contri-butions from all the conjugated-circuits. For all current models, the direction of the current contribution that corresponds to a cycle C is always the same. If C is a 4n-cycle, the current is clockwise and if C is (4n + 2)- cycle, the current is counterclockwise. This reflects the chemical and physical distinction between anti-aromatic and aromatic carbon monocycles, where the induced cur-rent flows in the diamagnetic sense for the aromatic system, and in the opposite paramagnetic sense for the anti-aromatic system. However conjugated-circuit currents are weighted differently for the different models. After summing the conjugated-circuit currents, because of different weighting schemes, one model can have positive current from u to v where another model may have positive current going from v to u.

Figure 5: Current c1 on G.

(16)

Figure 7: The current difference c1− c2.

2.4 Current Models

As mentioned in introduction we are considering conjugated-circuit models; the Randi´c model, the Ciesielski et. al. model and the Mandado model. This section provides details on their definitions.

2.4.1 The model proposed by Randi´c

To compute an approximation for the current, Randi´c [12] considers all ordered pairs of perfect matchings. For every cycle that results, there is assignment of one unit of current going clockwise if the cycle is a 4n-cycle and counterclockwise if it is a (4n + 2)-cycle. Then, all cycle currents are summed up and scaled by dividing by m(G)(m(G) − 1).

To illustrate the calculation of the Randi´c current, we have taken an example of a 3-hexagon zigzag fibonacene. Figure 8 shows all matchings for this molecular graph.

(17)

(a) Matching M1. (b) Matching M2.

(c) Matching M3.

(d) Matching M4. (e) Matching M5.

(18)

(a) Cycle C1. (b) Cycle C2.

(c) Cycle C3.

(d) Cycle C4. (e) Cycle C5.

(f) Cycle C6.

Figure 9: Cycles for the 3-hexagon zigzag fibonacene.

Figure 9 shows the six cycles C1, C2, C3, C4, C5 and C6 in the 3-hexagon

zigzag fibonacecne. The matrix of matchings in Table 1 is formed by considering the ordered pairs of perfect matchings. The current on an edge e is the sum of these cycle currents for all cycles that include edge e arising by pairing perfect

(19)

(a) Current for C1. (b) Current for C2.

(c) Current for C3.

(d) Current for C4. (e) Current for C5.

(f) Current for cycle C6.

Figure 10: Randi´c current weights for each cycle in the 3-hexagon zigzag fi-bonacene found by vector addition of bond currents.

(20)

M1 M2 M3 M4 M5 M1 × C3 C2 C1 + C3 C1 M2 C3 × C5 C1 C1+ C3 M3 C2 C5 × C6 C4 M4 C1+ C3 C1 C6 × C3 M5 C1 C1+ C3 C4 C3 ×

Table 1: Conjugated-circuits that result from each pair of matchings.

Figure 11: Sum of currents for all cycles on the 3-hexagon zigzag fibonacene, using the Randi´c current model.

Figure 12: Final currents on the 3-hexagon zigzag fibonacene, using the Randi´c current model.

(21)

matchings. For example, cycle C1 appears eight times in the matrix which

results in a current weight of eight units on cycle C1 and other cycles are also

weighted similarly as shown in Figure 10. Figure 11 shows sum of currents on each cycle. Figure 12 shows the final current on the zigzag fibonacene obtained by adding the currents on each edge, and then scaling.

2.4.2 The model proposed by Ciesielski et. al.

The Ciesielski et. al. model [5] also considers ordered pairs of perfect match-ings. The area of a cycle C in a benzenoid, Area(C), is the number of hexagons inside C. Thus, all face areas are specified in terms of a standard hexagon. The Ciesielski et. al. model puts Area(C) units of current on C for each pair of per-fect matchings that includes C. For every cycle, there is assignment of Area(C) units current going clockwise if the cycle is a 4n-cycle and counterclockwise if it is a (4n + 2)-cycle.

(22)

(a) Current for cycle C1. (b) Current for cycle C2.

(c) Current for cycle C3.

(d) Current for cycle C4. (e) Current for cycle C5.

Figure 13: Current weights in the Ciesielski et. al. model on each cycle in a 3-hexagon zigzag fibonacene.

Figure 14 shows sum of currents on each cycle. The final current is the sum of these currents over all the ordered pairs of perfect matchings scaled by dividing by m(G)(m(G) − 1). The final current computed using the Ciesielski et. al. current model for a 3-hexagon zigzag fibonacene is shown in Figure 15.

(23)

Figure 14: Sum of currents on each cycle of 3-hexagon zigzag fibonacene using the Ciesielski et. al. current model.

Figure 15: Final current on a 3-hexagon zigzag fibonacene using the Ciesielski et. al. current model.

2.4.3 The model proposed by Mandado

Mandado [9] considered all pairs of perfect matchings, as did Randi´c but discards the cases that have more than one cycle, see Table 2. For every cycle that results,

1

Area(C) units of currents is assigned to go clockwise if the cycle is a 4n-cycle

and counterclockwise if it is a (4n + 2)-cycle. The Area(C) is the number of hexagons inside C. Then, all cycle currents are summed up on edges in Figure 17 and are scaled by dividing by m(G) as shown in Figure 18.

(24)

M1 M2 M3 M4 M5 M1 × C3 C2 C1+ C3 C1 M2 C3 × C5 C1 C1 + C3 M3 C2 C5 × C6 C4 M4 C1 + C3 C1 C6 × C3 M5 C1 C1+ C3 C4 C3 ×

(25)

(a) Current for C1. (b) Current for C2.

(c) Current for C3.

(d) Current for C4. (e) Current for C5.

Figure 16: Mandado current weights on each cycle in the 3-hexagon zigzag fibonacene.

(26)

Figure 17: Sum of currents on each cycle of the 3-hexagon zigzag fibonacene, using the Mandado current model.

Figure 18: Final currents on the 3-hexagon zigzag fibonacene, using the Man-dado current model.

2.5 Framework for Combinatorial models

Myrvold and Fowler [6] provided a combinatorial framework that encompasses the conjugated-circuit models of Randi´c [12], Ciesielski et. al. [5] and Mandado [9]. A cycle C in the undirected graph G corresponds to two directed cycles in the corresponding directed graph ~G, one that is oriented clockwise and another that is oriented counterclockwise. The current is assigned to the clockwise directed cycle if C is a 4n-cycle and to the counterclockwise directed cycle if C is a (4n + 2)-cycle. For the cycles that are not conjugated-circuits, the

(27)

contribution is zero.

The models use different formulas for the weight of a circuit current. The Randi´c, Ciesieleski et. al. and Mandado models put a weight on each conjugated-circuit which can be formalized in the following framework:

w(C) = 1

S2m(G − C)

p_Area(C)k ₍₁₎

where the value of p is two for the Randić model and Ciesielski models and one for the Mandado model, the value of k is zero for the Randić model, one for the Ciesielski model and −1 for the Mandado model, and the value of S is equal to m(G)(m(G) − 1) for the Randić and Ciesielski models and m(G) for the Mandado model. The value of Area(C) is the number of hexagons inside C. Table 3 shows the models and their conjugated circuit current contributions.

Model Conjugated circuit current Randi´c 2[m(G − C)2_]

Ciesielski et. al. 2[m(G − C)2_]Area(C)

Mandado 2[m(G−C)]_Area(C)

Table 3: Conjugated circuit models and their currents

If G − C has m(G − C) perfect matchings, considering the same example of a 3-hexagon zigzag fibonacene as in Figure 8, we will again calculate Randi´c currents using equation (1). For every cycle that results, we compute the number of matchings for G − C where C can be any cycle C1, C2, C3, C4, C5 or C6 as

shown in Figure 19. We obtain the values of unscaled current on C shown in Table 4.

(28)

Cycle C m(G − C) Randi´c model Ciesielski et.al. model Mandado model C1 2 8 8 4 C2 1 2 2 2 C3 2 8 8 4 C4 1 2 4 1 C5 1 2 4 1 C6 1 2 6 2/3

Table 4: Unscaled current for each cycle for the Randi´c, Ciesielski et.al. and Mandado models

The Randi´c weights computed in Table 4 are the same as shown in Figure 10. Then, all cycle currents are summed up and final current is shown in Figure 12 by scaling by multiplying by 1

m(G)(m(G)−1), where m(G) = 5, is the number

(29)

(a) One perfect match-ing for G-C1. (b) A second perfect matching for G-C1. (c) A unique perfect matching for G-C2.

(d) One perfect match-ing for G-C3.

(e) A second perfect matching for G-C3. (f) A unique perfect matching for G-C4. (g) A unique perfect matching for G-C5. (h) A unique perfect matching for G-C6.

Figure 19: Perfect matchings for G − C graphs.

2.5.1 Two-term and Adjusted Two-term models proposed by Myr-vold and Fowler

The Two-term and Adjusted Two-term models [11] were inspired by Aihara’s formulas [1], [2], [8] for H¨uckel–London currents. These models are not published yet.

As mentioned in Section 2.1, the characteristic polynomial of a graph can be written as:

g0λ0+ g1λ1+ g2λ2+ ... + gnλn

(30)

benzenoids, g0 (the constant term) is equal to (−1)n/2m(G) m(G). Hence, g0

is always non-zero for graphs with perfect matchings. Let r0 be the minimum

value of i such that gi 6= 0.

The weight for the current for a Circuit C is computed as shown in Table 5.

Two-term model w(C) = 2 ∗ Area(C) ∗ [cr0/gr0+ cr0+2/gr0+2]

Adjusted Two-term model w(C) = 2 ∗ Area(C) ∗ [cr0/gr0+ 4 ∗ cr0+2/gr0+2]

Table 5: Two new models proposed by Myrvold and Fowler

If C is a (4n)-Cycle, w(C) will be positive and w(C) units are assigned to C in the clockwise direction. If C is a (4n+2)-cycle, w(C) will be less than or equal to zero and | w(C) | units are assigned to C in the counterclockwise direction. One approach to find eigenvalues and hence also characteristic polynomials is to use Jacobi’s method [10].

2.5.2 Examples of the Two-term and Adjusted Two-term models Consider the 3-hexagon linear polyacene with cycles C1, C2, C3, C4, C5 and C6

as shown in Figure 20. The coefficients of the characteristic polynomial of G are given in Table 6. The coefficients of the characteristic polynomial for each G−C graph are given in Table 7. Tables 8 and 9 shows the contribution of current for each cycle and finally Figures 21 and 22 show final currents corresponding to the Two-term and Adjusted Two-term models by summing up contribution of each cycle over every edge.

(31)

(a) Cycle C1. (b) Cycle C2.

(c) Cycle C3. (d) Cycle C4.

(e) Cycle C5. (f) Cycle C6.

Figure 20: Cycles in 3-hexagon linear polyacenes.

g0 g1 g2 g3 g4 g5 g6 g7

-16 0 148 0 -392 0 473 0

g8 g9 g10 g11 g12 g13 g14

-296 0 98 0 -16 -0 1

(32)

Graph ↓ ci→ c0 c1 c2 c3 c4 c5 c6 c7 c8 C1 1 0 -13 0 18 0 -8 0 1 C2 1 0 -6 0 11 0 -6 0 1 C3 1 0 -13 0 18 0 -8 0 1 C4 1 0 -3 0 1 C5 1 0 -3 0 1 C6 1

Table 7: The ci values for the G − C graphs for cycles C1 to C6.

Graph w(C) = 2 ∗ Area(C) ∗ [cr0/gr0+ cr0+2/gr0+2] computed w(C)

C1 2 ∗ 1 ∗ [(1/ − 16) + (−13/148)] -0.3006757 C2 2 ∗ 1 ∗ [(1/ − 16) + (−6/148)] -0.2060811 C3 2 ∗ 1 ∗ [(1/ − 16) + (−13/148)] -0.3006757 C4 2 ∗ 2 ∗ [(1/ − 16) + (−3/148)] -0.3310811 C5 2 ∗ 2 ∗ [(1/ − 16) + (−3/148)] -0.3310811 C6 2 ∗ 3 ∗ [(1/ − 16) + (0/148)] -0.375

Table 8: The current contributions for each cycle for the Two-term model.

(33)

Graph w(C) = 2 ∗ Area(C) ∗ [cr0/gr0+ 4 ∗ cr0+2/gr0+2] computed w(C) C1 2* 1 * [ (1 / -16) + 4 * (-13 / 148) ] 0.8277027 C2 2* 1 * [ (1 / -16) + 4 * (-6 / 148)] 0.4493243 C3 2* 1 * [ (1 / -16) + 4 * (-13 / 148) ] 0.8277027 C4 2* 2 * [ (1 / -16) + 4 * (-3 / 148)] 0.5743243 C5 2* 2 * [ (1 / -16) + 4 * (-3 / 148) ] 0.5743243 C6 2* 3 * [ (1 / -16) + 4 * (0 / 148) ] 0.375

Table 9: The current contributions for each cycle for the 3-hexagon linear poly-acene.

Figure 22: The Adjusted Two-term current for the 3-hexagon linear polyacene.

These two new models predict currents for molecules with no perfect match-ings where all existing CC models give zero current for those molecules. The 13-vertex benzenoid shown in Figure 23 is example of such a molecule.

Figure 23: No perfect matching for G − C7 graph.

(34)

For i = 1,2, . . .,7 , G–Ci has an odd number of vertices and hence m(G–Ci)

= 0. The Two-term model uses characteristic polynomials, the eigenvalues and the characteristic polynomials for G − C for each Cycle C to determine current. Using these terms in Two-term formula given in Table 5, the final current on the 13-vertex benzenoid is given in Figure 26. The resultant givalues are shown

in Table 10. Table 11 shows the resultant values of ci. The current contribution

for each cycle are given in Table 12 and the direction of current on each cycle is shown in Figure 25.

(35)

(a) Cycles C1, C2, C3. (b) Cycle C4.

(c) Cycle C5. (d) Cycle C6.

(e) Cycle C7.

(36)

(a) Two-term current for Cycle C1.

(b) Two-term current for Cycle C2.

(c) Two-term current for Cycle C3.

(d) Two-term current for Cycle C4.

(e) Two-term current for Cycle C5.

(f) Two-term current for Cycle C6.

(g) Two-term current for Cycle C7.

(37)

g0 g1 g2 g3 g4 g5 g6 g7

0 54 0 -207 0 309 0 -226

g8 g9 g10 g11 g12 g13

0 84 0 -15 0 1

Table 10: The gi values for the 13-vertex benzenoid from Figure 20.

Graph ↓ ci → c0 c1 c2 c3 c4 c5 c6 c7 G − C1 0 -4 0 10 0 -6 0 1 G − C2 0 -4 0 10 0 -6 0 1 G − C3 0 -4 0 10 0 -6 0 1 G − C4 0 -2 0 1 G − C5 0 -2 0 1 G − C6 0 -2 0 1 G − C7 0 1

Table 11: The civalues for the G−C graphs for cycles C1to C7for the 13-vertex

benzenoid.

Cycles w(C) = 2 ∗ Area(C) ∗ [cr0/gr0+ cr0+2/gr0+2] computed w(C)

C1, C2, C3 2 ∗ 1 ∗ [(−4/54) + (10/ − 207)] -0.2447665

C4, C5, C6 2 ∗ 2 ∗ [(−2/54) + (1/ − 207)] -0.1674718

C7 2 ∗ 3 ∗ [(1/54) + (0/ − 207)] 0.1111111

Table 12: The current contributions for each cycle for the Two-term model for the 13 vertex benzenoid.

(38)

Figure 26: The Final Two-term current for 13-vertex benzenoid.

3 Scaling Current Models

Suppose ~G is a current graph. Scaling ~G by a constant s results in a new current such that for each arc (u, v) with the current cu,v in ~G, the scaled current is s

cu,v. If the goal is to scale a current B to the best match a current A, then the

error on edge (u, v) is equal to | CA(u, v) − CB(u, v) |. The L1- error is defined

to be the sum of the errors on each edge of the graph G. The L∞- error is

defined to be maximum of the errors on each edge of the graph G.

4 Finding Scaling factors using Linear

Program-ming

Formally, a linear programming problem is the problem of maximizing (or minimizing) a linear function subject to a finite number of linear constraints. The standard form of a linear programming problem [4] is

Maximize cT_x

subject to Ax ≤ b x ≥ 0

(39)

where A is a n × m matrix of coefficients, c and b are vectors and x is the vector of variables to be found.

Suppose G1 and G2 are directed graphs with currents cA and cB flowing on

an arc (u, v), computed by models A and B respectively. Let us assume that cA _{is equal to the value of the current to match and c}B _{is the current we want}

to scale to get an optimal match to cA_{. The next two sections present the}

linear programming problems used to scale cB _{while minimizing the L} ∞-error

(section 4.1) and the L1-error (section 4.2).

4.1 Minimizing L

∞

-error

To minimize the L∞-error while scaling model B to match model A, using the

L∞approach:

The number of variables will be equal to two: an error r and a constant s which is the scaling factor. The direction of current is same for benzenoids considered in this paper so s is greater than zero or zero. The number of equations will be equal to twice number of edges of the graph.

The problem can be formulated as follows Minimize r

subject to

for each edge (u, v) in the undirected graph G, r ≥| cAu,v− s cBu,v| and

r, s ≥ 0.

The standard form of problem is: Maximize −r

subject to:

for each edge (u, v) in the undirected graph G, −r − scB

u,v≤ −cAu,v,

−r + scB

(40)

r, s ≥ 0.

4.2 Minimizing L

1

-error

To minimize the total absolute error while scaling model B to match model A by using the L1 approach:

The edges of G are numbered as e1, e2, e3...em. For each edge ei, there is a

variable ri that is the error for edge i. There is also the variable s (the scaling

factor). The direction of current is same for benzenoids considered in this paper so s is greater than zero or zero. The number of variables will be equal to the number of edges of the graph plus one. The number of equations will be equal to twice the number of edges of the graph.

The problem can be formulated as follows: Minimize r1 + r2+ r3+. . .+ rm

subject to:

for each edge (u, v) in the undirected graph G, ri ≥| cAu,v− scBu,v| and

r1, r2, r3. . ., rm, s ≥ 0.

The standard form for the L1-norm for scaling current is:

Maximize -(r1 + r2 + r3+ . . .+ rm)

subject to:

for each edge (u, v) in the undirected graph G −ri− scBu,v≤ −cAu,v,

(41)

4.3 Example

Suppose G1 and G2are directed graphs with currents cA and cB, computed

using the H¨uckel–London model and the Mandado model respectively as shown in Figures 27 (a) and (b) respectively.

The optimal scaling factor for the L1-norm is 1.183 and for the L∞-norm is

1.097. For each scaling factor, we multiply the values of current generated by Mandado model (cB_{) by the constant s (scaling factor) to get the scaled}

current (s × cB). Now, using the H¨uckel–London current, cA and the scaled Mandado current scB, the difference is computed as shown in Figures 27(e) and 27(f) for the L1-norm and L∞-norm. The L1-error is 0.605 and the

(42)

(a) H¨uckel–London current, cA_. _{(b) Mandado current, c}B_.

(c) Mandado current scaled by multiplying by 1.183 (optimal so-lution using the L1-norm).

(d) Mandado current scaled by multiplying by 1.097 (optimal so-lution using the L∞-norm).

(e) The difference between H¨uckel–London current and scaled the Mandado current from (c).

(f) The difference between H¨uckel–London current and scaled Mandado current from (d).

(43)

5 Results

For the Mandado, Randi´c, Ciesielski et. al., Two-term and Adjusted Two-term models, we have calculated scaling factors for the linear polyacenes and zigzag fibonacenes using the L∞ and L1 approaches and the results are shown in

Tables 13 and 16 respectively. The errors are shown in Tables 17-20. The cases which gave the best match to H¨uckel–London are highlighted in bold.

Number of hexagons → 1 2 3 4 5 6 7 8 9

Model ↓

Randi´c 1.000 1.639 1.919 2.176 2.438 2.712 2.927 3.085 3.172 Ciesielski et. al. 2.000 2.185 1.919 1.740 1.625 1.550 1.464 1.371 1.269 Mandado 1.000 1.093 1.097 1.119 1.158 1.211 1.254 1.281 1.290 Two-term 2.000 1.355 1.029 0.844 0.727 0.648 0.592 0.551 0.508 Adjusted Two-term 2.000 0.892 0.648 0.504 0.426 0.350 0.307 0.277 0.254

Table 13: Scaling factors for linear polyacenes (L∞).

Model ↓

(44)

Number of hexagons → 1 2 3 4 5 6 7 8 9 Model ↓

Table 15: Scaling factors for zigzag fibonacenes (L∞).

Model ↓

Randi´c 1.000 1.639 1.894 2.151 2.255 2.679 2.706 2.754 2.769 Ciesielski et. al. 2.000 2.185 2.526 2.689 2.618 2.670 2.720 2.604 2.618 Mandado 1.000 1.093 1.003 1.006 0.954 0.936 0.930 0.861 0.856 Two-term Model 2.000 1.355 1.154 1.038 0.936 0.888 0.856 0.833 0.789 Adjusted Two-term Model 2.000 0.892 0.638 0.526 0.466 0.424 0.395 0.376 0.358

Table 16: Scaling factors for zigzag fibonacenes (L1).

Model ( Error) ↓

Randi´c (L∞error) 0.000 0.000 0.125 0.198 0.243 0.274 0.312 0.356 0.405

Ciesielski et. al. (L∞error) 0.000 0.000 0.125 0.198 0.243 0.274 0.312 0.356 0.405

Mandado (L∞error) 0.000 0.000 0.079 0.136 0.174 0.201 0.231 0.268 0.309

Two-term (L∞error) 0.000 0.000 0.048 0.096 0.138 0.173 0.201 0.223 0.260

Adjusted Two-term (L∞) 0.000 0.000 0.068 0.045 0.046 0.055 0.099 0.136 0.166

(45)

Number of hexagons → 1 2 3 4 5 6 7 8 9 Model ( Error) ↓

Table 18: L1-errors for linear polyacenes.

Model ( Error) ↓

Table 19: L∞-errors for zigzag fibonacenes.

Model ↓

Randi´c 0.000 0.000 1.303 1.457 2.500 2.584 2.982 3.019 3.267 Ciesielski et. al. 0.000 0.000 0.213 0.848 1.325 1.621 1.970 2.271 2.374 Mandado 0.000 0.000 0.233 0.429 0.914 1.361 1.886 2.369 2.616 Two-term (L1error) 0.000 0.000 0.083 0.345 0.657 0.888 1.155 1.420 1.567 Adjusted Two-term 0.000 0.000 0.074 0.172 0.125 0.252 0.466 0.727 0.942

(46)

5.1 The best matches to H¨

uckel–London

Figures 28-31 gives the L1- errors and L∞- errors among all conjugated-circuit

models as the number of hexagons increases for both linear polyacenes and zigzag fibonacenes except for one case that is shown in Section 5.1.1 where the Two-term model is best. The Adjusted Two-term model has minimum L1 and

L∞-errors among all models discussed in this paper. Tables 19 and 20 show

that Mandado has minimum L1-error and L∞-error as the number of hexagon

increases for linear polyacenes. Tables 17 and 18 show, Mandado also gives minimum L1-error and L∞-error for zigzag fibonacenes as the number of

hexagon increases, although there are cases when Ciesielski gives the minimum or equivalent L1 or L∞-errors. For the 3-hexagon zigzag fibonacene, Ciesielski

gives minimum L1 and L∞-error and for the 9-hexagon zigzag fibonacene,

Mandado and Ciesielski gave equivalent L∞-error but Ciesielski gives

minimum L1-error.

(47)

Figure 29: L1-error for zigzag fibonacenes.

(48)

Figure 31: L∞-error for linear polyacenes.

5.1.1 A case where the Two-term model gives less L∞-error than

the Adjusted Two-term model

For the 3-hexagon linear polyacenes, the Two-term model gives less L∞-error

than the Adjusted Two-term model as shown in Figure 33. Figure 32(a) shows the H¨uckel–London current, Figure 21 shows the Two-term current and Figure 22 shows the Adjusted Two-term current. The scaling factor for the Two-term model is 1.029 and the scaling factor for Adjusted Two-term model is 0.648. The difference in both models is the factor 4 in the Adjusted Two-term model. Possibly there is too much weight on the second term of the Two-term model for this case.

(49)

(a) H¨uckel–London currents. (b) The scaled Two-term current.

(c) The scaled Adjusted Two-term current.

(50)

(a) The L∞−error for Two-term

cur-rent model.

(b) The L∞ − error for Adjusted

Two-term current model.

Figure 33: The Two-term model gives a smaller L∞-error than the Adjusted

Two-term model for 3-hexagon Linear Polyacene.

6 Conclusions and Future work

This paper work shows results which are computed using a set of software. For existing conjugated-circuit models currents were computed using Bill Bird’s project software [3]. Myrvold and Fowler provided their computations for the new models [11]. The Simplex method, generation of the benzenoid families, formulation of Linear Programming problems for families, scaled currents, errors between current models were programmed by the author and the linking of tasks was automated. The scaling factors and errors were checked

independently by Myrvold. These results can be verified using SCIP [7] for future work.

The results shows that the Adjusted Two-term model gives better

approximations to the H¨uckel–London model except for the one case discussed in Section 5.1.1. It is still debatable if the Two-term should be adjusted with a different value, possibly dependent on on some variable n or m(G). Do these models show consistent patterns of internal or external ring currents? Do they show common high or low currents in a certain pattern? Also, It is hard to say

(51)

from this research which error measure is better L1-error or L∞- error as the

model which gave the best or the worst results in L1approach was the same

for the L∞.

References

[1] Aihara. J. Journal of the American Chemical Society, 101:558, 1979.

[2] Aihara. J. Journal of the American Chemical Society, 101:5913, 1979.

[3] Bill Bird. Comparison of combinatorial models for ring currents. Honours Project, University of Victoria, pages 1–24, 2011.

[4] Vaˇsek Chv´atal. Linear programming. W.H. Freeman, New York, 1983.

[5] Arkadiusz Ciesielski, Tadeusz M. Krygowski, Micha l K. Cyra´nski,

Micha l A. Dobrowolski, and Jun-ichi Aihara. Graph-topological approach to magnetic properties of benzenoid hydrocarbons. Physical chemistry chemical physics : PCCP, 11(48):11447, 2009.

[6] Patrick W. Fowler and Wendy Myrvold. The “anthracene problem”: closed-form conjugated-circuit models of ring currents in linear

polyacenes. The journal of physical chemistry. A, 115(45):13191–13200, 2011.

[7] Gerald Gamrath, Tobias Fischer, Tristan Gally, Ambros M. Gleixner, Gregor Hendel, Thorsten Koch, Stephen J. Maher, Matthias

Miltenberger, Benjamin M¨uller, Marc E. Pfetsch, Christian Puchert, Daniel Rehfeldt, Sebastian Schenker, Robert Schwarz, Felipe Serrano, Yuji Shinano, Stefan Vigerske, Dieter Weninger, Michael Winkler, Jonas T. Witt, and Jakob Witzig. The scip optimization suite 3.2. Technical Report 15-60, ZIB, Takustr.7, 14195 Berlin, 2016.

[8] Horikawa J., Aihara. T. Bulletin of the Chemical Society of Japan, 56:1853, 1983.

(52)

[9] Marcos Mandado. Determination of London Susceptibilities and Ring Current Intensities using Conjugated Circuits. Journal of chemical theory and computation, 5(10):2694–2701, 2009.

[10] Gerard Meurant. Computer Solution of Large Linear Systems. 1999.

[11] Myrvold and Fowler. Personal Communication. 2017.

[12] Randi´c. Graph theoretical approach to π-electron currents in polycyclic conjugated hydrocarbons. Chemical Physics Letters,

115(1-3):13191–13200, 2011.

[13] Lionel Salem. The molecular orbital theory of conjugated systems. W.A. Benjamin, Reading, Mass, 1974.

[14] Douglas B. West. Introduction to graph theory. Prentice Hall, Upper Saddle River, NJ, 1996.

Scaling Factors for Models of Currents in Benzenoids

Contents

List of Tables

List of Figures

Acknowledgements

1

Introduction

2

Background

2.1

Basic Graph Theory Definitions

2.2

Conjugated-Circuits

2.3

Error between two current models

2.4

Current Models

2.5

Framework for Combinatorial models

3

Scaling Current Models

4

Finding Scaling factors using Linear

Program-ming

4.1

Minimizing L

-error

4.2

Minimizing L

-error

4.3

Example

5

Results

5.1

The best matches to H¨

uckel–London

6

Conclusions and Future work

References