Li and Na Adsorption on Graphene and Graphene Oxide Examined by Density Functional Theory, Quantum Theory of Atoms in Molecules, and Electron Localization Function

Adsorption of Li and Na on pristine and defective graphene and graphene oxide (GO) is studied using density functional theory (DFT) structural and electronic calculations, quantum theory of atoms in molecules (QTAIM), and electron localization function (ELF) analyses. DFT calculations show that Li and Na adsorptions on pristine graphene are not stable at all metal coverages examined here. However, the presence of defects on graphene support stabilizes both Li and Na adsorptions. Increased Li and Na coverages cause metal nucleation and weaken adsorption. Defective graphene is associated with the presence of band gaps and, thus, Li and Na adsorptions can be used to tune these gaps. Electronic calculations show that Li– and Na–graphene interactions are Coulombic: as Li and Na coverages increase, the metal valences partially hybridize with the graphene bands and weaken metal–graphene support interactions. However, for Li adsorption on single vacancy graphene, QTAIM, ELF, and overlap populations calculations show that the Li-C bond has some covalent character. The Li and Na adsorptions on GO are significantly stronger than on graphene and strengthen upon increased coverages. This is due to Li and Na forming bonds with both carbon and oxygen GO atoms. QTAIM and ELF are used to analyze the metal–C and metal–metal bonds (when metal nucleation is present). The Li and Na clusters may contain both covalent and metallic intra metal–metal bonds: This effect is related to the adsorption support selection. ELF bifurcation diagrams show individual metal–C and metal–metal interactions, as Li and Na are adsorbed on graphene and GO, at the metal coverages examined here.


Introduction
Rechargeable lithium ion batteries (LIB) serve as energy storage devices for a variety of applications, such as portable electronics and hybrid and electric vehicles [1][2][3][4]. The LIB performance depends on the structure and composition of the anode and cathode electrodes [5]. Graphite has been used as anode material for LIB [6], with graphite's theoretical capacity limit set at 372 mAh/g (intercalation compounds of LiC 6 configurations) [7]. Current and future needs dictate the search for LIB materials with capacities significantly higher than graphite to satisfy portable electronic devices' requirements [5]. The use of other carbon materials, such as graphene [8,9], C 60 , and carbon nanotubes, increases the LIB capacity to 540 mAh/g, 730 mAh/g, and 784 mAh/g, respectively [10]. Graphene is a monolayer honeycomb lattice of carbon atoms and is known for its high mechanical strength and high thermal and electrical conductivity [11][12][13]. Graphene B and N-doping further increase the LIB capacity to 1227 mAh/g and 872mAh/g, respectively [14].
Increased LIB graphene capacities relative to graphite seem indicative of increased Li atoms being adsorbed on graphene relative to graphite. However, the contrary has been observed via in situ Raman spectroscopy [15]. Although computational simulations support the adsorption of single Li atoms on graphene (one side [16][17][18] and both side adsorption [5,19]), this adsorption is less stable relative to bulk metallic Li. Lee et al. showed that Li adsorption on pristine graphene is not favorable relative to bulk metallic Li under any configuration and is only possible on graphene defects and edges [20]. Okamoto further supported Li adsorption on defects relative to edges [21]. This is due to the edges of the anode LIB materials being covered by low electron conductivity solid electrolyte interphase, thus affecting the reversibility of Li movement during charge and discharge processes. Li adsorption on defective graphene has been reported using computational calculations [21][22][23][24]. The presence of double carbon vacancies on N-doped graphene support can also improve LIB performance [25].
The main disadvantage of LIBs is Li cost, which has been rising steadily in recent years, with the exception of a price drop observed in the last quarter of 2018 [26]. Li is a very rare element (Li concentration is~35 ppm in the upper continental crust [27]), and thus there are concerns that Li may not be sufficient for meeting future energy demands. LIB, as well as sodium ion batteries (NIB) [26,28], suffer from dendrite growth due to metal nucleation (i.e., clustering), which lowers their performance [1,29]. Li forms clusters due to the small Li diffusion barrier on graphene support [30]. Li cluster formation on pristine graphene (i.e., graphene without defects) and defective graphene [21] has been explored using density functional theory (DFT) [30,31]. Defective graphene is produced during graphene oxide's (GO) thermal reduction process, as hydroxyl and epoxide O atoms are removed from the GO sheet, via oxidation and reduction [32,33]. Li 6 clusters collapse during adsorption on 24 vacancies defective graphene [21].
Sodium ion batteries could serve as alternative batteries due to the abundance of Na. However, the NIB capacity of graphite is only 35 mAh/g, substantially lower relative to LIB [26]. This is attributed to NaC 64 formation [34,35] due to the large Na + size relative to Li + (100 pm vs. 76 pm) and Na high ionization potential [21]. NIB capacities were recorded at 1450 mAh/g, for Na adsorption on defective graphene with divacancy [24]. DFT calculations on Na adsorbed on defective graphene have been reported [36]. Specifically, Na cluster formation for adsorption on pristine and defective graphene has been explored via DFT calculations by Liang et al. [37].
In this work, we examine Li and Na adsorptions on pristine and defective graphene, as well as on graphene oxide (GO), at various metal coverages, using periodic DFT, with and without long-range interactions. The metal-support interactions are elucidated the via changes in the adsorption and cohesive energies (E ads and E c , respectively) per metal, the metal-Carbon plane distances (height, h), and the electronic charge transfers from the metals to graphene and GO supports. The E ads and the E c for the adsorption of n-metals (n = 1, 2, . . . ) on graphene and GO are defined as where X is the metal adsorbate and the support is either graphene or GO. The E ads is the minimum energy required to remove a single Li or Na atom during adsorption on graphene and GO, whereas E c measures the cluster-support interaction. For n = 1, E ads = E c . Single metals and their clusters form a bond with the support if E c < 0. The metal-support interaction is considered stable if the |E ads | is larger than the metal binding energy of its bulk. The metal-support interactions are also examined via electronic structure analyses. The local metal-C and the metal-metal interactions are examined by two independent analyses: The quantum theory of atoms in molecules (QTAIM) and the electron localization function (ELF). QTAIM was developed by Bader and co-workers, who note that orbitals, much less wave functions, are unphysical in nature; thus, chemical bonds are not quantum observables [38,39]. However, the electron densities (ρ → r ) and their gradients are quantum observables and can be used to characterize chemical systems [39,40]. ELF [41] is a measure of the probability of finding an electron in the vicinity of another same-spin electron. It reveals core, binding, and lone electron pairs in molecular and crystalline systems. Results from both QTAIM and ELF analysis are basis-set and method-independent, as long as a minimally adequate basis set is used in these calculations [42,43].

Periodic Structure Modeling
Monolayer pristine graphene is modeled as a two-dimensional 4 × 4 hexagonal lattice with 32 carbon atoms in the unit cell. Figure 1 shows the DFT optimized structures for Li n and Na n (n = 1, 3, 5) adsorbed on pristine graphene, the Li and Na E c values, the heights h for the closest metal to the graphene plane, and the smaller of the metal-metal distances per configuration. The Li n /graphene and Na n /graphene (n =1, 3, 5) unit cells are built as follows: for n = 1, single metal atoms (i.e., Li and Na) are initially placed distant from graphene, followed by geometry optimizations; for n = 3, two metal atoms are placed distant from the previously optimized Li/graphene and Na/graphene structures, followed by geometry optimizations. A similar approach is followed for the Li 5 /graphene and Na 5 /graphene cases. Figure 2 shows the optimized structures for Li n and Na n (n =1, 3, 5) adsorbed on single and double vacancy defective graphene, as well as the E c and geometrical parameters, in a similar fashion as in Figure 1. Graphene defects are modeled by the use of "ghost" massless C atoms in the graphene lattice: These atoms share the same basis set as C, in agreement with our past report [44]. Here, the ghost atoms are placed in the center of the supercell and the defect is periodically repeated by the use of translational symmetry vectors [45]. Ghost atoms remain fixed in the locations occupied by the original C atoms. Figure 3 shows the optimized structures for Li n and Na n (n = 1, 3) adsorbed on graphene oxide (GO, O/C = 43.75%), with E c and geometrical parameters. Here, GO supports contain epoxides and hydroxyls at 25% and 18.75%, respectively (percentages relative to C). Our GO structure is in agreement with the proposed non-hydrated GO structure from Zhou and Bongiorno [46] (26% and 18% for epoxides and hydroxyls, respectively), found by using X-ray photoelectron spectroscopy and DFT.

DFT Parameters
CRYSTAL17 [47] is used to calculate the optimized geometries and electronic properties for Li and Na atoms adsorbed on graphene (with and without vacancies) and GO. An advantage of CRYSTAL is its treatment on two-dimensional slabs, without the need for artificial periodicity perpendicular to the surface, as is typically found in other DFT codes [48]. Periodic restricted DFT is used here, as in our past work [17,44,[49][50][51]. The PBE0 non-empirical/parameter-free functional is used [52,53]. All calculations are treated using the D2 semiempirical correction by Grimme [54]. This correction improves the DFT functionals description of the long-range electron correlations, which are responsible for van der Waals interactions between the adsorbed metals and graphene/GO supports. Moreover, for selective cases, our calculations are repeated without employing D2 corrections. All atoms are described by all-electron basis sets. The C and H atoms are treated by the 3-ζ valence plus polarization basis set (pob-TZVP) from Peintinger et al. [55] and the O atoms by the split-valence 8-411G(2d1f) basis set from Mahmoud et al. [56]. These basis sets are optimized for crystalline calculations. The Li and Na are described by the atomic basis sets (6s3p1d) and (9s5p1d), respectively, and were used in our past work [17]. Calculations are at 0 K, as expected. Temperature and pressure effects will be described in a future work. Figure 1. Optimized geometries, heights h (blue arrows) for the metal (Li, Na) closest to graphene plane, the closest metal-metal distances (red arrows), and E c for Li n (yellow spheres) (a-c) and Na n (purple spheres) (d-f) (n = 1, 3, 5) adsorbed on pristine graphene (C, gray spheres). Values in parentheses include D2 correction. The thin black line is the unit cell boundary.

Figure 2.
Optimized geometries, heights h (blue arrows) for the metal (Li, Na) closest to graphene plane, the closest metal-metal distances (red arrows), and E c for Li n (yellow spheres) (a-f) and Na n (purple spheres) (g-n) (n = 1, 3, 5) adsorbed on defective graphene (C, gray spheres, "ghost" C, blue spheres) with one and two C vacancies. All calculations include D2 correction. The thin black line is the unit cell boundary. The thick dashed black line is the carbon plane.
For geometry optimizations, Brillouin zone integrations (Monkhorst-Pack grid) [57], the Fermi energy, and the density matrix calculations (Gilat grid) [58,59] are performed on a 12 × 12 grid. The densities of states (DOS) and the electronic band structures are obtained using the denser 48 × 48 grid. The band structures are calculated on the reciprocal space Brillouin zone path Γ-M-K-M´-Γ, where Γ = (0, 0, 0), M = (1/2, 0, 0), K= (2/3, 1/3, 0), and M´= (0, 1/2, 0). The Fermi surface is smeared with a Gaussian of 0.005 Hartrees for convergence purposes. Moreover, the SCF energy convergence is achieved using Anderson quadratic mixing [60], coupled with additional mixing of the occupied with the virtual orbitals. The SCF energy threshold value for our calculations is set at 10 -9 Hartrees (default value is 10 -7 Hartrees). A large integration grid is used: this is a pruned grid with 75 radial and 974 angular points. The Crystal Orbital Overlap Populations (COOP) [61], for selected pairs of atoms, are calculated directly by CRYSTAL17. Charge transfers are calculated using Mulliken population analysis [62] and the iterative Hirshfeld population [63], the latter being based on the atoms in molecules (AIM) approach [38,39]. Hirshfeld populations are less basis-set depended than Mulliken [64] and smaller in magnitude (relative to Mulliken) [65]. The absolute magnitude of the calculated charges obtained through populations analysis has little physical meaning: These calculations may depend on the method and basis sets used [65], and thus we only focus on the relative changes of these values.
The stabilities of the final geometry conformations are secured via post-geometry optimizations. When Gaussian basis sets are used, E ads and E c calculations suffer from basis set superposition errors (BSSE) [66]. The BSSE errors are minimized using large basis sets and employing counterpoise corrections using "ghost" massless atoms in the fragment energy calculations of the metal-graphene structures. Here, we only report BSSE corrected E ads and E c values. Calculated E ads and E c values using the counterpoise correction are less negative in energy than corresponding uncorrected values [17,49,67].

QTAIM Parameters
The ρ → r topology is performed by the TOPOND program [68], which is incorporated in CRYSTAL17. TOPOND calculates ρ → r and its Laplacian (∇ 2 ρ → r ) at critical points, for periodic and non-periodic systems. The critical points in the ρ  [70,71]. Gatti stated that adoption of a single QTAIM criterion for bond assessment is challenging [70]. Therefore, we used ρ  [72]. For metallic bonds, the above conditions are modified (units of length 2 ) serves as an additional criterion to classify non-polar covalent and metallic bonds [73,74].

ELF Parameters
The electron localization function (ELF) is an alternative to QTAIM for analyzing chemical bonding [43,[75][76][77]. ELF measures the electron localization in atomic and molecular systems and it was originally described by Becke and Edgecombe [41] as follows: where χ → r is the relevant ELF kernel and 0 ≤ η → r ≤ 1 [70]. The η → r = 1 case corresponds to perfect localization, whereas the η → r = 0.5 to the electron gas like pair probability. The ELF function η → r reveals atomic shell structure and core, similar to −∇ 2 ρ → r . However, −∇ 2 ρ → r fails to reveal the shells structure for heavy atoms [78]. Kohout and Savin found that η → r can resolve all shells for all atoms from Li to Sr [79]. Savin et al. [80] interpreted ELF as a local measure of the Pauli repulsion on the kinetic energy ρ → r , and applied ELF to a variety of atoms, molecules, and solids.
Specifically, η → r ∼ 1 is indicative of weaker Pauli repulsion than in the uniform electron gas. ELF gradient analysis reveals attractors (local maxima) and corresponding basins [43]. In QTAIM, basins are regions in space bounded by zero flux surfaces in the gradient of ρ → r [81]. Basins are classified as core and valance and characterized by their synaptic order (i.e., the number of core basins that share a common boundary [82]). Only core basins contain a nucleus, whereas valence basins are always connected to at least one core basin. At low η → r values only one surface is visualized, whereas at higher η → r basins separate. Therefore, bifurcation (i.e., tree) diagrams [83], with the bifurcation points corresponding to η → r minima, can be used to examine bond types and regions in space, where lone pair electrons are located. The first bifurcation point corresponds the separation of the core and the valance basins from the same atom. The lower the bifurcation point the more localized are the corresponding basins [43]. The ELF 3D isosurfaces are plotted using Chimera [84].

Adsorption on Pristine Graphene
At 1/32 monolayer (ML) coverage, Li and Na are adsorbed on the hollow site, in agreement with past reports [5,[16][17][18]22]. At the higher coverages 3/32 ML and 5/32 ML, both metals form clusters (Li 3 , Li 5 , Na 3 , and Na 5 , Figure 1b,c,e,f), as expected [30,31,37]. Here, the effect of D2 corrections in our calculated E c values is small. However, disagreements between the DFT calculations with and without D2 corrections are observed in the h values and the metal-metal closest distances (e.g., for Na 5 /graphene, |∆h| = 0.11 Å, Figure 1f). For Li 3 and Li 5 adsorptions on pristine graphene, our closest calculated Li-Li distances are 2.85 Å and 2.77 Å, when D2 corrections are included in our calculations (Figure 1b,c). These values are smaller than the bulk Li closest Li-Li distance of 3.04 Å [85]. Similarly, for the Na case, our closest Na-Na distances during adsorption on pristine graphene are less than the bulk Na closest Na-Na distance of 3.72 Å [86]. These are indicative of the presence of stronger Li and Na intra-cluster bonds relative to the ones in their metallic bulks. Figure 4 shows the E ads for Li and Na adsorptions on graphene and GO. At the same metal coverages, Li adsorption on pristine graphene is significantly stronger relative to Na, as shown by E ads , E c , and h values. This statement is extended for adsorption on defective graphene and GO, with the exception of Na 3 /GO calculations without D2 corrections (vide infra). Our Li and Na |E ads | values are smaller than the binding energies of bulk Li (1.63 eV/Li [87]) and Na (1.13 eV/Na [88]): These are indicative of unstable Li and Na adsorptions on pristine graphene, at the coverages of this work, in agreement Lee at al. [20] and Liang et al. [37]. Although our calculations are at 0 K, the above metal stability statements are in agreement with experimental observations at room temperature. For example, Na adsorption on pristine graphene has only been detected under Na/C 64 configuration [89], where in this case Na |E ads | is larger than its bulk binding energy [37]. For adsorption on pristine graphene, the E c becomes less negative along with increased metal coverage, thus weakening the metal-support interaction, as expected [17,37]. This statement is extended for adsorption on defective graphene. Table 1 shows our E ads values and from past reports, for Li and Na adsorbed on pristine and single vacancy graphene: our E ads values are in good agreement with past work.   [94] PBE+D3 [91,93] −2.59 [94] PBE+vdW-DF2 [91,95] Na /Graphene 0 −0.59 (−0.59) −0.55 [90], −0.462 [16] PBE [91] −0.72 [92] a LDA −0.93 [90] PBE+D2 [54,91] −0.64 [90] PBE+D3 [91,93] −0.49 [90] PBE+vdW-DF2 [91,95] −0.88 [37] PBE+D2 [

Adsorption on Defective Graphene
Li and Na adsorptions on defective graphene are stable ( Figure 4). Here, Li and Na |E ads | values at all coverages are larger than the binding energies of their metallic bulk, in agreement with past reports [20,37]. The presence of defects on graphene significantly strengthens Li and Na adsorption relative to their adsorption on pristine graphene (Figures 1 and 2). Increased carbon vacancies decrease the metal-graphene interaction (Figure 2). For example, for Li/graphene at 1/32 ML coverage, increased C vacancies decrease |E c | from 3.44 eV to 2.14 eV. At 3/32 ML coverage, the presence of a single C vacancy causes the Li clusters to collapse, whereas in all other cases Li and Na form clusters, in a similar fashion as during adsorption on pristine graphene.
The strong adsorption of Na on defective graphene (and GO, vide infra) is indicative of its possible use in NIB. Although, Na adsorption on graphene is not as strong as Li, its abundancy is a significant advantage over Li.

Adsorption on Graphene Oxide
Li and Na adsorptions on GO are also stable ( Figure 4). Moreover, contrary to adsorption on pristine and defective graphene, increased metal coverage strengthens the metal-GO interaction. This statement does not apply for the Na 3 /GO geometry calculated without D2 corrections. Figure 3d,e show striking differences between the optimized geometries for Na 3 /GO with and without D2 corrections. Therefore, for Na adsorption on GO, at coverages > 1/32 ML, it is important to include long-range interactions due to Na atoms' large polarizations [37]. For Na 3 /GO, Na clustering is observed, in contrast to the Li 3 case. Table 2 shows the calculated Li and Na Hirshfeld and Mulliken charges and the metal-C COOP values for the smaller metal-C distances per configuration, for Li and Na absorbed on graphene and GO, under the metal coverages examined here. Chan et al. stated that alkalis are very close to ideal ionic bonding and, thus, Li-C and Na-C bonds are expected to be ionic [16]. There is a positive relationship between the average per metal charge transferred to the support and the metal-C bond strength. Both Hirshfeld and Mulliken calculations show that charge is transferred from Li and Na towards the support (graphene or GO). This charge transfer is decreased along with increased metal coverage, in agreement with decreases in the |E c | and, thus, the weakening of the metal-support interaction. Dobrota et al. showed that Na adsorption on graphene-OH compounds leads to Na being neutral [96]. This is in contrast with our work, which shows that charge is transferred from Li and Na to GO (Table 2). However, the support of Dobrota et al. represents a reduced GO structure (rGO) with a significantly smaller O/C ratio than in our case.

Charge transfers and COOP
The metal-C bond strengths can be estimated via changes in their COOP values. The majority of the metal-C COOP values in Table 1 are of the order expected for ionic bonding. However, there are some exceptions: The Li-C COOP value for Li adsorption on single C vacancy graphene (i.e., 0.138) is at the order expected for covalent bonding. Doan et al., using COOP, found that Li bonds covalently to the Nafion SO 3 group [97]. Our COOP calculations indicate that metal-C bonds, for adsorption on single vacancy graphene, are covalent (Table 1). COOP calculations are basis set dependent [98]. Therefore, these bonds will be also examined via the method and basis set independent QTAIM and ELF analyses in the next sections.

Adsorption on Pristine Graphene
The band structure calculations show that Li and Na adsorptions downshift the Dirac point (E D , DOS(E D ) = 0) below the Fermi energy ( Figures 5 and 6). Recall that pristine graphene E D coincides with the Fermi energy and appears in the band structure as a K-point, where the C pz valance (π) and conduction bands (π*) meet in a linear dispersion fashion [8]. Therefore, pristine graphene is a zero-gap semiconductor (i.e., semi-metal) [99]. The E D downshift, due to Li and Na adsorptions on pristine graphene, is indicative of n-doped graphene support via charge transfer from Li and Na. At low Li and Na coverages, the E D is intact (Figures 5d and 6d) and the Li-2s and Na-3s valance orbitals appear in the band structures as flat bands near and above the Fermi energies. Moreover, the DOS spectra below the Fermi energy are unaffected of the adsorption (Figures 5a and 6a), thus supporting a Columbic interaction between the metals and pristine graphene. The DOS and band structure calculations show that Li/graphene and Na/graphene are metallic.
At higher Li coverages, the E D is broken (Figure 5e,f), whereas for the Na case the E D is only broken at the 5/32 ML Na coverage (Figure 6f). However, only for the Li 3 /graphene case, a small band gap of 0.05 eV is observed in the E D region (Figure 5e,f and Figure 6f), turning Li 3 /graphene to a semiconductor. However, this effect may be an artifact of the calculations and it will be further examined in the future. All other metal-pristine graphene configurations are metallic. The flat bands below the Fermi energy in Figures 5e and 6e and the wavy bands in Figures 5f and 6f represent Li-2s, Na-3s, and Na-2p valence states. In the DOS spectrum, these bands appear as peaks near and below the Fermi energy (for Na 5 case this peak is broadened) and thus, partially occupied. The partial occupancy of these states decreases the Coulombic interaction between the metals and pristine graphene, as the metal coverage increases. This statement also applies for adsorption on defective graphene.

Adsorption on Defective Graphene.
Defects on graphene upshift the E D above the Fermi energy and create band gaps. For example, the double vacancy graphene is a semiconductor [37]. This E D upshift is indicative of a p-doped graphene support [37]. At low Li and Na coverages and for Na 3 adsorbed on double vacancy graphene, the E D appears above the Fermi energy (Figure 7d, Figure 8d, and Figure 10e). Liang et al. showed that doping defective graphene with Na ions downshifted E D , but still it remained above the Fermi energy. Therefore, Li and Na adsorptions on defective graphene decrease p-doping. Moreover, at low metal coverages, the observed band gaps further grow as vacancies increase. The DOS spectra (Figures 7a and 8a) show minimal hybridization of the metal valance states with the graphene bands, which supports metal-graphene Coulombic interaction. However, we will show in the next section that the Li-C bond for Li adsorbed on single vacancy graphene has some covalent character.
With the exception of Na 3 adsorbed on double vacancy graphene, increased Li and Na coverages on defective graphene push the E D below the Fermi energy, thus changing graphene to n-doped. For Li adsorption on single vacancy graphene the band gap grows due to Li clustering, whereas the opposite is observed for the Na case. Multiple band gaps are observed for select Li and Na clusters adsorbed on double vacancy graphene. Here, all metal-defective graphene configurations are semiconductors. Our Na band structure calculations are in agreement with Liang et al. [37].

Adsorption on GO
Graphene Oxide, at O/C = 43.75%, is an insulator. Its calculated band gap is 5.26 eV, when only epoxides are present in the GO structure (PBE0 calculations without D2 correction) [50]. This value is decreased to 4.40 eV, for GO that includes both epoxides and hydroxyls, under the same DFT calculations. The GO band structure calculations do not reveal E D (Figure 11). Li and Na adsorption on GO is associated with multiple band gaps below the Fermi energy. For example, the largest of these band gaps are 3.23 eV and 3.17 eV for Li/GO and Na/GO, respectively. The largest band gap width and the Li and Na coverages are negatively correlated ( Figure 11). The GO DOS spectra below the Fermi energy show minimal contributions from the metals valance states. However, the Li/GO and Na/GO DOS spectra below the Fermi energy are different, which implies structural changes on GO support due to metal adsorption. This statement also applies to adsorptions at higher metal coverages.    Tables 3 and 4 show the QTAIM properties ρ → r , ∇ 2 ρ → r , and (G/ρ) → r at the Li-C and C-Na bond critical points, respectively (closest metal-C distances). Moreover, for adsorption on GO, the above QTAIM properties are also evaluated at the metal-O bond critical points. Table 3 shows that some Li-C bond critical points for Li 3 adsorption on pristine graphene and GO are not revealed by TOPOND. Moreover, for the former case, the Li-C critical point at 2.28 Å appears as a (3, +1) point (i.e., charge minima). In all cases here, with the exception of Li adsorbed on single vacancy graphene, QTAIM shows that metal-C bonds are ionic. Recall that for ionic bonding ∇ 2 ρ → r cp > 0, (H/ρ) → r cp > 0, and (G/ρ) → r cp ≥ 1, whereas ρ → r cp is small. The Li-C bond, for Li adsorption on single vacancy graphene, is analyzed using the bond classification of Espinosa et al. [71].

The Metal-C Bond
Here (|V|/G) → r cp < 1 is indicative of ionic bonding and (|V|/G) → r cp = 1 is the lower limit of the transit region indicating incipient covalent bond formation. Recall that for this case, the Li-C COOP value is 0.138, and is the highest value observed in this work for metal-C bonds. Therefore, the presence of vacancies in the support may add a covalent character in the metal-C bonds.
QTAIM shows that for adsorption on GO, Li and Na form ionic bonds of similar strength with both C and O atoms. This explains the significantly stronger metal-GO interaction relative to adsorption on pristine graphene, at the same metal coverages. Figure 13 shows a negative correlation between (H/ρ) → r and ρ → r thus, verifying ionic metal-C bonding.    Tables 5 and 6 show the QTAIM properties ρ( r →), (H/ρ)( r →), (G/ρ)( r →), ξ( r →), and the sign of ∇ 2 ρ → r and ELF η → r values at the Li-Li and Na-Na bond critical points, respectively. Not all metal-metal bond critical points are revealed for both Li and Na clusters adsorbed on graphene and GO (Tables 4 and 5). Moreover, the Li-Li bond critical points are not revealed for Li 3 adsorbed on defective graphene and for Li 3 and Li 5 adsorbed on GO. Recall that for metallic bonds ξ > 1, whereas ξ < 0, due to ρ < 0, is indicative of non-metallic non-polar covalent bonds. For Li 3 adsorbed on pristine graphene, Na 3 adsorbed on pristine and double vacancy graphene, and for Li 5 and Na 5 adsorbed on defective graphene all intra metal-metal bonds are metallic (ξ > 1, Table 3). However, Li 5 , Na 3 , Na 5 show a mixture of covalent and metallic metal-metal bonds, when are adsorbed on other supports (e.g., Li 5 on pristine graphene, Table 5). Therefore, the support selection during metal cluster adsorption could alter the metal-metal bond types.

ELF Analysis
Ayers and Jenkins stated that metallic bonds rarely have ELF η → r values larger than 0.6 [73].
Our η → r values, for metallic metal-metal bonds at their bond critical points, are less than 0.43, which is within the Ayers and Jenkins ELF limit for metallic bonds. However, some Li-Li and Na-Na non-metallic covalent bonds have low η → r values comparable to the ones for metallic bonds.
For example, for Li 5 adsorbed on pristine graphene the Li-Li bond with bond length 2.91 Å has a negative metallicity but its η → r value is 0.39. A similar case appears for Na 3 /GO.  Figures 15 and 16 show the bifurcation diagrams and corresponding isosurfaces for Li and Na adsorbed on pristine and single vacancy graphene, respectively. The values shown in bifurcation diagrams correspond to the η → r value, where a separation between the basins occurs. We start our ELF analyzes with Li and Na adsorption on pristine graphene, at low metal coverages.   For Li/graphene, the bifurcation diagram reveals the Li and C core basins (C(Li) and C(C)) and the polysynaptic valence basin V(C,C) (Figure 15a). For Na/graphene, the bifurcation diagram additionally reveals a monosynaptic valance V(Na) (Figure 15g). For Li/graphene, the first bifurcation occurs at η → r = 0.04, where C(Li) separates from the disynaptic basin. Here, the second bifurcation occurs at η → r = 0.12, for the C(C) separation, whereas the V(C,C) separation appears at η → r = 0.69. Similar η → r separation values are found for Na/graphene. Our ELF separation values are in agreement with those found by Matito et al. [81], who studied the bonding in methylalkalimetals (CH 3 M) n (M = Li, Na, and K; n = 1, 4). The η → r C(Li) and C(Na) separation values are indicative of ionic bonding between Li and Na metals and C, in agreement with QTAIM and electronic structure analysis (vide supra). It is noteworthy that our C(Na) η → r separation values are never higher than C(Li) values, in agreement with the weaker Na-graphene interaction relative to Li-graphene (Figures 15 and 16).
At low metal coverages and for adsorption on single vacancy graphene, the first bifurcation point occurs at higher η → r values relative to adsorption on pristine graphene. These higher η → r separation values are indicative of stronger metal-C bonding and are in the covalent bond region (Figure 16a). COOP, QTAIM, and ELF agree that the Li-C bond, for adsorption on single vacancy graphene, has some covalent character. However, the Na-C covalent character, for Na adsorbed on single vacancy graphene, indicated by ELF analyses is not agreement with QTAIM, whereas it agrees with COOP. Disparagements between ELF and QTAIM have been reported [81].

Conclusions
Li and Na adsorptions on pristine and defective graphene and GO has been studied using periodic DFT structural and electronic calculations paired with COOP, QTAIM, and ELF analyzes. DFT optimized geometries show that increased metal coverage clusters Li and Na during adsorption on graphene, which lowers the efficiency of LIB and NIB. However, Li 3 clusters collapse during adsorption on single vacancy graphene and GO. The E ads calculations show that Li and Na adsorption on pristine graphene is unstable at all the metal coverages examined here, whereas the presence of defects in the graphene support stabilizes Li and Na adsorptions. Increased metal coverage weakens metal-graphene and strengthen metal-GO interactions. Adsorption on GO is strong due to simultaneous metal bonding with C and O GO atoms.
DOS calculations show that, at low metal coverages, Li and Na adsorptions on pristine does not affect the DOS spectra below the Fermi energy. This effect is indicative of Coulombic interactions between the adsorbates and the support. These interactions weaken as metal coverages increase, due to the downshift of the metal valances below the Fermi energy. The band calculations show that the E D is pushed below the Fermi energy due to Li and Na adsorptions on pristine graphene and may be associated with small band gaps at higher metal coverages. The presence of defects on graphene is associated with E D shifting and bang gap openings.
The metal-C bonds are analyzed using COOP, QTAIM, and ELF. The COOP metal-C COOP values at the order expected for ionic bonding, for all cases except for adsorption on single vacancy graphene, at low metal coverages. QTAIM shows that only the Li-C bond contains covalent character for adsorption on single vacancy graphene, whereas in all other cases the metal-C bond is ionic. ELF partially agrees with QTAIM and shows that both Li-C and Na-C bonds are covalent for adsorption on single vacancy graphene. QTAIM and EFT are also used to study Li-Li and Na-Na bonds: Li and Na clusters may contain metallic and covalent metal-metal bonds depending on adsorption support (i.e., pristine and graphene and GO).