Nitrogen Doping Enables Covalent-Like π–π Bonding between

For 0D PLY and its π-dimers, the size of the supercell was fixed to 25 Å in all three ... (70) For the single-point RPA calculations, we first repro...
0 downloads 0 Views 2MB Size
Download PDF
Letter pubs.acs.org/NanoLett

Nitrogen Doping Enables Covalent-Like π−π Bonding between Graphenes Yong-Hui Tian,*,† Jingsong Huang,‡,◇ Xiaolan Sheng,† Bobby G. Sumpter,‡,◇ Mina Yoon,‡ and Miklos Kertesz§ †

College of Life Sciences, Research Center of Analytical Instrumentation, Sichuan University, Chengdu, Sichuan 610064, P. R. China Center for Nanophase Materials Sciences, Oak Ridge National Laboratory, Oak Ridge, Tennessee 37831, United States ◇ Computer Science & Mathematics Division, Oak Ridge National Laboratory, Oak Ridge, Tennessee 37831, United States § Department of Chemistry, Georgetown University, 37th and O Streets, NW, Washington, D.C. 20057, United States ‡

S Supporting Information *

ABSTRACT: The neighboring layers in bilayer (and fewlayer) graphenes of both AA and AB stacking motifs are known to be separated at a distance corresponding to van der Waals (vdW) interactions. In this Letter, we present for the first time a new aspect of graphene chemistry in terms of a special chemical bonding between the giant graphene “molecules”. Through rigorous theoretical calculations, we demonstrate that the N-doped graphenes (NGPs) with various doping levels can form an unusual two-dimensional (2D) π−π bonding in bilayer NGPs bringing the neighboring NGPs to significantly reduced interlayer separations. The interlayer binding energies can be enhanced by up to 50% compared to the pristine graphene bilayers that are characterized by only vdW interactions. Such an unusual chemical bonding arises from the π−π overlap across the vdW gap while the individual layers maintain their in-plane π-conjugation and are accordingly planar. The existence of the resulting interlayer covalent-like bonding is corroborated by electronic structure calculations and crystal orbital overlap population (COOP) analyses. In NGP-based graphite with the optimal doping level, the NGP layers are uniformly stacked and the 3D bulk exhibits metallic characteristics both in the in-plane and along the stacking directions. KEYWORDS: graphene, nitrogen doping, interlayer separation, π−π bonding, pancake bonding, van der Waals while maintaining the in-plane π-conjugation of doped graphene framework. It is very intriguing to examine whether a π−π bonding can take place in stacked graphenes. To the best of our knowledge, there have been no discussions on this topic in the graphene literature. Our motivation comes from the unconventional intermolecular two-electron multicenter (2e/mc) π−π bonding, often termed pancake bonding that occurs in a broad range of neutral and charged radical species. These radicals can form πdimers with an intermolecular separation D shorter than the vdW contacts.20−22 Taking a substituted phenalenyl radical molecule (PLY) for instance (Figure S1a,b, Supporting Information), both the intermolecular separation of 3.20 Å23 and the binding energy of −9.50 kcal/mol (−0.41 eV) in the πdimer PLY2 characterize a π−π bonding instead of a pure vdW interaction.24 The π−π bonding interaction in the π-dimers of various substituted PLYs has been identified in several experimental21,24−27 and theoretical studies.28−34 Qualitatively, this form of π-dimerization can be ascribed to the bonding

ating back to the first demonstration of the preparation of graphene (GPH),1 there has been a continual growth in research and a broad interest in GPH and GPH-based materials.2 One important subject in this field is tuning the electronic properties of GPH by adjusting the thickness and stacking of few-layer GPH materials.3 Among them, bilayer GPH materials with AA,4−6 AB (or Bernal),7−9 or rotated stacking patterns10 have been experimentally prepared and studied, and their electronic properties have been broadly reviewed.11,12 Depending on the stacking motifs, GPH bilayers may exhibit many interesting properties such as band gap opening13−15 and electron localization.16,17 A common feature of those GPH bilayers is that the neighboring layers are separated at a distance around 3.4 Å corresponding to van der Waals (vdW) interactions as in the 3D bulk structure of hexagonal graphite. Thereby, the π-electron conjugation and accordingly the planarity in individual GPH layers are preserved despite the possible modification of the Fermi surface as a consequence of stacking effects.18,19 In this Letter, we propose a new aspect of graphene chemistry to interrogate whether a special chemical bonding other than vdW interactions can be created between stacked graphenes through chemical doping

D

© XXXX American Chemical Society

Received: May 18, 2015 Revised: July 6, 2015

A

DOI: 10.1021/acs.nanolett.5b01940 Nano Lett. XXXX, XXX, XXX−XXX

Letter

Nano Letters overlap interactions of the singly occupied π-orbitals of two constituent radical molecules featuring an in-plane π-electron conjugation. Graphene can be viewed as a giant conjugated “molecule” composed of sp2 hybridized carbon atoms.2,35 Of course, unlike the PLY radical, graphene cannot be considered open-shell and therefore all GPH bilayers prepared so far only display interlayer distances in the range of vdW contacts. We will discuss the possible reasons for the absence of π−π bonding in the pristine GPH bilayers and further develop strategies based on chemical doping to introduce unusual covalent-like π−π bonding between stacked layers. Nitrogen doping is an important approach to tune the electronic structure of graphene in order to endow it with potential applications as electromagnetic,36−38 molecular sensing,39 and catalytic materials.40 A variety of synthesis methods have been developed to make N-doped graphenes (NGPs) with different doping levels,41−45 and the electronic and geometric structures of these NGPs have been extensively characterized.46−51 According to the doping configurations, the nitrogen dopants in a graphene can be categorized into three types: pyridinic N, pyrrolic N, and graphitic N, all of which coexist in most of the synthesized NGP materials.48,52 Highly doped NGPs (up to ∼20 at. %) have been prepared in which the nitrogen atoms are mainly in the pyridinic and pyrrolic forms,42,53 whereas NGPs with predominantly graphitic Ndopants have been achieved only at relatively low doping concentrations (below 5 at. %).42,46,54,55 Nevertheless, theoretical investigations predict that doping concentration as high as 33.3−37.5 at. % is achievable for graphitic nitrogensubstituted graphene.56 Graphitic N-doping is advantageous over the other two types of doping configurations in that the πconjugation and planarity, together with lattice parameters, remain essentially the same as those of pristine graphene due to the similar C−C and C−N bond lengths. Herein, we demonstrate for the first time that at proper doping levels, the graphitic N-doping can induce covalent-like π−π bonding (pancake bonding) in NGP bilayers. The formation of π−π bonding is directly supported by short interlayer separations and stronger binding energies, and further by electronic structure calculations and crystal orbital overlap population (COOP) analyses. Because of the interlayer π−π bonding, an electrical conduction pathway is formed in the solid state along the stacking direction. This finding may open a way to assemble 2D NGPs into novel layered materials with distinctive properties associated with the novel covalent-like π−π bonding. Computational Methods. First-principles calculations were performed mainly with the Vienna Ab Initio Simulation Package (VASP, version 5.3.5).57 For density functional theory (DFT) calculations, the interactions between ions and electrons were described using the projector-augmented wave (PAW) method,58 and the exchange and correlation interactions of valence electrons were described with the Perdew−Burke− Ernzerhof (PBE) functional.59 The corrections for vdW interactions were addressed by four empirical schemes including Grimme’s D2 and D3,60,61 and Tkatchenko−Scheffler methods without (TS) and with self-consistent screening (TS +SCS).62−64 These methods are denoted as PBE-D2, PBE-D3, PBE-TS, and PBE-TS+SCS, respectively. In addition to these empirical corrections, five versions of nonlocal vdW-DF functionals of Langreth and Lundqvist and co-workers were also employed, ranging from the originally proposed method using revPBE exchange functional,65 the improved versions with optimized exchange functionals including optPBE, optB88,

and optB86b,66 and the latest vdW-DF2 method using the rPW86 exchange functional.67 Based on a careful set of validations of these vdW methods (Table S1, Supporting Information), the semiempirical PBE-TS and PBE-D3 methods, together with the self-consistent optB88 vdW-DF method, were selected for further computations. Specifically, the PBE-TS was used to optimize the structures (giving comparable results as the optB88 vdW-DF), whereas the PBE-D3 was used to calculate the single-point binding energies. This approach is denoted as PBE-D3//PBE-TS. As a further verification, binding energies were calculated for the pristine GPH bilayers by the adiabatic-connection fluctuation−dissipation theorem in the random phase approximation (ACFDT-RPA)68 using the PBEbased Kohn−Sham orbitals. For a few selected systems, we also performed PBE-based calculations including an accurate description of nonlocal many-body dispersion (MBD) interactions,63 which is denoted as PBE-MBD, using the FHIaims all-electron code.69 Both the RPA and MBD binding energies were also obtained from single-point calculations using the PBE-TS optimized geometries. Three types of structures were investigated including 3D bulk structures of hexagonal graphite in an AB stacking and of NGP-based graphite in an AA stacking, 2D single- and bilayer structures of NGPs, and 0D unsubstituted PLY C13H9 and its πdimers in both AA (eclipsed) and AB (staggered) conformations (Figure S1c,d, Supporting Information). For DFT calculations, the Brillouin-zone integrations were performed on Γ-centered 6 × 6 × 4, 6 × 6 × 1, and 1 × 1 × 1 k-point grids for 3D, 2D, and 0D structures, respectively. The kinetic energy cutoff for plane waves was set to 1000 eV and the convergence criterion for electronic self-consistency was set to 10−6 eV. During the structure optimizations of 2D single- and bilayer structures of NGPs, the vacuum regions were fixed to 25 Å to ensure the periodic images are well separated, whereas other lattice vectors were fully relaxed. For 0D PLY and its π-dimers, the size of the supercell was fixed to 25 Å in all three directions. All atoms were relaxed until the Hellmann−Feynman forces were smaller than 0.001 eV/Å unless otherwise noted. All of the optimized structures were further checked by spin-polarized DFT giving zero magnetization. The k-point grids, the kinetic energy cutoff, and the vacuum regions employed were based on systematic convergence studies (Tables S2 and S3, Supporting Information). The interlayer binding energies Ebind for bilayer and bulk NGPs were calculated by Ebind = (E0 − 2 × E0′)/N, where E0 and E0′ are the total free energy of the unit cell for the bilayer (or bulk) and monolayer, respectively, and N is the total number of carbon and nitrogen atoms. In the case of PLY πdimers, the intermolecular binding energies were calculated likewise but with the normalization based on only the 13 × 2 carbon atoms not including the hydrogen atoms. For the singlepoint PBE-MBD calculations, the vacuum size was increased to 400 Å. The k-point grids of 12 × 12 × 1 and 12 × 12 × 4 were used for the 2D systems and the pristine or NGP-based graphite, respectively. These conditions were carefully determined based on convergence studies, allowing us to reproduce the binding energies of −23 and −48 meV/C for AB-stacked GPH bilayer and hexagonal graphite, respectively, as reported in the literature.70 For the single-point RPA calculations, we first reproduced the interlayer binding energy of −48 meV/C for hexagonal graphite using the conditions reported in the literature,71 and then calculated the interlayer binding energies for the pristine GPH bilayers. The single-point COOP calculations were performed with the SIESTA packB

DOI: 10.1021/acs.nanolett.5b01940 Nano Lett. XXXX, XXX, XXX−XXX

Letter

Nano Letters age,72 where the PBE functional and Troullier−Martins pseudopotentials were used in combination with localized basis of pseudoatomic orbitals with double-ζ plus polarization functions.73 The energy cutoff for the real space grid is 400 Ry, whereas the k-point grids are the same as those given above. Validation of the Computational Methods. Table 1 shows the results of validation using a few selected vdW

nitrogen atoms. These structures are denoted as NGP-n, where n refers to the number of nitrogen atoms in the unit cells. Specifically, n takes the values of 0, 1, 2, 3, 4, 6, and 8 corresponding to doping concentrations in the range of 0−25 at. % where NGP-0 corresponds to the pristine GPH. Figure 1a illustrates the case of an NGP-4 monolayer. We note that at a given doping concentration a plethora of structures are possible depending on the doping sites.77 In this Letter, we do not attempt to cover all of the possible doping configurations but only studied the simplest scenario for each concentration where the nitrogen atoms are well separated in order to reach an even distribution of the doping atoms (Figure 1a and Figure S2, Supporting Information). Upon nitrogen substitution, the inplane covalent bond lengths were only slightly affected, and the hexagonal lattice parameters remained essentially the same as in the pristine GPH (Table S4, Supporting Information). For each NGP-n, bilayers are constructed by stacking NGP-n in an AA or AB motif, giving supercells with 64 atoms per cell. The cases for AA- and AB-stacked NGP-4 bilayers are illustrated in Figure 1b,c, respectively. Interlayer Separations and Binding Energies in NGP Bilayers. In Figure 2, the optimized interlayer distances and calculated binding energies were plotted as a function of the doping levels for NGP-n bilayers in both AA- and AB-stackings. For the pristine GPH bilayers in both AA- and AB-stackings, depending on the theoretical methods used for calculations, the interlayer distances (D) in the literature vary in a wide range.78−81 According to our calculations for the AB-stacked GPH bilayer, both PBE-TS and optB88-vdW methods predict D value 3.37 Å in good agreement with the experimental value of 3.4 Å.6 For the pristine AA-stacked GPH bilayer, the PBE-TS and optB88-vdW methods predicted D values of 3.51 and 3.55 Å, respectively, close to the observed interlayer spacing of 3.55 Å in the bulk structure of AA graphite.82 Figure 2a shows a consistent trend from the two different theories that the intermolecular separation D for the AB-stacking motif is slightly and monotonically shortened with increasing doping level, whereas that for the AA-stacking shows a significant variation as a function of the doping level. Figure 2b shows essentially the same trend of interlayer binding energies as the interlayer separations. The binding energies calculated by PBE-TS and optB88-vdW are close to each other, whereas those from the single-point PBE-D3 calculations have smaller magnitude. Here, we will focus our attention for the binding energies on the PBE-D3 curves, as the results for GPH bilayers were confirmed by the MBD and high-level RPA calculations as indicated by the isolated markers for GPH, NGP-4, and NGP-8, for which the 4a × 4b supercells can be reduced to 2a × 2b cells by symmetry. In comparison, the RPA calculations are prohibitively expensive and therefore can only be applied to

Table 1. Interlayer Separations D and Binding Energies Ebind Normalized to Per Carbon Atom of Hexagonal Graphite Calculated by Various Methods and Compared with Experiments methods

D (Å)a

Ebind (meV/atom)b

PBE-TS optB88-vdW PBE-D3 PBE-MBD RPA experiment

3.35 3.35 3.48 3.35e 3.34g 3.35h

−83c −70c −48c,d −48d,f −48d,g −43i, −52j

a

The interlayer distances between carbons that are in direct contact. Ebind = (E0 − 2 × E0′)/Nc, where E0 and E0′ are the total free energies of graphite and single-layer GPH, and Nc = 4 is the number of carbon atoms per unit cell. cBinding energies based on the structures optimized at the respective theory levels. dSingle-point binding energies based on the structures of graphite and GPH optimized by PBE-TS. ePBE-TS geometry. fResult from ref 70 and reproduced in this work. gResult from ref 71 and reproduced in this work. hRef 74. i Ref 75. jRef 76. b

methods against the available experimental data for hexagonal graphite (the full list of validation results are provided in Table S1, Supporting Information). It can be seen that PBE-TS and optB88-vdW perform the best with respect to the experimental interlayer separation of D = 3.35 Å,74 whereas PBE-D3, RPA, and PBE-MBD based on the PBE-TS geometries yield consistent binding energies that fit well in the experimental range.75,76 As will be further shown below, the binding energies obtained by the PBE-D3 method agree well with the RPA and MBD results for the pristine GPH bilayers and selected NGP bilayers. All of these results indicate that PBE-TS and optB88vdW should be suited for geometry optimizations giving reasonable interlayer separations while the PBE-D3 single-point energy calculations based on the PBE-TS geometries may provide trustworthy interlayer binding energies for the NGP materials. Structure Models of N-Doped Graphenes (NGPs). Single-layer NGPs with different doping levels were constructed by starting from the 4a × 4b supercell of the pristine GPH consisting of 32 atoms and by varying the number of

Figure 1. Structural models for N-doped graphene (NGP). (a) NGP-4 monolayer consisting of 28 carbon and 4 nitrogen atoms per unit cell (see Figure S2 of the Supporting Information for the entire set of monolayer structures). (b) and (c) Perspective views of the AA- and AB-stacked NGP-4 bilayers. C

DOI: 10.1021/acs.nanolett.5b01940 Nano Lett. XXXX, XXX, XXX−XXX

Letter

Nano Letters

Figure 2. Interlayer distances D (a) and binding energies per atom (b) as a function of doping levels for NGP bilayers in both AA and AB stackings. The geometries and the respective energies were calculated by the PBE-TS (TS) and optB88-vdW (BO) methods. The binding energies for GPH bilayers were also predicted at the PBE-D3 (D3), PBE-MBD (MBD), and ACFDT-RPA (RPA) level based on the TS geometries. The interlayer separations were averaged over all of the atom pairs in direct contact. For AB-stacked bilayers of NGP-3 (9.4 at. %) and NGP-4 (12.5 at. %) as denoted by the hollow markers and dashed line segments, the convergence of structural relaxations was reached with the Hellmann−Feynman forces set to a threshold of 0.01 eV/Å.

the GPH bilayers for which the 2a × 2b supercells can be further reduced to the 1a × 1b primitive cells by symmetry. As can be seen from Figure 2a,b, at both small and high doping levels, the bilayer NGPs including GPH, NGP-1, NGP2, and NGP-8 in the AA-stacking motif all have larger interlayer spacing and accordingly weaker binding than the AB-stacking motif. However, the trend is inverted for intermediate doping levels, thus creating two crossovers. It is noteworthy that other available vdW methods predicted similar trends in terms of interlayer separations and binding energies for the entire set of bilayer NGPs (Figure S3, Supporting Information). For NGP-4 with an optimal doping level of 12.5 at. %, the D of the AAstacking motif is reduced to the shortest value of 3.1−3.2 Å depending on the applied theory, significantly shorter than those of the AA-stacked (3.51−3.55 Å) and AB-stacked (3.37 Å) GPH bilayers. Correspondingly, its binding energy is enhanced by ca. 50% (corresponding to a change from −19 meV/atom for GPH to −30 meV/atom for NGP-4, both in AA-stacks). This change amounts to a −11 meV/atom enhancement, or −11 × 64 = −704 meV of additional binding energy per unit cell. Despite the good performance of pairwise PBE-D3 method for graphene,61 we are aware of the fact the inclusion of a repulsive three-body (Axilrod−Teller−Muto) term would underestimate the dispersion interactions in the absence of higher order dispersion corrections.83,84 However, the D3 and MBD methods predict similar binding energies for the bilayer GPH, NGP-4, and NGP-8. These results indicate that significant binding enhancement occurs at the doping level of 12.5 at. % corresponding to NGP-4. Next, we discuss Figure 2 in more detail regarding the origins for this unexpected behavior. At relatively low doping levels for NGP-1 (3.1 at. %) and NGP-2 (6.2 at. %), the AB-stacked bilayers are very close to the pristine GPH bilayer in terms of both interlayer spacing and binding energies, which implies predominant vdW interactions in the N-doped AB-bilayers. Meanwhile, similar to the case of pristine GPH, AB-stacks are energetically favored over the AA-stacks, in good agreement with the experimentally observed AB-stacked NGP bilayers at the relatively low doping concentration around 4 at. %.42 On the other hand, as shown in Figure 2, the NGP-1 and NGP-2 bilayers in the AA stacking motif start to exhibit shorter interlayer separations and meanwhile stronger interlayer

bonding than the pristine GPH bilayer with the same packing motif. Both these two indices, although still not significant, already suggest the presence of some special bonding other than merely vdW interactions in the N-doped AA-stacks. It is plausible to assume π−π bonding formation in the two NGP bilayers, as good interlayer orbital overlap is accessible in this eclipsed AA-stacking motif. More solid evidence for the π−π bonding formation is shown below. At higher doping levels in NGP-3 and NGP-4, the interlayer spacings for the AA-stacked bilayers were significantly shortened, whereas their binding energies were significantly enhanced when compared with the pristine GPH or doped NGP bilayers at lower doping levels. At NGP-4, the binding strength reached the maximum corresponding to the shortest interlayer separation among the set of N-doped AA-bilayers. However, for the AB-stacked motifs of both NGP-3 and NGP4, the interlayer bindings are still close to the pristine GPH bilayer. These differences clearly indicate the existence of a bonding effect in the NGP-3 and NGP-4 AA-bilayers that is different from the vdW interaction in the AB-bilayers. Meanwhile, the AA-stacked NGP-3 and NGP-4 bilayers became more stable than the corresponding AB stacks. This is in contrast to the preference of an AB-stacking in the pristine GPH or NGP-1 and NGP-2 bilayers as a result of the reduced exchange repulsion in the AB-stacks.85 These observations imply the existence of a competition between π−π bonding and reduced exchange repulsion favoring AA and AB stacking configurations, respectively. It is noteworthy that the preference of AB- over AA-stacked configuration has also been observed for PLY π-dimers in the solid state,23 in line with theoretical calculations.33,86 For the bilayer NGP-3 and NGP-4, we suggest that the preference of AA stacking could be an indication that the predominant π−π bonding effect outweighs the exchange repulsion. Starting with bilayer NGP-6, the interlayer spacing of the AA stacking motif begins to increase, whereas its binding energy starts to be weakened. However, it can be seen in Figure 2 that the π−π bonding still exists, as reflected by the slightly preferred AA-stacking motif. Further increase of the doping level to NGP-8 results in the AA-stacked bilayer again having a longer interlayer spacing than the AB motif. Accordingly, it became energetically less favored over the AB-stack, similar to D

DOI: 10.1021/acs.nanolett.5b01940 Nano Lett. XXXX, XXX, XXX−XXX

Letter

Nano Letters

Electronic Structures of NGP bilayers. In Figure 4, we compare the energy band structure and DOS for the single- and bilayer sheets of the pristine GPH and doped NGP-4. For the pristine GPH, the π and π* bands from two identical monolayers (Figure 4a) split into πA ∓ πA′ and π*A ∓ π*A′ bands (Figure 4b) upon AA-stacking, where the negative and the positive signs stand for bonding and antibonding combinations, respectively, and the labels A and A′ refer to the two GPH layers. The energy level crossing between πA + πA′ and π*A − π*A′ results in a Fermi surface in the shape of a circle with a small radius around the K point in the Brillouin zone. As a consequence of the small DOS at Ef, the π−π bonding effect is negligible, as reflected by the long interlayer distances corresponding to vdW interactions (Figure 2). To verify the π−π bonding formation in the N-doped systems, we examined the electronic structure for the AA-stacked NGP-4 bilayer that shows the strongest interlayer binding among the set of NGPs. Upon doping, an energy gap of ca. 0.7 eV at the K point is created in the NGP-4 monolayer (Figure 4c), which is consistent with the previous calculations for N-doped graphene.87 Meanwhile, Ef is raised up due to the extra electrons brought into the system by N. Consequently, the initially empty bands in GPH become partially filled until Ef, which sits right on a peak of the DOS. The considerable DOS around Ef facilitates to create an effective π−π bonding between two stacked layers. As exemplified by the band structure in Figure 4d for the NGP-4 bilayer, two of the π* bands near the Fermi level in monolayers split pairwise into parallel bonding and antibonding bands that are separated by energy gaps about 0.8 and 1.0 eV for π*1 and π*2 bands, respectively. Notably along K−Γ−M−K (Figure 4d), a few π* bands are involved in the energy stabilization around Ef and hence in the formation of π−π bonding, which is consistent with the significant bonding energies and short bonding distances for the AA-stacked NGP4 bilayer. The reason for NGP-4 instead of other doping levels to display the strongest π−π bonding should be ascribed to the fact that when n = 4, Ef is raised to the highest DOS peak (Figure 4c), as verified by a comparison across different doping levels (Figure S4, Supporting Information). From another aspect, the preference of AA- over AB-stacking for NGP-4 is the outcome of the same orbital symmetries and efficient orbital overlaps between the two layers ensured by the face-to-face stacking favoring π−π bonding formation. Crystal Orbital Overlap Population Analysis (COOP) Analysis. COOP analyses can provide helpful information about chemical bonding in extended systems.88 A relevant index is the integrated COOP (INTCOOP). According to the COOP scheme, a positive value of INTCOOP up to Ef is indicative of overall bonding orbital overlap of two chemical fragments, whereas a negative value points toward antibonding interactions. In addition, the magnitude of the INTCOOPs is qualitatively proportional to the bonding strength. For the NGP bilayers under studies, the π−π bonding formation is the outcome of pz orbital overlap along the stacking directions, and hence, only the pz orbital contribution to the INTCOOPs is of interest here. As can be seen from Figure 5a, for the entire set of AAstacked NGP-n bilayers, the INTCOOPs up to Ef are found to increase with increasing doping levels until NGP-4 and then start to decrease thereafter. This trend appears to be a mirror image of the dependence of interlayer separations D on the doping levels. The maximum of INTCOOP and the minimum of D are consistent with the maximum binding energies at the

their low doping analogues (NGP-0, NGP-1, and NPG-2). For all of the AB-stacked NGP bilayers, the interlayer binding appears to be only vdW interactions in nature in the entire range of doping level of 0−25 at. %. The slight and monotonic shortening of interlayer separations may be ascribed to the smaller vdW radius for N than C in combination with the gradually increased doping levels. Qualitative Description of the π−π Bonding Formation in 2D Systems. Before presenting the results of the electronic structure calculations, we illustrate the necessary conditions for the π−π bonding formation in 2D-bilayers by using the schematic energy diagram shown in Figure 3. For a

Figure 3. Schematic diagram showing 2D covalent-like π−π bond formation. (a) Energy level splittings around the Fermi level Ef (dashed line) lead to level crossings (red circles) as a result of the interlayer π−π orbital overlaps. The occupied and empty levels in the monomers (in this case the monolayer) are artificially separated by a gap to guide the eye. (b) Relative positions of Ef between the pristine and the N-doped graphenes. Ef is shifted from the Dirac point where the density of states (DOS) is low to a place with enhanced DOS, making it possible for a significant π−π bonding effect.

given system with small or zero band gaps, the near energy degeneracy of the occupied and the empty energy levels around the Fermi level, Ef, would be lifted due to interlayer orbital overlaps when two “monomers” are stacked together (Figure 3a). In a certain energy range around Ef, the unoccupied states in the “monomers” would become occupied in the “dimers” in the form of interlayer π−π bondings as a result of an effective in-phase orbital overlap between two stacked layers. On the other hand, due to an out-of-phase orbital overlap, some of the occupied states near Ef in the “monomers” would become empty in the “dimers” in the form of interlayer π−π antibonding interactions. The net outcome from this interlayer orbital overlap is the level crossing, as highlighted by the red circles in Figure 3a. Therefore, for 2D extended systems, the small or zero band gap in conjunction with the large density of states (DOS) around Ef are essential for the π−π bonding formation when two identical layers are stacked in an eclipsed fashion. For pristine GPH, although the band gap is zero at the Dirac points (Figure 3b left and Figure 3a where a gap is deliberately created for illustration), the π−π bonding effect is negligible due to the small DOS near Ef. In contrast, for Ndoped graphene, the energy level crossing would occur at the Fermi surface that is rich in electronic states (Figure 3b right), and thereby significant amount of electrons would be dumped from the originally occupied states in the “monomers” to the interlayer π−π bonding states in the “dimers”, half of which are originally occupied while the other half are originally empty in the “monomers”. Accordingly, a strong π−π bonding is expected in the NGP bilayers at a proper doping level. This pictorial description will be exemplified by the following electronic structure calculations. E

DOI: 10.1021/acs.nanolett.5b01940 Nano Lett. XXXX, XXX, XXX−XXX

Letter

Nano Letters

Figure 4. Band structure and DOS of GPH and NGP-4. (a) Pristine GPH monolayer in a 4a × 4b supercell consisting of 32 atoms. (b) AA-stacked GPH bilayer in a 4a × 4b supercell consisting of 64 atoms. (c) Single layer of NGP-4 with the structure shown in Figure 1a. (d) AA-stacked NGP-4 bilayer with the structure shown in Figure 1b. All band structures were calculated by PBE-TS based on the geometries optimized by the same theory. In the bilayer band structures, some of the π* bands near the Fermi level crucial for π−π bonding formation are labeled (see the text).

Next, we show further evidence of π−π bonding formations from the COOP curve together with integrated density of states (INTDOS). As displayed in Figure 5b, for the AA-stacked NGP-4 bilayer, the COOP curve shows an oscillating feature. At the critical point Ec as marked by the red circle where the COOP values change from negative to positive, the INTCOOP reaches a local minimum (0.027). Beyond this point until Ef, the COOPs are positive indicating the π−π bonding states in the energy range from Ec to Ef. This point of view is corroborated by INTDOS represented by the black dashed line in Figure 5b. The INTDOS up to Ec corresponds to 256 electrons amounting to the total number of electrons in each unit cell for the pristine GPH bilayer. Therefore, the extra 8 electrons from the eight N-dopants in the AA-stacked NGP-4 bilayer are responsible for the π−π bonding formation. Additional evidence supporting this argument can be obtained by comparing the COOP curves for the AA-bilayers of NGP-4 and GPH. Noticeably, for a GPH bilayer with D intentionally reduced to 3.18 Å, its INTCOOP up to Ef is 0.028 (Figure 5c), nearly equal to the value of 0.027 for NGP-4 up to the energy level of Ec. In addition, the INTDOS corresponds to 256 valence electrons for both the GPH bilayer (up to Ef) and the NGP-4 bilayer (up to Ec). Therefore, the nearly identical

doping level of NGP-4. Therefore, the COOP analysis confirms the occurrence of π−π bonding interactions in the AA-stacked NGP bilayers. In contrast, for the AB-stacked NGP-n bilayers, no obvious INTCOOP changes were observed with the change of doping concentrations (Figure S5, Supporting Information), consistent with the predominant vdW interactions in the Ndoped AB-bilayers. To further understand the π−π bonding nature in the AA-stacked NGP bilayers, a comparison is directly made between the AA-stacked pristine GPH and NGP-4 bilayers. As shown in Figure 5a, for the AA-stacked pristine GPH bilayer at the equilibrium distance of D = 3.51 Å, the INTCOOP of 0.037 is significantly smaller than 0.115 for the AA-stacked NGP-4 bilayer with an interlayer separation of D = 3.18 Å. A critical question is whether the much smaller INTCOOP of the AA-stacked GPH bilayer is mainly caused by the lack of N-doping or it could be directly reduced by its large interlayer separation. As can be seen from Figure 5a, even when the interlayer separation for the AA-stacked GPH bilayer is reduced from D = 3.51 to 3.18 Å, the INTCOOP is only slightly changed from the initial 0.037 to the final 0.028 reflecting the nonbonding nature between GPH layers. The subtle change clearly demonstrates the effect of N-doping on the π−π bonding formation in NGP bilayers. F

DOI: 10.1021/acs.nanolett.5b01940 Nano Lett. XXXX, XXX, XXX−XXX

Letter

Nano Letters

Figure 5. Crystal orbital overlap population (COOP) analyses based on the pz orbital contributions for the AA-stacked NGP bilayers and PLY πdimer. (a) Integrated COOP up to Ef and interlayer separation D as a function of doping levels. The INTCOOP values are also shown for the AAstacked PLY π-dimer at D = 3.33 Å and the pristine GPH bilayer at D = 3.18 Å. (b) COOP, INTCOOP, and INTDOS of the AA-stacked NGP-4 bilayer. The red circle indicates the critical point Ec (see the text for description). (c) COOP and INTCOOP of the AA-stacked pristine GPH bilayer calculated for D = 3.18 Å. (d) COOP and INTCOOP of the AA-stacked PLY π-dimer at the equilibrium D = 3.33 Å. For comparability with the NGP bilayers consisting of 32 atoms pairs per unit cell, the COOP and INTCOOP for the AA-stacked PLY π-dimer containing 13 carbon pairs are scaled by a factor of 32/13.

INTCOOP implies that the binding in the AA-stacked NGP-4 bilayer would be basically vdW interactions without the extra 8 electrons from the N-dopants. For comparative purpose, the COOP analysis was also carried out for the AA-stacked PLY π-dimer, which serves as a prototype of the intermolecular π−π bonding. As can be seen from Figure 5d, an oscillating feature is also observed in its COOP curve. Similarly to the NGP-4 bilayer, the top peak right below the Femi level is a positive COOP corresponding to the π−π bonding involving two radical electrons. This feature suggests that the bonding in NGP-4 bilayer is analogous in nature to the prototypical π−π bonding in the PLY2 π-dimer. Furthermore, the AA-stacked NGP-4 bilayer and the PLY2 πdimer show comparable INTCOOP values of 0.115 (Figure 5a) and 0.093 (Figure 5a,d). Note that the INTCOOP for the PLY2 π-dimer consisting of 13 carbon pairs was scaled by a factor of 32/13 to ensure comparability with the NGP-4 bilayer containing 32 atom pairs per unit cell. Thus, the high INTCOOP value for the AA-stacked NGP-4 bilayer is solid evidence for the novel π−π bonding discussed in this Letter. Three-Dimensional Material Based on NGP-4. The unit cell used for the electronic structure calculations for a 3D bulk material contains two NGP-4 layers in order to investigate a possible symmetry breaking between the layers by allowing the interlayer distance to alternate between a larger and a smaller value. However, this symmetry breaking was not found to lower

the total energy. As shown in Figure S6 of the Supporting Information, the NGP-4 layers are uniformly stacked with an equilibrium interlayer spacing of D = 3.24 Å, which is slightly longer than that of the isolated AA-stacked bilayer (D = 3.18 Å). In the solid state, the binding energy is found to be −60 meV/atom at the PBE-D3//PBE-TS level, which is twice the value of −30 meV/atom for the isolated bilayer. The enhanced binding energy in the bulk can be attributed to the formation of two equivalent π−π bondings for each individual layer, one on each side. In addition, the binding energy of −60 meV/atom for the NGP-4-based graphite is significantly enhanced compared to that of −48 meV/atom for the pristine hexagonal graphite, in line with the shortening of interlayer separations from 3.35 to 3.24 Å. These binding energies at the PBE-D3 level are further confirmed by PBE-MBD calculations giving −63 and −48 meV/atom, respectively. The π−π bonding between layers makes it possible to form an electrical conducting pathway along the stacking z direction. As can be seen from the band structure in Figure 6 for the AAstacked NGP-4 bulk material, a pair of back-folded bands that flank the Fermi level can be found along the K−H or Γ−A line corresponding to the stacking direction. The back-folding results from the two equivalent layers in the unit cell. By comparing the bands of the single-layer NGP-4 (Figure 4c) and its bulk at the Γ point (Figure 6), this pair of bands originated from the π*2 band. From the large bandwidth W of 1.84 eV at G

DOI: 10.1021/acs.nanolett.5b01940 Nano Lett. XXXX, XXX, XXX−XXX

Letter

Nano Letters

expected that the π−π pancake bonding may cause changes to spectroscopic signatures such as the interlayer shear and breathing modes of the bilayer structures. On the basis of the interlayer π−π bonding, the NGP-4 bulk material can in principle be obtained with an interlayer separation of 3.24 Å. For this solid-state NGP-4 material, the predicted binding energy is twice that of the isolated AA-stacked bilayer, indicating a strong stabilization effect arising from π−π bonding delocalization between the uniformly packed layers. The strong π−π bonding interactions also lead to dispersive π* bands crossing the Fermi level suggesting that the NGP-4 bulk material should exhibit metallic characteristics both in the graphene plane and in the normal direction.

■

ASSOCIATED CONTENT

* Supporting Information S

Figure 6. Band structure and DOS for the AA-stacked NGP-4 bulk material, C7N, calculated at the PBE level following the PBE-TS optimization. The special reciprocal k-points are K(2/3, 1/3, 0), Γ(0, 0, 0), M(1/2, 1/2, 0), H(2/3, 1/3, 1/2), A(0, 0, 1/2). The band splitting along K−Γ−M−K or the bandwidth (W) along K−H is twice that of the isolated AA-stacked bilayers as indicated in Figure 4d. The π1* and π2* bands around the Fermi level were labeled for illustration.

Validations and convergence studies of computational methods, comparison of C−C and C−N bond lengths of single-layer GPH and NGPs, structures of PLY and its π-dimers, optimized structure models and DOS of single-layer NGPs, interlayer separations and binding energies as a function of doping levels calculated by additional vdW-DF methods, COOP analyses and interlayer separations of the AB-stacked NGP bilayers, and optimized structure model of the 3D bulk material of NGP-4. The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.nanolett.5b01940.

the Γ point, we estimate the transfer integral to be t = 0.46 eV according to the approximate relationship of W = 4t for the bulk. In comparison, for the AA-stacked NGP-4 bilayer, the corresponding transfer integral is estimated to be t′ = 0.51 eV based on the dimer level splitting at the Γ point according to W′ = 2t′ for the bilayer structure. The slight difference between t for the bulk and t′ for the bilayer is simply due to the difference of the interlayer distance. The rather large interlayer transfer integral t serves to provide another quantification of the strong π−π bonding formation between NGP-4 layers. As can be further seen from Figure 6, along the stacking direction, the zero band gap and the large DOS at the Fermi level suggest high electrical conductivity in the c direction in the NGP-4 bulk material, in sharp contrast to hexagonal graphite. Meanwhile, the band structure along K−Γ−M−K is basically the same as the NGP-4 monolayer (Figure 4c) except for the energy level splitting as a result of π−π bonding effect. Due to the n-type doping, solid state NGP-4 material also shows metallic characteristics in the plane as reflected by the partially filled bands around the Fermi level. Consequently, this solidstate material is predicted to be a 3D metallic conductor of a rather low density of 2.4 g/cm3 with conducting pathways both in the graphene plane and in the normal direction. Conclusions. In summary, we have demonstrated that delocalized 2D π−π bonding can be realized with N-doped graphenes (NGPs). Compared to pristine graphene, the interlayer binding energies for the AA-stacked NGP bilayers are enhanced by up to 50%, whereas the interlayer separations between neighboring NGP layers are significantly reduced to as short as 3.18 Å, clearly indicating covalent-like π−π bonding formations and not simply vdW interactions. The strongest binding interaction occurs in the AA-stacked NGP-4 bilayer corresponding to a doping concentration of 12.5 at. %. Our electronic structure calculations and crystal orbital overlap population (COOP) analyses indicate that the enhanced bindings upon N-doping are attributed to the π−π pancake bonding formations. We suggest that the absence of π−π bonding in the pristine graphene bilayers is due to the negligible DOS near the Fermi level, whereas the significant DOS for the NGPs makes the π−π bonding possible. It can be

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS This research was supported by the National Science Foundation of China (Grant No. 21443012) and the Faculty Startup Grant of Sichuan University, by the Center for Nanophase Materials Sciences, a U.S. DOE Office of Science user facility, and by the U.S. National Science Foundation support at Georgetown University (Grant No. CHE-1006702). Y.T. and X.S. thank the National Supercomputing Center in Shenzhen for providing the computational resources and HPC engineers Boyang Li and Zhifang Song at Beijing Paratera Technology Co., Ltd. for their technical support. J.H., B.G.S., and M.Y. are indebted to the computational resource of the National Energy Research Scientific Computing Center, which is supported by the U.S. DOE Office of Science under Contract No. DE-AC02-05CH11231. M.K. is member of the Georgetown Institute of Soft Matter.

■

REFERENCES

(1) Novoselov, K. S.; Geim, A. K.; Morozov, S. V.; Jiang, D.; Zhang, Y.; Dubonos, S. V.; Grigorieva, I. V.; Firsov, A. A. Science 2004, 306, 666−669. (2) Geim, A. K. Science 2009, 324, 1530−1534. (3) Latil, S.; Henrard, L. Phys. Rev. Lett. 2006, 97, 036803. (4) Liu, Z.; Suenaga, K.; Harris, P. J.F.; Iijima, S. Phys. Rev. Lett. 2009, 102, 015501. (5) Papagno, M.; Pacilé, D.; Topwal, D.; Moras, P.; Sheverdyaeva, P. M.; Natterer, F. D.; Lehnert, A.; Rusponi, S.; Dubout, Q.; Calleja, F. ACS Nano 2012, 6, 9299−9304. H

DOI: 10.1021/acs.nanolett.5b01940 Nano Lett. XXXX, XXX, XXX−XXX

Letter

Nano Letters (6) Borysiuk, J.; Soltys, J.; Piechota, J. J. Appl. Phys. 2011, 109, 093523. (7) Feldman, B. E.; Martin, J.; Yacoby, A. Nat. Phys. 2009, 5, 889− 893. (8) Mayorov, A. S.; Elias, D. C.; Mucha-Kruczynski, M.; Gorbachev, R. V.; Tudorovskiy, T.; Zhukov, A.; Morozov, S. V.; Katsnelson, M. I.; Fal’ko, V. I.; Geim, A. K.; Novoselov, K. S. Science 2011, 333, 860− 863. (9) Yan, Z.; Peng, Z.; Sun, Z.; Yao, J.; Zhu, Y.; Liu, Z.; Ajayan, P. M.; Tour, J. M. ACS Nano 2011, 5, 8187−8192. (10) Hass, J.; Varchon, F.; Millan-Otoya, J. E.; Sprinkle, M.; Sharma, N.; de Heer, W. A.; Berger, C.; First, P. N.; Magaud, L.; Conrad, E. H. Phys. Rev. Lett. 2008, 100, 125504. (11) McCann, E.; Koshino, M. Rep. Prog. Phys. 2013, 76, 056503− 056531. (12) Castro Neto, A. H.; Guinea, F.; Peres, N. M. R.; Novoselov, K. S.; Geim, A. K. Rev. Mod. Phys. 2009, 81, 109−162. (13) Ohta, T.; Bostwick, A.; Seyller, T.; Horn, K.; Rotenberg, E. Science 2006, 313, 951−953. (14) Sboychakov, A.; Rakhmanov, A. L.; Rozhkov, A. V.; Nori, F. Phys. Rev. B: Condens. Matter Mater. Phys. 2013, 87, 121401. (15) Rutter, G. M.; Jung, S.; Klimov, N. N.; Newell, D. B.; Zhitenev, N. B.; Stroscio, J. A. Nat. Phys. 2011, 7, 649−655. (16) de Laissardiere, G. T.; Mayou, D.; Maguad, L. Nano Lett. 2010, 10, 804−808. (17) Lopes Dos Santos, J. M. B.; Peres, N. M. R.; Castro Neto, A. H. Phys. Rev. Lett. 2007, 99, 256802. (18) Rakhmanov, A. L.; Rozhkov, A. V.; Sboychakov, A. O.; Nori, F. Phys. Rev. Lett. 2012, 109, 206801. (19) de Andres, P. L.; Ramírez, R.; Vergés, J. A. Phys. Rev. B: Condens. Matter Mater. Phys. 2008, 77, 045403. (20) Novoa, J. J.; Lafuente, P.; Del Sesto, R. E.; Miller, J. S. Angew. Chem., Int. Ed. 2001, 40, 2540−2545. (21) Itkis, M. E.; Chi, X.; Cordes, A. W.; Haddon, R. C. Science 2002, 296, 1443−1445. (22) Huang, J.; Kingsbury, S.; Kertesz, M. Phys. Chem. Chem. Phys. 2008, 10, 2625−2635. (23) Goto, K.; Kubo, T.; Yamamoto, K.; Nakasuji, K.; Sato, K.; Shiomi, D.; Takui, T.; Kubota, M.; Kobayashi, T.; Yakusi, K.; Ouyang, J. J. Am. Chem. Soc. 1999, 121, 1619−1620. (24) Small, D.; Zaitsev, V.; Jung, Y.; Rosokha, S. V.; Head-Gordon, M.; Kochi, J. K. J. Am. Chem. Soc. 2004, 126, 13850−13858. (25) Suzuki, S.; Morita, Y.; Fukui, K.; Sato, K.; Shiomi, D.; Takui, T.; Nakasuji, K. J. Am. Chem. Soc. 2006, 128, 2530−2531. (26) Pal, S. K.; Itkis, M. E.; Tham, F. S.; Reed, R. W.; Oakley, R. T.; Haddon, R. C. Science 2005, 309, 281−284. (27) Kubo, T.; Shimizu, A.; Sakamoto, M.; Uruichi, M.; Yakushi, K.; Nakano, M.; Shiomi, D.; Sato, K.; Takui, T.; Morita, Y.; Nakasuji, K. Angew. Chem., Int. Ed. 2005, 44, 6564−6868. (28) Takano, Y.; Taniguchi, T.; Isobe, H.; Kubo, T.; Morita, Y.; Yamamoto, K.; Nakasuji, K.; Takui, T.; Yamaguchi, K. J. Am. Chem. Soc. 2002, 124, 11122−11130. (29) Mota, F.; Miller, J. S.; Novoa, J. J. J. Am. Chem. Soc. 2009, 131, 7699−7707. (30) Huang, J.; Kertesz, M. J. Am. Chem. Soc. 2007, 129, 1634−1643. (31) Tian, Y.-H.; Kertesz, M. J. Am. Chem. Soc. 2010, 132, 10648− 10649. (32) Tian, Y.-H.; Kertesz, M. J. Phys. Chem. A 2011, 115, 13942− 13949. (33) Cui, Z.-H.; Lischka, H.; Beneberu, H. Z.; Kertesz, M. J. Am. Chem. Soc. 2014, 136, 5539−5542. (34) Tian, Y.-H.; Sumpter, B. G.; Du, S.; Huang, J. J. Phys. Chem. Lett. 2015, 6, 2318−2325. (35) Morita, Y.; Suzuki, S.; Sato, K.; Takui, T. Nat. Chem. 2011, 3, 197−204. (36) Sawada, K.; Ishii, F.; Sito, M.; Okada, S.; Kawai, T. Nano Lett. 2009, 9, 269−272.

(37) Muchharla, B.; Pathak, A.; Liu, Z.; Song, L.; Jayasekera, T.; Kar, S.; Vajtai, R.; Balicas, L.; Ajayan, P. M.; Talapatra, S.; Ali, N. Nano Lett. 2013, 13, 3476−3481. (38) Ponomarenko, L. A.; Geim, A. K.; Zhukov, A. A.; Jalil, R.; Morozov, S. V.; Novoselov, K. S.; Grigorieva, I. V.; Hill, E. H.; Cheianov, V. V.; Fal’ko, V. I.; Watanabe, K.; Taniguchi, T.; Gorbachev, R. V. Nat. Phys. 2011, 7, 958−961. (39) Lv, R. T.; Li, Q.; Botello-Mendez, A. R.; Hayashi, T.; Wang, B.; Berkdemir, A.; Hao, Q. Z.; Elias, A. L.; Cruz-Silva, R.; Gurierrez, H. R.; Kim, Y. A.; Muramatsu, H.; Zhu, J.; Endo, M.; Terrones, H.; Charlier, J. C.; Pan, M. H.; Terrones, M. Sci. Rep. 2012, 2, 586. (40) Jeon, I. Y.; Shin, Y. R.; Sohn, G. J.; Choi, H. J.; Bae, S. Y.; Mahmood, J.; Jung, S. M.; Seo, J. M.; Kim, M. J.; Chang, D. W.; Dai, L. M.; Baek, J. B. Proc. Natl. Acad. Sci. U. S. A. 2012, 109, 5588−5593. (41) Some, S.; Kim, J.; Lee, K.; Kulkarni, A.; Yoon, Y.; Lee, S.; Kim, T.; Lee, H. Adv. Mater. 2012, 24, 5481−5486. (42) Deng, D.; Pan, X.; Yu, L.; Cui, Y.; Jiang, Y.; Qi, J.; Li, W.-X.; Fu, Q.; Ma, X.; Xue, Q.; Sun, G.; Bao, X. Chem. Mater. 2011, 23, 1188− 1193. (43) Wang, X.; Li, X.; Zhang, L.; Yoon, Y.; Weber, P. K.; Wang, H.; Guo, J.; Dai, H. Science 2009, 324, 768−771. (44) Wang, H.; Wang, Q.; Cheng, Y.; Li, K.; Yao, Y.; Zhang, Q.; Dong, C.; Wang, P.; Schwingenschlögl, U.; Yang, W.; Zhang, X. X. Nano Lett. 2012, 12, 141−144. (45) Sun, Z.; Yan, Z.; Yao, J.; Beitler, E.; Zhu, Y.; Tour, J. M. Nature 2010, 468, 549−552. (46) Fujimoto, Y.; Saito, S. Phys. Rev. B: Condens. Matter Mater. Phys. 2011, 84, 245446. (47) Zhao, L.; He, R.; Rim, K. T.; Schiros, T.; Kim, K. S.; Zhou, H.; Gutiérrez, C.; Chockalingam, S. P.; Arguello, C. J.; Pálová, L.; Nordlund, D.; Hybertsen, M. S.; Reichman, D. R.; Heinz, T. F.; Kim, P.; Pinczuk, A.; Flynn, G. W.; Pasupathy, A. N. Science 2011, 333, 999−1003. (48) Schiros, T.; Nordlund, D.; Pálová, L.; Prezzi, D.; Zhao, L.; Kim, K. S.; Wurstbauer, U.; Gutiérrez, C.; Delongchamp, D.; Jaye, C.; Fischer, D.; Ogasawara, H.; Pettersson, L. G. M.; Reichman, D. R.; Kim, P.; Hybertsen, M. S.; Pasupathy, A. N. Nano Lett. 2012, 12, 4025−4031. (49) Meyer, J. C.; Kurasch, S.; Park, H. J.; Skakalova, V.; Künzel, D.; Groß, A.; Chuvilin, A.; Algara-Siller, G.; Roth, S.; Iwasaki, T.; Starke, U.; Smet, J. H.; Kaiser, U. Nat. Mater. 2011, 10, 209−214. (50) Nath, P.; Chowdhury, S.; Sanyal, D.; Jana, D. Carbon 2014, 73, 275−282. (51) Nath, P.; Sanyal, D.; Jana, D. Phys. E 2014, 56, 64−68. (52) Lu, Y.-F.; Lo, S.-T.; Lin, J.-C.; Zhang, W.; Lu, J.-Y.; Liu, F.-H.; Tseng, C.-M.; Lee, Y.-H.; Liang, C.-T.; Li, L.-J. ACS Nano 2013, 7, 6522−6532. (53) Park, S. H.; Chae, J.; Cho, M.-H.; Kim, J. H.; Yoo, K.-H.; Cho, S. W.; Kimde, T. G.; Kim, J. W. J. Mater. Chem. C 2014, 2, 933−939. (54) Usachov, D.; Vilkov, O.; Gr€uneis, A.; Haberer, D.; Fedorov, A.; Adamchuk, V. K.; Preobrajenski, A. B.; Dudin, P.; Barinov, A.; Oehzelt, M.; Laubschat, C.; Vyalikh, D. V. Nano Lett. 2011, 11, 5401−5407. (55) Wei, D.; Liu, Y.; Wang, Y.; Zhang, H.; Huang, L.; Yu, G. Nano Lett. 2009, 9, 1752−1758. (56) Shi, Z.; Kutana, A.; Yakobson, B. I. J. Phys. Chem. Lett. 2015, 6, 106−112. (57) Kresse, G.; Furthmüller, J. Comput. Mater. Sci. 1996, 6, 15−50. (58) Blöchl, P. E. Phys. Rev. B: Condens. Matter Mater. Phys. 1994, 50, 17953−17979. (59) Perdew, J. P.; Burke, K.; Ernzerhof, M. Phys. Rev. Lett. 1996, 77, 3865−3868. (60) Grimme, S. J. Comput. Chem. 2006, 27, 1787−1799. (61) Grimme, S.; Antony, J.; Ehrlich, S.; Krieg, S. J. Chem. Phys. 2010, 132, 154104. (62) Tkatchenko, A.; Scheffler, M. Phys. Rev. Lett. 2009, 102, 073005. (63) Tkatchenko, A.; Di Stasio, R. A.; Car, R.; Scheffler, M. Phys. Rev. Lett. 2012, 108, 236402. (64) Bučko, T.; Lebègue, S.; Hafner, J.; Á ngyán, J. G. Phys. Rev. B: Condens. Matter Mater. Phys. 2013, 87, 064110. I

DOI: 10.1021/acs.nanolett.5b01940 Nano Lett. XXXX, XXX, XXX−XXX

Letter

Nano Letters (65) Dion, M.; Rydberg, H.; Schröder, E.; Langreth, D. C.; Lundqvist, B. I. Phys. Rev. Lett. 2004, 92, 246401. (66) Klimeš, J.; Bowler, D. R.; Michaelides, A. J. Phys.: Condens. Matter 2010, 22, 022201. (67) Lee, K.; Murray, E. D.; Kong, L.; Lundqvist, B. I.; Langreth, D. C. Phys. Rev. B: Condens. Matter Mater. Phys. 2010, 82, 081101. (68) Harl, J.; Kresse, G. Phys. Rev. Lett. 2009, 103, 056401. (69) Blum, V.; Gehrke, R.; Hanke, F.; Havu, P.; Havu, V.; Ren, X.; Reuter, K.; Scheffler, M. Comput. Phys. Commun. 2009, 180, 2175− 2196. (70) Gao, W.; Tkatchenko, A. Phys. Rev. Lett. 2015, 114, 096101. (71) Lebegue, S.; Harl, J.; Gould, T.; Angyan, J. G.; Kresse, G.; Dobson, J. F. Phys. Rev. Lett. 2010, 105, 196401. (72) Soler, J. M.; Artacho, E.; Gale, J. D.; García, A.; Junquera, J.; Ordejón, P.; Sánchez-Portal, D. J. Phys.: Condens. Matter 2002, 14, 2745. (73) Troullier, T.; Martins, J. L. Phys. Rev. B: Condens. Matter Mater. Phys. 1991, 43, 1993−2006. (74) Zhao, Y. X.; Spain, I. L. Phys. Rev. B: Condens. Matter Mater. Phys. 1989, 40, 993−997. (75) Girifalco, L. A.; Lad, R. A. J. Chem. Phys. 1956, 25, 693−697. (76) Zacharia, R.; Ulbricht, H.; Hertel, T. Phys. Rev. B: Condens. Matter Mater. Phys. 2004, 69, 155406. (77) Lee, Y. S.; Kertesz, M. J. Chem. Soc., Chem. Commun. 1988, 75− 75. (78) Gould, T.; Lebègue, S.; Dobson, J. F. J. Phys.: Condens. Matter 2013, 25, 445010. (79) Dappe, Y. J.; Bolcatto, P. G.; Ortega, J.; Flores, F. J. Phys.: Condens. Matter 2012, 24, 424208. (80) Lebedeva, I. V.; Knizhnik, A. A.; Popov, A. M.; Lozovikda, Y. E.; Potapkin, B. V. Phys. Chem. Chem. Phys. 2011, 13, 5687−5695. (81) Chakarova-Käck, S. D.; Schröder, E.; Lundqvist, B. I.; Langreth, D. C. Phys. Rev. Lett. 2006, 96, 146107. (82) Lee, J.-K.; Lee, S.-C.; Ahn, J.-P.; Kim, S.-C.; Wilson, J. I. B.; John, P. J. Chem. Phys. 2008, 129, 234709. (83) von Lilienfeld, O. A.; Tkatchenko, A. J. Chem. Phys. 2010, 132, 234109. (84) Ambrosetti, A.; Reilly, A. M.; DiStasio, R. A., Jr.; Tkatchenko, A. J. Chem. Phys. 2014, 140, 18A508. (85) Charlier, J.-C.; Gonze, X.; Michenau, J.-P. Europhys. Lett. 1994, 28, 403−408. (86) Our PBE-D3//PBE-TS calculations indicate that the PLY πdimer in AB-stack configuration is 9 meV/C more stable than that of AA-stacked motif, which is in accordance with MR-AQCC calculations in ref 33 (9 meV/C). (87) Rani, P.; Jindal, V. K. RSC Adv. 2013, 3, 802−812. (88) Hoffmann, R. Solids and Surfaces: A Chemist’s View of Bonding in Extended Structures; Wiley-VCH: New York, 1988.

J

DOI: 10.1021/acs.nanolett.5b01940 Nano Lett. XXXX, XXX, XXX−XXX