Electronic Properties of Triangle-Shaped Graphene Nanoflakes from TAO-DFT.

Reliable prediction of the properties of nanosystems with radical nature has been tremendously challenging for common computational approaches. Aiming to overcome this, we employ thermally-assisted-occupation density functional theory (TAO-DFT) to investigate various electronic properties (e.g., singlet-triplet energy gaps, vertical ionization potentials, vertical electron affinities, fundamental gaps, symmetrized von Neumann entropy, active orbital occupation numbers, and visualization of active orbitals) associated with a series of triangle-shaped graphene nanoflakes with n fused benzene rings at each side (denoted as n-triangulenes), which can be extended from triangulene. According to our TAO-DFT results, the ground states of n-triangulenes are singlets for all the values of n studied (n = 3, 5, 7, 9, ..., and 21). Moreover, the larger the values of n, the more significant the polyradical nature of n-triangulenes. There are approximately (n - 1) unpaired electrons for the ground state of n-triangulene. The increasing polyradical nature of the larger n-triangulenes should be closely related to the fact that the active orbitals tend to be mainly concentrated at the periphery of n-triangulenes, apparently increasing with the molecular size.

Introduction

Because of its promising properties and potential applications, graphene has been extensively studied by several researchers in recent years.[1−6] For instance, the high carrier mobility, saturation velocity, and long spin diffusion length of graphene have yielded the possibility of developing fascinating electronics and spintronics based on graphene.[1,2,6] Nonetheless, owing to the lack of an energy band gap, graphene is not suitable for transistor applications. A method to create a finite band gap in graphene is to cut it into graphene nanoflakes.[7−9] However, the properties of graphene nanoflakes can be highly influenced by their shapes, sizes, and periphery conditions.[10,11] Consequently, it is essentially important to investigate how these factors affect the properties of graphene nanoflakes to fully utilize the promising potential of graphene in electronics and spintronics.

Among graphene nanoflakes with several possible shapes, graphene nanoribbons (GNRs), which are narrow strips of graphene, have attracted considerable attention from the research community and have been regarded as fascinating quasi-one-dimensional materials for electronic nanodevices.[7−9,12−24] Owing to their low dimensionality and quantum size effect, GNRs can possess energy band gaps for transistor operation with excellent switching speed and high carrier mobility. In particular, zigzag GNRs are anticipated to possess edge-localized states, which could be important ingredients for electronics and spintronics based on GNRs.

On the other hand, other shapes of graphene nanoflakes may deserve our investigation as well. Very recently, the on-surface synthesis of triangulene (a molecule with two unpaired electrons) has been successfully achieved by Pavliček et al. using a needle-like scanning tunneling microscopy tip (i.e., an unconventional synthesis).[25] Besides, they found that triangulene did not bond to the underlying metal surface, which was highly unexpected. Note, however, that the microscopy tip used was not magnetic, and hence unable to directly identify the spin polarization of the ground state of triangulene.[26] Nonetheless, it can be anticipated that triangle-shaped graphene nanoflakes, which can be extended from triangulene, may be created using this new type of chemical synthesis in the near future. Consequently, in this study, we focus on a series of triangle-shaped graphene nanoflakes with n fused benzene rings at each side (denoted as n-triangulenes), as illustrated in Figure . Note that n-triangulenes (CH3), which belong to the category of polycyclic aromatic hydrocarbons, are delocalized π-conjugated systems. A systematic study of the properties of n-triangulenes is crucially important to fully understand their potential applications.

Figure 1

Geometry of n-triangulene, consisting of n fused benzene rings at each side.

To date, the properties of n-triangulenes have been investigated mainly via computational approaches.[27−34] However, Kohn–Sham density functional theory (KS-DFT)[35] with commonly used (e.g., the local density approximation (LDA), generalized-gradient approximation (GGA), hybrid, and double-hybrid) exchange–correlation (XC) density functionals cannot adequately describe the ground-state properties of systems with radical nature,[36−38] such as n-triangulenes (as will be shown and discussed later). On the other hand, it is common to carry out ab initio multireference electronic structure calculations[13,16,39−41] for an accurate prediction of the ground-state properties of systems with radical nature. Nonetheless, as the system size (n) grows, the number of electrons in n-triangulene, N = 6n2 + 27n + 9, can quickly increase. Therefore, it is computationally intractable to carry out sufficiently accurate multireference electronic structure calculations for n-triangulenes, especially for those with the larger n values. Accordingly, it remains tremendously challenging for common computation approaches to reliably describe the ground-state properties of the larger n-triangulenes.

With the aim to study the ground-state properties of nanosystems with radical nature at low computational cost with reasonable accuracy, we have recently developed thermally-assisted-occupation density functional theory (TAO-DFT),[15] a density functional theory with fractional orbital occupations produced by the Fermi–Dirac distribution function. In TAO-DFT, the entropy contribution (e.g., see eq 26 of ref (15)), which are explicitly dependent on the fictitious temperature (θ) and orbital occupation numbers, can approximately describe the static correlation energy of a system (e.g., see Section III.E of ref (15) for the physical arguments and Section IV of ref (15) for the numerical investigations), even when the simplest LDA XC density functional is adopted. More sophisticated XC density functionals, such as the GGA,[19] global hybrid,[24] and range-separated hybrid[24,42] XC density functionals, can be adopted in TAO-DFT as well. Besides, an approach that determines the fictitious temperature in TAO-DFT in a self-consistent manner[43] has been recently developed to improve the overall accuracy of TAO-DFT for general applications. As TAO-DFT is computationally efficient (i.e., similar to KS-DFT), we have employed TAO-DFT for the study of the electronic properties of various nanosystems with radical nature in recent years.[20,22,44−49] Very recently, TAO-DFT and related approaches have also been employed for the study of the properties of several nanosystems of radical nature (e.g., cyclic nanorings and carbon nanotubes) by other researchers.[50,51]

Despite attracting great interest recently, graphene nanoflakes with pronounced radical nature can be wrongly treated with standard electronic structure methods (e.g., KS-DFT with commonly used XC density functionals); hence, more investigations and new conceptions are really needed at the nanoscale. Owing to an excellent trade-off between computational cost and accuracy, employing TAO-DFT for a comprehensive study on graphene nanoflakes with different shapes, edges, and sizes is well-justified. Consequently, in this work, TAO-DFT is adopted for the study of various electronic properties associated with n-triangulenes (n = 3, 5, 7, 9, ..., and 21).

Computational Details

All calculations are carried out with Q-Chem 4.0,[52] using the 6-31G(d) basis set and the numerical grid consisting of 75 radial points in the Euler–Maclaurin quadrature and 302 angular points in the Lebedev grid. Results are computed using TAO-LDA[15] (i.e., TAO-DFT with the LDA XC functional[53,54] and the LDA θ-dependent functional (e.g., see eq 41 of ref (15))) with the fictitious temperature θ = 7 mhartree (e.g., see Section V of ref (15) for this choice of θ). For comparison, some results are also obtained with KS-LDA (i.e., KS-DFT with the LDA XC functional[53,54]). As KS-LDA is simply TAO-LDA with θ = 0, the performance of KS-LDA is shown here to assess the significance of TAO-LDA.

Although the XC density functionals at the upper rungs of Jacob’s ladder, such as the GGA[19] and global hybrid[24] XC density functionals, can also be adopted in TAO-DFT, they improve upon TAO-LDA mainly for the properties governed by short-range XC effects (e.g., the atomization energies and barrier heights of systems with nonradical character) not for the properties governed by strong static correlation effects (e.g., the singlet–triplet energy gaps and fundamental gaps of systems with radical character).[15,19,24] As shown in our previous studies, the GGA and global hybrid XC density functionals in TAO-DFT have similar performance as TAO-LDA for the electronic properties of linear acenes (i.e., systems with polyradical character).[15,19,24] Accordingly, the electronic properties of n-triangulenes obtained with TAO-LDA are expected to be qualitatively similar to those obtained with the GGA and global hybrid XC density functionals in TAO-DFT.

It is worth mentioning that if the fundamental gap is estimated by the highest occupied molecular orbital (HOMO)–lowest unoccupied molecular orbital (LUMO) gap (i.e., the energy gap between the highest occupied molecular orbital (HOMO) and the lowest unoccupied molecular orbital (LUMO)) in KS-DFT, the LDA XC functional tends to underestimate the fundamental gap due to the lack of the derivative discontinuity (DD) of the LDA XC functional.[55] Therefore, in this work, the fundamental gap is obtained with multiple energy-difference calculations (see eqs –4), circumventing the above DD issue and showing a reasonably good accuracy for a wide variety of atoms and molecules (see, e.g., Table 5 of ref (56)). For TAO-LDA, the fundamental gap cannot be directly estimated by the HOMO–LUMO gap due to the possibility of fractional orbital occupations (see, e.g., Section III of ref (43)); hence, the fundamental gap is obtained with the aforementioned multiple energy-difference calculations.[19]

Besides, the expectation value of the total spin-squared operator, ⟨Ŝ2⟩, is computed for the assessment of the degree of spin contamination in KS-DFT. For a system with significant radical nature, the ⟨Ŝ2⟩ value calculated using commonly used XC functionals in KS-DFT can differ considerably (e.g., 10% difference or more) from the exact value S(S + 1),[57] where S can be 0 (singlet), 1/2 (doublet), 1 (triplet), 3/2 (quartet), 2 (quintet), and so forth. For such a system, the results obtained with the commonly used XC functionals in KS-DFT can be unreliable due to spin contamination (i.e., the artificial mixing of different electronic spin states).[15,19,20,44,48,49,58−60]

Similar to the case of spin-unrestricted Hartree–Fock (HF) theory, the expectation value of the total spin-squared operator of spin-unrestricted KS-DFT is evaluated with the Kohn–Sham (KS) wavefunction (i.e., a single Slater determinant constructed from a set of KS orbitals) associated with the ground-state density (e.g., see eq 3 in the Supporting Information (SI) of ref (21)). However, in TAO-DFT,[15] the ground-state density of a physical system at zero (physical) temperature is represented by the thermal equilibrium density of a noninteracting system at the fictitious temperature θ; hence, the wavefunction associated with this ground-state density remains unknown. Therefore, the expectation value of the total spin-squared operator of spin-unrestricted TAO-DFT remains undefined due to the lack of wavefunction information. Therefore, we also investigate the consequences of spin contamination by comparing the difference between the spin-restricted and spin-unrestricted energies for the lowest singlet state of n-triangulene, which should be zero for the exact theory due to the symmetry constraint.[15,17,19,20,44,48,49]

Results and Discussion

Singlet–Triplet Energy Gap

To obtain the ground state (i.e., the energetically preferred spin state) of n-triangulene (n = 3, 5, 7, 9, ..., and 21), calculations based on spin-unrestricted TAO-LDA and KS-LDA are carried out to determine the lowest singlet and lowest triplet states of n-triangulene, with the respective geometries being fully optimized. Afterward, the singlet–triplet energy gap of n-triangulene is calculated aswhere ES and ET are the energies of the lowest singlet and lowest triplet states, respectively, of n-triangulene. Owing to severe spin contamination (see the discussion below), it is extremely difficult to obtain the converged KS-LDA results for the larger n-triangulenes (e.g., n > 9). Therefore, the KS-LDA results are only reported for the smaller n-triangulenes (e.g., up to n = 9) due to convergence problems.

Figure plots the EST value of n-triangulene as a function of the side length, computed using spin-unrestricted TAO-LDA and KS-LDA. Based on TAO-LDA, EST decreases monotonically with the increase of n, and n-triangulenes possess singlet ground states for all the values of n studied (n = 3, 5, 7, 9, ..., and 21). By contrast, from the KS-LDA results, the triplet states of n-triangulenes (e.g., n = 3, 5, 7, and 9) have lower energies than the corresponding singlet states (see Table S1 in the Supporting Information (SI) for detailed numerical results).

Figure 2

Singlet–triplet energy gap of n-triangulene computed using spin-unrestricted TAO-LDA and KS-LDA.

Here, we investigate the reasons of discrepancies by computing the ⟨Ŝ2⟩ values for the lowest singlet and lowest triplet states of n-triangulene using spin-unrestricted KS-LDA. On the basis of the computed ⟨Ŝ2⟩ values (see Table ), the lowest singlet and lowest triplet states of n-triangulene, computed using spin-unrestricted KS-LDA, are strongly influenced by spin contamination, with the sole exception of the lowest triplet state of 3-triangulene. Moreover, the degree of spin contamination in the lowest singlet and lowest triplet states of n-triangulene rapidly increases with increasing n, implying that the larger n-triangulenes may have an increasing polyradical nature in the lowest singlet and lowest triplet states. Therefore, the contradictory results obtained with KS-LDA should be artifacts associated with spin contamination.[15,19,20,44,48,49,58−60] It should be mentioned that the errors caused by spin contamination are not systematic errors, which can adversely affect the difference in energy between the spin-states of n-triangulene (e.g., EST).

Table 1

Expectation Value of the Total Spin-Squared Operator ⟨Ŝ2⟩ for the Lowest Singlet and Lowest Triplet States of n-Triangulene Computed Using Spin-Unrestricted KS-LDAa

n	3	5	7	9
singlet	1.0049	2.0033	3.0085	4.0202
triplet	2.0137	3.0121	3.9960	5.0067

Note that the exact value S(S + 1) is 0 for the singlet state and 2 for the triplet state.

On the other hand, owing to the symmetry constraint,[15,17,19,20,44,48,49] for the lowest singlet state of n-triangulene, the spin-restricted energy obtained with the exact theory should be identical to the respective spin-unrestricted energy. Nevertheless, KS-DFT with commonly used XC functionals can fail to fulfill this condition due to the aforementioned spin contamination. Here, we examine if the symmetry-breaking effects occur by additionally performing spin-restricted TAO-LDA and KS-LDA calculations for the lowest singlet states of n-triangulenes, with the respective geometries being fully optimized. Because of the spin contamination issue mentioned above, the difference between the spin-restricted and spin-unrestricted energies, obtained with KS-LDA, for the lowest singlet state of n-triangulene is considerably large (e.g., the energy difference is 10.78 kcal/mol for n = 3, 8.70 kcal/mol for n = 5, and 12.26 kcal/mol for n = 7) when compared with the magnitude of the corresponding EST value, obtained with spin-unrestricted KS-LDA (e.g., see Figure ). By contrast, the spin-restricted energy obtained with TAO-LDA is essentially the same as the respective spin-unrestricted energy (i.e., within the numerical precision considered in the present work) for the lowest singlet state of n-triangulene, indicating that unphysical symmetry-breaking solutions are not generated by our spin-unrestricted TAO-LDA calculations.

Vertical Ionization Potential, Vertical Electron Affinity, and Fundamental Gap

Here, we examine whether n-triangulenes are promising for photovoltaic applications. Calculations based on spin-unrestricted TAO-LDA are carried out at the ground-state (i.e., the lowest singlet state) geometry of n-triangulene to obtain the vertical ionization potentialvertical electron affinityand fundamental gapwhere E(X), E(X+), and E(X–) are the total energies of n-triangulene in the neutral, cationic, and anionic states, respectively.

As presented in Figure , 4, 5, with the increase of system size, IPv decreases monotonically, EAv increases monotonically, and thus, Eg decreases monotonically. It is worth mentioning that the Eg values of n-triangulenes (n = 5, 7, 9, ..., and 19) range from 1 to 3 eV, lying in the ideal regime for photovoltaic applications (see Table S2 in the SI for detailed numerical results).

Figure 3

Vertical ionization potential for the ground state of n-triangulene computed using spin-unrestricted TAO-LDA.

Figure 4

Vertical electron affinity for the ground state of n-triangulene computed using spin-unrestricted TAO-LDA.

Figure 5

Fundamental gap for the ground state of n-triangulene computed using spin-unrestricted TAO-LDA.

Symmetrized von Neumann Entropy

To assess the radical nature of the ground state of n-triangulene, we compute the symmetrized von Neumann entropy[17,19,20,44,49]using spin-unrestricted TAO-LDA. Here, f (i.e., a value between zero and one) is the ith σ-spin (σ = α or β) orbital occupation number obtained with spin-unrestricted TAO-LDA, approximately yielding the ith σ-spin natural orbital occupation number.[15,19,22,24] For a system of nonradical nature, {f} take values in the vicinity of zero or one; hence, the corresponding SvN is very small. However, for a system with significant radical nature, {f} can differ greatly from either zero or one for active spin orbitals (i.e., the spin orbitals that are fractionally occupied) and take the values in the vicinity of zero or one for other spin orbitals; hence, the corresponding SvN can grow rapidly with the number of spin orbitals that are fractionally occupied.

As shown in Figure , SvN grows monotonically with increasing n. Therefore, the larger n-triangulenes are expected to possess increasing polyradical nature in their ground states (see Table S2 in SI for detailed numerical results).

Figure 6

Symmetrized von Neumann entropy for the ground state of n-triangulene computed using spin-unrestricted TAO-LDA.

Active Orbital Occupation Numbers

To understand why SvN grows rapidly with n, we plot the occupation numbers of active orbitals for the ground state of n-triangulene, computed using spin-restricted TAO-LDA. For the ground state of n-triangulene, the (N/2)th orbital is referred to as the HOMO, the (N/2 + 1)th orbital is referred to as the LUMO, and so forth, with N being the number of electrons in n-triangulene.[15,20,22,44,49]

As presented in Figure , with the increase of the side length of n-triangulene, there are growing number of orbitals with occupation number in the vicinity of one (i.e., there are growing number of spin orbitals with occupation number in the vicinity of 0.5). Interestingly, the occupation numbers of active orbitals are found to be very close to one (e.g., between 0.9 and 1.1), indicating that the radical nature of n-triangulenes should be very significant. On the basis of our TAO-LDA results, 3-triangulene exhibits a pronounced diradical nature (e.g., the active orbitals are HOMO and LUMO), 5-triangulene exhibits a pronounced tetraradical nature (e.g., the active orbitals are HOMO – 1, HOMO, LUMO, and LUMO + 1), 7-triangulene exhibits a pronounced hexaradical nature (e.g., the active orbitals are HOMO – 2, HOMO – 1, HOMO, LUMO, LUMO + 1, and LUMO + 2), 9-triangulene exhibits a pronounced octaradical nature (e.g., the active orbitals are HOMO – 3, HOMO – 2, HOMO – 1, HOMO, LUMO, LUMO + 1, LUMO + 2, and LUMO + 3), and so forth. In other words, there are approximately (n – 1) unpaired electrons for the ground state of n-triangulene.[22,33] Accordingly, the larger n-triangulenes should indeed exhibit increasing polyradical nature in their ground states.

Figure 7

Occupation numbers of active orbitals for the ground state of n-triangulene computed using spin-restricted TAO-LDA. Here, for simplicity, HOMO and LUMO are referred to as H and L, respectively.

Similar to our previous TAO-DFT studies on hexagon-shaped graphene nanoflakes[48] and zigzag GNRs,[20] the larger n-triangulenes also possess a polyradical nature in their ground states. However, the electronic properties of n-triangulenes, hexagon-shaped graphene nanoflakes, and zigzag GNRs are distinctly different. For example, n-triangulenes possess a very significant polyradical nature (e.g., the occupation numbers of active spin orbitals are all very close to 0.5), yielding approximately (n – 1) unpaired electrons in their ground states. By contrast, with increasing system size, there is a monotonic transition from the nonradical nature of the smaller hexagon-shaped graphene nanoflakes to the increasing polyradical nature of the larger hexagon-shaped graphene nanoflakes (see Figures 5–9 of ref (48)). Besides, the zigzag GNRs exhibit an oscillatory polyradical nature with increasing ribbon length, and the polyradical nature is also highly dependent on the ribbon width (see Figures 11–13 of ref (20)). On the other hand, the narrowest armchair GNRs (i.e., a series of planar poly(p-phenylene) oligomers) possess a nonradical nature in their ground states,[21] regardless of the ribbon length. These examples clearly support the statement that the radical nature of graphene nanoflakes (e.g., triangle-shaped graphene nanoflakes (n-triangulenes), hexagon-shaped graphene nanoflakes, zigzag GNRs, and armchair GNRs) is intimately correlated with the shapes, edges, and sizes of graphene nanoflakes. Accordingly, it is essential to conduct a comprehensive study on graphene nanoflakes of different shapes, edges, and sizes for guiding the modulation of electronic and magnetic properties.

Visualization of Active Orbitals

For the ground states of some representative n-triangulenes (e.g., n = 3, 5, 7, and 9), we explore the visualization of active orbitals (i.e., the orbitals with an occupation number in the vicinity of one), obtained with spin-restricted TAO-LDA. As presented in Figures , 9, 10, 11, the highest two active orbitals (e.g., HOMO and LUMO for n = 3, LUMO and LUMO + 1 for n = 5, LUMO + 1 and LUMO + 2 for n = 7, LUMO + 2 and LUMO + 3 for n = 9, and so forth) are delocalized over the entire n-triangulene molecule. Besides, the lower the energies of active orbitals, the more concentrated at the periphery of n-triangulene the active orbitals. This indicates that these edge states, when available, are energetically more favorable than others. For larger n-triangulenes, there are more edge states available; hence, there are more unpaired electrons near the edges of n-triangulenes. However, the reason of this fact (i.e., that the active orbitals with lower energies (and hence, with larger occupation numbers) have an increased tendency to localize at the edges of n-triangulenes) may not be obvious; hence, it will be interesting to build a simple model to explain this fact in the near future.

Figure 8

Visualization of the active orbitals (HOMO and LUMO) for the ground state of 3-triangulene computed using spin-restricted TAO-LDA at isovalue = 0.02 e/Å3. Here, for simplicity, HOMO and LUMO are referred to as H and L, respectively. The orbital occupation numbers are shown in parentheses.

Figure 9

Visualization of the active orbitals (HOMO – 1, HOMO, LUMO, and LUMO + 1) for the ground state of 5-triangulene computed using spin-restricted TAO-LDA, at isovalue = 0.02 e/Å3. Here, for simplicity, HOMO and LUMO are referred to as H and L, respectively. The orbital occupation numbers are shown in parentheses.

Figure 10

Visualization of the active orbitals (HOMO – 2, HOMO – 1, HOMO, LUMO, LUMO + 1, and LUMO + 2) for the ground state of 7-triangulene computed using spin-restricted TAO-LDA at isovalue = 0.02 e/Å3. Here, for simplicity, HOMO and LUMO are referred to as H and L, respectively. The orbital occupation numbers are shown in parentheses.

Figure 11

Visualization of the active orbitals (HOMO – 3, HOMO – 2, HOMO – 1, HOMO, LUMO, LUMO + 1, LUMO + 2, and LUMO + 3) for the ground state of 9-triangulene computed using spin-restricted TAO-LDA at isovalue = 0.02 e/Å3. Here, for simplicity, HOMO and LUMO are referred to as H and L, respectively. The orbital occupation numbers are shown in parentheses.

Similar to previous findings for zigzag GNRs,[13,16,17,20] the increasing polyradical nature of the larger n-triangulenes should be closely related to the fact that active orbitals tend to be mainly concentrated at the periphery of n-triangulenes, apparently increasing with the molecular size.

Conclusions

We have investigated the electronic properties, such as EST, IPv, EAv, Eg, SvN, active orbital occupation numbers, and visualization of active orbitals, of n-triangulenes (n = 3, 5, 7, 9, ..., and 21) using TAO-DFT due to its decent balance between cost and performance for studying nanosystems with radical nature. As the larger n-triangulenes exhibit an increasing polyradical nature in their ground states, the properties of n-triangulenes may not be reliably predicted by KS-DFT with the commonly used XC functionals. On the other hand, it is computationally intractable to study the properties of the larger n-triangulenes using a sufficiently accurate multireference electronic structure method. Accordingly, employing TAO-DFT to study the electronic properties of n-triangulenes in this work is certainly justified.

According to our TAO-DFT results, n-triangulenes have singlet ground states for all the values of n studied. Relative to the smaller n-triangulenes, the larger n-triangulenes possess smaller EST values, smaller Eg values, larger SvN values, and more significant polyradical nature. In particular, 3-triangulene exhibits a pronounced diradical nature, 5-triangulene exhibits a pronounced tetraradical nature, 7-triangulene exhibits a pronounced hexaradical nature, 9-triangulene exhibits a pronounced octaradical nature, and so forth. Accordingly, there are approximately (n – 1) unpaired electrons for the ground state of n-triangulene. The increasing polyradical nature of the larger n-triangulenes should be closely related to the fact that the active orbitals tend to be mainly concentrated at the periphery of n-triangulenes, apparently increasing with the molecular size. On the other hand, in view of the approximate nature of the method adopted in this work, to examine the accuracy of our TAO-DFT results, the electronic properties of n-triangulenes obtained with relatively inexpensive and reliably accurate ab initio multireference electronic structure methods are called for.

Since the polyradical nature of n-triangulenes (i.e., triangle-shaped graphene nanoflakes) is closely related to the shapes, edges, and sizes of graphene nanoflakes, several possible directions can be further explored to establish the structure–, edge–, and size–property relationships guiding the modulation of electronic and magnetic properties. For example, if two graphene nanoflakes are covalently bonded together, it will be interesting to examine if the polyradical nature remains or not. Besides, the role of structural defects (e.g., atomic vacancies, pentagon–heptagon defects, adatoms, and dopants) in graphene nanoflakes can be further investigated. Moreover, to passivate the polyradical nature of graphene nanoflakes, we may consider different edge types and chemical termination of the edges in graphene nanoflakes. Furthermore, we may switch from graphene nanoflakes to the corresponding nanoflakes with two or more heteroatoms (e.g., boron nitride nanoflakes) to see if the polyradical nature persists or not. We plan to address some of these questions in the near future.