Computational study of adsorption of cobalt on benzene and coronene

We investigate the binding of the cobalt atom on small aromatic model systems as a proxy for interaction with graphene, using density functional theory, coupled-cluster theory, and combinations of them using projector-based quantum embedding. We set out in some detail the electronic structure of the cobalt atom alone, because some nuances of atomic structure appear to have been overlooked in previous studies. Two states of the complex in particular are studied: those formed from the a 4F ground state of the atom; and from c 2D, the lowest doublet state with configuration 3d94s0. We highlight the difficulties in extracting reliable results from typical approximate density functionals, and demonstrate that embedding calculations using the coupled-cluster theory in an active subsystem greatly reduce functional dependence, and produce a picture more consistent with the available experimental information. Our results cast doubt on previous calculations that have predicted strong chemisorptive binding between graphene and the c 2D state of cobalt.


Introduction
Since its first isolation in 2004, graphene has been the subject of intense scientific interest [1,2]. The peculiar structure, characterised by a regular hexagonal sheet of carbon atoms, each with three sp 2 bonds to its neighbours, gives the material exceptional electronic and physical properties [3][4][5][6]. The carbon p z orbitals give rise to a valence band of π and a conduction band of π * orbitals. These bands touch at the Brillouin-zone corners, making graphene a zero-band-gap semiconductor, and leading to high chargecarrier mobility, a quantum Hall effect, and massless relativistic carriers (see for example [7]).
Many groups are actively researching ways to tune electronic and magnetic properties of graphene-based materials to meet requirements for the design of new devices. Introduction of adsorbates (or dopants in general) [7] can help to modulate such properties, for example, through the alteration of the effective dielectric constant [8], shifting of the chemical potential [9], and opening the band gap [10]. The last can be achieved in at least three ways: by confining graphene in one dimension to form nanoribbons, by forming bilayers, or by applying strain to graphene [10].
Adsorption of transition metals (TMs) on graphene surfaces is relevant for the development of various technological applications (e.g. materials for spintronics [11][12][13]), because the TMs induce localised magnetic moments in graphene [14]. In particular, cobalt has been studied extensively in this context. For example, both theoretical [15] and experimental studies have shown that an appropriate tuning * Corresponding author. Email: fred.manby@bristol.ac.uk of the electronic structure by cobalt on graphene can give rise to physical phenomena such as a Kondo effect [16][17][18]. This also opens the prospect of the development of novel analytical devices such as chemical sensors [19].
Theoretical modelling of such metal-graphene systems can provide insight to the nature of the binding between them and the electronic states involved. Many ab initio calculations have been performed on various metals approaching graphene, in both finite models and periodic boundary conditions [20][21][22][23][24][25][26][27]. Kohn-Sham density functional theory (KS-DFT), in particular, has often been used because of its favourable scaling with system size.
However, modelling the open-shell electronic structure of TMs adsorbed on a surface can be a challenging task [28]. These systems typically have many neardegenerate states, and may show strong electron-correlation effects. Approximate DFT descriptions of cobalt approaching a graphene model can be strongly functional dependent [29,30], and correlated wavefunction-based methods are probably needed to accurately analyse the energetics of such systems. Although wavefunction-based methods can provide a more consistent description [30,31], the computational costs typically rise very quickly with system size.

Electronic structure of the cobalt atom
In order to set out a clear picture of the electronic structure of cobalt on finite models of graphene, we first rehearse the key facts about the electronic structure of the atom alone.  The ground state a 4 F is dominated by the configuration 3d 7 4s 2 , and a second quartet, b 4 F(3d 8 4s 1 ), lies around 0.51 eV higher in energy [32]. These two configurations also give rise to 4 P states that appear at energies 1.7-2.0 eV above the ground state. The 3d 7 4s 2 and 3d 8 4s 1 configurations also give rise to a succession of seven doublet terms, lying in an energy range 2.0-2.8 eV above the ground state. However, previous DFT studies of cobalt adsorbed on graphene have focused on the much higher energy c 2 D doublet state that arises from the 3d 9 4s 0 configuration, because this state is found to bind very strongly to the surface using certain approximate functionals. A key purpose of this paper is to determine whether the chemisorption of the c 2 D state predicted using DFT is feasible.
The electronic structure of the complex is intricate, because each term leads to a multitude of different states. In the relatively simple, but related example of the CoH molecule, multi-reference configuration interaction calculations find the ground-state (a 4 F) term of the cobalt atom to form 3 , 3 , 3 , 3 − , 5 , 5 − , 5 , and 5 states of the diatom [33]. A similar profusion of states can be expected to exist in the interaction of cobalt with graphene.
The difficulties in mapping this out using the DFT are threefold: first, one is often interested in states that are not the lowest of a given symmetry; second, single-determinant methods can describe some, but not all of the multiplets that arise; and third, there is the obvious issue of accuracy of the approximate exchange-correlation functional. Wavefunction methods aimed at ground states -including coupled cluster (CC) theory -suffer from at least the first of these problems.
Using conventional ground-state mean-field theories, the accessible energies are those of the lowest states for each term symbol. The a 4 F ground state therefore presents no problem. Higher states with the same spin multiplicity and spatial symmetry are sometimes accessible in meanfield calculations by fixing orbital occupancies, in what are essentially delta self-consistent field ( -SCF) calculations [34].
Even when a -SCF calculation is possible, the excited state of the complex does not necessarily correspond to a specific term of the cobalt atom in the dissociation limit. For example, quartet states arising from the configuration 3d 8 4s 1 include the b 4 F state. For the m = ±3 components of this term, there are valid single-determinantal wavefunctions. But for other components, single determinants are unphysical mixtures of the b 4 F and b 4 P terms. These same considerations apply to the doublet states lying between 1.73 and 2.87 eV. They are found very close in energy and with low-spin configurations of either 3d 7 4s 2 or 3d 8 4s 1 . Even imposing a specific occupancy, we cannot confidently discriminate among these states with single-determinant approaches.
As noted above, the b 4 F (3d 8 4s 1 ) state is experimentally observed to lie 0.51 eV above the ground state, a 4 F. In approximate DFT calculations, the ordering of the states is typically reversed, with the lowest quartet state arising from 3d 8 4s 1 being 0.13 and 0.60 eV below a 4 F with B3LYP and PBE, respectively.
Previous studies have implicated the c 2 D state of the atom in strong binding to graphene [20,21,29,35]. This state is observed 3.45 eV above the atomic ground state, is the lowest doublet arising from the 3d 9 configuration, and is well separated from other doublet states. Approximate DFT gives a qualitatively reasonable splitting (2.80 and 1.96 eV with B3LYP and PBE, respectively). Although a Hartree-Fock (HF) calculation can be converged on this state (by fixing the correct occupation), there is no a priori reason why a CCSD(T) calculation based on this excited-state reference should behave reasonably. However, what we find is that the CCSD(T) splitting (3.47 eV) is in excellent agreement with experiment. The T 1 diagnostic is found to be around 0.02 for both states, just within the bounds where single-reference methods can be expected to perform well [36].

Projector-based embedding method
As shown in the previous section and elsewhere [31], an accurate description of the electronic structure of cobalt requires methods beyond the accuracy of commonly used approximations in DFT. On the other hand, because of their computational cost, the use of accurate wavefunction-based electronic structure methods is limited to small system sizes. Quantum embedding schemes [37] can be used to bridge this gap, providing better accuracy than DFT alone at lower computational cost than wavefunction-based quantum chemistry on the whole system.
Here, we will assess the use of projector-based CCSD(T)-in-DFT embedding calculations [38][39][40]. DFT-based embedding schemes are constructed by partitioning the electronic density into separate contributions that define subsystems A and B [41][42][43] The energy then partitions into contributions for each subsystem and a coupling term that contains all of the nonadditive effects All terms in Equation (2) can easily be evaluated in the framework of conventional, approximate KS-DFT, except for the kinetic-energy contribution to the nonadditivity term, This either has to be treated with approximate kinetic-energy functionals [44][45][46][47], or through potential inversion methods [48][49][50]. The former is not yet sufficiently accurate for partitions across covalent bonds, and the latter presents extreme numerical challenges.
Projector-based embedding simplifies things by focusing on partitions of the density for which the nonadditive kinetic energy, T s [ρ A , ρ B ], is exactly zero. In the present implementation, this is done by first performing a KS calculation on the whole system, localising the occupied orbitals (here using the Pipek-Mezey scheme [51]), and partitioning the local orbitals into two subsets: The densities of the two subsystems are defined in terms of these orbitals, and the kinetic energy of the whole system is then just the sum of kinetic energies for each subsystem.
Next, we form a HamiltonianĤ A in B for subsystem A that includes all interactions between subsystem A and B, and a projection operator that level shifts the orbitals in subsystem B to high energy. This ensures that any calculation on the electrons in subsystem A satisfies the Pauli principle, by blocking access to the occupied orbitals of subsystem B.
A KS calculation usingĤ A in B reproduces exactly the orbitals {φ i } A , but instead a high-level wavefunction-based calculation can be performed in subsystem A; the wavefunction (WF)-in-DFT energy is defined as where E B is the DFT energy of subsystem B. Further details on the method can be found in reference [38]. The code is implemented in the development version of the Molpro package [52,53], which was used for all calculations presented in this work.

Interaction of cobalt and benzene
We start by computing binding curves of C 6v complexes of cobalt and benzene, a system for which CCSD(T) on Distance Co-benzene (Å)  Figure 2. Binding curves of cobalt on benzene as function of distance measured relative to the centre of the ring, using DFT (top), CCSD(T) (middle), and CCSD(T)-in-DFT (bottom). The curves approaching zero at the dissociation limit refer to the a 4 F3d 7 4s 2 state of cobalt, and those tending to a higher dissociation energy refer to c 2 D3d 9 4s 0 . the entire complex is easily feasible. Previous DFT and WF-based computational studies [30,31,54] have identified the hollow site as being a stable binding site for cobalt on the substrate, although some experimental studies show evidence of binding at the top site as well [55]. Here, we Table 1. Binding energies of cobalt and benzene for the a 4 F and c 2 D states of the atom, calculated at 3.5 and 1.5Å, respectively. The energies (in eV) are computed using CCSD(T)-in-DFT with the specified functional. The integer N A specifies the number of electrons in subsystem A, which is treated at the CCSD(T) level; N A = 0 denotes the DFT calculation on the whole complex. The experimental binding energy is −0.34 eV, attributed to the bound c 2 D state [60].  [56]) and DFT, in particular using PBE [57] and B3LYP [58] as representative examples of generalised gradient (GGA) and hybrid functionals. All calculations are run using the cc-pVDZ basis for the substrate and cc-pVTZ for the cobalt atom. In the CCSD(T) calculations, cobalt core orbitals are kept frozen. We use a fixed geometry for the cobalt-benzene complex with C-C and C-H bond lengths of 1.42 and 1.09Å, respectively. Binding energies, computed relative to the groundstate energies of the fragments, are shown in the top two panels of Figure 2.
The top panel of Figure 2 illustrates the considerable functional dependence of approximate DFT calculations, and shows the very strong binding of the c 2 D state found using the PBE GGA and reported elsewhere for PBE and other GGAs [20,21,35,54,59]. In the middle panel, it can be seen that CCSD(T) also predicts binding of the c 2 D state, albeit to a much lesser degree. The B3LYP curves are qualitatively similar to CCSD(T).
We now consider CCSD(T)-in-DFT embedding calculations with the number of electrons in subsystem A (N A ) ranging from 33 (27 on cobalt and 6 π electrons) through to 63 (the total number of electrons is 69). The bottom panel of Figure 2 shows CCSD(T)-in-DFT calculations with 33 electrons in subsystem A. It can be seen that the difference between the two functionals is almost completely eliminated, and the behaviour of the curves is qualitatively similar to CCSD(T) on the whole complex.
The behaviour of the binding energy as a function of the number of electrons in subsystem A is shown in Table 1. The DFT binding energies of the complex formed with the a 4 F state of the atom differ considerably, with PBE in good agreement with CCSD(T), but B3LYP predicting essentially no binding. The CCSD(T)-in-DFT energies are much less functional dependent, even with small N A , and converge to within a few hundredths of an eV of the CCSD(T) value.
The DFT predictions for the c 2 D state differ more widely, this time with B3LYP giving a reasonable value and PBE strongly overbinding the complex, as noted above. Again, the embedding calculations almost completely eliminate the functional dependence, but for this state the binding energy converges less quickly, so that even with N A = 63 (out of 69) electrons there is a discrepancy of a few tenths of an eV relative to the CCSD(T) reference calculation. The reason for the slow convergence for the doublet case is not entirely clear and merits further investigation on different systems.
This preliminary set of calculations shows that the embedding scheme greatly reduces the functional dependence of DFT calculations on the cobalt-benzene complex. For the complex formed from the atomic ground state, we find the embedding scheme converges smoothly. However, quantitative accuracy is not achieved for the c 2 D state until CCSD(T) is performed on the entire cobalt-benzene complex. This is not particularly surprising because in this small system all parts of the substrate are in close proximity to the cobalt atom, and much closer for the chemisorbed c 2 D state than for the physisorbed a 4 F state.

Cobalt on coronene
We now consider binding of the cobalt atom above the centre of coronene (C 24 H 12 ). We use this molecular system as a prototype for a graphene sheet, where on the hollow site of the central ring, a cobalt atom is introduced as a local defect. Metal-coronene complexes are currently of experimental [61][62][63][64] and theoretical [65,66] interest, because the delocalised π electrons in coronene make it a suitable model system for studying the interaction of metals with extended systems. These papers show that many different metals (TMs and otherwise) can form bound ionic and neutral complexes with coronene.
Despite the extensive experimental literature on metalcoronene complexes, not much is specifically known about the cobalt-coronene system itself. Kandalam et al. [67] reported a joint theoretical and experimental investigation on Co m [coronene] and Co m [coronene] − complexes (with m = 1, 2), where the cobalt is found either bridging a C-C bond or above the hollow site of one of the outer benzene rings. They used anion photoelectron spectroscopy to obtain the spectra of the anionic species, and combinations of these spectra and GGA calculations to obtain information about the neutral. According to their work, for the singleatom anion complex the preferred binding is to a bridging site with cobalt in a triplet state. The single-atom neutral complex appears to prefer the hollow site with the TM in a doublet spin state, although as we will see, this may be a consequence of relying on GGA calculations.  Figure 3. DFT binding energy curves of cobalt along the hollow site of coronene as a function of the distance from the centre of the central ring. The curves approaching zero at the dissociation limit refer to the a 4 F3d 7 4s 2 cobalt configuration. The curves tending to a higher dissociation energies refer to c 2 D3d 9 4s 0 cobalt configuration. The curves are computed with PBE (dashed line) and B3LYP (solid line) functionals, employing cc-pVTZ basis set for Co and cc-VDZ for the benzene ring. Akima spline method [68] is used to fit the curves. We used a fixed geometry for the calculations with C-C and C-H bond lengths of respectively 1.42 and 1.09Å.
The approximate DFT potentials for the cobaltcoronene system (shown in Figure 3) are broadly similar to those found for benzene, but with some notable differences. First, the binding of the a 4 F ground state becomes even weaker (PBE) or disappears altogether (B3LYP). Second, the strength of binding to the c 2 D state is reduced (relative to benzene) in both functionals, leaving PBE predicting a complex bound by 1.4 eV, and B3LYP predicting no binding.
The two example functionals we are using give an inconclusive picture of binding in the cobalt-coronene complex. In the cobalt-benzene case, we found CCSD(T)-in-DFT embedding calculations to eliminate most of the functional dependence; therefore, we now attempt to model binding in the cobalt-coronene complex using this method.
Building on the benzene calculations in the previous section, we explore different embedding partitions, with N A = 57 and N A = 69. The electronic density for subsystem A (with N A = 69) is shown in Figure 4. The active regions are chosen to formally include the same orbitals as the ones included in the active regions for the benzene calculations.
In Table 2, we report the binding energies for the two states considered with the various methods used. The distance between the atom and the ring centre is 1.6Å for the doublet (corresponding approximately to the minimumenergy configuration), and to 3.5Å for the quartet (approximately the optimmum distance for the quartet of the cobalt-benzene complex).
As seen above, the DFT calculations give a picture of little or no interaction with the a 4 F state, and either In particular, neither functional exhibits binding, and difference between the predicted binding energies is reduced by almost an order of magnitude.
By treating the N A electrons of subsystem A at the CCSD(T) level, a key part of the dispersive interaction between the atom and the substrate molecule is captured. But the interaction between cobalt and atoms in subsystem B is treated at the approximate DFT level, and therefore does not account for long-range dispersion. As a simple test of whether this omission could affect our conclusions, we have added Grimme-type pairwise dispersion corrections Table 2. Binding energies of a 4 F and c 2 D states of cobalt on coronene. In the complex, the atom is a distance 3.5Å from the coronene centre for the quartet, and 1.6Å for the doublet. The integer N A specifies the number of electrons in subsystem A, which is treated at the CCSD(T) level; N A = 0 denotes the DFT calculation on the whole complex. The reference coupled-cluster calculation has not yet been possible. The final row (69 * ) denotes the 69-electron embedded calculation corrected with Grimmetype dispersion interactions between the cobalt and the noncentral atoms of coronene [69].  [69] between the cobalt and all coronene atoms apart from the central ring. This approach can only be considered indicative of the magnitude of the effect, because the dispersion interactions that are added are not perfectly neglected in the embedding calculation due to imperfect orbital localisation; and also because it is assumed that the C 6 parameters for interaction with cobalt are appropriate not only for the ground state but also for c 2 D. The corrections amount to a few tenths of an eV, and the corrected 69-electron result is shown in the final row of Table 2.
The corrected data give a broadly similar conclusion. It appears from these calculations that the a 4 F ground state of cobalt weakly binds to coronene with an atom-ring-centre distance of around 3.5Å. (It should be noted that we have performed only nonrelativistic calculations, and the spinorbit effect is of a comparable magnitude to the binding effect observed for a 4 F.) The c 2 D state does not form a bound complex with coronene. Hence, the prediction made elsewhere [20,25,59] of strong binding of the doublet state is almost certainly an artefact of the PBE functional (and other GGA functionals [21,29,30,35,67,70]). We have not used counterpoise corrections in this study. Even though these would be expected to slightly reduce the computed binding energies, they would be unlikely to change our overall conclusions.

Conclusion
The computational modelling of open-shell TM complexes still represents a very challenging task. Previous work has shown that in order to describe these systems correctly, accurate wavefunction-based methods are required [29][30][31]. Often the most affordable and popular approach -approximate DFT -struggles to give a reliable description of the energetics, and results can be profoundly dependent on the approximate exchange-correlation functional used.
In this work, we have applied the projector-embedding approach [38] to describe such systems with the accuracy required at affordable computational cost. The calculations presented underline a key benefit of the CCSD(T)-in-DFT framework: dependence on the choice of density functional is greatly reduced.
CCSD(T)-in-DFT calculations on the binding in the cobalt-benzene complex are almost independent of exchange-correlation functional even with rather small numbers of electrons in the subsystem treated with the coupled-cluster theory. However, the onset of binding as a function of N A is rather slow, appearing only when 57 (out of 69) electrons are treated at the CCSD(T) level.
The cobalt-coronene complex is quite large and we have not been able to perform CCSD(T) on the whole system. We have used CCSD(T)-in-DFT projector embedding calculations to reduce functional dependence and increase accuracy in computationally affordable calculations.
The embedded calculations clearly indicate that the c 2 D state of cobalt does not bind to coronene, contrary to the situation for the cobalt-benzene complex, and contrary to the prediction of PBE (and other GGA functionals). On the other hand, the a 4 F state appears weakly bound using dispersion-corrected embedding calculations with either functional. These results are in accord with very recent experimental work [55] which finds that, on highly ordered pyrolitic graphite, early TMs (i.e. Fe and Co) prefer to be weakly bound (physisorbed) with high-spin configurations.