Electronic Properties of Graphene Nano-Parallelograms: A TAO-DFT Computational Study

1. Introduction

Carbon nanostructures are ubiquitous. This versatility comes from the ability of carbon to form bonds through sp, sp², and sp³ hybridization. The different types of hybridization can lead to various allotropes and complexes. Carbon is also available in nature as diamond and graphite. It forms all organic matters and is the basis of life. It is also an essential element in many drugs [1,2]. The nanoforms of allotropes were first introduced by the discovery of C₆₀ molecules [3,4], followed by the observation of carbon nanotubes [5], then followed by the isolation and many groundbreaking experiments in graphene [6,7]. All these nanostructures have various special properties that are interesting and can be advantageous for specific applications. The series of carbon nanoallotropes have created a fresh search for new 0D, 1D, and 2D nanomaterials [2]. Among carbon nanomaterials, graphene, which is a 2D nanomaterial, stands out, due to its fascinating properties and vast potential applications [7,8].

As graphene lacks a band gap, its finite nanostructures, such as graphene nanoribbons (GNRs), are also of great interest because of a wide range of possibilities to tune their properties for specific applications. GNRs can be thought of quasi-1D nanostructures carved from graphene. GNRs can exhibit intriguing electronic, optical, and magnetic characteristics [9,10,11,12,13,14,15,16], which can be useful for many applications [16,17,18,19]. In general, GNRs can be classified into two main categories: zigzag GNRs and armchair GNRs, based on their edge patterns [12,13,14]. Both zigzag and armchair GNRs of various widths have been successfully synthesized [12]. There are also other types of GNRs with different kinds of edge structures, such as the chevron, fjord, chiral, junction, cove, and gulf types of GNRs [20]. Although there is a rich literature on these electronic systems, ongoing theoretical and experimental studies aim to further explore the properties of GNRs.

Although zigzag and armchair GNRs have been extensively documented in the literature, reports on the parallelogram-shaped GNRs with zigzag edges (see Figure 1 for illustration), denoted as graphene nano-parallelograms (GNPs), are relatively scarce. In this study, the narrowest GNPs, which are n-acenes (i.e., acenes with n linearly fused benzene rings), are labeled as GNP[1,n]; the second narrowest GNPs are labeled as GNP[2,n], containing two parallelly fused n-acenes forming parallelograms; the third narrowest GNPs are labeled as GNP[3,n], containing three parallelly fused n-acenes forming parallelograms, and so on. For GNP[m,n], the m value specifies the GNP width, and the n value specifies the GNP length. Some GNPs, such as GNP[2,4] and GNP[2,5] (also called [n,2]peri-acenoacenes with n = 4 and 5, respectively) have been recently synthesized in molecular and crystalline forms [21,22]. Besides, for GNP[2,n] (with n = 3–5), their properties and lasing applications in the near-infrared (NIR) region have been recently reported [23]. In addition, a computational study has recently reported a few electronic properties of GNP[m,n] (with $m=n$ = 2–4) [24]. However, several properties of GNP[m,n] with different widths and lengths remain unavailable. Investigating the finite-size effects of GNP[m,n] is essential to gain a comprehensive understanding of the width- and length-dependent electronic properties of GNP[m,n]. Therefore, this computational study aims at exploring the electronic properties of GNP[m,n] with various values of m and n. Our previous study has revealed the presence of strong static correlation effects and intriguing properties in the longer zigzag GNRs [25]. Accordingly, it can be anticipated that GNPs, which have similar shapes as zigzag GNRs, can also possess strong static correlation effects in their ground states (as will be shown later).

Among existing electronic structure methods, Kohn-Sham density functional theory (KS-DFT) [26,27] is popular and valuable for investigating the ground-state properties of electronic systems, especially for single-reference (SR) systems (i.e., electronic systems possessing SR character in their ground states). In fact, KS-DFT has served as the backbone of theoretical condensed matter physics, quantum chemistry, and computational materials science [28]. For the sake of computational efficiency, KS-DFT calculations are commonly carried out using the local density approximation (LDA) [29,30] and generalized gradient approximation (GGA) [31] exchange-correlation (xc) energy functionals. However, KS-DFT with the LDA and GGA xc energy functionals can encounter severe problems for issues related to the self-interaction error, non-covalent interaction error, and static correlation error [32,33,34,35]. Among these qualitative errors, the self-interaction error can be reduced by the use of hybrid xc energy functionals [36,37,38], incorporating a fraction of Hartree-Fock exchange into the parent LDA or GGA xc energy functionals. Besides, a number of computationally efficient dispersion correction schemes [39,40], which can be directly added to the parent LDA or GGA xc energy functionals, are readily available for improving the description of non-covalent interactions. Nevertheless, the presence of strong static correlation effects in multi-reference (MR) systems (i.e., electronic systems possessing MR character in their ground states) has posed a formidable challenge, particularly for KS-DFT with the conventional LDA, GGA, and hybrid xc energy functionals. Generally, SR electronic structure methods (e.g., KS-DFT with the conventional LDA, GGA, and hybrid xc energy functionals as well as the Hartree-Fock theory) are unreliable for studying the ground-state properties of MR systems.

Typically, to explore the ground-state properties of MR systems, ab initio MR electronic structure methods [41,42,43,44,45,46,47], such as the MR configuration interaction (MRCI) methods and density-matrix renormalization group (DMRG) algorithms, are essential [48,49]. However, these methods can become impractical for large electronic systems due to the prohibitively high cost of performing ab initio MR electronic structure calculations. Accordingly, the demand for a reliable and efficient electronic structure method for large MR systems is very high.

Recently, thermally-assisted-occupation density functional theory (TAO-DFT) [50] has emerged as an effective solution for addressing the challenges posed by large MR systems. Generally, TAO-DFT focuses on improving the representability of ground-state electron density by incorporating fractional orbital occupation numbers, which can be efficiently computed using the Fermi-Dirac (FD) distribution function with some fictitious temperature $\theta$. The LDA [50], GGA [51], global hybrid [52], and range-separated [52,53] exchange-correlation-$\theta$ (xc$\theta$) energy functionals (i.e., the combined xc and $\theta$-dependent energy functionals [54]) can be incorporated in TAO-DFT. Besides, simple models for defining the optimal system-independent [55] and system-dependent [56] fictitious temperatures of an energy functional in TAO-DFT have been recently proposed. Note also that the difference among KS-DFT [26,27], TAO-DFT [50], and finite-temperature density functional theory (FT-DFT) [27,57] (i.e., three generally different electronic structure methods) has been properly discussed in a recent work [54].

Very recently, various TAO-DFT-related extensions, such as TAO-DFT-based ab initio molecular dynamics (TAO-AIMD) [58], TAO-DFT with the polarizable continuum model (TAO-PCM) [59], and a real-time extension of TAO-DFT (RT-TAO-DFT) [54], have been developed, expanding the capabilities of TAO-DFT to handle a more diverse range of applications. Moreover, TAO-DFT has also been adopted to investigate the various properties (e.g., electronic properties [25,60,61,62,63,64,65,66,67,68], hydrogen storage properties [61,63], spectroscopic properties [58,69,70], and equilibrium thermodynamic properties [58]) of MR systems at the nanoscale.

Consequently, in this computational study, we utilize TAO-DFT to obtain the electronic properties of GNP[m,n] (with m = 1–4 and n = 2–30), including the singlet-triplet energy gaps, vertical electron affinities / ionization potentials, fundamental gaps, symmetrized von Neumann entropy, active orbital occupation numbers, and real-space representation of active orbitals.

2. Computational Details

We perform all calculations with Q-Chem 4.4 [71], using the 6-31G(d) basis set. The electronic properties of GNP[m,n] (with m = 1–4 and n = 2–30) are computed using TAO-LDA (i.e., TAO-DFT with the LDA xc$\theta$ energy functional) with the recommended fictitious temperature $\theta$ = 7 mhartree [50].

3. Results and Discussion

3.1. Singlet-Triplet Energy Gap

For a neutral molecule, the singlet-triplet energy gap (ST gap) provides an understanding on the ground-state nature and how far the lowest singlet and triplet states are separated in the energy landscape. The ST gap can also be useful for understanding the radical nature of molecules [72,73,74,75,76] as well as for providing valuable information for photovoltaic applications [77].

Here, the ST gap ($E_{\mathrm{ST}}$) of GNP[m,n] is calculated as the energy difference between the lowest singlet and triplet states of GNP[m,n] (evaluated on the respective optimized geometries):

$$E_{\mathrm{ST}}=E_{\mathrm{UT}}-E_{\mathrm{US}},$$

(1)

where $E_{\mathrm{US}}$ / $E_{\mathrm{UT}}$ is the energy of the lowest singlet / triplet state of GNP[m,n], obtained with spin-unrestricted TAO-LDA.

As presented in Figure 2, GNP[m,n] possess singlet ground states (also see Table S1 in Supplementary Information (SI)) for all the cases examined (i.e., m = 1–4 and n = 2–30). The ST gap decreases monotonically, as the GNP length increases. Since MR systems typically possess very small ST gaps, the longer GNPs can have MR character in their ground states.

It is well known that for the lowest singlet state of an MR system, the spin-symmetry constraint, which must be satisfied by an exact theory, can be violated by SR electronic structure methods (e.g., KS-DFT with the conventional LDA, GGA, and hybrid xc energy functionals as well as the Hartree-Fock theory) [25,33,50,51,52,54,58,59,78]. In other words, spin-unrestricted and spin-restricted SR electronic structure calculations can yield different energies for the lowest singlet state of an MR system, leading to unphysical spin-symmetry breaking effects. Here, to examine if this spin-symmetry constraint can be satisfied by TAO-LDA, spin-restricted TAO-LDA calculations are additionally carried out for the lowest singlet energies of GNP[m,n] (evaluated on the respective optimized geometries). Within the numerical accuracy considered in the present study, spin-unrestricted and spin-restricted TAO-LDA calculations essentially yield the same energies for the lowest singlet states of GNP[m,n] (with m = 1–4 and n = 2–30), indicating that essentially no unphysical spin-symmetry breaking effects occur in our spin-unrestricted TAO-LDA solutions for all the cases examined.

3.2. Vertical Ionization Potential, Vertical Electron Affinity, and Fundamental Gap

The vertical ionization potential, vertical electron affinity, and fundamental gap of a neutral molecule in the ground state are important electronic properties. The vertical ionization potential is the energy change when an electron is removed from the neutral molecule (without changing its geometry), and the vertical electron affinity is the energy change when an electron is added to the neutral molecule (without changing its geometry). The fundamental gap is the difference between the vertical ionization potential and vertical electron affinity.

According to their definitions, on the spin-unrestricted TAO-LDA optimized geometry of ground-state GNP[m,n], we compute the vertical ionization potential

$${\mathrm{IP}}_v=E_{N-1}-E_N,$$

(2)

vertical electron affinity

$${\mathrm{EA}}_v=E_N-E_{N+1},$$

(3)

and fundamental gap

$$E_g={\mathrm{IP}}_v-{\mathrm{EA}}_v$$

(4)

of ground-state GNP[m,n], using multiple energy-difference calculations, with $E_N$ being the total energy of the N-electron molecule (i.e., GNP[m,n]) obtained with spin-unrestricted TAO-LDA.

The vertical ionization potential (${\mathrm{IP}}_v$), vertical electron affinity (${\mathrm{EA}}_v$), and fundamental gap ($E_g$) of ground-state GNP[m,n] are shown in Figure 3, Figure 4, and Figure 5, respectively (also see Tables S2 to S5 in SI). With increasing GNP length, the ${\mathrm{IP}}_v$ decreases monotonically, and the ${\mathrm{EA}}_v$ increases monotonically, leading to a monotonically decreasing $E_g$. These quantities are essential in materials science, spectroscopy, catalyst selection, etc., and are also important for the selection of candidate materials for electronic devices, solar cells, etc. Generally, electronic systems with the fundamental gaps in the range of 1–3 eV are desirable for photovoltaic applications. According to our TAO-LDA results, for all the cases studied in this work, GNP[1,n] (n = 10–30), GNP[2,n] (n = 8–30), GNP[3,n] (n = 7–30), and GNP[4,n] (n = 5–30) can be suitable candidates for photovoltaic applications.

3.3. Symmetrized von Neumann Entropy

To estimate the degree of MR character associated with ground-state GNP[m,n], we compute the symmetrized von Neumann entropy [51,52,64,78]

$$S_{\mathrm{vN}}=-\frac{1}{2}\sum _{\sigma =\alpha ,\beta }\sum _{i=1}^{\infty }\left\{f_{i,\sigma }\ln \left(f_{i,\sigma }\right)+(1-f_{i,\sigma })\ln (1-f_{i,\sigma })\right\},$$

(5)

where $f_{i,\sigma }$ (a number between 0 and 1) is the occupation number of the $i^{\mathrm{th}}$$\sigma$-spin (i.e., $\alpha$-spin or $\beta$-spin) orbital, calculated by spin-unrestricted TAO-LDA [50], approximately yielding the $i^{\mathrm{th}}$$\sigma$-spin natural orbital occupation number [60,79]. As can be seen in Eq. (5), a term with $f_{i,\sigma }$ = 0 or 1 has no contribution to the $S_{\mathrm{vN}}$, while a term with $f_{i,\sigma }$ significantly differing from 0 and 1 can cause a large increase in the $S_{\mathrm{vN}}$. As a result, for a SR system, all $f_{i,\sigma }$ values should be very close to 0 or 1, and hence, the $S_{\mathrm{vN}}$ should be vanishingly small. By contrast, for an MR system, some $f_{i,\sigma }$ values can significantly differ from 0 and 1, and hence, the $S_{\mathrm{vN}}$ can become very large.

As shown in Figure 6, the symmetrized von Neumann entropy ($S_{\mathrm{vN}}$) of ground-state GNP[m,n] increases monotonically with increasing GNP length (also see Tables S2 to S5 in SI). This indicates that the degree of MR character associated with ground-state GNP[m,n] should generally increase, as the GNP length increases.

3.4. Active Orbital Occupation Numbers

With increasing GNP length, to justify the increasing trend of $S_{\mathrm{vN}}$ of ground-state GNP[m,n], in spin-unrestricted TAO-DFT, there should be more active spin-orbitals [i.e., spin-orbitals with considerable fractional occupation numbers (e.g., in the range of 0.1–0.9)] and/or the occupation numbers of active spin-orbitals are closer to 0.5. In spin-restricted TAO-DFT, this indicates that there should be more active orbitals (i.e., orbitals with occupation numbers in the range of 0.2–1.8) and/or the occupation numbers of active orbitals are closer to 1.

To examine this, we plot the occupation numbers of active orbitals for the ground state of GNP[m,n], calculated using spin-restricted TAO-LDA [see Figure 7 (for m = 1), Figure 8 (for m = 2), Figure 9 (for m = 3), and Figure 10 (for m = 4)]. For the ground state of GNP[m,n] (containing N electrons), the HOMO (highest occupied molecular orbital) is the ${(N/2)}^{\mathrm{th}}$ orbital, the LUMO (lowest unoccupied molecular orbital) is the ${(N/2+1)}^{\mathrm{th}}$ orbital, and so on. In the present study, for brevity, HOMO and LUMO are denoted as H and L, respectively.

A common observation is that the shorter GNPs possess nonradical nature (i.e., all the orbital occupation numbers are very close to 0 or 2), and the longer GNPs possess increasing polyradical nature (i.e., as the GNP length increases, there are more active orbitals and/or the occupation numbers of active orbitals are closer to 1). Therefore, as the GNP length increases, the general feature of a transition from the nonradical nature of the shorter GNPs to the increasing polyradical nature of the longer GNPs is observed for each GNP width examined (i.e., m = 1–4). For the wider GNPs, the evolution of polyradical nature is faster.

3.5. Real-Space Representation of Active Orbitals

Here, we report the real-space representation of active orbitals (HOMO and LUMO) for the ground state of some representative GNP[m,n], including GNP[m,5] with m = 1–4 (see Figure 11), GNP[m,10] with m = 1–4 (see Figure 12), GNP[m,15] with m = 1–4 (see Figure 13), and GNP[m,20] with m = 1–4 (see Figure 14), obtained with spin-restricted TAO-LDA.

As shown, for a very narrow and short GNP (e.g., GNP[1,n] with n = 2–5), the HOMO and LUMO are delocalized over the entire GNP. However, as the width and length of GNP (i.e., m and n) increase, the HOMO and LUMO of ground-state GNP[m,n] have an increasing tendency to localize at the longer edges as well as the corners of smaller interior angles. In contrast to the edge localization of active orbitals found for the wider and longer zigzag GNRs [25,78], the edge/corner localization of active orbitals for the wider and longer GNPs is observed in this study, possibly due to the parallelogram shapes of GNPs.

4. Conclusions

In summary, we have reported the electronic properties of GNP[m,n] (with m = 1–4 and n = 2–30) obtained with TAO-DFT. As the longer GNPs exhibit significant MR character (e.g., polyradical nature), KS-DFT with the conventional LDA, GGA, and hybrid xc energy functionals can lead to unreliable results. Besides, because of the large electronic systems considered, performing ab initio MR electronic structure calculations can be impractical due to the prohibitively high cost. Owing to its nice compromise between efficiency and accuracy, TAO-DFT appears to be an ideal electronic structure method for this computational study.

According to TAO-DFT, the ST gap decreases monotonically with increasing GNP length, and GNPs possess singlet ground states for all the cases examined. As the GNP length increases, the vertical ionization potential and fundamental gap decrease monotonically, whereas the vertical electron affinity and symmetrized von Neumann entropy increase monotonically. Similar to the findings of previous studies on zigzag GNRs [25,78], the shorter GNPs possess nonradical nature, and the longer GNPs possess increasing polyradical nature. Therefore, there is a transition from the nonradical nature of the shorter GNPs to the increasing polyradical nature of the longer GNPs. Besides, the evolution of polyradical nature is faster for the wider GNPs. However, in contrast to the edge localization of active orbitals found for the wider and longer zigzag GNRs [25,78], we observe the edge/corner localization of active orbitals for the wider and longer GNPs in this study, possibly due to the parallelogram shapes of GNPs.