Corrugated graphene exposes the limits of a widely used ab initio van der Waals DFT functional

Theoretical formulations capable of modeling chemical interactions over 3–4 orders of magnitude of bond strength, from covalent to van der Waals (vdW) forces, are one of the primary goals in materials physics, and chemistry. Development of vdW corrections for density-functional theory has thus been a major research ﬁeld for two decades. While many of these corrections are semiempirical, more theoretically rigorous ab initio functionals have been developed. The ab initio functional vdW-DF2, when coupled with the reoptimized B86 exchange function (vdW-DF2-rB86), has typically performed as well, if not better than most semiempirical formulations. Here we present a system, Co intercalation of graphene on Ir(111), for which a semiempirical correction predicts local corrugation maxima in locations at which the vdW-DF2-rB86 functional predicts global minima. Sub-angstrom precision quantitative structural measurements show better agreement with the semiempirical correction. We posit that it is balancing the weak vdW interaction with the stronger, even covalent, interactions that proves a challenge for the vdW-DF2-rB86 functional.

(e.g., the interplay of covalent, polar, and van der Waals interactions determining the shape of proteins). Thus, theoretical formulations that can span these 3-4 orders of magnitude of bond strength, modeling interactions with disparate origins, are one of the primary goals in theoretical materials physics and chemistry. To this end, the last two decades have seen a rapid development of van der Waals corrections for densityfunctional theory (DFT) formalizations.
Two primary methodologies have ultimately emerged out of this work. The first, typified by the DFT-D2 [3] formalization, is a semiempirical correction that is added to the Kohn-Sham energy. This correction broadly mirrors the first term in London's Taylor series expansion for the London dispersion force, with a 1 /R 6 relationship, summed over all atoms with a semiempirical dispersion coefficient, scaling factors, and a cutoff radius below which the dispersion force is damped. The other methodology, typified by DFT vdW-DF2, attempts a more ab initio approximation of the dispersion force, effectively taking the exchange energy from the generalized gradient approximation (GGA), the local correlation at local density approximation (LDA) level, and calculating the nonlocal correlation energy by an integration over electron densities on different space locations that interact by a given response function [4]. Both methods have seen considerable success in accurately modeling experimental data [5][6][7][8][9][10][11][12][13][14][15] and have become regarded as essential in any DFT calculation; calculations that lack such corrections, especially for systems in which molecular π orbitals play a significant role, are generally regarded as inadequate.
In most of the direct comparisons of the two methodologies, e.g., Refs. [15][16][17][18][19] (although with some exceptions: e.g., Refs. [20][21][22][23]) only minute differences, especially in the structure, have been observed. Here we present a system, cobalt intercalation of graphene (Gr) on Ir(111), in which the DFT vdW-DF2 functional, utilizing the reoptimized Becke (rB86) exchange functional [24] (vdW-DF2-rB86), predicts a significantly different structure from that of the semiempirical DFT-D2 formulation; this difference would seem to be due to the coexistence of the strongly varying bonding strengths present in the system. The results of these calculations are benchmarked against quantitative structural measurements, utilizing the normal incidence x-ray standing-wave (NIXSW) method, leading to the conclusion that it is the semiempirical DFT-D2 formulation that is more consistent with the experimental data. We do not draw general conclusions as to the efficacy, or otherwise, of the DFT-D2 functional, as there are a multitude of successful semiempirical corrections for dispersion forces (e.g., vdW-surf [25], DFT-D3 [26], H06-L [27], etc.). Instead, we use the good agreement with experiment found in the DFT-D2 calculations to understand the potential causes for the challenge that this system proves to present for the DFT vdW-DF2-rB86 functional, which has provided highly accurate results for other adsorbed graphene and graphenelike (e.g., boron nitride) layers in the past [6,9,10,28,29] (including, notably, comparison to NIXSW measurements [10]), by analyzing the level of covalent interaction between the C and Co atoms.
Graphene, as is well known, exhibits unique electronic and mechanical properties, leading to an extraordinarily diverse range of suggested applications from the comparatively mundane: e.g., water filtration/desalination [3,30], transparent electrodes [31], and high-speed electronics [32,33], to the more exotic: e.g., supercapacitors [34], and "valleytronics" [35,36]. However, many of these applications, especially those in the field of semiconductors or magnetism, require interfacing graphene with other materials, and interfacial interactions can have wildly varying effects on the electronic properties of the graphene [37]. Intercalating a layer of metal atoms, between the Gr layer and its substrate, is a popular method to tailor the electronic properties of Gr; for example, it can induce ferromagnetism [38][39][40], or even induce and tune a band gap [41]. However, it is clear that such significant modification of the electronic structure of Gr may disrupt some of its desired properties, such as the remarkably high mobility of charge carriers. More significantly, it may result in a lifting of the sp 2 -hybridized character that defines graphene, resulting in a more sp 3 -hybridized, graphanelike, layer. Of course, graphane, as a two-dimensional semiconductor, is itself also the focus of much research [42][43][44][45], so the conversion of a graphenelike layer to a graphanelike layer is not necessarily deleterious.
If the interaction between graphene and the intercalant is indeed strong enough to drastically perturb the sp 2hybridization this would present a fundamentally different system from any previously explored using the vdW-DF2-rB86 functional. If the intercalant even forms a direct covalent bond between the intercalant and C atoms, it would present a system where both van der Waals and covalent interactions play an important role, and thus require the theoretical construct to span a large difference in bonding strength. As we demonstrate in this paper, it is this perturbation of the sp 2 hybridization and the necessity of balancing the covalent and van der Waals interactions that most probably cause the significant disagreement between the predictions of the vdW-DF2-rB86 functional and the experimental measurements.

A. Experimental details
The measurements were performed at the I09 beam line of the Diamond Light Source in the permanent end station, base pressure ∼4 × 10 −10 mbar, which is equipped with the standard equipment for in situ sample cleaning by sputtering and annealing and surface characterization. The I09 beam line provides both soft (110-1100 eV) and hard (2200-15 000 eV) x-ray radiation in parallel. The NIXSW technique exploits the standing wave created by the interference between the incident and reflected x-ray beams when a Bragg condition is satisfied in a crystal. This standing wave extends into and out of the bulk of the material and has a periodicity equal to the spacing between the Bragg planes, d hkl . The standing wave exists over a narrow range of incident photon energies, within which its phase shifts relative to the atomic scattering planes as a function of the incident photon energy, causing the nodes and antinodes of the standing wave to move relative to the Bragg planes. The photoemission intensity at these photon energies from an atom situated anywhere in or on the surface is therefore dependent on its location relative to the Bragg planes. The NIXSW measurements reported here were acquired, at almost normal incidence, by stepping the photon energy through the (111) Bragg reflection of Ir (∼2800 eV at 300 K). The relative intensity of the x-ray reflectivity as a function of photon energy was measured optically from a fluorescent screen mounted on the port through which the incident photons passed and was used to define the energy scale of the NIXSW measurements relative to the Bragg energy. Fitting of this reflectivity required the convolution of the theoretically calculated Darwin curve from the Ir(111) and from the Si(111) double-crystal monochromator with a Gaussian line shape with a width of between 0.20 and 0.25 eV. This Gaussian broadening models imperfections in the monochromator and the mosaicity of the Ir single-crystal substrate. The NIXSW absorption profiles of the overlayer constituent atoms were acquired by monitoring the variation in the intensity of the synchrotron x-ray photoelectron spectra (SXPS) from the C 1s and Co 2p 3/2 core levels using a VG Scienta EW4000 HAXPES hemispherical electron analyzer with an angular acceptance range of 56°. The analyzer was mounted such that the angle between the center of its angular acceptance and the incident light was 90°, in the horizontal plane of the photon polarization. For the NIXSW measurements the integrated intensity over all the collected emission angles between the photon polarization direction and the photon incidence direction was used. The backward-forward asymmetry parameter Q, due to nondipole effects in the photoelectron angular distribution (as defined in Ref. [46]), which is required to analyze the NIXSW data, was calculated using an intensity-weighted average value for this angle of 18°. Higher-resolution Co 2p, C 1s, and Ir 4 f SXPS data, used to characterize the surface, were acquired at photon energies of 950, 400, and 176 eV, respectively, and the absolute energy calibration performed by comparison against a Fermi energy measurement performed at the same photon energies. The Co 2p 3/2 and C 1s spectral peaks recorded in the NIXSW measurements were fitted with a convolution of a Gaussian and a Doniach-Šunjić line shape [47]. The resulting NIXSW profiles can be uniquely fitted by two parameters, the coherent fraction ( f ) and coherent position (p) [48]. The coherent position corresponds to the height above the extended bulk-scattering plane immediately below the adsorbing atoms in units of the spacing of these bulk-scattering planes. For a perfectly ordered structure, with all absorbing atoms at the same height and no thermal vibrations, the value of the coherent fraction is unity. Lower values can be assigned to static or dynamic disorder or, for values less than ∼0.6-0.7, to a multiplicity of co-occupied heights. The coherent position is a measure of the average height of the absorbing atom relative to the extended bulk crystal scattering planes.
A clean Ir(111) crystal was prepared by repeated cycles of Ar + sputtering and annealing to 1375 K, monitored by a pyrometer (CellaTemp PA 35, λ = 0.82-0.93 μm, emissivity of 0.3), with heating achieved by e-beam bombardment (V = 1 kV, I em = 19 mA). A monolayer of graphene was prepared in situ by briefly exposing the clean Ir(111) sample at room temperature to 5 × 10 −8 mbar of ethylene via a high-precision leak valve mounted on the chamber, before annealing the sample, in vacuum, to 1300 K. Once at this temperature the pressure of ethylene was increased to 2 × 10 −7 mbar and maintained in this condition for 8 min. Finally, the sample was annealed in vacuum for a further 5 min at 1400 K. This resulted in a sharp low-energy electron diffraction pattern of a moiré structure, similar to that already presented in the literature (e.g., Ref. [49]). Co was deposited onto the sample using an Omicron EFM-3 e-beam evaporator. Approximately 1.7-Å average thickness of Co [equivalent to approximately 80% of a single layer of Co(0001)] was deposited, using a water-cooled quartz microbalance to monitor the evaporation. Intercalation of the Co was achieved by annealing the sample, mounted on a pyrolytic boron nitride heater, to the relevant temperatures (see below) in UHV.

B. Computational details
Our ab initio spin-polarized electronic structure calculations were performed using the DFT [50,51] and the projector augmented-wave method [52] as implemented in the VASP code [53][54][55]. To describe the van der Waals interactions present between the graphene and Co/Ir(111) substrate we employed (1) the semiempirical method DFT-D2 [3] with the C6 coefficient for Ir as evaluated in Refs. [5,56] and (2) the nonlocal correlation-energy functional vdW-DF2 [57] together with a reoptimized [24] Becke (B86b) exchangeenergy functional [58]. The DFT-D2 calculations are the same as those presented in our prior work [56]. The moiré pattern of the graphene/Co/Ir(111) system was modeled by a (10 × 10) graphene layer adsorbed onto a (9 × 9) Co/Ir(111) with an in-plane unit cell of 21.93 × 21.93 Å 2 containing one layer of Co and three layers of Ir with a vacuum region of ∼20 Å. The structural relaxations were performed with a cutoff energy of 500 eV by relaxing the atomic positions of the graphene, Co monolayer, and first Ir surface layer with a threshold value of the calculated forces of ∼3 meV/Å using gammapoint sampling of the Brillouin zone, while the charge-density differences were obtained for 10 k points in its irreducible part. The spin polarization above the Co/Ir(111) surface and the out-of-plane magnetization direction is reproduced when using a U = 4 eV in our calculations [56]. The total binding energy per C atom was found to be 102 meV for the DFT-D2 correction and 105 meV for the DFT-vdW2-rB86 functional.

A. Synchrotron x-ray photoelectron spectroscopy
The fabrication of the Co intercalated graphene was monitored using synchrotron x-ray photoelectron spectroscopy and normal incidence X-ray standing waves in the following steps: after pristine graphene formation on Ir [Gr/Ir(111)] after room-temperature deposition of Co [Co/Gr/Ir(111)], after annealing to 600 K [Co/Gr/Ir(111) + 600 K] and after annealing to 750 K [Co/Gr/Ir(111) + 750 K]. C 1s SXP spectra from all four steps are shown in Fig. 1(a). The first and last of these spectra agree well with the results of Vita et al. [38], Pacilé et al. [59], and Avvisati et al. [49]. The pristine graphene film exhibits a single C 1s peak at a binding energy of 284.1 eV (C graphene ); after intercalation two clearly resolved peaks are observed at 284.4 (C weak ) and 284.9 eV (C strong ). These two clearly resolved peaks exhibited an intensity ratio of (2.1 ± 0.3):1 throughout the annealing process. In the work of Pacilé et al. [59] these peaks were assigned to, respectively, areas of the Gr layer that are weakly bound to the intercalated Co and areas that are strongly bound to the intercalated Co. This assignment was supported by model DFT calculations of core-level shifts for graphene on Co(0001) at different heights [49]. The origin of the two distinct bonding regimes stem from the different lattice parameters of the Gr and the intercalated Co monolayer resulting in a moiré superstructure, as is observed in the majority of systems where graphene is adsorbed on a metal substrate (e.g., Refs. [49,[60][61][62]). In such a superstructure, three distinct registry subregions can be identified, generally labeled as hcp, fcc, and atop [60]. They are identified by the lateral adsorption sites in that subregion that an individual graphene ring is above (or, alternately, the lateral adsorption site in that subregion above which no C atom is found). Typically, in theoretical modeling of such systems, it is the atop subregion that is predicted to be most distant from the surface [56,61,62]. The SXPS from the Ir 4 f and Co 2p core levels, after intercalation, can be found Fig. S2 in the electronic Supplemental Material (ESI) [63]. No changes are observed in the binding energy of the Co 2p XP spectra, although as the Co intercalates the intensity does increase. We interpret this counterintuitive behavior as indicating that the Co above the graphene, prior to intercalation, is in the form of three-dimensional clusters that self-attenuate the Co 2p signal, whereas after intercalation the Co is in the form of a single two-dimensional layer so the Co 2p signal is only attenuated by the graphene layer.  (3) 3.40 (4) 4.07 (2) 4.77 (4) 2.06 (3) 2.76 (5) (c) C graphene C strong C weak Co (b) NIXSW from the three components shown in the C 1s SXPS data, as well as that from the Co 2p 3/2 core level (Co 2p SXPS shown in ESI, Fig. S2 [63]). The C graphene NIXSW was obtained prior to Co deposition; the C strong , C weak , and Co data were obtained after Co intercalation at 600 K. (c) shows schematically the resulting mean adsorption heights of all three components, indicating the topmost Ir surface atoms (grey spheres), the graphene layer (solid black line), and the intercalated Co atoms (purple spheres).

B. Experimental structure determination-normal incidence x-ray standing waves
The NIXSW absorption profiles from the Ir(111) bulk reflection, obtained from the three different components of the C 1s and the single Co 2p 3/2 component, are shown in Fig. 1(b), and the resulting layer structure, indicated by the NIXSW fitting parameters (Table I), in Fig. 1(c). The corrugated graphene layer is represented by two distinct layer spacings and the average heights are given by the coherent positions. Note that coherent positions in NIXSW are a measure of the height of an absorbed species relative to the scattering planes defined by the extension of the underlying bulk lattice and not relative to the outermost atomic layer of the substrate; the technique is "blind" to surface-layer relaxations. As a result, the height of adsorbed atoms above the Ir surface may include a systematic error due to this effect [although in general fcc (111) surface relaxations are no more than a few hundredths of an ångstrom]. However, this effect can have no impact on the measured difference in height between the intercalated Co and the graphene layer above it.
Above the projected bulklike termination of the substrate, the intercalated Co layer has an adsorption height of 2.01 ± 0.03 Å, which is 0.21 Å less than the Ir(111) bulklayer spacing. The graphene, prior to intercalation, has an adsorption height of 3.40 ± 0.04 Å, a result very similar to that found by Busse et al. [5]. The two C species, C strong and C weak , are found at average adsorption heights of 4.07 ± 0.02 Å and 4.77 ± 0.04 Å, yielding differences of height, with respect to the intercalated Co layer, of 2.06 ± 0.03 Å and 2.76 ± 0.05 Å. The coherent fraction values were 0.85 ± 0.10 and 0.40 ± 0.10 for C strong and C weak , respectively, implying that C strong is associated with a much narrower range of heights than C weak .
The NIXSW measurements were taken as a function of the annealing temperature used to intercalate the Co. A subsection of the absorption profiles is shown in Fig. S1 in the ESI [63]. They show that the structure of the C graphene , C strong , and C weak is effectively invariant throughout the intercalation process, with only the relative areas covered by the three components varying [as seen in the C 1s SXPS in Fig. 1(a)]. However, there is a significant variation in the Co 2p NIXSW profile. This variation is due to the coherent fraction of the Co species increasing with annealing temperature, with the lowest coherent fraction found after room-temperature deposition, when only a small amount of Co intercalation is observed in the C 1s SXPS. Figure 2 shows the variation of the coherent fraction of the Co atoms as a function of the annealing temperature, juxtaposed against the relative area of the C 1s SXP spectra   2. (Red) Coherent fractions obtained from fitting the XSW profiles from the Co 2p 3/2 core level as a function of the annealing temperature, compared to (blue) the off-Bragg intensity of the two peaks in the C 1s XP spectra (Fig. 1) related to the graphene above the intercalated Co normalized to the total C 1s off-Bragg intensity.
that corresponds to graphene intercalated by Co. This further supports our interpretation that Co atoms originally above the Gr layer are intercalating below it and that the Co atoms above the Gr layer are clustered, resulting in a low coherent fraction. Interestingly, as the intensity of the intercalated graphene component in the C 1s SXPS approaches saturation, the Co 2p coherent fraction continues to increase, suggesting that an excess of Co was deposited on the surface, and annealing to elevated temperatures is driving some of the excess Co subsurface [64].

C. Calculated structure-density-functional theory
In Figs. 3(a) and 3(b), maps of the predicted adsorption heights of graphene above intercalated Co can be seen for a monolayer (1 ML) of intercalated Co with both the semiempirical DFT-D2 [ Fig. 3(a)] and the vdW-DF2-rB86 [ Fig. 3(b)] functionals. A similar map for graphene above Ir(111) without intercalated Co can be found in the ESI (Fig.  S3) [63]. Both calculations predict that the majority of the carbon atoms are bound strongly to the Co (adsorption height <2.2 Å) in the fcc and hcp subregions, interspersed with weakly bound regions (adsorption height >3 Å) in the atop subregion. However, there is a striking difference between the semiempirical DFT-D2 and the vdW-DF2-rB86 functional in the predicted structure of this layer in the regions that lie in the shortest distances between two weakly bound regions, i.e., in the boundary area between the fcc and hcp subregions. In the semiempirical DFT-D2 calculation these regions exhibit a local maximum in adsorption height, whereas in the calculations utilizing the vdW-DF2-rB86 functional, these regions are predicted to have the overall lowest adsorption heights. This effect is further illustrated in Figs. 3(c) and 3(d) that show the regions of the minimum adsorption heights for both dispersion corrections in an enhanced color range. Note that this difference in structure is not due to local energetic minimum, as the vdW-DF2-rB86 calculations were started in the same configuration as the energy-minimized DFT-D2 structure, and the energy minimization of the vdW-DF2-rB86 calculation resulted in this difference. These rescaled regions also highlight another significant feature of both predicted structures of the Co intercalated graphene, that the most strongly bound C species exhibit notable nearest-neighbor bucking (∼0.05 Å), in contrast to the more smoothly varying corrugation in the weakly bound areas observed in Figs. 3(a) and 3(b) and highlighted in Fig. S4 in the ESI [63].
To compare the NIXSW results with these structures predicted by the DFT calculations, we must assign the C strong and C weak species to a particular range of heights in the structural models. Figure 4(a) shows histograms of the predicted carbon atom height distributions obtained from the two different DFT   Table I. In both cases the adsorption height is relative to that of the underlying Co layer.
calculations. At first glance these seem difficult to reconcile with the XP spectra that are apparently dominated by two components, C strong and C weak , corresponding to C atoms in the corrugated layer have short and long Co-C bonding distances. Each histogram shows a clear group of closely spaced peaks for C heights less than ∼2.2 Å that might be expected to correspond to the C strong atoms, but at larger C-Co spacings the distributions show only long tails with no distinct peaks. With hindsight, however, this seems less surprising. It is clear that at large separations, between the C atoms and the Co layer, the photoelectron binding energy shift must be weakly dependent on that separation (at sufficiently large separation there is no dependence), whereas at small separations this dependence must be much stronger. Qualitatively, at least, the assignment of C strong to a relatively narrow distribution of heights, and C weak to a wide distribution of heights, is therefore to be expected. Of course, in reality, there must be many different values of the XP chemical shift, and thus many component peaks, but the clustering of these into two main groups leads to a two-peak fit to the spectra. There is no meaningful way of achieving a unique multipeak fit to the experimental spectra without the use of constraints based on other information. The ideal solution to this problem would be to use DFT calculations to predict the expected chemical shift for each C atom in the unit mesh, but for this very large unitmesh structure containing many symmetrically inequivalent C atoms, such a set of calculations is unrealistic in terms of the computational resources required.
However, to aid this assignment of the different C atoms in these histograms to C strong and C weak we have one key experimental datum, namely that the XP intensities of these two components are in the approximate ratio of 2:1. Moreover, XP peak intensities are known to be broadly linear in coverage, especially at the high photoelectron kinetic energies that were used in this study, a fact that is the basis of the extensive use of XPS for surface analysis in a wide range of materials problems. We have therefore assumed that 66% of the C atoms with the lowest adsorption heights must be assigned to C strong and the remaining 34% to C weak ; the resulting cutoff heights separating the two species are superimposed on the histogram of Fig. 4(a). These assignments allow us to determine, for each of the two structural models, the predicted NIXSW coherent positions and coherent fractions for C strong , C weak , and the intercalated Co with respect to the projected bulklayer spacing of the underlying crystal for all the calculated models; these are shown in Table I. Clearly the best match to the experimental values is found for the DFT-D2 calculations. Moreover, no possible alternate occupation ratio assignment of C strong and C weak assignment (from 1%:99% to 99%:1%), when applied to the structural results of the vdW-DF2-rB86 calculations, can match the experimentally measured adsorption heights. This is demonstrated by the results shown in Fig. 4(b), with the experimental results superimposed. For the DFT-D2 calculations the adsorption heights for both species agree well with experiment when the fraction of C atoms attributed to C strong is between 60 and 70%. For the vdW-DF2-rB86 calculations, no value of this fraction is consistent with the measured height of the C strong atoms.

D. Covalent bonding and origin of the difference between the semiempirical and ab initio van der Waals methods
All of the systems listed above, in which the vdW-DF2-rB86 functional performed as well as, or better than, semiempirical DFT-D2 corrections, are typified by much weaker bonding than what is observed in this system. Thus, the bonding motif is likely to be dominated either by van der Waals interactions or by an interplay between van der Waals interactions, electrostatic forces, and weak metalligand bond formation. By contrast, the interaction between the Co intercalation layer and the graphene in the C strong region is far stronger and the adsorption height much smaller (2.06 ± 0.03 Å). This could be indicative of a covalent bond formation between the C strong region and the underlying Co layer. The predicted C-C buckling in this region (see Fig. 3) may also indicate that the interaction with the substrate is sufficiently strong to partially lift the graphene sp 2 hybridization. To probe the possibility of a covalent interaction between the Co and the C atoms, additional calculations were undertaken to determine the charge redistribution between the Co and graphene layer and the projected density of states (PDOS) of specific C and Co atoms. The results of the DFT-D2 calculations of this system, focusing on charge transfer from the Ir(111) into the graphene layer, the induction of ferrimagnetism into the graphene layer, and a net charge transfer from the Co to the strongly bound graphene C atoms, have been previously published by some of us [56]. However, this work did not focus on the local charge redistribution and, specifically, did not discuss the behavior in the areas between the strongly bound regions, where we find the main differences between the results of the DFT-D2 and vdW-DF2-rB86 calculations to occur. Figure 5 displays a series of charge redistribution isosurfaces calculated using the semiempirical DFT-D2 and the vdW-DF2-rB86 functionals. Figures 5(a) and 5(d) display the charge redistribution specifically in the graphene layer, showing a depletion between the carbon atoms, slightly below the carbon plane, and a significant accumulation above and below alternate carbon atoms that are directly atop a cobalt atom. Figures 5(b) and 5(e) display the charge redistribution between the C and the Co atoms, for the C atom in an atop site, showing charge depletion around both the C and the Co atoms, and accumulation between the C and the Co and around a set of Co orbitals that seem to have an overlap with the p z orbitals of the C atoms. Thus, the relaxed geometries obtained using both van der Waals correction schemes predict a sharing of charge density between the C and Co atoms, at the expense of the bonding between the C atoms, consistent with a C-Co covalent interaction. In the area of the corrugated Gr where the two methods differ significantly, shown in Fig. 3, the local adsorption site of the C atoms has the Co atom lying directly below a C-C bond. The charge redistribution isosurface in this region, between the two C atoms and the Co atom, is shown in Figs , with a similar depletion in the electron density in the C-C bonding, and an accumulation between the Co and the C atoms. However, in the semiempirical DFT-D2 calculations for this same region [ Fig. 5(c)], a rather more asymmetrical redistribution is observed. The depletion around the Co atom remains broadly the same, but the depletion between the C atoms is much weaker and there is only significant accumulation between one of the C atoms and the Co atom. Looking at the PDOS in this region (Figs. 6 and 7) we see that there is little difference in the energetic distribution of the density of states between the two calculations for the Co atoms [Figs. 6(a) and 6(b)], but that the C p orbitals [Figs. 6(c) and 6(d)] differ significantly between the DFT-D2 and vdW-DF2-rB86 calculations. Looking at the overlap between the C p z and Co p orbitals [ Fig. 7(a)], it can be seen that the differences in the C p z orbital result in the vdW-DF2-rB86 calculations predicting a greater overlap in the density of states of the Co p orbitals below the Fermi edge (FE) at −0.4 [emphasized in Fig. 7(b)] and −3 eV. Though the density of states around the FE of the C p x orbital is around a factor of 10 weaker, there is also an improvement in overlap between this orbital and the Co p orbitals in the vdW-DF2-rB86 calculations [ Fig. 7(c)]. In the fcc, hcp, and atop subregions there are only minor differences between the DFT-D2 and vdW-DF2-rB86 calculations (see Figs. S7-S9). The full spin-polarized PDOS for the Co p-and d orbitals as well as the C s-and p orbitals are shown in Figs. S10-S21 in the ESI [63].
Both the charge-redistribution maps and the PDOS maps suggests that the vdW-DF2-rB86 calculations predict a greater covalency in the boundary area between the fcc and hcp subregions. We should stress that we do not suggest that the 124001-7 charge redistribution leads to the structural difference, nor that the structural difference leads to the charge redistribution; the two are inextricably linked. However, differences in the charge redistribution do give some insight into differences in bonding character that accompany the differences in structure. Specifically, we posit from these results that the difference in the predicted structures from the semiempirical DFT-D2 and vdW-DF2-rB86 calculations would appear to be related to how well the two dispersion corrections describe the balance between the local covalent and nonlocal dispersion interactions between the graphene and cobalt layer.

IV. DISCUSSION AND CONCLUSIONS
The generalized gradient approximation-Perdew Burke-Ernzerhof (GGA-PBE) functional typically gives an accurate description of covalent chemical bonds [65][66][67][68][69], and, within the semiempirical DFT-D2 method, the dispersion correction is added directly to the GGA-PBE Kohn-Sham energy. However, in the case of the vdW-DF2-rB86 calculation, to avoid a double counting of the semilocal contributions present in the GGA approximation, the exchange-correlation energy is calculated in a different way. Specifically, the correlation energy is calculated by a sum of the local-only correlation energy at the LDA level and the vdW-DF2 nonlocal contribution. Additionally, in the vdW-DF2-rB86 calculations, the exchange energy is calculated from a reoptimized Becke (B86b) exchange functional [24], which gave excellent agreement with the previously benchmarked (by NIXSW measurements) DFT-D2 calculations for graphene on Ir(111) [5]. This approach should give a very good description of the nonlocal contributions (e.g., the van der Waals forces), but these differences in the vdW-DF2-rB86 functional may be the origin of the inaccuracy in describing a system, such as Co intercalated graphene, that has such a coexistence of covalent and van der Waals dom-inated regions. Conversely, utilizing a pure GGA functional would completely fail to describe this system, as the GGA functional would clearly be unable to model the large van der Waals dominated regions. Yet, this interpretation is still too simple, as even in the covalent regions of the graphene layer, the van der Waals forces play a role; excluding them results in a repulsive interaction between the Co and the graphene [56]. Instead, it is the contributions from the two interactions, covalent and van der Waals, that must be balanced correctly, and failing to model either one properly will result in the failure to model the whole system accurately.
We are not aware of any other system, published thus far in the literature with stringent benchmark experimental parameters, that presents the same sort of challenge and has been modeled with the vdW-DF2 family of functionals. Perhaps the closest comparison would be the adsorption of h-BN on Ir(111) [10]. In this system there is also a coexistence of strongly and weakly bound regions, although with much smaller areas of strongly bound regions than in the C-Co-Ir system of our study. However, the h-BN does exhibit a short adsorption height on Ir(111), and a slight buckling between the B and N atoms. On the other hand, the charge redistribution map in the strongly bound region, despite being far closer to the substrate, is remarkably similar to that of graphene on Ir(111) [5], which, as discussed previously, occupies an adsorption height even greater than the weakly bound region of the C-Co-Ir system. Specifically, comparison of Fig. 2(e) of Ref. [5] and the inset of Fig. 2(c) of Ref. [10] yields a similar oscillation in charge redistribution between the surface Ir atom and the C/N atom directly above it, with both studies exhibiting a node in that oscillation approximately halfway between the Ir and C N atom. This is in complete contrast to the results of the calculations presented in this study [see Figs. 3(b) and 3(e)].
In the specific case of the C-Co-Ir system, the strong bonding between the majority of the carbon atoms to the cobalt atoms, and the apparent covalent nature of this interaction, suggests that the sp 2 hybridization of the graphene layer is at least partially lifted. This may suggest that these C strong regions may have some graphanelike characteristic. The measured adsorption height, to the Co layer, of 2.06 ± 0.03 Å, is even shorter than that found between the sp 3 regions of hydrogenated graphene and the Ir(111) surface [70]. However, the predicted neighbor-neighbor buckling of hydrogenated graphene (∼0.4 Å) [71] is much greater than that predicted here (∼0.05 Å). This suggests that, although the sp 2 character of the graphene may be weakened by the Co intercalation, the hybridization of the C orbitals is certainly not fully converted to sp 3 character. However, the present covalent interaction between the Co and C atoms suggests that Co intercalated graphene may be a good target for hydrogenation, and might result in larger hydrogenated regions than were observed for hydrogenation of graphene on Ir(111), without intercalant [42,43].
The system studied here exhibits both a local and a global balancing of covalent and van der Waals interactions. On the global scale we have a periodic array of strongly and weakly bound regions in which covalent and van der Waals forces, respectively, dominate. This periodic array will result in a balance between the strongly bound region being as close to the substrate as possible, and the weakly bound region as distant from the surface as possible, probably putting significant strain on the graphene network. Locally, the adsorption of each individual C atom is determined by a balance between the attractive covalent forces, the repulsive van der Waals forces, and the attractive van der Waals forces. We posit that it is this mixture of strong covalent interactions and weak van der Waals interaction, both local and global, which is proving to be a challenge for the vdW-DF2-rB86 functional. Therefore, this system could present a class of systems that vdW-DF2-rB86 calculations fail to model accurately, demonstrating the pressing need for further reliable quantitative experimental structural benchmarks against which to test such computational methods.