From Zero to Infinity: Customized Atomistic Calculations for Crystalline Solids — Applications to Graphene and Diamond

methodology has been very successfully applied and extrapolated to Si, Be, BeH, CdSe, MgH, crystals and nanocrystals, with almost chemical accuracy in most cases. Here, after a pedagogical and critical review of the earlier results, we introduce a new combined and expanded approach to comparatively describe the electronic and cohesive properties of diamond and graphene. For the later a drastically enlarged sequence of “nanocrystals” of well-chosen geometries and sizes up to 1440 atoms or 8190 electrons is used to verify earlier predictions and results. We have obtained in a simple and fast way the bandgap (5.4 eV) and the cohesive energy (7.34 eV/atom) of diamond with almost chemical accuracy; and we have fully rationalized (in a different perspective and prospective) the electronic and cohesive properties of graphene, with a tentative value of cohesive energy of 7.52 eV/atom. Strangely enough this value is larger than the one for diamond and is currently under investigation. Finally, we suggest that this methodology in its current simple and transparent form can be a first-line diagnostic, functional, and inexpensive computational tool. This is particularly true for quick assessments and comparative estimates, size-dependence studies, or cases where standard k-space methods or other advanced techniques either fail or demand unavailable computational resources. We review and combine two different atomistic-calculation approaches for macroscopic solids, applying them successfully to 2D graphene, in comparison to the 3D diamond with a dual target: 1) to gain novel physical insight about the Dirac


Introduction
Although the "natural" and customary way to study crystalline and nanocrystalline solids is the traditional band theory developed in k-space taking advantage of the periodicity of the system in one, two, or three dimensions, appropriate real-space models can give surprisingly good results for selective macroscopic properties of the infinite solids. Historically [1,2], besides the studies of atomic clusters per se, such methods were initially restricted to local or "localized" properties [1][2][3][4][5][6][7], such as impurity states or inner-core electronic properties [5][6][7] (in the bulk); and chemisorption (or adsorption) at the surfaces and interfaces [3,4]. However, "the cluster method" was occasionally (and rather unconventionally) extended to nonlocal "macroscopic" structural and electronic properties [8][9][10]; especially for relative comparisons, or when limited information was needed, or even when computational resources were scarce. The biggest problem in such calculations, which involve "an act of truncation" from a larger system, separating only a small part of it was the effect (i.e., the interaction) of the environment (rest of the system) with the cluster, which has to be, or assumed to be "small''. Otherwise, the validity of the results could be doubtful. Here these types of atomistic calculations, which for reasons which will become clear below we could call We demonstrate that a suitable atomistic method with judicially selected nanoclusters/ nanocrystals (in real space) supplemented with general symmetry and dimensionality arguments, can give surprisingly good results for macroscopic properties of the infinite crystalline solid, such as bandgaps, cohesive energies, as well as aromaticity (if any), at minimal computational cost and maximum physical insight. For graphene on top of these properties the present approach can successfully describe in real space and illuminate many of its exotic properties, which are usually introduced in k-space, such as Dirac points or topological insulators. An early version of this methodology has been very successfully applied and extrapolated to Si, Be, BeH, CdSe, MgH, crystals and nanocrystals, with almost chemical accuracy in most cases. Here, after a pedagogical and critical review of the earlier results, we introduce a new combined and expanded approach to comparatively describe the electronic and cohesive properties of diamond and graphene. For the later a drastically enlarged sequence of "nanocrystals" of well-chosen geometries and sizes up to 1440 atoms or 8190 electrons is used to verify earlier predictions and results. We have obtained in a simple and fast way the bandgap (5.4 eV) and the cohesive energy (7.34 eV/atom) of diamond with almost chemical accuracy; and we have fully rationalized (in a different perspective and prospective) the electronic and cohesive properties of graphene, with a tentative value of cohesive energy of 7.52 eV/atom. Strangely enough this value is larger than the one for diamond and is currently under investigation. Finally, we suggest that this methodology in its current simple and transparent form can be a first-line diagnostic, functional, and inexpensive computational tool. This is particularly true for quick assessments and comparative estimates, size-dependence studies, or cases where standard k-space methods or other advanced techniques either fail or demand unavailable computational resources. ___________________________________________________________________________________ "static" or "single", are not the subject of the present work, although a relatively large part of it could be thought of as a review article. The real focus of the present investigation is a conceptually different type of atomistic calculations, which for distinction we can call "dynamic", although such name could be misleading. In reality, such types of cluster calculation involve a judicially chosen sequence of clusters (n in number) of successively increased sizes, representing the molecular evolution to the infinite system approached in the limit of n→∞. Clearly such evolution is spatial not temporal. However, the name could be somehow justified under some conditions or assumptions, in which either a hidden velocity is implied, or the sequence converges to a superposition of states, which could be interpreted as a dynamic interchange. The corresponding author, who has work for a long time with real space cluster calculation and is familiar with both types of atomistic calculations, "static" [3][4][5][6][7][8][9][10], and "dynamic" [11][12][13][14][15][16][17][18][19][20][21][22][23][24][25][26][27], has recently introduced a special methodology in the last category of atomistic calculations based in guided extrapolation in accord with general dimensionality arguments, which in many cases can lead to very accurate (sometimes within "chemical accuracy") results [11][12][13][14][15][16][17][18]. Moreover, lately a different approach within the same "dynamic" methodology has been proposed by the same author (A. D. Z.) and collaborators for 2-dimensional (2D) graphene-based structures and graphene [19][20][21][22][23][24][25][26][27], which is based in a simple "geometrical" notion of aromaticity, general symmetry arguments, and the "shell model", [20] which underlines the properties of the polycyclic aromatic hydrocarbons (PAHs) constituting graphene nanocrystals or "graphene dots". It should be mentioned at this point that in modern terms the atomistic approach encompasses a large variety of molecular-type real space methods, which have the advantage over traditional band structure methods in kspace to be able to include, if needed, high level correlations by quantum chemistry methods, such as manybody perturbation theory (MBPT), coupled clusters, etc. In addition, current atomistic calculations can also include periodic boundary conditions, at various levels of correlation, most of the times surpassing in accuracy related band structure calculations. However, at the same time such calculations become more complex or cumbersome and more computationally demanding, losing the prime advantage of simplicity and transparency, which are the key features of the present work. Here we will focus on the recent "dynamic" approach mentioned above. It has been shown that such approach for 2D graphene-based structures is robust, powerful, and transparent, reproducing "advanced" many-body effects and peculiarities of graphene [21][22][23][24][25][26][27][28], such as Dirac cones and topological insulators, at a minimum amount of mathematical and computational cost. The same quantities for diamond, as will be shown below, have been performed in the literature with a large variety of techniques of usually much higher computational cost and (most of the times) lower accuracy [29][30][31][32][33][34][35][36]. A short summary of the methodology used especially in the two "dynamic" atomistic approaches, which are relatively new will be given below in the next two sections, which describe both the computational and conceptual framework of this work. This summary having a special pedagogical character is dedicated to help the familiarization and the working knowledge of the reader with relatively new and not widely known "techniques", in which the present authors happen to be very familiar. As such, this kind of summary unavoidably describes at the same time an account of the authors' recent work, with the risk to be considered somewhat lopsided. After the summary of previous results new results (in comparison to earlier results) are presented for graphene and larger hexagonal nanographenes (NGRs) and PAHs in the subsequent two sections. These new results expand and verify earlier predictions and results [23][24][25][26][27][28] leading also to new interpretations and predictions. We present two "dynamic" atomistic approaches for graphene, which we combine to derive additional new results and insight for its electronic and cohesive properties. These are subsequently compared with similar results for 3D diamond.

Computational details
The computational and technical details have been described before [23][24][25][26][27]. All geometrical structures have been optimized (with tight convergence criteria), using the hybrid PBE0 [37] functional, employing the 6-31G(d) basis set. This level of theory is considered adequate for the purpose of this work. The complete set of calculations was performed with the GAUSSIAN program package [38].

Conceptual and computational framework The numerical extrapolation technique
Historically and conceptually the cluster model of a solid (crystal) was developed and used for more or less "localized" properties (impurity levels, chemisorption, etc.) [1][2][3][4][5][6][7][8][9][10], as was mentioned earlier, using a single isolated ("small") cluster to model the solid (i.e., the region of interest of the solid). This approach, having the added advantage of the possibility to use (when needed, and possible) advanced correlation techniques of quantum chemistry [10,11,13] was quite successful [1][2][3][4][5][6][7][8][9][10][11], even when was used in borderline cases of questionable validity [8,9]. With the advent of computers (and the resulting increase of computer power) and the discovery of nanocrystals and nanocomposed materials, as porous silicon, the study of the variation of properties, such as the "bandgap" or the "optical gap" in terms of size of the nanoparticles was inevitable and very fruitful [10][11][12], culminated with the optical gap (and resulting photoluminescence) of silicon (and other [16][17][18] related) nanoparticles [10][11][12]. The "bandgap" in real space corresponds to the energy difference between the highest occupied molecular orbital (HOMO) and the lowest unoccupied molecular orbital (LUMO), i.e., the HOMO-LUMO gap. For better comparisons with (photoluminescence) experiments instead of the HOMO-LUMO gap, the "optical gap" is more appropriate [10][11][12], which includes electronic correlation through "many-body" effects. With time, besides HOMO-LUMO and optical gaps, which were very important for the observed photoluminescence from silicon nanoparticles despite the very low (optical) gap of silicon (≈1 eV), other properties, such as desorption and/or cohesive energies (for other nanoparticles or nano composed materials, used for hydrogen storage or similar applications) emerged [13][14][15] to be (nanoparticle) size-dependent. For hydrogen storage applications our group has worked with Mg, Be, as well as MgH and BeH nanocrystals and nonoclusters [13,15]. During these studies, a new feature emerged for the size dependence of (at least) these properties, which include HOMO-LUMO and optical gap, as well as cohesive or binding energies (to be defined below). It was found that the variation of these properties with size was more or less uniform and could be approximated by a simple algebraic relation of the form: where is the energy under question (HOMO-LUMO gap, optical gap, binding energy or cohesive energy), D is the "diameter" of the nanoparticle and N the number of (heavy) atoms of the nanoparticle. In the above expression A, B, C, F, m, and n are free parameters (constants) to be determined by the fit of the (experimental or theoretical) points (E, D) or (E, N). In the present study only theoretical (not experimental) results have been considered. Such semiempirical relations, which are in accord with quantum confinement [12], were widely used in the literature for the optical gap of Si nanoparticles [10][11][12][13], with quoted values for the exponents m or n varying between 0.76 and 1.3, whereas our initial value was 13 = 0.89 ± 0.15, close to unity. In the above relation for ( ) or ( ) observe that the value of the parameter C or A corresponds to the energy gap (band gap) of the infinite crystal, since as → ∞ (or → ∞), A (or C) becomes equal to (∞), e.g., The value we initially obtained was = 0.89 ± 0.15 and A = 1.02 ±0.25 eV, almost identical to the bandgap of crystalline Si. Subsequently, when the fits were further improved by including larger size Si nanocrystals [16], with higher levels of electronic correlation [10] (and better fit to the experimental measurements [10,16]) we have obtained a better value of the exponent m, m= 1.01 ± 0.05, which was practically equal to unity with a predicted optical gap [16,17] of Eg = E(∞) ≈ A = 1.08±0.1 eV. Therefore, it was highly suggestive to us that m should be equal to unity (m=1), and the gap would vary inversely proportional to the diameter D of the Si nanoparticles. In this case for spherical (or nearly spherical) nanoparticles the exponent n should be n = 1/3, e.g., the gap would vary inversely proportional to the cubic root of the number of nanoparticles N. Thus, by keeping the values of the exponents m and n fixed eq. (1) becomes: Obviously, for 1-dimensional (1D) or 2-dimensional systems the denominator of the exponent n would be 1 and 2 respectively representing the ("effective") dimensionality of the nanoparticles. Moreover, since the energy (HOMO-LUMO) gap is a zero th order measure of chemical hardness or (kinetic) chemical stability, it could be expected that binding energy/cohesive energy would also vary with the same exponents (-1 or -1/3) in terms of D or N, and furthermore, if relations (1) and (2), were a general "geometrical" relationship, it should be expected to be valid for other nanoparticles (besides Si) as well. Indeed, Vanithakumari, and Nanda [29][30] have illustrated, based on the liquid drop model and other general arguments and experimental results that that the cohesive energy [29] and other related quantities depend inversely proportional on sized D and that this is a rather "universal" relationship, independent of the particular material [30]. The binding (or cohesive) energy is the energy required to break a solid in isolated atomic species from which the solid was built, and as such is related with melting temperature, boiling temperature, Debye temperature, and desorption energy (in case of metal hydrides MgH, BeH, etc., used for hydrogen storage), etc. For a silicon or other nanoparticle (NP) or "quantum dot" (QD) terminated with hydrogen, consisting of NSi Si and NH H atoms, , the binding energy BENP or BEQD is defined as: where, ( ), ( ), and [ ] are the total energies of Si atom, H atom, and QD, respectively. Alternatively, we can introduce the cohesive energy ℎ, in which the interaction energy between the Si and H at the surface has been effectively removed: where is the chemical potential of hydrogen. The value of is obtained by DFT calculations (here PBE0/6-31g*) of the total energy differences before and after removing one or more equivalent H atoms. In several cases (not here) for simplicity the values of obtained from methane are used. Usually, we are interested in the binding or cohesive energies per (heavy) atom (here Nsi), which correspond to the cohesive energy of the infinite solid for Nsi→∞ This quantity is obtained by dividing the expressions (3, 4) by NSi. Using these expressions for the binding and cohesive energies and the ansatz of inverse D dependence we have reproduced the experimental binding energy of Be crystal (3.32 eV) with unexpectedly very good agreement (3.26 ± 0.06 eV), and we have predicted a value of 7.85 eV ± 0.02 eV for the binding energy of for [BeH2]∞ infinite system [15]. We have also predicted that the majority of the lowest energy stoichiometric BenH2n nanoparticles are chains or chain-like structures (1D). With the same methodology we have found [13] the experimental desorption energy of 75.5 kJ/mol for the infinite MgH system with remarkable accuracy (76.5 ± 1.5 kJ/mol). The absolute (not normalized) desorption or dissociation energy, ΔΕdtot(MgnH2n), for a MgnH2n nanocluster is defined as usual by the relation: whereas the normalized per H2 molecule desorption energy, ΔΕd (MgnH2n), is given as: In the above relations (5, 6) E(H2), E(Mgn), and E(MgnH2n) are the total energies of the H2 molecule, and of the Mgn and MgnH2n clusters, respectively. Normally these values include the zero-point energy (ZPE) corrections for all relevant structures (H2, Mgn, and MgnH2n). Finally, with the same methodology, as is illustrated in in Fig. 1, we have obtained with excellent accuracy the (optical) gap of crystalline Si (1.08 ± 0.01 eV) and with relatively good accuracy its cohesive energy, 5.01 ± 0.40 eV, whereas the experimental value is 4.63 -4.68 eV/atom [31,32]. The 5.01 eV value was obtained [17] with the M06 DFT functional, compared to 4.13 eV/atom with the B3LYP functional. Needless to say, that within the statistical error of the fit (±0.40 eV) the experimental result is bracketed by these values. It is important at this point to emphasize that the structures of the nanoparticles used in these calculations are nanocrystals (not nanoclusters), and also make the distinction between nanocrystal and nanoclusters. For a given number of Si (or other) atoms, the structure of the nanocluster corresponds to the lowest energy structure of the free Si (or other) cluster. In contrast, the structure of the nanocrystal corresponds to the lowest energy structure having the symmetry (and geometrical properties of the infinite crystal. These nanocrystals are terminated by hydrogens to saturate the surface dangling bonds (similarly to their experimental counterparts). Thus, the two could be totally different. For analogous isoelectronic (and "isoatomic") Ge nanoparticles of the same geometry and symmetry (Td), the calculated optical gap of Ge crystal was similarly obtained [16,18] as Eg= 0.95 ± 010 eV, well within the experimental value. Results for carbon (diamond) isoelectronic QDs are given in Section 3. Finally, it should be mentioned that in all or most of the nonlinear fits involved in these calculations attempts to improve the fits by letting the exponents m or n in (1) vary, fail. In most cases the quality of the fit worsens. On the other hand, in cases where the effective dimensionality of the samples is not obvious (as in Be or graphene nanocrystals), freeing the exponent n can be very illuminating, vide infra.

The numerical extrapolation technique
Besides the "classical" or general ("dynamical") atomistic approach described above, a conceptually different approach has emerged recently for graphene and its "nanocrystals" (NGRs and PAHs), suggested by the corresponding author and collaborators [20][21][22][23][24][25][26][27][28]. In this approach the series of nanocrystals not only have the same symmetry and dimensionality as the infinite system, but they also form a well-defined and well-coordinated (with mutually interrelated members) sequence of hexagonal (D6h symmetric) PAHs, named the "main sequence" of PAHs, which is proven to precisely (we could risk to call it "mathematically") converge to the (physical, chemical, and topological) properties of graphene, including the exotic ones (Dirac points, topological insulators, etc.) [23][24][25][26]. The members of the main sequence are not just model systems but (many of them) exist and are naturally abundant polycyclic aromatic molecules with well-defined number of carbon and hydrogen atoms (at their edges). Each PAH contains in a form of a babushka Russian doll all the previous PAHs (one inside the other) defining a Shell structure [24] with a benzene "nucleus", which renders this sequence of PAHs, a real molecular bridge from benzene to graphene. [25] The critical properties of the n th member of this sequence (electronic, aromatic, and topological properties) are uniquely determined by their sequence number, n, which is also called the shell number, for obvious reasons. Such properties include the morphology, symmetry and parity of both HOMO and LUMO (which in the limit of n→∞ converge to the valence and conduction bands, respectively), as well as the precise distribution of aromatic and non-aromatic ("empty") rings which defines the aromaticity pattern and the symmetry of the electronhole pair) [24,25]. These properties alternate as the shell number n increases, alternating between odd and even numbers. Thus, the full properties of the nth member of the sequence are naturally linked with those of the (n − 1) th and the (n + 1) th members. This provides not only a description of the end-system (graphene) but also of the intermediate members of the sequence, many of which could be considered as equally important as the end-system, defining a route or a spatial (and time?) evolution at the same time. As a result, in the limit of n→∞, where the properties of the (n − 1)th and the (n+1)th PAHs coincide with the properties of the nth, we are lead to a resonance (superposition) of aromaticity patterns, HOMO-LUMO orbitals, as well as full and empty rings [23][24][25] since the HOMO and full rings of the nth PAH are identical with the LUMO and empty rings of the (n−1) th and (n+1) th PAHs. This renders graphene a flat aromatic resonance structure, pretty much like benzene, justifying the name "super benzene" [24,25]. At the same time the interchange of HOMO and LUMO could naturally describe (in the limit n→∞) the touching of the valence and conduction bands (at the Dirac points) and the resulting valence-conduction [25,26] band on, as well as the associated dynamical breakdown of parity (which is a many-body effect) [31]. Moreover, as in the pair of nth and (n±1) th , one of them obeys both fundamental aromatic rules, Hückels rule of 4m + 2π electrons and Clar's rule of sextets (6k), whereas the other does not, we can deduce that graphene is "dynamically" aromatic by superposition (or resonance) of the two forms [20,24,25]. All this is illustrated schematically in Fig. 2. In the bottom part of Fig. 2, in Fig. 2(b) is illustrated that each of these PAHs always contains a "hard core" of characteristic benzene (occupied) orbitals, verifying the shell structure (geometrical and electronic) [24]. The same shell structure is responsible for the periodic variation of the electronic and aromatic properties of armchair graphene nanoribbons (AGNRs) with width [25,26] (W) and their resulting classification according to the number N of carbon atoms across the width in three categories: W = 3N, W = 3N−1, and W = 3N + 1, which in fact is classification according to their aromaticity patterns [26]. This is true not only for infinite length AGNRs, but also for finite length AGNRs, for which the study of the length dependence of their properties can reveal very interesting results, including phase transitions [27]. In general, the properties of graphene and graphene nanostructures, including PAHs are very sensitive to geometry, edge morphology, and, last but not least, aromaticity, as the present approach has revealed recently. Thus, the great importance of aromaticity for the electronic (and other) properties of graphene itself and of the large variety of NGRs is one of the important lessons learned here [27]. We should also realize that in many cases, where the variation of the properties with width or length is a real physical (Chemical), and not a model question, the atomistic or molecular description is the most appropriate and "natural" method to use, which may be considered as a borderline case where the initial and final system are just the same.

Combining the two approaches: The numerical extrapolation for graphene
It would be interesting at this point, after this discussion of the "aromatic" extrapolation method, to apply the numerical extrapolation method to graphene. According to our earlier conclusions we should try to fit the HOMO-LUMO energy gap in terms of N (the number of carbon atoms) to a form like: Fig. 3 describes the results of this attempt. As we can see in Fig. 3(a), both attempts (with C = 0 and C ≠ 0) with exponent n = -1/2 fail. Strangely enough the best fit is obtained for non-zero C, and exponent n = -1/4; but what is even more strange is that the parameter C, which as we have seen earlier should be equal to the band gap of graphene, C= E(∞), is negative: C ≈ −2.60 eV. Before we try to access or rationalize this result, we have to make sure of its essential correctness, if any. To this end, in the next section we have extended and expanded our calculations to include more and much larger PAHs from the main sequence.

Expanding and extending the combined approach
Enlarging the molecular bridge to graphene In Fig. 3(b), we show the analogue of Fig. 3(a), for the PAHs with shell number from 1 to 15. The 15th PAH consists of 1350 C atoms and 90 H atoms in the periphery containing a total 8190 electrons.
As we can see in Fig. 3(b), the fit with exponent 1/4 is by far the best. But again, as before, Fig. 3(a), the parameter C = E(∞), is still negative C ≈ −2.60 eV. However, before we venture a plausible interpretation of these results, we should inspect the morphology and symmetry of both HOMO and LUMO orbitals and their interrelation, for the larger PAHs. This is done in Fig. 4. As we can see in Fig. 4, the alteration of symmetry and parity of HOMO, LUMO gaps described above for the smaller members of the main sequence is preserved and remains the same. Moreover, we can observe in Fig. 4 that as the size is getting larger, the largest values of the electronic density in the frontier orbitals are concentrated in the periphery of the PAHs. This is related with the shell structure combined with the zigzag morphology of the edges [24,28]. The same effect in rectangular AGNRs leads to topological edge (or end) states [27]. The trend of parity alteration, dictated by inversion symmetry and reflecting the electron-hole symmetry [24,25], also demands that the two degenerate HOMOs of even parity, e1g (with 2-dimendional (2D) symmetry group representations) obtained for odd shell numbers, will be associated and linked with the corresponding (2D) LUMOs of odd parity e2u. This is complemented and corroborated also by the linking of the full and empty rings between successive PAHs [25]. In addition, the same even parity HOMO pair will be also linked with the odd-parity HOMOs (e2u) of the previous and next PAHs in the sequence, with even shell numbers. Obviously the same is true for the odd parity (e2u) HOMOs, which are linked to the e1g LUMOs, and the even parity e1g HO-MOs of the neighboring PAHs. This is emphasized in Figure S2, which shows in better resolution the morphology and symmetry of the two components of the 2D degenerate molecular orbitals, HOMO (1, 2) and LUMO (1, 2), for one odd-shell PAH (with shell number n = 9) and one PAH with even shell number (n = 14). Moreover, to illustrate that the effective dimensionality of 4, and the accompanying negative value at infinity are purely electronic characteristics (involving the electrons in the frontier orbitals), we show in Fig. 5 the corresponding fits and predictions for the cohesive and binding energies. These are given by Equations (3) and (4) above, divided by the number of heavy atoms, i.e. the number of carbon atoms N, in place of NSi in eqs. (3)(4).  Moreover, to illustrate that the effective dimensionality of 4 (and the accompanying negative value at infinity) is a purely electronic characteristic (in particular of the electrons in the frontier orbitals), we show in Fig. 5 the corresponding fits and predictions for the cohesive and binding energies given by Equations (3) and (4) above, divided by the number of heavy atoms, i.e. the number of carbon atoms N. The results for cohesive and binding energies are practically identical.
As we can see in Fig. 5, the variation with the number of carbon atoms for both of them correctly and fully is described by the exponent in agreement with the twodimensional nature of graphene, as expected. In addition, both curves predict a (positive) cohesive energy of graphene of about 7.7 eV/atom, which appears to be larger than the cohesive energies of both diamond and graphite [31]. This is a very significant result, but we will discuss and compare this value later in comparison to analogous results for diamond.

Possible interpretations and plausible suggestions
The 4th dimensional variation of the HOMO LOMO gap with the number atoms should be fully related with inversion symmetry conflict between the full molecular group (D6h) and the (sub)lattice group (D3h), which is also responsible for all or most of the exotic properties of graphene [25], including the existence of two aromaticity patterns with trigonal (D3h) and hexagonal (D6h) symmetries, respectively [25,28]. In the present case inversion symmetry conflict demands the alteration and interlinking between HOMO and LUMO orbitals (as well as with LUMO and HOMO orbitals of neighboring PAHs) as well as of aromaticity patterns, as is shown in Fig. 2(a). In addition, this trend reflecting the electron-hole symmetry [24,25], is also responsible for the interlinking of the full and empty rings between successive PAHs [25]. In all PAHs (and, by extrapolation, graphene) The HOMO and LUMO molecular orbitals are 2 × 2. The 2D (doubly) degenerate orbitals, holding a total of 4 electrons and 4 holes. This can be interpreted as an effective dimensionality of 4. Quantum mechanically this is also in full accord with the 4-dimensional nature of the effective relativistic Dirac equation, which (effectively) describes the electronic properties of graphene.
The negative value of the fitted parameter C = Eg(∞) ≈ −2.60 eV is not independent of 4th dimensional variation of the HOMO LOMO, if we recall valence-conduction band inversion and related "dynamical" breakdown of parity, the latter been a many body effect related with inversion symmetry breaking [33]. Moreover, the value of −2.60 eV, which is practically equal to the hoping integral -t entering the tight binding Hamiltonian (and the Hubbard model Hamiltonian) [28] (where are the creation, annihilation operators, respectively, which create and annihilate an electron at site i with spin σ) is highly suggestive that C in this case, with coupled HOMO-LUMO orbitals could be equal to t. After all, the HOMO-LUMO gap is normally associated with the transfer of an electron from the ground (HOMO) to the 1st excited state (LUMO). In the case of hoping it could be argued that this is rather related to the transfer of an electron from one cite A to the nearby cite B. Otherwise, the mixing of HOMO and LUMO implied for large values of the shell number n, dictated by inversion symmetry could be considered to lead to open shell singlet state, within the one-body DFT framework, which loosely speaking, could be linked with an effectively negative value of the gap, equal to the hoping (or resonance) integral. Indeed, initial estimations (based on the fits of the first 7 PAHs of the main sequence) for the value of the shell number n at which band-crossing or band inversion could occur suggested values around n ≈ 13 with a margin of 2-3 units. This could be considered as appositive prediction since the lowest energy state of the n = 10 PAH is found to be open singlet, with small distortions of symmetry. Nevertheless, all these suggestions are at least plausibility arguments which need further studies and work, which we plan to explore in future work. The 4th dimensionality and the negative value of Eg(∞) discussed above are clearly plausibility arguments and interpretations which are not formulated mathematically or proven by stringent facts. This would be completely outside the scope and the principles of the present investigation, which emphasizes simplicity and physical insight. However, it is natural (and welcomed) these results to generate possible doubts and questions concerning the suitability of the present approach based on DFT calculation, and in particular the choice of (only) one specific functional, PBE0. Obviously, the full answer involves Dirac's equation and/or many-body-theory of graphene, which are well beyond the scope of the present work. However, we can examine the same results with a different hybrid functional such as B3LYP. This functional is very popular and well-known, and includes "exact exchange", which is absolutely necessary for this problem, in view of inversion symmetry conflict mentioned earlier. This is because the exchange interaction is more sensitive to inversion symmetry. More elaborated or advanced methods are no needed at this point since we are only interested in plausibility arguments and tentative suggestions, focusing on trends, which sometimes can be considered equally or even more important than "exact" numbers. Fig. 6 shows the B3LYP fits in comparison to the PBE0 ones. As we can see in this figure, both the "4th dimensionality and the negative Eg(∞) are fully reproduced. In addition, the value of Eg(∞) is the same, within statistical uncertainty, Eg(∞) = t ≈ -2.50 eV, which is also well within the uncertainty of the value of t in the literature. Thus, the results revealed here, although plausibility arguments and tentative suggestions, could be considered rather universal in the sense that they are well reproduced by other (hybrid) functionals. Moreover, we should also take into account that similar preliminary calculations for silicene (in the buckled D3d geometry) produced very good quality fits with the same key characteristics i.e. 4th -dimensionality and negative parameter C = Eg(∞) = t. The value of the parameter t in this case is equal to ~ -1.2 eV, in excellent agreement with the quoted results in the literature for silicene. Fig. 7(a) shows the HOMO-LUMO gap variation of carbon nanocrystals. As we can see in Fig. 7, the N −1/3 fits are not so good as the N -1/4 and N -1/3 for graphene in Figs. 3 and 5. Nevertheless, these fits are better than any other choice of exponent, and even better than the ones with free exponents (to be determined from the fit). In addition, the N −1/3 is the only one which reproduces with almost chemical accuracy the band gap of diamond 5.42 eV, compared to the experimental value of 5.46 eV. This is unexpectedly good result. For obvious reasons we do not have here artificial (or real) negative values, as in graphene. Of very good accuracy is also the fit of cohesive and binding energies shown in Fig. 7(b), which predict a 7.50 eV/atom cohesive energy which in full agreement with the results obtained by the advance Quantum Monte Carlo (QMC) [31] method. (7.503 eV/atom).

Comparisons to diamond
Both quoted values do not include vibrational zeropoint energy (ZPE) which can be corrected by a uniform shift of 0.176 eV/atom [31] which brings both values (ours 7.34 eV/atom and the QMC 7.327 eV/atom) in excellent agreement with the experimental value of 7.346 eV/atom. Interestingly enough, the latest quoted values in the literature for the cohesive energy of diamond using periodic second-order Moller-Plesset perturbation theory (MP2) is 7.91 eV/atom, whereas the results of periodic couple clusters theory with single and double excitations (CCSD) give 7.04 eV/atom respectively [32]. Both values are quite far away from the present results, the Monte Carlo results, and the experimental values. It is also important to emphasize that the CCSD results took 10000 CPU-hours in a cluster of 512 GB memory, while the present results were obtained in a few CPU-hours using a single 4-processors workstation with 16 MB memory. Obviously, not all results can be obtained with the same accuracy using our present approach, although the CPU requirements would be expected to be of about the same level. Corresponding periodic DFT calculations are clearly much less computationally demanding compared to the CCSD calculations, but undoubtedly more demanding, compared to the present approach, with varying accuracy. For example, He et. al., [34] using periodic DFT calculation at the local density, and generalized density approximations (LDA, and GGA, respectively), have obtained cohesive energies of 8.997 eV/atom (LDA), and 7.693 eV/atom (GGA) respectively (about 4% off the experimental for the GGA). However, the corresponding bandgaps were found at 4.648 eV, and 4.635 eV respectively, almost 1 eV smaller than the experimental value (more than 20% off). Earlier periodic LDA results employing norm conserving pseudopotentials [35] have obtained the cohesive energy of diamond within a range [-16%, +6%] depending on the quality of the basis set, at a very low computational cost, but with bandgaps underestimated by about 20%. To get better bandgaps one has to use hybrid functionals. Muscat et. al., [36] have used hybrid B3LYP and obtained a bandgap of 5.8 eV (+0.3 eV above the experimental value). Thus, the present method, on the basis of the computational cost and the accuracy of both electronic (bandgap) and cohesive (cohesive energy) properties it can clearly be considered superior, compared to the methods mentioned above. Coming back to graphene now, for which cohesive property calculations are rare, and experimental results, as far as we know, are not yet available, our present calculations should be considered theoretical predictions, which could further test this approach. For graphene the ZPE correction amounts to 0.166 eV/atom [31], which renders the predicted value for the cohesive energy of graphene (obtained from the fit in Fig. 6) equal to 7.52 eV/atom, still larger than the one for diamond, whereas the value obtained by QMC is 7.298 eV/atom, smaller than both diamond and graphite. This could be settled in future work.

Conclusion
We have successfully applied and critically reviewed, analyzed in a rather tutorial way a modern version of a relatively old technique. This technique originally known as the "cluster" model, was later generalized to "atomistic calculations", especially by solid state physicists. In modern times the atomistic approach has been evolved into a full-fledged successful field of materials science of very good accuracy, which includes a large variety of real space methods, incorporating high level correlation by quantum chemistry methods, as well as periodic boundary conditions, allowing super shell calculations; at the cost of becoming more complex, less transparent, and more computationally demanding. We have seen computational cost of 10000 CPU-hours in clusters with 512 GB memory, [32] in comparison to a few CPU hours in single 4processor workstations with 16 MB memory, using the present method, with comparable or better accuracy. Such kind of accuracy was rather unexpected, since the original target was simplicity and transparency. Nevertheless, such result clearly illustrates the hidden power of "simple" and general symmetry and dimensionality arguments which are largely misunderstood or undervalued. This is the major point of the present work. The viewpoint presented here gives emphasis in the physical content, simplicity, and consistency, demanding a critical selection of the various parameters and methods. This approach was initially invoked as an order-of-magnitude estimate for quick comparisons, seeking maximum information at minimum computational cost, relying in simplicity and transparency. However, it was finally evolved to be a simple, useful, and powerful tool, without losing its initial target. We have illustrated and reviewed the high potential of this "strategy" in a large variety of applications from our own work ranging from comparative studies, quick estimates, sizedependence studies, all the way up to full and accurate studies, in a rather large range of applications. In the present work we have applied this approach to a comparative study of 3D diamond and 2D graphene, in which we have obtained with almost chemical accuracy the bandgap (5.4 eV) and the cohesive energy (7.34 eV/atom) of diamond, and we have fully rationalized in a different perspective (and prospective) the electronic and cohesive properties of graphene. For the cohesive energy of graphene, we have obtained a tentative value of 7.52 eV/atom, which is currently under further investigation since it is larger than the cohesive energy of both diamond and graphite. Finally, we suggest that this methodology can be very useful and successful in a large variety of future applications in computational material science. This methodology should be particularly useful for length or size-dependence studies, relative comparisons, and semi-infinite systems with mixed nano, micro, or macro components. Obviously, not all results nor all properties can be obtained with the same accuracy for all materials, using this approach. We anticipate, however, its expanded future applications for more new results and applications with more rich physical insight and understanding, at minimum computational cost. The present authors share the view that trends (stressing the physical content and insight) should be more important than numbers. Nevertheless, sometimes (as in the present study) one could practically have both.

Supporting information
Supporting informations are available online at journal website.