The Localization–Delocalization Matrix and the Electron-Density

Oct 24, 2014 - Department of Chemistry, Saint Mary's University, Halifax, Nova Scotia ... Hunter College and the Graduate School, City University of N...
0 downloads 0 Views 3MB Size
Subscriber access provided by NORTH DAKOTA STATE UNIV

Article

The Localization-Delocalization Matrix and the Electron Density-Weighted Connectivity Matrix of a Finite Graphene Nanoribbon Reconstructed from Kernel Fragments Matthew J Timm, Cherif F. Matta, Lou Massa, and Lulu Huang J. Phys. Chem. A, Just Accepted Manuscript • DOI: 10.1021/jp508490p • Publication Date (Web): 24 Oct 2014 Downloaded from http://pubs.acs.org on October 28, 2014

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

The Journal of Physical Chemistry A is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 35

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

The Localization-Delocalization Matrix and the Electron Density-Weighted Connectivity Matrix of a Finite Graphene Nanoribbon Reconstructed from Kernel Fragments Matthew J. Timm,a Chérif F. Matta,a,b,c,* Lou Massa,d Lulu Huange (a)

Department of Chemistry and Physics, Mount Saint Vincent University, Halifax, Nova Scotia, Canada B3M2J6. (b) Department of Chemistry, Dalhousie University, Halifax, Nova Scotia, Canada B3H4J3. (c) Department of Chemistry, Saint Mary's University, Halifax, Nova Scotia, Canada B3H3C3. (d) Hunter College and the Graduate School, City University of New York, New York, NY 10065, USA. (e) Center for Computational Materials Science, Naval Research Laboratory, Washington, DC 20375-5341, USA. * Corresponding author: Tel.: +1-(902)-457-6142, Fax: +1-(902)-457-6134, E-mail: [email protected].

Abstract: Bader’s quantum theory of atoms in molecules (QTAIM) and chemical graph theory, merged in the localization-delocalization matrices (LDMs) and the electron density-weighted connectivity matrices (EDWCM), are shown to benefit in computational speed from the kernel energy method (KEM). The quantum chemical graph matrices (LDM and EDWCM) of a 66 atom C46H20 hydrogen-terminated armchair graphene nanoribbon, in 14 (2×7) rings of C2v symmetry, are accurately reconstructed from kernel fragments (this includes the full sets of electron densities at 84 bond critical points (BCPs) and 19 ring critical points (RCPs), and the full sets of 66 localization and 4290 delocalization indices (LIs and DIs)). The average absolute deviations between KEM and directly-calculated atomic electron populations, obtained from the sum of the LIs and ½ of the DIs of an atom, are 0.0012±0.0018 e− (~ 0.02±0.03%) for carbon atoms and are 0.0007±0.0003 e− (~ 0.01±0.01%) for hydrogen atoms. The integration errors in the total electron population (296 electrons) are +0.0003 e− for the direct calculation (+0.0001%) and +0.0022 e− for KEM (+0.0007%). The accuracy of the KEM matrix elements is, thus, probably of the order of magnitude of the combined precision of the electronic structure calculation and the atomic integrations. KEM appears capable of delivering not only the total energies with chemical accuracy (which is well documented), but also local, and non-local properties accurately including the DIs between the fragments (crossing fragmentation lines). Matrices of the intact ribbon, the kernels, the KEM-reconstructed ribbon, and of errors are provided as Supporting Information. Keywords: Graphene, kernel energy method (KEM), chemical graph theory, molecular graph, quantum theory of atoms in molecules (QTAIM), electron localization index, electron delocalization index, bond critical points, fragment calculations, linear scaling.

1 ACS Paragon Plus Environment

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

1.

Introduction

Electronic fingerprinting using the quantum theory of atoms in molecules' (QTAIM)1-3 descriptors of bonding and of electron delocalization (sharing) have recently been shown to be useful in the prediction of a number of molecular properties.4-7 The organized listing of bonding and of electron delocalization descriptors derived from QTAIM in matrix format allows for their analysis with the mathematical apparatus of chemical graph theory.8-15 In chemical graph theory, matrix representatives of the molecular graph, whether hydrogen-suppressed or not, are manipulated mathematically to extract graph invariants that do not depend on atom numbering for the use in predictive modeling of physicochemical properties. On the other hand, QTAIM defines local molecular descriptors of bonding such as the bond critical point (BCP) data between every pair of chemically-bonded atoms, and non-local descriptors such as the delocalization (or electron sharing) index16-21 that counts the number of electrons delocalized between every pair of atomic basins in a molecule (whether bonded or not) or the inter-atomic components of “interacting quantum atoms” (IQA) energies.22-26 Atom-atom pairwise connectivity which yields the molecular graph is established within QTAIM1-3 by the presence of a topological feature of the electron density [ρ(r)], namely, the bond path.1,27-29 The bond path is a line of maximum electron density linking the nuclei that share an interatomic zero-flux surface and hence are chemically bonded. There exists one, and only one, bond path between any two chemically bonded atoms regardless of the nature or strength of this bonding, a strength reflected in the BCP local properties. The bond path is an allor-none topological feature, and hence when an associated bond property (e.g. the electron density at the bond critical point (BCP)) is organized in a matrix format, the matrix has non-zero entries only in those non-diagonal elements that belong to bonded pairs of atoms and are zero otherwise. Thus, an electron-density weighted connectivity matrix (EDWCM) can be defined similar to the chemical graph theory connectivity matrix but instead of the ones used to indicate bonding the EDWCM multiplies the ones by the electron density at the BCP (ρb).4 A description of electron sharing between atomic basins Ωi and Ωj is the QTAIM delocalization index (or DI, δ(Ωi,Ωj)),16 a non-local property since it is calculated by integrating the shared density over each of the two atomic basins. The DI measures, thus, the number of electrons delocalized or shared between two atomic basins Ωi and Ωj, and for a closed-shell 2 ACS Paragon Plus Environment

Page 2 of 35

Page 3 of 35

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

molecule is defined within Hartree-Fock theory as:16

δ ( Ωi , Ω j ) = 2 F α ( Ωi , Ω j ) + 2 F β (Ω i , Ω j ) ,

(1)

where occ occ

F σ (Ωi , Ω j ) = −∑∑ ∫

∫ ϕ (r )ϕ (r )ϕ (r )ϕ (r )dr dr * k

1

l

1

* l

2

k

2

1

2

(2)

l Ωi Ω j

k

occ occ

= −∑∑ Skl (Ω i ) Slk (Ω j ) k

l

(3) is the Fermi correlation, and Skl(Ωi) (which is equal to Slk(Ωi)) is the overlap integral of two spin orbitals φk and φl over the volume of the atomic basin Ωi, and σ can refer to either α or β spins. For post-Hartree-Fock methods the φ’s are natural orbitals and the double sums must, therefore, also include occupation numbers. In the equations above, when i = j, both integrals are over the same atomic basin, say the th

i , giving the Fermi correlation for the electrons in Ωi which, at the limit of total localization (the limit of no delocalization to other atomic basins), approaches the negative of the σ-spin population of atom Ωi, −Nσ(Ωi). This limit is never reached exactly in any system since electrons in Ωi always exchange to some extent with electrons in other atoms in the molecule, and hence generally |Fα(Ωi, Ωi)| ≤ Nα(Ωi) where the equality is only realized for atoms at the infinite separation limit. Thus, a localization index (or LI, Λ(Ω)) can be defined as:16 Λ (Ω i , Ω i ) = F α (Ωi , Ωi ) + F β (Ω i , Ω i ) .

(4)

Since electrons are either localized within the ith atomic basin or delocalized between that basin and every other basin in the molecule, the LI of that atom plus half of the sum of all DIs for this atom (totalling n − 1, where n is the number of atoms in the molecule), yields the total electron population of that atom, N(Ωi), a bookkeeping of electrons concisely expressed by:16

N (Ωi ) = Λ(Ωi ) +

1 n ∑ δ(Ωi , Ω j ) . 2 j ≠i

(5)

The atomic electron population N(Ωi), which can be obtained in this manner or through the direct integration of the electron density over the volume of the atomic basin (usually bound from the exterior by the ρ = 0.001 atomic units (a.u.) isodensity envelope), determines the atomic charge, which is expressed in a.u.: 3 ACS Paragon Plus Environment

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

q(Ωi ) = Z Ωi − N (Ωi )

Page 4 of 35

(6)

where Z Ωi is the atomic number. The total number of electrons in the molecule N can, hence, be expressed as the sum of two separate populations: A localized electron population (Nloc) and a delocalized electron population (Ndeloc) since the sum of all N(Ωi) equals N, that is: n

n

i =1

i =1

N = ∑ N (Ωi ) = ∑ Λ(Ωi ) +

1 n n ∑∑ δ(Ωi , Ω j ) = Nloc + Ndeloc , 2 i =1 j ≠i

(7)

where n

N loc ≡ ∑ Λ (Ω i ) ,

(8)

i =1

and

Ndeloc ≡

1 n n ∑∑ δ(Ωi , Ω j ) = N − tr (ζ) = N − Nloc . 2 i =1 j ≠i

(9)

This paper demonstrates that one can predict a local property such as the electron density at the BCP (ρ(rBCP)), integrated atomic properties such as the atomic charges and the localization indices, and even a non-local property such as delocalization index of a large molecule from fragments. The test presented in this paper is stringent for a fragmentation method since it is applied to a delocalized two-dimensional electronic system (graphene) to predict delocalization indices and bond critical point electron densities including those that cross from one fragment to another. A number of soft-scaling electronic structure calculations are known to circumvent the steep numerical scaling of electronic structure calculations with a quasi-linear scaling.30-56 A particularly accurate approximation is known as the kernel energy method (KEM). This method is well tested and established for the fast and accurate prediction of ab initio energies (total energies, binding energies, interaction energies, etc.) of aperiodic macromolecules such as proteins and DNA.57-62 More recently, KEM has been shown to deliver the energies63 and the external electric field-induced dipole moments, changes in the energies, and polarizabilities of a delocalized two-dimensional (2D) electronic systems, namely, a finite hydrogen-terminated graphene armchair nanoribbon.64 A strength of the method is that of delivering these properties for aperiodic finite systems that can include local defects and/or impurities or dopants, in contrast to regular lattice systems which can be addressed by periodic calculations. 4 ACS Paragon Plus Environment

Page 5 of 35

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

KEM uses molecular geometries (e.g. obtained from crystallography or from a molecular dynamics optimization) to estimate the total energy by breaking a large molecule into “kernel” fragments. Taken in pairs all “double kernels” account for interactions between single kernels however remote from one another. The total energy is calculated by summing over that of the double kernels reduced by that of the single kernels which have been over counted:57-62 m −1

EKEM = ∑

m



a =1 b = a +1

m

Eab − ( m − 2)∑ Ec

(10)

c =1

where EKEM is the KEM energy of the full system, Eab is the energy of the abth double kernel, Ec is the energy of the cth single kernel, and m the number of single kernels. The computational difficulty of the KEM is only quadratic in the number of single kernels representing a molecule (See Ref. 64).

2.

The hypothesis: Pairwise-local and non-local properties reconstructed

from kernel fragments An hypothesis verified numerically in this paper is that the density of the intact full molecule can be reconstructed from kernel fragments at the critical points of the density where the gradient vanishes (other than at the nuclear critical points). These points include the bond and ring critical points, respectively abbreviated as BCPs and RCPs. The reconstruction is achieved through the application of the following approximation: m −1

m

ρ KEM (rBCP( Ω |Ω ) ) = ∑ ∑  ρ ab (rBCP( Ω |Ω ) )  i

j

a =1 b = a +1

i

j

m

Ωi ∧Ω j∈Kab

− (m − 2)∑  ρ c (rBCP( Ωi |Ω j ) )  c =1

, Ωi ∧Ω j∈Kc

(11) where ρ KEM (rBCP( Ωi |Ω j ) ) is the KEM-reconstructed approximate electron density at a given BCP in the full molecule between atoms Ωi and Ωj separated by (and sharing) the interatomic surface of zero-flux in the gradient field of the electron density symbolized by the vertical bar “|”, ρab and ρc are the electron densities at the BCPs of the abth double kernel (Kab) and of the cth single kernel (Kc), and m is the number of single kernels. The condition for including a kernel’s contributions to a sum in the R.H.S. of Eq. 11 is, of course, that the given BCP exists in that kernel, otherwise the kernel’s contribution is zero. This latter point distinguishes local properties 5 ACS Paragon Plus Environment

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 6 of 35

that depend on two particular atoms such as critical point densities from global properties such as the total energy, whereby all fragments have non-vanishing contributions. This paper provides numerical evidence that Eq. 11 is physically reasonable by accurately reconstructing the bond critical point densities (BCP) and ring critical point (RCP) densities of a graphene ribbon. The electron density at the BCP is a pairwise property since it involves two atoms, but is “local” since it is evaluated at a precise point in three-dimensional space located on the bond path connecting the nuclei of these two atoms. In contrast, the electron delocalization index is a non-local property since it results from a double integral involving any two atoms in the molecule whether in direct contact (bonded, share a bond path) or not. This part of the tested hypothesis is particularly challenging considering that the DI, a pairwise non-local property, is defined between any pair of atoms in the molecule – no matter how distant – and can involve two atoms separated by the fission fragmentation lines separating the kernel fragments. The KEM reconstruction of QTAIM’s DIs is expressed mathematically as: m−1

δ KEM (Ωi , Ω j ) = ∑

m



a =1 b = a +1

δ ab (Ωi , Ω j ) 

m

Ωi ∧Ω j∈Kab

− (m − 2)∑ δ c (Ωi , Ω j )  c =1

,

(12)

Ωi ∧Ω j∈Kc

while that of the LIs as: m−1

Λ KEM (Ω i ) = ∑

m

∑ [Λ

a =1 b = a +1

3.

m

ab

(Ω i ) ]Ω ∈K − ( m − 2)∑ [ Λ ab (Ω i ) ]Ω ∈K . i

ab

c =1

i

c

(13)

The definition of QTAIM-quantum chemical graph (QCG) matrices

The QTAIM counterpart of chemical graph theory’s (CGT) edges (or bonds) are the bond paths. The connection between QTAIM and chemical graph theory is brought about by the molecular graph as stressed by Bader1 and others,14,65 and is exploited by Jenkins, Kirk and coworkers to create predictive QTAIM phase diagrams.65-69 While most QTAIM molecular graphs reproduce the bonding patterns assigned on the basis of classical chemical arguments, discrepancies are known especially in borderline cases of weak and non-classical modes of bonding. Molecular graphs of QTAIM and CGT thus generally coincide but occasionally differ. An example of such differences occur in systems exhibiting

6 ACS Paragon Plus Environment

Page 7 of 35

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

hydrogen-hydrogen weak interactions between closed-shell hydrogen atoms70,71 as the adjacent bay hydrogen atoms of the hydrogen-terminated graphene ribbon. Therefore, there can be topological differences between connectivity matrices based on bond paths and those based on classical chemical structures. QTAIM non-nuclear attractors72-76 and groups condensed in “super-atoms”5 can also be treated like an atom, a difference with CGT connectivity matrices. There are four possible critical points (CP) (where the ∇ρ = 0) in the 3D electron density distribution: Nuclear critical points (NCP) exhibiting three negative curvatures since ρ is a local maximum (includes nuclei and non-nuclear attractors); bond critical points (BCP) characterized by two negative curvatures as ρ is a maximum in the plane of two eigenvectors perpendicular to one another and to the bond path but is a minimum along the third direction parallel to the bond path; ring critical points (RCP) with two positive curvatures given that ρ is a minimum in the ring plane but a maximum along the axis perpendicular to the ring plane; and finally cage critical points (CCP) when the three curvatures are positive and ρ is a local minimum as in the centroid of cubane. The number and types of CPs are interrelated through the equalities: 1 (Isolated molecules) nNCP - nBCP + nRCP - nCCP =  0 (Infinite crystals)

,

(14) where nx is the number of the corresponding xCP, x = “N”, “B”, “R”, or “C”. The first equality, the Poincaré-Hopf relationship (PH),1 applies for isolated molecules while the second, known as the Morse equation, applies in the case of an infinite periodic system.77 The complete set of BCP and RCP electron densities of the intact graphene ribbon are displayed in Fig. 1 along with the electron density contours, the associated gradient field lines, and QTAIM integrated atomic charges. There are 66 NCP, 84 BCPs (which includes the hydrogen−hydrogen BCPs), 19 RCP, and no CCP, a topological set that satisfies Eq. 14. The complete set of BCP electron densities can be cast in matrix format to form what we term an electron-density weighted connectivity matrix (EDWCM) representative of the molecule. An EDWCM is not unique since there exists n! possible labeling schemes for a set of n atoms and hence n! equivalent matrices. This problem is irrelevant to the present study since we reproduce a given matrix representative of a model graphene ribbon after properly re-labelling the atoms of the kernels to correspond to their counterparts in the full intact ribbon. The 7 ACS Paragon Plus Environment

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 8 of 35

EDWCM, like any connectivity matrix, contains entries that are either finite or zero. QTAIM localization (LIs) and a delocalization indices (DIs) count, respectively, the number of electrons that on average are localized within an atomic basin or delocalized from that basin to another. The complete set of n LIs and n(n – 1) DIs can be organized together into a localization-delocalization matrix (LDM, or ζ-matrix), after dividing the DIs each by 2. Thus, an LDM represents the molecule by placing the set of LIs, {Λ(Ωi)}, along the diagonal according to a given atom labeling scheme and then filling the remaining matrix elements with the set of halfDIs, {δ(Ωi,Ωj)/2}:4

∑ row

δ(Ω1 , Ω 2 ) 2 L δ(Ω1 , Ω n ) 2   Λ (Ω1 ) δ(Ω , Ω ) 2 Λ (Ω 2 ) L δ(Ω 2 , Ω n ) 2  2 1  ζ≡   M M O M   Λ (Ω n )  δ(Ωn , Ω1 ) 2 δ(Ω n , Ω 2 ) 2 L

∑ column

= N (Ω1 ) = N (Ω 2 ) = N (Ω n ) 1444444424444444 3

n×n

= N (Ω1 )  = N (Ω 2 )  n  ∑ N (Ω i ) = N M  i =1 . = N (Ω n ) 

(15)

tr (ζ ) = N loc

n

∑ N (Ω ) = N i

i =1

In this manner, the sum of the matrix elements in a given row, say row i, equals the sum of the ith column and both sums equal the total number of electrons in the Ωith atomic basin, N(Ωi), a number that can be obtained by the direct integration of the electron density over the volume of that atom bounded by the external van der Waals (0.001 a.u.) isodensity envelope: N (Ωi ) = Λ (Ωi ) +

1 n ∑ δ(Ωi , Ω j ) = ∫ ρ (r )dr . 2 j ≠i Ωi

(16)

From the populations, the atomic charge are readily obtained from Eq. 6, and the total number of electrons in the molecule N is the sum of all N(Ωi). Note that the atomic charges in

Fig. 1 have been obtained by direct integration of the electron density over each atomic basin, that is, through the second equality of Eq. (16), while those tabulated and discussed otherwise are obtained from the localization and delocalization indices (the first equality), the two sets of calculated atomic populations typically agreeing to the 4th decimal.

4.

Computational details 8 ACS Paragon Plus Environment

Page 9 of 35

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

A finite, electrically neutral, hydrogen-terminated, armchair graphene nanoribbon is the test target molecule to be reconstructed. The ribbon is terminated with hydrogen atoms and has the stoichiometry C46H20 consisting of 14 (2×7) aromatic 6-membered rings and a total of 296 electrons. A ball-and-stick representation of this planar ribbon of C2v symmetry is displayed in the uppermost panels of Fig. 2. The diagram of the ribbon in these figures is demarked with inverted hues (dark) regions that indicates the fission lines at which the molecule is partitioned into the computational fragments as described elsewhere.63,64 The depicted fission lines do not sever any aromatic ring, preserving the number of aromatic sextets78 and hydrogen–hydrogen bonding interactions.70,71 These latter contacts, not graphically displayed for simplicity in the figures despite the presence of the associated bond paths and BCPs, occurring between H54−H55, H56−H57, H58−H59, H62−H63, and H64−H65. KEM has been repeatedly shown not to depend on the underlying model chemistry, this conjecture has been tested under the stringent conditions of strong external electric fields.64 Given the size of the systems considered in this study, the chemical model used in all calculations is the HF/6-31+G(d,p) as implemented in the Gaussian 09 programme.79 The geometry of the intact graphene ribbon is first optimized without constraints then fissioned63 at dark regions of the top panel in Fig. 2 to form the double and single kernels, respectively. The dangling bonds resulting from the fissioning are saturated with hydrogens at the positions of the missing carbon atoms. Three single kernels (K1, K2, and K3) and three double kernels (K12, K13, and K23) were obtained in this manner, the properties of the latter representing the first double summation terms in Eqs. 10-13 while the former deliver the correction terms (the second single summation with a leading minus sign). The AIMAll programme80 has been used to obtain the bond paths, the critical point data, and the atomic integration files of the full intact ribbon (which is used as the standard for comparison) and the six hydrogen-saturated kernels for reconstructing the KEM approximation to the total density. The KEM approximations of the electron density at the bond critical points and at the ring critical points are then reconstructed from Eq. (11) and compared with the directly calculated values. The delocalization indices were obtained from a slightly updated version of AIMDELOC81 since it lists the LIs and DIs together conveniently as a lower triangular matrix. The reconstructed DIs and LIs obtained from Eqs. (12) and (13) are then compared with those

9 ACS Paragon Plus Environment

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 10 of 35

obtained directly for the full intact graphene ribbon.

5.

Results and discussion

5.1

KEM Reconstruction of the Electron Density Weighted Connectivity Matrix (EDWCM).

The molecular EDWCM is the complete set of all the values of the electron density at the bond critical point, {ρ(rBCP)} (Fig. 1), in matrix format. While the full 66×66 matrices of ρ(rBCP), directly calculated and reconstructed from KEM, are available in the Supporting Information along with the difference matrix and the sub-matrices of the six kernels, Table 1 lists the unique non-zero elements of the directly calculated ρ(rBCP) matrix for the intact ribbon and those obtained from KEM via Eq. 11. The ribbon has C2v symmetry and hence consists of two mirror image halves. There are 84 chemical bonds in total, 79 C−C bonds and 5 H−H bonds. The two mirror images are connected by three C−C bonds (C21−C26, C23−C28, and C19−C24) and one H−H bond (H56−H57), hence there exists (84 − 4)/2 + 4 = 44 unique bonds which are listed in

Table 1. (The table should be consulted simultaneously with the numbering scheme in Fig. 2). It is customary to list bond critical point data in a.u. up to a precision of three or four decimals, but the entries in Table 1 are given to the fifth decimals to avoid concealing several small differences between directly calculated values and their KEM approximations. The average absolute deviations are in the fifth decimal place and so are the maximum absolute deviations which are marginally in the fourth decimal place. The non-classical H−H bond density is reproduced with no detectable discrepancy up to the fifth decimal place despite that these weak hydrogen-hydrogen interactions are severed in the fissioning procedure. The electron density at the BCPs of strong inter-kernel C−C bond that are split in the fissioning are slightly less accurately reproduced in two cases, C13−C18 and C17−C22, both exhibiting an error of 0.00017 a.u. There are two artifactual bond critical points, listed at the bottom of Table 1, that arise from the double kernel K13 between pairs of inner capping hydrogen atoms, respectively, 18−29 and 22−27, which are hydrogens that replace carbon atoms with these numbering labels in the intact ribbon (fourth panel in Fig. 2). This artifact is not a limitation of the method since, apart from it, all other bond critical point densities are well reproduced.

10 ACS Paragon Plus Environment

Page 11 of 35

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

5.2

KEM Reconstruction of the Ring Critical Points Vector (RCPV). As with EDWCMs and

LDMs, two molecular fingerprinting descriptors which were defined above, a vector can also be defined and sorted in any desired but consistent way to represent the molecular electron density of polycyclic systems sampled at the ring critical points (RCPs). The ordered set of the RCPs may be termed a “ring critical point vector” or RCPV and may be used, in the future, as an additional inexpensive fingerprinting tool. The electron density (and other scalar field densities) determined at the RCPs of aromatic molecules have been shown to correlate with some measures of aromaticity82-85 which lends an added importance to these critical points in the molecular electron density.

Table 2 lists the symmetry-unique part of the RCPV along with the labels of the atoms forming the given ring. Two types of rings exist in this hydrogen-terminated graphene ribbon: The conventional aromatic rings (arom) and also rings formed by the closure of the bay-region due to the formation of a hydrogen−hydrogen (HH) bridge. The agreement between the KEMreconstructed densities at these points and the exact values is remarkable given the cancellation of contributions from the capping (inner) hydrogens for example (see Fig. 2). The values of the density at the RCP has to be listed to the fifth decimal place, of doubtful or marginal significance, in order that any difference may become discernable. The agreement between the direct and calculated values of ρ(rRCP) data that can be gleaned from Table 2 obviates the need for a statistical analysis. It is noted in passing that the pairs of capping hydrogens pointing inwards head-to-head in K13 are linked by bond paths that crosse the gap between the two parts of this kernel, these are: H22−H27, and H18−H29 (Fig. 2, fourth panel). These long H−H bond path are artifactual and

result

in

the

closure

of

a

10-membered

ring

composed

of

the

atoms

C12−C13−H18−H29−C34−C33−C32−H27−H22−C17. The RCP is characterized by of ρ(rRCP) = 6.38×10−4 a.u., two orders of magnitude lower than any of the “legitimate” RCP densities tabulated in Table 2. This RCP is ignored in KM reconstructions.

5.3

Reconstruction of the Localization-Delocalization Matrix (LDM) and the Complete set

of Atomic Charges. The LDM lists the complete set of the LIs and ½DIs in matrix format. The full 66×66 matrix of the intact ribbon, the KEM reconstructed approximate matrix, their difference matrix, and the sub-matrices of the six kernels are available in the Supporting 11 ACS Paragon Plus Environment

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 12 of 35

Information. Some key QTAIM atomic properties which were directly calculated for the intact ribbon are listed beside their counterparts estimated from the KEM approximation for the 33 unique atoms of the graphene ribbon in Table 3. These properties are: The atomic electron populations (N(Ω)) and their corresponding atomic charges q(Ω), along with the average number of electrons localized within the atomic basin, that is, the LI or Λ(Ω). The results in the table indicate that KEM reproduces these integrated atomic properties with errors typically in the fourth decimal place and occasionally at the third (and generally marginally). The atomic populations listed in Table 3 are obtained for the directly and the KEMreconstructed ribbons from Eq. 5, that is, the sum of the rows or columns of the LDM. These values are identical to at least the third decimal (and typical up to the fourth) to the values displayed in Fig. 1 and which have been obtained by the direct integration of the electron density over the respective atomic basins. The sum of sums defined in Eqs. 7 or 15 yields N, the total molecular population which is given at the bottom of Table 3. The accuracy of the numerical QTAIM basin integration is evident in the very low deviations from the integer number of electrons (296 e−), in the fourth decimal place for the direct calculation and barely the third decimal place for the KEM sum of sums. The calculated average absolute deviation of the atomic population/charge is 0.0012 e−/a.u. for the carbon atoms and 0.0007 e−/a.u. for the hydrogen atoms, values that are probably borderline-insignificant and likely within the usual combined precisions of the electronic structure calculations and the following numerical QTAIM integrations. An even higher accuracy is achieved by KEM in estimating the LIs, which are also listed in Table 3. The average absolute deviation of Λ(Ω) is 0.0007 e− for the carbons and 0.0006 e− for the hydrogen atoms. The molecular Nloc (defined in Eq. 8), that is the total average number of electrons that is localized within their atomic basins, is 188.6698 (direct calculation) and 188.6698 (KEM), differing by a mere − 0.0162 e−, KEM being in error by less than − 0.01%. The LIs and DIs are not independent as can be readily seen from Eqs. 7-9. The localization of ~ 188.7 electrons implies that there exists ~ 107.7 electrons that are delocalized in this system, that is, about 36.3% of the electron population of this ribbon are delocalized between pairs of atomic basins on average at this level of theory. The accurate recovery of the atomic localization indices Λ(Ωi), the total population of localized electrons Nloc, and of the total

12 ACS Paragon Plus Environment

Page 13 of 35

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

molecular electron population N (which includes sums of the delocalization indices), Eq. 7-9, leads one to anticipate the accurate reproduction, by KEM, of the delocalization indices as well, a stringent test for this approximation given the non-local character of the DI.

Table 4 lists a sample of 58 delocalization indices (out of a total of 66 × 65 / 2 = 2145, the total number of distinct atom pairs in the ribbon). This sample includes the DIs with the largest magnitudes (1.585 > δ > 0.956) in the set which are those between atoms sharing a covalent bond path: 31 unique C−C bond paths (28 in each half of the C2v molecule and three crossing the σv mirror plane bisecting the molecule in two halves) and 10 unique C−H bond paths. The table also lists a sampling of significantly lower DIs: Those associated with the closed-shell H−H bond paths (two unique and one crossing the σv plane), and a sampling of some DI’s between non-bonded pairs of carbon atoms and non-bonded pairs of carbon and hydrogen atoms with (0.081 > δ > 0.023).

5.4

Global quality indicators.

Table 5 lists a number of final quantitative comparison criteria to evaluate the overall performance of KEM in recovering the corresponding full matrices from the direct calculations. First, the sum of absolute deviations of all 4356 matrix elements is listed and found to be a mere 0.218 in the case of the LDM and, even much smaller for the EDWCM (0.0127) as can be expected since the latter matrix is much more sparse than the former and has zero entries along the diagonal elements as opposed to the localization indices the numerical values of which can approach the atomic number of the atom in question. Further, the off-diagonal elements of these two matrices, when ρ(rBCP) is expressed in a.u. have comparable magnitudes. The last two remarks concerning the general numerical differences between the LDM and the EDWCM is also reflected on the sum of all the entries in these matrices, 296 electrons delivered by the former and an accumulated density at all bond critical points in the molecule of ~ 48.8 a.u. The % deviations of the KEM-reconstructed matrices from the reference directly calculated counterparts is at most 0.07%. Next, Table 5 lists the Frobenius distance of the KEM-reconstructed and the corresponding directly calculated matrix or vector. This distance is a widely used measure of the distance between matrices. The distance between two n×n matrices A and B is measured by the "Frobenius norm" of the difference matrix between matrices A and B, that is: 13 ACS Paragon Plus Environment

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

d ( A, B) ≡ A − B ≡

∑α

2

ij

− β ij ,

Page 14 of 35

(17)

i, j

where αij and βij are two corresponding matrix elements in A and B. The Frobenius distance of the KEM matrix from the full matrix directly-calculated is 0.014 while the distance in the case of the EDWCM is 0.003. These distances acquire meaning when compared with a reference taken here as the Frobenius norm of the directly calculated matrix, that is, the square-root of the sum of the squares of all matrix elements. The Frobenius norms of the direct LDM, and EDWCM are, respectively, 27.590 and 3.874, which implies corresponding overall discrepancies with regards to the Frobenius distance criterion of only 0.05% and 0.09%.

6.

Concluding remarks

KEM has been amply shown to recover total, interaction, and binding energies of giant molecules (up to a record of 33,000 atoms).61,62 Recently, KEM has been shown to deliver electric-field induced response properties with accuracy. In the present work, the domain of application of KEM is extended to reach the local properties of the electron density and the full set of QTAIM localization delocalization indices. The atomic populations and charge (Eq. 6), normally obtained via numerical integration of the electron density over the volume of the atom (Eq. 16, second equality), is well predicted from the sum of the localization index of the atom and its share in the delocalization index with each remaining atom in the molecule (Eq. 16, first equality). In recovering the full set of localization and delocalization indices, KEM may be useful in studies of the aromaticity of very large aromatic systems, especially those with defects, irregularities, or very large unit cells precluding periodic calculations, given the known relation between the DIs and aromaticity.86-95 The ability of KEM to closely reproduce the directly calculated electron density has only been verified at the critical points (both BCPs and RCPs), and previously in the ability of KEM to reproduce the field-induced dipole moment.64 We are currently investigating the reconstruction of the entire electron density scalar field from the densities of kernel fragment in a manner equivalent to Eq. 11 but applied at every point in three-dimensional space. The KEM approximate density can then be used to generate all one-electron properties that are represented 14 ACS Paragon Plus Environment

Page 15 of 35

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

by multiplicative operators. Further, the molecular electrostatic potential (MESP) can readily be obtained from the KEM density and the known charges Zi (in a.u.) and positions (Ri) of the nuclei:96 VKEM (r ) =

M



i =1 R i ≠r

Zi ρ (r') − ∫ KEM dr ' , Ri − r r' − r

(18)

where the ith term in the first sum must be eliminated when Ri = r. Alternatively, and instead of using Eq. 18, the VKEM(r) can be reconstructed from the MESP calculated for the kernel fragments and manipulated, in its turn, on a point-by-point basis using an equation of the same form of Eq. 11. A point which was not investigated in this paper is the ability of KEM to reproduce the position vectors (rCP) of the bond, ring, and cage critical points in space. While the nuclear critical points are exactly those of the intact molecule since KEM calculations are based on the frozen geometry of the intact system, the position of a given critical point is expected to be slightly different from kernel to kernel and from the reconstructed KEM total density gauged against the position of this critical point obtained from the direct calculation on the intact molecule. The accuracies of KEM’s rCP’s are anticipated to be high given the high accuracy by which it delivers the critical point densities demonstrated in this work. To verify this point requires the reconstruction of the full density scalar field (using an equation of the same form of Eq. 11 applied at every point of space) then determining the positions where the gradient of the KEM density vanish. We conclude by presenting two thoughts related to an understanding of the accuracy of the KEM, which is evident from the numerical results of this paper. First, KEM represents a numerical method which displays what Bader and his school have called “transferability” of the electron density.97-103 The assertion is that the electron density of a group of atoms is transferable from one molecule to another, and that quantum properties which are functionals of that density are also transferable.103 The KEM may be conceived of as a procedure in which groups of atoms are extracted from a full molecule transferring thereby their density and their quantum mechanical properties into single and double kernel fragments. Insofar as the kernels by transferability replicate the density of the fragments they mimic within the full molecule, they also mimic the quantum mechanical properties of those fragments. Thus it is that an appropriate summation over the kernel quantum properties should indeed replicate the corresponding 15 ACS Paragon Plus Environment

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 16 of 35

quantum properties of the full molecule. Second, numerical experience suggests that the usual quantum mechanical properties, obtained from the usual chemical models, can be replicated by the KEM with that accuracy appropriate to a given chemical model. Consider for example the total molecular energy which is obtained with accuracy by KEM.57 This suggests that KEM also determines with accuracy the density, the one body density matrix, and the diagonal elements of the two body density matrix, since these determine respectively, the external potential energy, the kinetic energy, and the electronic repulsion energy. But given the density and these density matrices, all of the usual quantum mechanical properties follow, applying the quantum mechanical rules for calculation of expectation values. In closing, this work may be seen as a resulting from a synergy of three branches of theoretical chemistry: Bader’s quantum theory of atoms in molecules (QTAIM), chemical graph theory (CGT), and the Kernel Energy method (KEM). While we consider this paper to conclusively show that the QTAIM results discussed here are inherently compatible with KEM, we understand this is only one molecular example, which provides a proof of principle argument. There are obviously a wide variety of molecules of differing shapes and sizes, and differing types of chemical bonds against which to test the accuracy and QTAIM-CGT/KEM compatibility that may be achievable. It has not been within the scope of this paper to undertake the wider task of calculating the QTAIM properties discussed here for the wide variety of molecules to which KEM has been applied. However, it is encouraging that KEM itself has proven to be highly accurate for molecules large and small of a very wide variety of geometrical shapes across the whole spectrum of normally used chemical models as documented in the literature cited above. That is at least consistent with the notion that the KEM density related to such a wide variety of results is accurately calculated.

Supporting information A text file listing all matrices in the following order: EDWCM of the full ribbon from direct calculation; EDWCM of the full ribbon reconstructed from KEM calculation; ∆EDWCM (direct calculation matrix element minus the corresponding KEM matrix element - notice that this is the reverse convention than the text where the difference is taken as (KEM – direct)); then the EDWCMs of K1, K2, K3, K12, K13, and K23. These are followed by the LDMs in the same 16 ACS Paragon Plus Environment

Page 17 of 35

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

corresponding order as for the EDWCM (the difference is opposite in its convention than that in the text as well).

Acknowledgments The authors thank the two anonymous reviewers of this paper for their helpful corrections and suggestions. Funding for this project was provided by the Natural Sciences and Engineering Research Council of Canada (NSERC), Canada Foundation for Innovation (CFI), and Mount Saint Vincent University - (C. F. M.), the U.S. Naval Research Laboratory (project # 4720300 01) and by a PSC CUNY Award (project # 63842-00 41) - (L. M.); and the Office of Naval Research (ONR) through the Naval Research Laboratory's Basic Research Program - (L. H.)

17 ACS Paragon Plus Environment

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

References (1) Bader, R. F. W. Atoms in Molecules: A Quantum Theory; Oxford University Press: Oxford, U.K., 1990. (2) Popelier, P. L. A. Atoms in Molecules: An Introduction; Prentice Hall: London, 2000. (3) Matta C. F.; Boyd, R. J. (Eds.) The Quantum Theory of Atoms in Molecules: From Solid State to DNA and Drug Design; Wiley-VCH: Weinheim, 2007. (4) Matta, C. F. Modeling biophysical and biological properties from the characteristics of the molecular electron density, electron localization and delocalization matrices, and the electrostatic potential. J. Comput. Chem. 2014, 35, 1165-1198. (5) Sumar, I.; Ayers, P. W.; Matta, C. F. Electron localization and delocalization matrices in the prediction of pKa's and UV-wavelengths of maximum absorbance of p-benzoic acids and the definition of super-atoms in molecules. Chem. Phys. Lett. 2014, 612, 190-197. (6) Dittrich, B.; Matta, C. F. Contributions of charge-density research to medicinal chemistry. Int. U. Cryst. J. (IUCrJ) 2014, accepted, in press. (7) Matta, C. F. Localization-delocalization matrices and electron density-weighted adjacency matrices: New electronic fingerprinting tools for medicinal computational chemistry. Future Med. Chem. 2014, 6, accepted, in press. (8) Hall, L. H.; Kier, L. B. Molecular Connectivity in Chemistry and Drug Research; Academic Press: Boston, 1976. (9) Balaban, A. T. Chemical Applications of Graph Theory; Academic Press: New York, 1976. (10) Diudea, M. D.; Gutman, I.; Lorentz, J. Molecular Topology; Nova Science Publishers, Inc.: Hauppauge NY, 1999. (11) Janezic, D.; Milicevic, A.; Nikolic, S.; Trinajstic, N. Graph Theoretical Matrices in Chemistry (Mathematical Chemistry Monographs, Vol. 3); University of Kragujevac: Kragujevac, 2007. (12) Bonchev, D.; Rouvray, D. H. Complexity in Chemistry: Introductions and Fundamentals; Taylor and Francis: London, 2003. (13) Bonchev, D.; Rouvray, D. H. Chemical Graph Theory: Introduction and Fundamentals; OPA: Amsterdam, 1991. (14) Dmitriev, I. S. Molecules without Chemical Bonds (English Translation); Mir Publishers: Moscow, 1981. (15) Todeschini, R.; Consonni, V. Molecular Descriptors for Chemoinformatics (Second Edition; Vols. I and II); Wiley-VCH Weinheim: Weinheim, 2009. (16) Fradera, X.; Austen, M. A.; Bader, R. F. W. The Lewis model and beyond. J. Phys. Chem. A 1999, 103, 304-314. (17) Poater, J.; Solà, M.; Duran, M.; Fradera, X. The calculation of electron localization and delocalization indices at the Hartree-Fock, density functional and post-Hartree-Fock levels of theory. Theor. Chem. Acc. 2002, 107, 362-371. (18) Fradera, X.; Poater, J.; Simon, S.; Duran, M.; Solà, M. Electron-pairing analysis from localization and delocalization indices in the framework of the atoms-in-molecules theory. Theor. Chem. Acc. 2002, 108, 214-224. (19) Matito, E.; Duran, M.; Solà, M. A novel exploration of the Hartree-Fock homolytic bond dissociation problem in the hydrogen molecule by means of electron localization measures. J.

18 ACS Paragon Plus Environment

Page 18 of 35

Page 19 of 35

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Chem. Edu. 2006, 83, 1243-1248. (20) Matta, C. F.; Hernández-Trujillo, J.; Bader, R. F. W. Proton spin-spin coupling and electron delocalisation. J. Phys. Chem. A 2002, 106, 7369-7375. (21) Wang, Y.-G.; Matta, C. F.; Werstiuk, N. H. Comparison of localization and delocalization indices obtained with Hartree-Fock and conventional correlated methods: Effect of Coulomb correlation. J. Comput. Chem. 2003, 24, 1720-1729. (22) Blanco M. A.; Martin Pendas, A.; Francisco, E. Interacting quantum atoms: A correlated energy decomposition scheme based on the quantum theory of atoms in molecules. J. Chem. Theor. Comput. 2005, 1, 1096-1109. (23) Martin Pendas, A.; Francisco, E.; Blanco, M. A.; Gatti, C. Bond paths as privileged exchange channels. Chem. Eur. J. 2007, 13, 9362-9371. (24) Francisco, E.; Martin Pendas, A.; Blanco, M. A. A molecular energy decomposition scheme for atoms in molecules. J. Chem. Theory Comput. 2006, 2, 90-102. (25) Martin Pendas, A.; Francisco, E.; Blanco, M. A. Binding energies of first row diatomics in the light of the interacting quantum atoms approach. J. Phys. Chem. A 2006, 110, 12864-12869. (26) Martin Pendas, A.; Blanco, M. A.; Francisco, E. The nature of the hydrogen bond: A synthesis from the interacting quantum atoms picture. J. Chem. Phys. 2006, 125, 184112. (27) Bader, R. F. W. Bond paths are not chemical bonds. J. Phys. Chem. A 2009, 113, 10391-10396. (28) Bader, R. F. W. A bond path: A universal indicator of bonded interactions. J. Phys. Chem. A 1998, 102, 7314-7323. (29) Runtz, G. R.; Bader, R. F. W.; Messer, R. R. Definition of bond paths and bond directions in terms of the molecular charge distribution. Can. J. Chem. 1977, 55, 3040-3045. (30) Yang, W.; Lee, T.-S. A density-matrix divide-and-conquer approach for the electronic structure calculations of large molecules. J. Chem. Phys. 1995, 103, 5674-5678. (31) Lee, T.-S.; Lewis, J. P.; Yang, W. Linear-scaling quantum mechanical calculations of biological molecules: The divide-and-conquer approach. Comput. Mater. Sci. Volume 12, 1998, Pages 259277 1998, 12, 259-277. (32) Kitaura, K.; Ikeo, E.; Asada, T.; Nakano T.; Uebayasi, M. Fragment molecular orbital method: an approximate computational method for large molecules. Chem. Phys. Lett. 1999, 313, 701-706. (33) Fedorov, D.; Kitaura, K. The Fragment Molecular Orbital Method: Practical Applications to Large Molecular Systems; CRC Press: Boca Raton, Florida, 2009. (34) Zhang, D. W.; Zhang, J. Z. H. Molecular fractionation with conjugate caps for full quantum mechanical calculation of protein-molecule interaction energy. J. Chem. Phys. 2003, 119, 35993605. (35) Chen, X. H.; Zhang, D. W.; Zhang, J. Z. H. Fractionation of peptide with disulfide bond for quantum mechanical calculation of interaction energy with molecules. J. Chem. Phys. 2004, 120, 839-844. (36) He, X.; Zhang, J. Z. H. A new method for direct calculation of total energy of protein. J. Chem. Phys. 2005, 122, 031103-1-031103-4. (37) Li, S.; Li, W.; Fang, T. An efficient fragment-based approach for predicting the ground-state energies and structures of large molecules. J. Am. Chem. Soc. 2005, 127, 7215-7226. (38) Li, W.; Fang, T.; Li, S. A fragment energy assembler method for Hartree-Fock calculations of large molecules. J. Chem. Phys. 2006, 124, 154102-1-15402-6.

19 ACS Paragon Plus Environment

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(39) Jiang, N.; Ma, J.; Jiang, Y. Electrostatic field-adapted molecular fractionation with conjugated caps for energy calculations of charged biomolecules. J. Chem. Phys. 2006, 124, 114112-1114112-9. (40) Hua, S.; Hua, W.; Li, S. An efficient implementation of the generalized energy-based fragmentation approach for general large molecules. J. Phys. Chem. A 2010, 114, 8126-8134. (41) Deev, V.; Collins, M. A. Approximate ab initio energies by systematic molecular fragmentation. J. Chem. Phys. 2005, 122, 154102-1-154102-12. (42) Collins, M. A.; Deev, V. A. Accuracy and efficiency of electronic energies from systematic molecular fragmentation. J. Chem. Phys. 2006, 125, 104104-1-104104-15. (43) Netzloff, H. M.; Collins, M. A. Ab initio energies of non-conducting crystals by systematic fragmentation. J. Chem. Phys. 2007, 127, 134113-1-134113-13. (44) Bettens, R. P. A.; Lee, A. M. A new algorithm for molecular fragmentation in quantum chemical calculations. J. Phys. Chem. A 2006, 110, 8777-8785. (45) Lee, A. M.; Bettens, R. P. A. First principles NMR calculations by fragmentation. J. Phys. Chem. A 2007, 111, 5111-5115. (46) Le H.-A.; Lee, A. M.; Bettens, R. P. A. Accurately reproducing ab initio electrostatic potentials with multipoles and fragmentation. J. Phys. Chem. A 2009, 113, 10527-10533. (47) Ganesh, V.; Dongare, R. K.; Balanarayan, P.; Gadre, S. R. Molecular tailoring approach for geometry optimization of large molecules: Energy evaluation and parallelization strategies. J. Chem. Phys. 2006, 125, 104109. (48) Hung, L.; Carter, E. A. Accurate simulations of metals at the mesoscale: Explicit treatment of 1 million atoms with quantum mechanics. Chem. Phys. Lett. 2009, 475, 163-170. (49) Mueller, C.; Paulus, B. Wavefunction-based electron correlation methods for solids. Phys. Chem. Chem. Phys. (PCCP) 2012, 14, 7605-7614. (50) Voloshina, E.; Paulus, B. Development of a wavefunction-based ab Initio method for metals applying the method of increments. Z. Physik. Chem. 2010, 224, 369-381. (51) Paulus, B.; Rosciszewski, K. Application of the method of increments to the adsorption of H2S on graphene. Int. J. Quantum Chem. 2009, 109, 3055-3062. (52) Paulus, B. The method of increments - A wavefunction-based ab-initio correlation method for solids. Int. J. Mod. Phys. B 2007, 21, 2204-2214. (53) Paulus, B. Towards an ab initio incremental correlation treatment for metals. Chem. Phys. Lett. 2003, 371, 7-14. (54) Paulus, B.; Rosciszewski, K.; Stoll, H.; Birkenheuer, U. Ab initio incremental correlation treatment with non-orthogonal localized orbitals. Phys. Chem. Chem. Phys. (PCCP) 2003, 5, 5523-5529. (55) Matta, C. F. Theoretical reconstruction of the electron density of large molecules from fragments determined as proper open quantum systems: the properties of the oripavine PEO, enkephalins, and morphine. J. Phys. Chem. A 2001, 105, 11088-11101. (56) Bader R. F. W.; Matta, C. F.; Martín, F. J. In: Medicinal Quantum Chemistry; Alber, F.; Carloni P. (Eds.), Wiley-VCH: Weinheim, 2003. (57) Huang, L.; Massa, L.; Karle, J. Kernel energy method: Application to insulin. Proc. Natl. Acad. Sci. USA 2005, 102, 12690-12693. (58) Huang, L.; Massa, L.; Karle, J. Kernel energy method illustrated with peptides. Int. J. Quantum Chem. 2005, 103, 808-817.

20 ACS Paragon Plus Environment

Page 20 of 35

Page 21 of 35

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

(59) Huang, L.; Massa, L.; Karle, J. Kernel energy method: Application to DNA. Biochemistry 2005, 44, 16747-16752. (60) Huang, L.; Massa, L.; Karle, J. The Kernel Energy Method: Application to a tRNA. Proc. Natl. Acad. Sci. USA 2006, 103, 1233-1237. (61) Huang, L.; Massa, L.; Karle, J. Kernel energy method applied to vesicular stomatitis virus nucleoprotein. Proc. Natl. Acad. Sci. USA 2009, 106, 1731-1736. (62) Huang L.; Massa, L.; Karle, J. In: Quantum Biochemistry: Electronic Structure and Biological Activity; Matta, C. F. (Ed.), Wiley-VCH: Weinheim, 2010. (63) Huang, L.; Bohorquez, H.; Matta, C. F.; Massa, L. The Kernel Energy Method: Application to Graphene and Extended Aromatics. Int. J. Quantum Chem. 2011, 111, 4150-4157. (64) Huang, L.; Massa, L.; Matta, C. F. A graphene flake under external electric fields reconstructed from field-perturbed kernels. Carbon 2014, 76, 310-320. (65) Balasubramanian, K. Integration of graph theory and quantum chemistry for structure-activity relationship. SAR & QSAR Eviron. Res. 1994, 2, 59-77. (66) Jenkins, S. Quantum topology phase diagrams for molecules, clusters, and solids. Int. J. Quantum Chem. 2013, 113, 1603-1608. (67) Xu, T.; Jenkins, S.; Xiao, C.-X.; Maza, J. R.; Kirk, S. R. The Pt site reactivity of the molecular graphs of Au6Pt isomers. Chem. Phys. Lett. 2013, 590, 41-45. (68) Jenkins, S.; Restrepo, A.; David, J.; Dulin Yina; Kirk, S. R. Spanning QTAIM topology phase diagrams of water isomers W4, W5 and W6. Phys. Chem. Chem. Phys. 2011, 13, 11644-11656. (69) Jenkins, S.; Rong, C.; Kirk, S. R.; Yin, D.; Liu, S. Spanning set of silica cluster isomer topologies from QTAIM. J. Phys. Chem. A 2011, 115, 12503-12511. (70) Matta, C. F.; Hernández-Trujillo, J.; Tang, T. H.; Bader, R. F. W. Hydrogen-hydrogen bonding: a stabilizing interaction in molecules and crystals. Chem. Eur. J. 2003, 9, 1940-1951. (71) Hernández-Trujillo, J.; Matta, C. F. Hydrogen-hydrogen bonding in biphenyl revisited. Struct. Chem. 2007, 18, 849-857. (72) Cao, W. L.; Gatti, C.; MacDougall, P. J.; Bader, R. F. W. On the presence of non-nuclear attractors in the charge distributions of Li and Na clusters. Chem. Phys. Lett. 1987, 141, 380-385. (73) de Vries, R. Y.; Briels, W. J.; Feil, D.; te Velde, G.; Baerends, E. J. Charge density study with maximum entropy method on model data of silicon. A search for non-nuclear attractors. Can. J. Chem. 1996, 74, 1054-1058. (74) Martin-Pendas, A.; Blanco, M. A.; Costales, A.; Mori-Sanchez, P.; Luana, V. Non-nuclear maxima of the electron density. Phys. Rev. Lett. 1999, 83, 1930-1933. (75) Bader, R. F. W.; Platts, J. A. Characterization of an F-center in an alkali halide cluster. J. Chem. Phys. 1997, 107, 8545-8553. (76) Taylor, A.; Matta, C. F.; Boyd, R. J. The hydrated electron as a pseudo-atom in cavity-bound water clusters. J. Chem. Theor. Comput. 2007, 3, 1054-1063. (77) Coppens, P. X-ray Charge Densities and Chemical Bonding; Oxford University Press, Inc.: New York, 1997. (78) Clar, E. The Aromatic Sextet; John Wiley and Sons Ltd.: London, 1972. (79)

Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Mennucci, B.; Petersson, G. A.; Nakatsuji, H.; Caricato, M.; Li, X.; Hratchian, H. P.; Izmaylov, A. F.; Bloino, J.; Zheng, G.; Sonnenberg, J. L.; Hada, M.; Ehara, M.;

21 ACS Paragon Plus Environment

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Vreven, T.; Montgomery Jr, J. A.; Peralta, J. E.; Ogliaro, F.; Bearpark, M.; Heyd, J. J.; Brothers, E.; Kudin, K. N.; Staroverov, V. N.; Keith, T.; Kobayashi, R.; Normand, J.; Raghavachari, K.; Rendell, A.; Burant, J. C.; Iyengar, S. S.; Tomasi, J.; Cossi, M.; Rega, N.; Millam, J. M.; Klene, M.; Knox, J. E.; Cross, J. B.; Bakken, V.; Adamo, C.; Jaramillo, J.; Gomperts, R.; Stratmann, R. E.; Yazyev, O.; Austin, A. J.; Cammi, R.; Pomelli, C.; Ochterski, J. W.; Martin, R. L.; Morokuma, K.; Zakrzewski, V. G.; Voth, G. A.; Salvador, P.; Dannenberg, J. J.; Dapprich, S.; Daniels, A. D.; Farkas, O.; Foresman, J. B.; Ortiz, J. V.; Cioslowski, J.; Fox, D. J., Gaussian 09, Revision B.01 (Gaussian Inc., Wallingford CT, 2010). (80) Keith, T. A. AIMAll. http://aim.tkgristmill.com/ 2014. (81)

Matta, C. F., AIMDELOC: Program to calculate AIM localization and delocalization indices (QCPE0802), Quantum Chemistry Program Exchange, Indiana University, IN, 2001.

(82) Howard, S. T.; Krygowski, T. M. Benzenoid hydrocarbon aromaticity in terms of charge density descriptors. Can. J. Chem. 1997, 75, 1174-1181. (83) Palusiak, M. and Krygowski, T. M. Application of AIM parameters at ring critical points for estimation of π-electron delocalization in six-membered aromatic and quasi-aromatic rings. Chem. Eur. J. 2007, 13, 7996-8006. (84) Ebrahimi, A. A., Ghiasi, R., and Foroutan-Nejad, C. Topological characteristics of the ring critical points and the aromaticity of groups IIIA to VIA hetero-benzenes. J. Mol. Struct. (THEOCHEM) 2010, 941, 47-52. (85) Firme, C. I., Galembeck, S. E., Antunes, O. A. C., and Esteves, P. M. Density, degeneracy, delocalization-based index of aromaticity (D3BIA). J. Braz. Chem. Soc. 2007, 18, 1397-1404. (86) Schleyer P.v-R. (Guest Ed.). Aromaticity (Special Issue). Chem. Rev. 2001, 101. (87) Mandado, M.; Gonzalez Moa, M. J.; Mosquera, R. A. Aromaticity: Exploring Basic Chemical Concepts with the Quantum Theory of Atoms in Molecules; Nova Science Publishers, Inc.: New York, 2008. (88) Wannere, C. S.; Corminboeuf, C.; Wang, Z. X.; Wodrich, M. D.; King, R. B.; Schleyer, P. v. R. Evidence for d orbital aromaticity in square planar coinage metal clusters. J. Am. Chem. Soc. 2005, 127, 5701-5705. (89) Poater, J.; Fradera, X.; Duran, M.; Solà, M. The delocalization index as an electronic aromaticity criterion: application to a series of planar polycyclic aromatic hydrocarbons. Chem. Eur. J. 2003, 9, 400-406. (90) Poater, J.; Duran, M.; Solà, M.; Silvi, B. Theoretical evaluation of electron delocalization in aromatic molecules by means of atoms in molecules (AIM) and electron localization function (ELF) topological approaches. Chem. Rev. 2005, 105, 3911-3947. (91) Matito, E.; Duran, M.; Solà, M. The aromatic fluctuation index (FLU): A new aromaticity index based on electron delocalization. J. Chem. Phys. 2005, 122, 014109. (92) Matito E.; Poater, J.; Solà, M. In: The Quantum Theory of Atoms in Molecules: From Solid State to DNA and Drug Design; Matta, C. F.; Boyd, R. J. (Eds.),Wiley-VCH: Weinheim, 2007. (93) Cyranski, M. K. Energetic aspects of cyclic σ-electron delocalization: Evaluation of the methods of estimating aromatic stabilization energies. Chem. Rev. 2005, 105, 3773-3811. (94) Matta, C. F.; Hernández-Trujillo, J. Bonding in polycyclic aromatic hydrocarbons in terms of the the electron density and of electron delocalization. J. Phys. Chem. A 2003, 107, 7496-7504 (Correction: J. Phys. Chem A, 2005, 109, 10798). (95) Wolstenholme, D.; Matta, C. F.; Cameron, T. S. Experimental and theoretical electron density

22 ACS Paragon Plus Environment

Page 22 of 35

Page 23 of 35

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

study of a highly twisted polycyclic aromatic hydrocarbon: 4 Methyl-[4]helicene. J. Phys. Chem. A 2007, 111, 8803-8813. (96) Bonaccorsi, R.; Scrocco, E.; Tomasi, J. Molecular SCF calculations for the ground state of some three-membered ring molecules: (CH2)3, (CH2)2NH, (CH2)2NH2+, (CH2)2O, (CH2)2S, (CH)2CH2, and N2CH2. J. Chem. Phys. 1970, 52, 5270-5284. (97) Bader R. F. W.; Cortés-Guzmán, F. In: Quantum Biochemistry: Electronic Structure and Biological Activity (Volume 1); Matta, C. F. (Ed.), Wiley-VCH: Weinheim, 2010. (98) Rykounov, A. A.; Tsirelson, V. G. Quantitative estimates of transferability of the QTAIM descriptors. Case study of the substituted hydropyrimidines. J. Mol. Struct. (Theochem) 2009, 906, 11-24. (99) Turovtsev, V. V.; Orlov, Y. D.; Lebedev, Y. A. The inductive effect of a radical center and transferability of the properties of the functional groups in n-alkyl radicals. Russ. J. Phys. Chem. A 2009, 83, 245-253. (100) Luger P.; Dittrich, B. In: The Quantum Theory of Atoms in Molecules: From Solid State to DNA and Drug Design; Matta, C. F.; Boyd, R. J. (Eds.),Wiley-VCH: Weinheim, 2007. (101) Cortés-Guzmán, F.; Bader, R. F. W. Transferability of group energies and satisfaction of the virial theorem. Chem. Phys. Lett. 2003, 379, 183-192. (102) Bader, R. F. W.; Keith, T. A.; Gough, K. M.; Laidig, K. E. Properties of atoms in molecules: additivity and transferability of group polarizabilities. Mol. Phys. 1992, 75, 1167-1189. (103) Bader, R. F. W.; Becker, P. Transferability of atomic properties and the theorem of Hohenberg and Kohn. Chem. Phys. Lett. 1988, 148, 452-458.

23 ACS Paragon Plus Environment

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 24 of 35

Table 1 The values of the electron density at the 44 unique bond critical points of the intact graphene ribbon obtained from the topological analysis of the electron density of the intact ribbon directly (ρdirect) and estimated from the KEM approximation (ρKEM), in atomic units (a.u.) Bond(a) C−C bonds C1−C2 C1−C5 C2−C3 C2−C6 C3−C4 C4−C8 C5−C10 C6−C7 C6−C11 C7−C8 C7−C12 C8−C9 C9−C14 C10−C11 C11−C16 C12−C13 C12−C17 C13−C14 C13−C18 C15−C16 C15−C20 C16−C17 C17−C22 C18−C19 C18−C23 C19−C24 C20−C21 C21−C22 C21−C26 C22−C23 C23−C28 C−H bonds C1−H52 C3−H47 C4−H48 C5−H53 C9−H61 C10−H54 C14−H62

ρdirect

ρKEM

0.31843 0.32623 0.30716 0.30447 0.33761 0.30801 0.32186 0.30958 0.30844 0.30858 0.30569 0.31608 0.33195 0.31765 0.29322 0.31123 0.30521 0.31142 0.29632 0.31596 0.33166 0.30939 0.30312 0.31545 0.30943 0.33201 0.31503 0.31069 0.29662 0.30489 0.30275

0.31845 0.32626 0.30717 0.30448 0.33762 0.30803 0.32188 0.30958 0.30843 0.30859 0.30563 0.31611 0.33193 0.31765 0.29319 0.31115 0.30518 0.31137 0.29615 0.31598 0.33162 0.30935 0.30295 0.31548 0.30939 0.33201 0.31502 0.31066 0.29662 0.30492 0.30273

0.29251 0.29262 0.29260 0.29401 0.29255 0.29547 0.29598

0.29260 0.29269 0.29263 0.29404 0.29260 0.29554 0.29606

24 ACS Paragon Plus Environment

Page 25 of 35

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Bond(a) C15−H55 C19−H63 C20−H56 H−H bonds H54−H55 H62−H63 H56−H57 Artifacts C/H18−C/H29 C/H22−C/H27

ρdirect ρKEM 0.29577 0.29576 0.29610 0.29610 0.29605 0.29602 0.01369 0.01369 0.01455 0.01455 0.01455 0.01455 0 0.00173 0 0.00165

(a) See Fig. 2 for atom numbering. The average absolute deviation ( ρ KEM − ρ direct ) , and the maximum absolute deviation (max( ρ KEM − ρ direct )) in a.u. are, respectively: (0.00003, 0.00017 for C−C bonds), and (0.00005, 0.00009 for C−H bonds).

25 ACS Paragon Plus Environment

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 26 of 35

Table 2 The values of the electron density at the 8 unique aromatic (arom) ring critical points (RCP), ρ(rRCP), and 3 unique rings closed by a hydrogen−hydrogen (HH) bond path obtained directly(ρdirect) and estimated from the KEM approximation (ρKEM), in atomic units (a.u.) Type arom arom arom arom arom arom arom arom HH HH HH

Ring members(a) C1−C2−C6−C11−C10−C5 C2−C3−C4−C8−C7−C6 C6−C7−C12−C17−C16−C11 C7−C8−C9−C14−C13−C12 C12−C13−C18−C23−C22−C17 C15−C16−C17−C22−C21−C20 C18−C19−C24−C29−C28−C23 C21−C22−C23−C28−C27−C26 C10−C11−C16−C15−H55−H54 C13−C14−H62−H63−C19−C18 C20−C21−C26−C25−H57−H56

ρdirect 0.01939 0.01891 0.01806 0.01913 0.01802 0.01919 0.01917 0.01803 0.01175 0.01222 0.01224

(a) See Fig. 2 for atom numbering.

26 ACS Paragon Plus Environment

ρKEM 0.01938 0.01891 0.01806 0.01913 0.01802 0.01920 0.01917 0.01805 0.01174 0.01222 0.01224

Page 27 of 35

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Table 3 Average atomic electron populations N(Ω), the corresponding atomic charges q(Ω) (in atomic units (a.u.)), and the average number of electrons localized within the atomic basin (the atomic localization index, LI, or Λ(Ω)) directly calculated and estimated from the KEM approximation for the 33 unique atoms of the graphene ribbon. Ω(a) N(Ω)direct N(Ω)KEM(b) q(Ω)direct q(Ω)KEM Λ(Ω)direct Λ(Ω)KEM(b) C1 5.9714 5.9711 0.0286 0.0289 3.9130 3.9127 C2 5.9984 5.9991 0.0016 0.0009 3.8877 3.8877 C3 5.9680 5.9674 0.0320 0.0326 3.9112 3.9107 C4 5.9677 5.9668 0.0323 0.0332 3.9106 3.9101 C5 5.9576 5.9571 0.0424 0.0429 3.8999 3.8996 C6 6.0117 6.0114 -0.0117 -0.0114 3.9023 3.9020 C7 6.0096 6.0098 -0.0096 -0.0098 3.8991 3.8994 C8 5.9979 5.9976 0.0021 0.0024 3.8871 3.8869 C9 5.9653 5.9640 0.0347 0.0360 3.9059 3.9050 C10 5.9664 5.9654 0.0336 0.0346 3.9062 3.9055 C11 6.0022 6.0030 -0.0022 -0.0030 3.8929 3.8939 C12 6.0076 6.0068 -0.0076 -0.0068 3.8969 3.8962 C13 6.0021 6.0004 -0.0021 -0.0004 3.8924 3.8902 C14 5.9619 5.9627 0.0381 0.0373 3.9013 3.9019 C15 5.9572 5.9567 0.0428 0.0433 3.8958 3.8952 C16 6.0029 6.0026 -0.0029 -0.0026 3.8940 3.8936 C17 6.0077 6.0059 -0.0077 -0.0059 3.8977 3.8955 C18 6.0016 6.0086 -0.0016 -0.0086 3.8916 3.8938 C19 5.9592 5.9596 0.0408 0.0404 3.8975 3.8977 C20 5.9591 5.9594 0.0409 0.0406 3.8977 3.8979 C21 6.0014 6.0012 -0.0014 -0.0012 3.8914 3.8913 C22 6.0067 6.0133 -0.0067 -0.0133 3.8969 3.8991 C23 6.0071 6.0071 -0.0071 -0.0071 3.8975 3.8974 H47 1.0258 1.0251 -0.0258 -0.0251 0.4675 0.4669 H48 1.0264 1.0255 -0.0264 -0.0255 0.4681 0.4675 H52 1.0294 1.0286 -0.0294 -0.0286 0.4701 0.4694 H53 1.0271 1.0262 -0.0271 -0.0262 0.4714 0.4707 H54 1.0357 1.0349 -0.0357 -0.0349 0.4672 0.4667 H55 1.0333 1.0337 -0.0333 -0.0337 0.4655 0.4658 H56 1.0345 1.0348 -0.0345 -0.0348 0.4658 0.4659 H61 1.0292 1.0280 -0.0292 -0.0280 0.4704 0.4695 H62 1.0338 1.0325 -0.0338 -0.0325 0.4650 0.4639 H63 1.0344 1.0346 -0.0344 -0.0346 0.4655 0.4656 (c) Σ 296.0003 296.0022 -0.0003 -0.0022 188.6860 188.6698 (a) (b)

See Fig. 2 for atom numbering. The average absolute deviation

(c)

(max( P (Ω ) KEM − P (Ω)direct )) are, respectively: (0.0010, 0.0070 for P(Ω) = N(Ω)), and (0.0007, 0.0022 for P(Ω) = Λ(Ω)). The sums include the other unlisted mirror image molecular half.

( P (Ω ) KEM − P (Ω ) direct ) , and the maximum absolute deviation

27 ACS Paragon Plus Environment

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Table 4 A sampling of delocalization indices δ(Ωi,Ωj) directly calculated and estimated from the KEM approximation between one of the two sets of unique pairs of strongly (covalently) bonded atoms, weakly bonded pairs of hydrogen atoms, and a few non-bonded delocalization indices in the graphene ribbon. (The full set of 66×66-66= 4290 DIs/2 and the 66 LIs for the direct, KEM, and kernel systems are all available in the Supporting Information). (Note that the listed DIs in this table are not divided by 2 as in the matrix tabulation in the Supporting Information). Ωi−Ωj(a) C−C Bonds C1−C2 C1−C5 C2−C3 C2−C6 C3−C4 C4−C8 C5−C10 C6−C7 C6−C11 C7−C8 C7−C12 C8−C9 C9−C14 C10−C11 C11−C16 C12−C13 C12−C17 C13−C14 C13−C18 C15−C16 C15−C20 C16−C17 C17−C22 C18−C19 C18−C23 C19−C24 C20−C21 C21−C22 C21−C26 C22−C23 C23−C28 C−H Bonds C1−H52 C3−H47 C4−H48 C5−H53 C9−H61 C10−H54 C14−H62 C15−H55 C19−H63 C20−H56 H−H Bonds H54−H55

δ(Ωi,Ωj)direct

δ(Ωi,Ωj)KEM(b)

1.3095 1.4123 1.1718 1.2431 1.5853 1.1765 1.3626 1.1912 1.2478 1.2704 1.2208 1.2678 1.4716 1.3584 1.1039 1.3053 1.1582 1.2708 1.1185 1.2972 1.4444 1.2937 1.2319 1.2886 1.2892 1.4507 1.2933 1.2887 1.1153 1.1626 1.2304

1.3100 1.4124 1.1722 1.2431 1.5852 1.1762 1.3624 1.1911 1.2480 1.2706 1.2204 1.2674 1.4719 1.3586 1.1035 1.3070 1.1569 1.2704 1.1143 1.2962 1.4449 1.2948 1.2290 1.2878 1.2909 1.4512 1.2923 1.2901 1.1146 1.1611 1.2303

0.9702 0.9685 0.9686 0.9746 0.9686 0.9589 0.9569 0.9564 0.9561 0.9563

0.9699 0.9682 0.9684 0.9745 0.9683 0.9586 0.9567 0.9565 0.9563 0.9564

0.0281

0.0280

28 ACS Paragon Plus Environment

Page 28 of 35

Page 29 of 35

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Ωi−Ωj(a) H62−H63 H56−H57 Non-bonded (C,C) C1,C3 C1,C6 C1,C10 C1,C11 C2,C5 C2,C10 C2,C11 C5,C11 C6,C10 C20,C25 Non-bonded (C,H) C1,H53 C2,H52 C5,H54 C10,H53

δ(Ωi,Ωj)direct 0.0296 0.0296

δ(Ωi,Ωj)KEM(b) 0.0295 0.0296

0.0579 0.0591 0.0695 0.0812 0.0669 0.0723 0.0557 0.0665 0.0618 0.0226

0.0578 0.0592 0.0695 0.0813 0.0669 0.0722 0.0556 0.0665 0.0618 0.0215

0.0501 0.0462 0.0473 0.0485

0.0501 0.0465 0.0473 0.0485

(a) See Fig. 2 for atom numbering. (b) For the bonded interactions, the average absolute deviation ( δ (Ω, Ω ')KEM − δ (Ω, Ω ')direct ) , the maximum absolute deviation (max( δ (Ω, Ω ') KEM − δ (Ω, Ω ')direct )) , the average % absolute error, and the maximum % absolute error are, respectively: (0.0008, 0.0042, 0.07%, 0.38% for C−C bonds), (0.0002, 0.0003, 0.02%, 0.03% for C−H bonds).

29 ACS Paragon Plus Environment

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 30 of 35

Table 5 Overall quality indicators of the KEM results. Comparison criterion(a) Sum of absolute deviations of all matrix elements Sum of all matrix elements of the direct calculation(b) Overall discrepancy of the matrix elements Frobenius distance of the KEM matrix from the matrix calculated directly Frobenius norm of the direct matrix Overall discrepancy of Frobenius distance with respect to the norm of the direct matrix (a) There are 4356 matrix elements in total. (b) The direct calculation is taken as the reference.

30 ACS Paragon Plus Environment

LDM 0.2182 296.0003 0.07% 0.0142 27.5986 0.05%

EDWCM 0.0127 48.7984 0.03% 0.0034 3.8741 0.09%

Page 31 of 35

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

31 ACS Paragon Plus Environment

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 32 of 35

Fig. 1 Top: The electron density contour maps in the σ v' plane along with the molecular graph. The lines joining the nuclei are the calculated bond paths, all solid lines except the H−H bond path in broken lines. The little dots between the nuclei on the bond path indicate the location of the bond critical points (BCP), each labeled by the value of the electron density at that point in a.u. The little dot inside rings are the ring critical points RCP also labelled by the density at these points in a.u. The contours from the outside inward are in a.u.: 0.001, 0.002, 0.004, 0.008, 0.02, 0.04, 0.08, 0.2, 0.4, 0.8, ... Bottom: The same molecular graph but this time displayed with the gradient vector field associated with the density showing the interatomic surfaces and the shapes of the atomic bains. The numerical labels near the nuclei are the integrated QTAIM atomic charges q(Ω). (All results are for the direct calculation (intact ribbon) at the HF/6-31+G(d,p) level of theory).

32 ACS Paragon Plus Environment

Page 33 of 35

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

33 ACS Paragon Plus Environment

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Fig. 2 Top: The numbering scheme of the intact hydrogen-terminated armchair graphene ribbon (C46H20) (2×7 aromatic rings in C2v symmetry). The two dark regions are the fission lines where the molecule is broken into kernel fragments. Second panel: Three single kernels (Kc, c = 1, 2, 3), each capped with hydrogen atoms to saturate the dangling bonds created by the fission. The capping hydrogens are placed at the location of the original carbon atoms of the neighboring kernel they replace. This results in C−Hcapping bonds that are much longer than the normal C−H bonds the reason for which they are displayed disconnected from the carbon (capping hydrogen atoms are enclosed in a thin circle to be easily spotted). The single kernels give rise to the second term in Eqs. 10-13. Third-to-fifth panels: Three hydrogen capped double kernels K12, K13, and K23 (Kab, a, b = 1, 2, 3; a < b). Enough atoms in the kernels are numbered to allow a unique mapping of every atom in a given kernel to the corresponding atom in the other kernels or the full intact ribbon. The double kernels give rise to the first term (the double sum) in Eqs. 10-13.

34 ACS Paragon Plus Environment

Page 34 of 35

Page 35 of 35

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 ACS Paragon Plus Environment