An analysis of van der Waals density functional components: Binding and corrugation of benzene and C60 on boron nitride and graphene
The adsorption of benzene and C60 on graphene and boron nitride (BN) is studied using density functional theory with the non-local correlation functional vdW-DF. By comparing these systems we can systematically investigate their adsorption nature and differences between the two functional versions vdW-DF1 and vdW-DF2. The bigger size of the C60 molecule makes it bind stronger to the surface than benzene, yet the interface between the molecules and the sheets are similar in nature. The binding separation is more sensitive to the exchange variant used in vdW-DF than to the correlation version. This result is related to the exchange and correlation components of the potential energy curve (PEC). We show that a moderate dipole forms for C60 on graphene, unlike for the other adsorption systems. We find that the corrugation is very sensitive to the variant or version of vdW-DF used, in particular the exchange. Further, we show that this sensitivity arise indirectly through the shift in binding separation caused by changing vdW-DF variant. Based on our results, we suggest a concerted theory-experiment approach to assess the exchange and correlation contributions to physisorption. Using DFT calculations, the corrugation can be linked to the optimal separation, allowing us to extract the exchange-correlation part of the adsorption energy. Molecules with same interfaces to the surface, but different geometries, can in turn cast light on the role of van der Waals forces.
The development of methods that include van der Waals forces in density functional theory (DFT) has made it possible to describe an abundance of important material systems — like water,Kelkkanen et al. (2009); Møgelhøj et al. (2011) DNA,Cooper et al. (2008) molecular crystals,Nabok et al. (2008a); Berland and Hyldgaard (2010); Nabok et al. (2008b); Berland et al. (2011a) and metal-organic frameworksKong et al. (2009); Li and Thonhauser (2012) — at the electronic level even without use of empirical parameters. Such materials can be categorized as sparse matter, because of the essential role of the low (electron) density regions or “voids” in these materials. These regions make it paramount to go beyond semi-local approximations like the generalized-gradient approximation (GGA)Langreth and Perdew (1977); Gunnarsson and Lundqvist (1976); Perdew and Yue (1986); Becke (1986); Perdew and Wang (1992); Perdew et al. (1996); Armiento and Mattsson (2005) and properly account for van der Waals interactions. A good description of sparse matter is important for developing new nanotechnology. The softer van der Waals bonds makes it easy for fragments to move, like in the self assembly of molecules on surfaces Wyrick et al. (2011); Berland et al. (2009) or for moving pieces of nano-mechanical devices.Chan et al. (2001)
Sparse matter constitutes a research-frontier for DFT and DFT-based modelling.Klimeš and Michaelides (2012); Langreth et al. (2009) Practical and robust schemes to describe sparse matter are all fairly recently developed and several competing and complimentary methods exist. Two prominent approaches are vdW-corrected DFT and the use of non-local correlation functionals. The first approach adds a semi-empirical pair-potential correction accounting for van der Waals forces on top of semi-local DFT.Wu and Yang (2002); Grimme (2004, 2006); Grimme et al. (2010) There exist many variation thereof and refinements, like adjusting the pair-potentials by using charge-density input.Tkatchenko and Scheffler (2009) Non-local correlation functionals aim to stay fully within DFT keeping the density as the key variable. The van der Waals density functional (vdW-DF) methodRydberg et al. (2000, 2003); Dion et al. (2004); Langreth et al. (2005); Thonhauser et al. (2007); Langreth et al. (2009) for approximating the exchange-correlation energy has led to a series of such functionals that lack empirical input. The vdW-DF1,Dion et al. (2004) introduced in 2004, and vdW-DF2,Lee et al. (2010a) introduced in 2012, emphasise a constraint-based design of the density fluctuation propagator,Lundqvist (1967a, b, c) while the precursor for layered systemsRydberg et al. (2000, 2003) emphasise an anisotropic screening account. The vdW-DFs vanish seamlessly in the uniform limitDion et al. (2004) and exhibits the correct asymptotic limit between molecular dimers.Dion et al. (2004) The prediction of various sparse matter methods can differ quite much.
Hand in hand with the development and refinements of sparse matter methods, it is essential to identify experiments and quantum-chemistry calculations that serves to test and elucidate these methods for various physical effects and across length scales.Sernelius and Björk (1998); Perdew et al. (2012a); Johnson and Becke (2006); Berland and Hyldgaard (2010); Barask (1988); Perdew et al. (2012b) One approach is to construct a representative data set that account for many physical effects,Björkman et al. (2012a); Jurečka et al. (2006); GraÌfovaÌ et al. (2010); ŘezaÌč et al. (2011) like the S22 data setJurečka et al. (2006) or the more recent and larger S66.ŘezaÌč et al. (2011) The average performance of methods can then be compared to each other. Another approach is to identify particular systems and associated experiments that reveal several properties of sparse matter binding.Dobson and Gould (2012) An example is to calculate potential energy curves (PEC) of H on Cu(111) and compare with the experimental ones constructed from backscattering experiments.Lee et al. (2011, 2012) Since the balance between the different contributions to PEC change with separation, much information can be extracted.
In this paper, we show that the details of adsorption, especially the predicted corrugation for adsorbate dynamics is a discriminator for assessing the balance between exchange and correlation components in a non-local functional design. We substantiate our conclusion by comparing the results of calculations using vdW-DF1 and vdW-DF2 and associated exchange choices for a group of adsorbate systems that are closely related, yet the nature of the binding differs in a well-defined manner.
We study the adsorption of benzene and C60 on graphene and BN. We show how these systems allow us to span, in a systematic manner, a wide range of different contributions to sparse matter binding. First, the interface between a benzene and graphene (or BN) is very similar to that of C60 on graphene (or BN), both in terms of contact area and character of the local bonds. Yet C60 is bigger and the difference in adsorption is mostly induced by the increased van der Waals attraction of the C60 molecule. Second, boron nitride and graphene are both flat one-atom layer thick sheets with eight valence electrons per unit cell, but whereas graphene is a semi-metal, BN is an insulator. The BN sheet also has a somewhat different charge-density landscape than the graphene sheet.
We show that the corrugation is very sensitive to the binding separation and therefore very sensitive to the version of correlation and exchange in particular used in vdW-DF. Yet at a given separation to the surface, the corrugation is largely independent of the functional version used. The shorter binding separation for C60 adsorbants compared to benzene increase corrugation in the same indirect manner.
Based on our results and analysis, we suggest a combined theory-experiment approach that can lead to insight into the nature of sparse matter binding and the magnitude of the different terms in the exchange-correlation functional.
The plan of the paper is as follows: After the introduction we summarize vdW-DF method and functional versions and variants and outline computational details specific to this study. Section III presents our results for adsorption of benzene and C60 on graphene and boron nitride. In section IV we analyze the various components to this binding, both the importance of GGA exchange choice and vdW-DF correlation version. We also analyze the role of different density separations and density regimes in the non-local correlation of vdW-DF. Section V compares the corrugation of the different systems and for the different functionals. Section VI discuss the possibility of using the selected systems to assess the quality of non-local correlation functionals and sparse-matter methods in general. Finally, we summarize our conclusions.
Ii vdW-DF calculations
The vdW-DF method is designed to handle both short-range covalent bonds and long-range van der Waals interactions without introducing damping functions or empiricism. Its exchange-correlation functional combines exchange at the GGA level, with the combination of correlation in the local density approximation (LDA) and non-local correlation: . The non-local correlation is responsible for the van der Waals interaction.
The non-local correlation term vanishes in the uniform limit and takes the form of a six dimensional integral
The kernel can be expressed in terms of a universal kernel that depends on two dimensionless length scales and where and is a modulation of the local Fermi vector that account for the local response. By depending on gradients as swell as local density, this local response reflects a broader dependence on the density variation through a plasmon-pole approximation.Dion et al. (2004)
The universal kernel kernel can also be tabulated in terms of and . Fig. 2 shows the kernel as a function of for three selected . The parameter is a measure of how different the local response of the two density regions are.
The GGA exchange energy can be expressed as a modulation of the LDA exchange energy , as follows
Here is a measure of how rapidly the density is changing compared to length scale set by the local Fermi wave vector. In the small limit, the modulations of the GGAs go like , while in the large limit they differ significantly both quantitatively and qualitatively.
ii.1 Functional versions
Because the nature of the adsorption systems considered here differ in a systematic manner, it becomes interesting to study how the predictions of vdW-DF changes with correlation version and choice of exchange partner. For vdW-DF1 correlationDion et al. (2004) (the version of 2004), we will in, in addition to the canonical choice of revPBE, use C09,Cooper (2010) PW86r,Murray et al. (2009); Perdew and Yue (1986) and optPBEKlimeš et al. (2010) as exchange partner. For vdW-DF2 correlation,Lee et al. (2010a) we use C09, in addition to the canonical choice of PW86r.
Standard vdW-DF1 with the revPBE Zhang and Yang (1998) exchange functional as companion to its non-local correlation have shown an impressive robustness, giving good results for binding between molecules of various sizes and properties,Dion et al. (2004); Thonhauser et al. (2007); Langreth et al. (2009); Puzder et al. (2006); Cooper et al. (2008); Chakarova-Käck et al. (2010); Kelkkanen et al. (2009); Møgelhøj et al. (2011); Hamada (2010); Berland and Hyldgaard (2010); Berland et al. (2011a) adsorption onto insulating surfaces,Viitala et al. (2012); Moses et al. (2009); Lee et al. (2010b); Johnston et al. (2008) graphene,Chakarova-Käck et al. (2006); Berland et al. (2011b); Londero et al. (2012) and metals,Sony et al. (2007); Berland et al. (2009); Lee et al. (2011, 2012); Wellendorff et al. (2010); Sun et al. (2010); Li et al. (2012) and between various layered compounds.Björkman et al. (2012b); Londero and Schröder (2010) Yet, its accuracy for sparse matter does not provide the full accuracy of traditional DFT calculations for typical covalently bonded syste, where the generalized-gradient approximation excels. Some of this inaccuracy can be attributed to the overly repulsive revPBEZhang and Yang (1998) exchange functional causing overestimated separations. This exchange functional was chosen because many other standard choices like PBE Perdew et al. (1996) tend to produce spurious binding arising from exchange effects.Wu et al. (2001); Kannemann and Becke (2009); Murray et al. (2009)
Some sparse matter systems are also particularly challenging, like those with weak chemisorption. These systems are highly sensitive to the balance between the different terms in the exchange-correlation functional.Kelkkanen et al. (2011); Hamada and Tsukada (2011); Hamada and Otani (2010); Peköz et al. (2012); Atodiresei et al. (2009); Busse et al. (2011); Brako et al. (2010) The overestimation of binding separations in vdW-DF can make it miss a van der Waals-induced charge transfer, a signature of weak chemisorption.
Refinements of the vdW-DF1 functional have been suggested, in particular it has been suggested to replace its exchange functional with a less repulsive one. Cooper Cooper (2010) designed an exchange functional C09 based on different fundamental constraints than revPBE, following the principles that underpins the introduction of the PBEsol functional.Perdew et al. (2008) C09 is less repulsive than revPBE, but becomes similar in the large -limit, thereby minimizing contributions to unphyscial exchange binding.
Klimeš and coworkersKlimeš et al. (2010) reparameterized a set of exchange functionals by optimizing their performance to the S22 data set.Jurečka et al. (2006) One of these, the optPBE functional, corresponds to a reparameterization of PBE and is therefore somewhat similar to revPBE. This functional will be included in our study.
Recently, a new version of vdW-DF called vdW-DF2 has been developed.Lee et al. (2010a) This functional sets the parameter by using a different description of the plasmon that is suggested by an analysis of the exchange-variation in the high-density limit.Schwinger (1980, 1981); Elliott and Burke (2009) The updated plasmon description makes vdW-DF2 less sensitive to low-density regions and generally reduces the non-local correlation energy. In addition to adjusting the version of correlation, this functional updates the account of exchange by relying on a refitted version of PW86, here labeled PW86r.Perdew and Yue (1986); Murray et al. (2009) This exchange functional agree well with the exchange of Hartre-Fock calculations and avoids spurious exchange binding.Harris (1985); Murray et al. (2009); Wu et al. (2001) This functional is also less repulsive than revPBE at binding separations.
In almost every case, the modifications of vdW-DF description reduces binding separations compared to the results of vdW-DF1. For the S22 data set,Jurečka et al. (2006) vdW-DF1 with C09 exchange and vdW-DF2 produces good results; unsurprisingly so does using vdW-DF1 with optPBE, since it is was fitted to this data set. When comparing the versions for a broader class of systems, we find that no version clearly outperforms the others.
To exemplify this observation, vdW-DF2 stands out by predicting a potential energy curve for H on Cu(111) that compares well with that generated from experimental data,Lee et al. (2011, 2012) both in terms of optimal separation and adsorption energy. However it underestimates the asymptotic form. For C60 crystals,Berland et al. (2011a) vdW-DF1 predicts a cohesive energy 1.6 eV in good agreement with experimental measured interval 1.6-1.9. vdW-DF1 using C09 exchange functional predicts an excellent lattice constant and a cohesive energy of about 2 eV. In contrast vdW-DF2 predicts a cohesive energy of merely 1.3 eV. For adsorption of C60 on Au(111), the combination of vdW-DF2 and C09 has been found to yield the best results.Hamada and Tsukada (2011) For this system, neither vdW-DF1 nor vdW-DF2 are able to predict the right balance between repulsive and attractive terms that results in the experimental finding of a weak chemisorption.Kelkkanen et al. (2011) We find that vdW-DF1 predicts the best binding energies for both benzene and C60 on graphene.
Björkman and coworkers recently asked the question “are we van der Waals ready?” in a large survey of the ability of electronic structure methods, including vdW-DF1 and vdW-DF2, to accurately describe binding of layered systems.Björkman et al. (2012a) None of the DFT functional methods they tested are able to consistently predict accurate binding energies, separations, and elastic constants. For instance both vdW-DF1 and vdW-DF2 in general overestimates interlayer separations.
To get a better understanding of why different functionals prevail for different systems, a more systematical analysis of their binding mechanisms are needed. A suggestion to test performance by comparing with experiments that reveal many interaction properties for the same one system was launched in Refs. Lee et al., 2011, 2012 Here, we seek to map out how both the exchange and correlation functionals contribute to PECs for related molecules of different sizes and based on this analysis propose an additional adsorption-based strategy to refine functional testing. Some of the analysis approaches taken here can be extended to other van der Waals bonded systems, or directly applied to other non-local correlation functionals like those of Vydrov and Voorhis.Vydrov and Voorhis (2009a, b, 2010)
ii.2 Computational method
The vdW-DF adsorption study presented here relies on a post-GGA procedure similar to many previous vdW-DF studies: the charge density is first generated in a self-consistent PBE calculation with ultrasoft pseudopotentials using the Dacapo software package;DAC () next, in a post-processing procedure, we calculate the semi-local exchange at the GGA level and non-local vdW-DF correlation. These terms replace those of the self-consistent calculation. The non-local correlation part is calculated with an in-house real-space code presented and discussed in Ref. Berland et al., 2011a.
This study is similar to previous vdW-DF adsorption studies detailed in Refs. Berland et al., 2009; Berland et al., 2011b; Rohrer and Hyldgaard, 2011 and resolves computational issues like grid sensitivity and handling of supercells in the same manner. Here, the plane-wave cutoff is set to 500 eV and a -sampling of is used. The supercell consists of graphene (boron nitride) primitive unit cells. This size makes inter-adsorbate interaction negligible even for the case of C60 adsorbates (see lower left of Fig. 1), since the non-local interaction between adsorbates is explicitly corrected for.
It is useful to consider the part of the total energy in vdW-DF that contains one-particle kinetic, electrostatic, local correlation, and exchange separately from the remaining non-local correlation . In a purely van der Waals bonded system, the former part typically gives rise to a repulsive contribution to the potential energy curve (PEC), while the non-local correlation part, an attractive. In calculating the contributions to the PEC between a molecule and surface, we here calculate the former by subtracting off the energy off a reference system where the molecule is far from the surface
in practice Å is sufficient. While for the contribution arising from non-local correlation, which have longer range, we subtract off the energy of isolated fragmentsKleis et al. (2008); Berland et al. (2011a)
Iii Adsorption results
We describe the adsorption results of our vdW-DF study of adsorption of benzene and C60 on graphene and boron nitride. We present the binding energies and sheet-to-molecule separation for the different systems and for the different functionals versions. We also calculate the dipole induced by adsorption and show that only C60-on-graphene induces a significant dipole at binding separations.
iii.1 Binding energy and separation
To determine the optimal site and corresponding binding energy and separation, we optimize the separation to the sheet for each high-symmetry site. Fig. 3 shows the six high symmetry sites considered for benzene on graphene. The molecule is assumed to lay flat on the sheet. On boron nitride, we consider eight high symmetry sites, since the center ring can be on-top of a nitrogen (top-N/top-N-rot) and on-top of a boron (top-B/top-B-rot). For C60 we consider both configurations where the carbon hexagon or pentagon face the sheet and find that the hexagon is the preferred one. For benzene, the optimal site is the top on graphene and top-N on boron nitride. The optimal hexagonal configuration are rotated degrees about the surface normal for C60 compared to the optimal one for benzene (top-rot/top-N-rot). The same optimal sites were predicted by the considered versions and variants of vdW-DF.
Figure 4 shows PECs for benzene and C60 on graphene as calculated in vdW-DF1 and vdW-DF2. Table 2 displays the adsorption energies and the optimal molecule-to-sheet separations for the optimal sites including adsorption on BN and use of non-canonical exchange partners.
The benzene adsorption energies on graphene almost match those on BN, while C60 binds slightly stronger to graphene. The optimal benzene-on-graphene separation is slightly larger than the benzene-on-BN, while the C60-on-graphene and C60-on-BN separations are the same. That benzene and C60 does not follow identical trends could stem from a moderate charge transfer between graphene and C60, compared to a much smaller one for the C60 on BN. The difference in charge transfer is discussed in the next subsection.
vdW-DF1 predicts adsorption energies very close to the experimental for both benzene and C60 on graphene. vdW-DF2 somewhat underestimates the experimental value. vdW-DF1 correlation combined with the non-canonical exchange variants considered here all result in overestimated adsorption energies, while vdW-DF2 with C09 exchange underestimates the adsorption energies.
Both vdW-DF1 and vdW-DF2 describe well the expected increased van der Waals attraction for the larger system. vdW-DF1 predicts that C60 binds 73% stronger than benzene on graphene, vdW-DF2 predicts 67%, and the experimental number is 70%. Note though, the specific factor change if a different exchange functional is used.
We are unaware of explicit experimental adsorption energies for benzene and C60 on BN. Reinke and coworkers did find that those C60 molecules not associated with defects desorbed at slightly lower temperature on boron nitride than on graphene. Reinke et al. (2003) This observation is consistent with our results. Binding separations are unavailable for any of the considered systems, complicating a clear assessment of the functionals. vdW-DF1 typically overestimates binding separations by about 0.2-0.3 Å for purely van der Waals bonded systems.Berland and Hyldgaard (2010); Berland et al. (2011a); Langreth et al. (2009) We find that the other variants reduce the separations, in particular those using C09 exchange.
Caciuc and coworkers Caciuc et al. (2012) also calculated benzene adsorption on graphene and BN both for vdW-DF1 and vdW-DF2. Our results agree well with theirs, with the exception their separations are about 0.1 Å shorter for benzene on BN than ours. Using with semi-empirical DFT-D2,Grimme (2006) they found the same optimal sites as us. The fact that their study is based on a projector-augemented waves implementation in VASP,Kresse and Joubert (1999) while ours rely on ultasoft pseudopotential and DACAPO,DAC () strengthens our confidence in the overall precision of both our and their study. The adsorption results for benzene on graphene also closely match the results in the seminal study by Chakarova-Käck and coworkers.Chakarova-Käck et al. (2006)
iii.2 Induced dipole
The upper panel of Fig. 5 visualizes the charge transfer density induced by the adsorption of C60 onto graphene, as given by . In the DFT calculations, grid sensitivity was minimized by using the same the same atomic coordinates relative to the underlying grid for the full systems as for the separate fragments. The induced charge transfer density has a rich dipolar structure. The magnitude of this dipole depends on the separation between the molecule and the sheet.
The lower panel of Fig. 5 shows the induced dipole as function of the benzene-to-sheet (dashed) and C60-to-sheet separation (full) for both graphene (black) and BN (blue) sheets. At typical binding separations and beyond, only C60-on-graphene produces a non-negligible dipole. The dipole that forms between C60 on BN is tiny, yet if the molecule is pushed towards the surface, the magnitude of the dipole increases; while for benzene adsorbants, the magnitude of the dipole barely changes with separation.
To get a crude estimate for the total charge transfer from C60 to graphene, we first project the charge transfer density onto the z-axis . We next identify the position where that lies between the adsorbant and the sheet . Finally, we integrate over the projected charge transfer density smaller than : . The estimated charge transfer, of about , is quite insensitive to binding separation, just like the induced polarization (as shown in the lower panel).
The estimated charge transfer is moderate compared to 0.3 e found for tetrafluoro-tetracyanoquinodimethane on graphene.Pinto et al. (2009) This molecule, basically a hexagonal carbon ring with two arms of (NC) and four fluorine (F) in place of H, has been proposed as a p-dopant on graphene. That this systems results in a significant charge transfer is not surprising considering the electronegativity of fluorine. Since we find some charge transfer for C60 and none for benzene, we speculate that functionalized C60 might be a better candidate than functionalized benzene for doping graphene.
Iv Analysis of exchange and correlation part of vdW-DF binding
We analyze the role of exchange and correlation in the vdW-DF account of adsorption of benzene and C60 on graphene. This analysis benefits from comparing these two molecules, since they have a similar interface to the sheets, but different sizes. This property is first used to illustrate the inherit failure of local approximations to the exchange-correlation functional, as they are unable to capture the increased attraction of C60. Next, we examine the role of different exchange and correlation components to the potential energy curves of benzene and analyze the adsorption energy trends in light of these. We also confirm that the bulk of the difference between the PECs curves of benzene and C60 adsorption arise from non-local correlation. A more detailed analysis of the non-local correlation of vdW-DF follows including the role of density separations and density regimes.
iv.1 Failure of local approximations
The fact that LDA does not include van de Waals forces is well established, yet one can still find recent LDA-based studies for sparse matter. Perhaps this is due to lack of awareness of DFT-D or non-local correlation functional methods. Perhaps a pragmatic point of view is taken: The underlying cause of the binding can be irrelevant for a given investigation, often the case for band structure calculations.Berland et al. (2013)
For some van der Waals bonded systems, an LDA account of exchange-correlation may predict binding energies and separations that compare well with experimental. This occurs because the LDA exchange (as well as correlation) part is inadequate for sparse matter;Murray et al. (2009) once replaced with a more sophisticated GGA functional, most or all this exchange binding vanish. We illustrate how local approximations to exchange-correlation functional is incapable of capturing the larger binding energy of C60 compared to benzene on graphene.
Figure 7 compares the potential energy curves for C60 (black full curves) and benzene (orange full curves). The full thick curves are obtained using vdW-DF1 while the two upper thin curves are with LDA. It produces binding separations that are within the range of what the different version of vdW-DF predict. The binding energy is 0.25 eV, about half of the experimental. Most revealing, however, is the puny 9% increase in LDA binding energy when comparing with C60 adsorption.
As Fig. 7 confirms, LDA does not capture the long-range van der Waals forces. Neither would of course any other local approximation to the exchange-correlation functional. To illustrate how a fitted local approximation might fail, we consider the following toy exchange-correlation functional, in the spirit of Slater (1963),
Here is taken, for the sake of the argument, as a fitting parameter, which we set to to produce the experimental adsorption energy of benzene on graphene, with corresponding PEC given by the dashed orange curve in Fig. 7. The optimal separation of about Å is likely an underestimation, but it is not that different from Å, as predicted by vdW-DF1(C09). However, when using this fitted description for a system of different size, the limitations, as also discussed in Ref. Rydberg et al., 2003, become evident. Applied to C60 on graphene, the same toy functional produces the black dashed curve of Fig. 7. The corresponding adsorption energy is merely 8% larger than that for benzene, completely off the experimental increase of 70%.
Fig. 7 shows density isosurfaces for a benzene-to-graphene separation of 3.0 Å. The well separated inner isosurfaces encapsulate 92 % of the total density density, while the outer, which is connected in some regions at the interface, encapsulate all but 1.6% of the total density. This small density overlap, even at such a short separation, gives a clear indication why it is futile to capture long van der Waals interactions with local density functionals, in which the binding must be described alone by the local density magnitude.
The presented LDA analysis also indicates the problem of combining accounts of van der Waals forces with GGA exchange functionals that induce some exchange binding, in particular if fitting to unrepresentative training sets are used. For bigger systems, in the interface-normal direction, the balance between repulsive and attractive terms can be shifted because a GGA exchange attraction only scales roughly with the size of the interface and it cannot keep up with the increased non-local attraction.
We note that the popular PBE functionalWu et al. (2001); Kannemann and Becke (2009); Murray et al. (2009) has a certain unphysical exchange attraction and that it is often used together with pair-potential accounts of van der Waals forces. The optPBE functional may also induce a small exchange attraction at moderate separations, and is in this sense not optimal. Formally, so does revPBE, but this exchange functional is so repulsive that this is not an issue in practice.
iv.2 Role of exchange and correlation version
Figure 8 presents the exchange and correlation components of the PEC (upper panel) and their derivates (lower panel) for benzene on graphene as obtained with different versions of vdW-DF. The values presented earlier in table 2 are discussed in terms of these curves. The upper blue curve in the upper panel (circular dots) shows the revPBE exchange part of the vdW-DF1 PEC, hereafter referred to as the revPBE curve; and in the same manner for the other exchange curves. The revPBE curve is above the other exchange curves. Hence, in agreement with table 2, canonical vdW-DF1 (with revPBE) gives the smallest binding energy. The revPBE curve is also the least steep among the exchange curves, as shown by its derivative in the lower panel. This result is related to revPBE producing the largest optimal separations.
The C09 curve (green with square markers) is the steepest explaining why vdW-DF1(C09) and vdW-DF2(C09) produce the shortest separations. The C09 curve becomes similar to the revPBE curve at 3.8 Å. This similarity reflects the similar enhancement factors of revPBE and C09 in the large -limit, as mentioned in section II.1.
The optPBE curve has about the same depth as the C09 (upper panel), but is less steep (lower panel). Hence the adsorption energy of vdW-DF1(optPBE) is similar to that of vdW-DF1(C09), while the optimal separation is larger.
The revPBE, optPBE, and C09 curves have a positive value at a separation of 3.8 Å. Since these curves inevitably will decline towards zero as the overlap between the fragments vanish, these exchange variants have the potential to induce spurious binding at large separations; a result related to the saturation of their enhancement factor at large separations.Murray et al. (2009) Whether spurious binding indeed arise for a given binding curve depend on the balance between the different terms contributing to a given PEC.
The PW86r curve stays below at all separations, thus this exchange choice does not induce binding for this system. For small molecular dimers, this non-binding property has been shown to be a consequence of the fact that the PW86r enhancement factor always increases with .Murray et al. (2009) Our result indicates that this non-binding property carries over to to extended systems, which further motivates choosing this exchange partner.
The full curves in the bottom part of the upper panel of Fig. 8 show the non-local correlation contributions to the PEC. The dashed include contributions from LDA correlation.
The vdW-DF1 correlation curve is significantly deeper than the one for vdW-DF2. This result relates to vdW-DF1’s larger adsorption energies. The steepness of these two curves, shown in the lower panel, however differ much less from each other than that the steepness of the exchange curves. Thus, within the selected range of functional variants considered here, the binding separation is primary set by which exchange variant is chosen.
Table 2 shows that for a given exchange variant, vdW-DF1 predicts merely about 0.05 Å shorter binding separations than vdW-DF2 despite the much larger binding energies. On the one hand, vdW-DF1 and vdW-DF2 correlation greatly affects the binding energy, while on the other hand, the correlation version matters less for the binding separation. This observation makes the binding separation a better parameter to assess the quality of an exchange partner to vdW-DF than the energy. In the same vein, if a strategy of fitting exchange to high-accuracy data sets is chosen, we expect the binding separation to be the better target.
The extremely steep LDA exchange curve (dotted curve in upper pane) seems to explain why LDA for some systems produce binding energies comparable to a van der Waals bond but for others severely underestimate it. At larger separations, the curve goes rapidly to zero as the overlap between the fragments vanish. At short, the curve is even deeper than the vdW-DF1 correlation curve. A slight shift in the kinetic-energy repulsion or electrostatic interaction could move the molecule slightly inwards or outwards, adjusting the adsorption energy greatly.
iv.3 Benzene versus C60 adsorption
Figure 10 confirms our assumption that the reason why C60 binds more and at shorter separations than benzene arise from the increased van der Waals interaction. The upper curves gives on graphene. has a mostly-repulsive contribution the PEC of C60 and benzene. The lower gives the same for the non-local correlation part . Since the two lower curves are far deeper than the two upper, the increased binding of C60 over benzene on graphene arise primarily from non-local correlation. Similar results are found for benzene and C60 on BN.
The figure further shows that the non-local correlation of vdW-DF1 deepens the PEC more than that of vdW-DF2 does, while the small shift from is similar for the two functionals. This result is related to the bigger role of non-local correlations in vdW-DF1. That the binding energies of the two functionals do not differ more is related to that revPBE exchange functional is more repulsive than PW86r. As we go from benzene to C60, the repulsive component do not change much. Therefore, the larger van der Waals component of vdW-DF1 compared to vdW-DF2, reduces the binding separation more, and increases the binding energy more.
The result that van der Waals forces is the primary cause of the increased binding of C60 compared to benzene on graphene (and BN) shows that these systems serve as a good starting point to explore properties these forces. Results for these systems can highlight the role played by different length scales or density separations in the non-local correlation.
iv.4 Role of different density separations
To elucidate the role size and geometry have on the van der Waals attraction, we also study how much various density separations contribute to the vdW-DF interaction energy. We perform this analysis for benzene and C60 on graphene and focus on the non-local correlation. For this purpose, we introduce a lower and upper separation cutoff in the vdW-DF kernel:
where is the Heaviside function. Just like for a full calculation, we subtract off, for a given kernel, the non-local correlation energy of the separate fragments (A,B),
This subtraction allows us to obtain a Separation-decomposed measure of the Non-Local correlation contribution to the Interaction energy (SNLI). This approach is similar in nature to the non-local correlation energy density analysis used by Lazić and coworkers.Lazić et al. (2010); Caciuc et al. (2012); Busse et al. (2011) However, while their analysis has emphasis on which spatial regions dominate the response, ours focus on the role of density separations .
The columns of Fig, 10 show the SNLI for separation between and , with . The upper panel is for benzene on graphene and the lower is for C60 (and benzene once more for comparison). The curves are guides to the eye. In all calculations, the adsorbant has an atomic separation of from the graphene sheet, the optimal one for benzene with vdW-DF1. At separations smaller than about , the SNLI is positive, while it is negative beyond. Its oscillatory shape reflects the shape of vdW-DF kernel, displayed in Fig. 2. The SNLI curves in Fig. 10 are dominated by negative values because we have subtracted off the interaction within isolated fragments; short separations only contribute at the interface between the two fragments.
The SNLI of vdW-DF1 and vdW-DF2 are very similar for small density separations in these systems. Beyond the repulsive peak at 0.7 Å, they start to deviate, but first at about 1.7 Å do this deviation become substantial. This result is linked to the stronger asymptote of vdW-DF1 compared to vdW-DF2.Lee et al. (2011); Vydrov and Van Voorhis (2010); Björkman et al. (2012a) The similarity for small reflects the fact that their description of the plasmon dispersion [underlying the non-local correlation of expression (1)] becomes similar for large .Dion et al. (2004); Lee et al. (2010a)
Comparing the benzene (full) and C60 (dashed) curves in the lower panel, we see that they are almost identical for short separations: The benzene-graphene and C60-graphene interface is essentially the same as far as the non-local correlation is concerned. Beyond about 2 Å, the benzene and C60 curves start to differ noticeably, and at 8 Å, the contribution to the non-local correlation interaction is more than four times larger for C60 than for benzene, both in the case of vdW-DF1 and vdW-DF2. Thus, if we gradually consider larger systems (compared to the size of interface), the difference between the non-local correlation of vdW-DF1 and vdW-DF2 increases. This observation is clearly also related to the stronger asymptote of vdW-DF1.
The separation analysis of Fig. 10 also highlights the strikingly different nature of a vdW-DF and a DFT-D account of van der Waals forces. A similar histogram for DFT-D, or strictly pair-potential account for that matter, would begin at the atomic separation, which in this case is 3.6 Å. In contrast, and even despite the extended nature of the systems considered, most of the non-local correlation interaction energy in vdW-DF is already accounted for at separation shorter than 3.6 Å. It is not surprising that vdW-DF is enhanced at moderate separations compared to a corresponding asymptotic pair potential account as discussed in Refs. Berland and Hyldgaard, 2010; Berland et al., 2011a. This enhancement is a reflection of the plasmon-nature of the vdW-DF design.
iv.5 Analysis of density regimes
Our analysis shows that most of the van der Waals interaction in vdW-DF arise from separations shorter than the ionic separation between the molecules. Because low-density regions are closer to each other, this result suggests that the low-density regions contribute the most to the interaction energy arising from the non-local correlation energy of vdW-DF. To quantify the role of low-density regions, we introduce a low-density cutoff in evaluating the non-local correlation of Eq. (1) (rather than separation cutoffs). By gradually increasing this cutoff, we systematically gauge the role of various-density regions.
Figure 11 shows the result of this density analysis for benzene on graphene at two representative separation of Å and Å. The vertical axis indicates the fraction of the non-local correlation interaction energy that remains after a given low-density cutoff has been introduced. The same density cutoff is used in the full calculation as in the reference calculations. The horizontal axis indicates the fraction of the total number of electrons that has been removed in the full system for a given density cutoff. The insert relates the density cutoff [in atomic units] (vertical axis) to the fraction of electrons removed (horizontal axis).
The figure reveals that the vdW-DF interaction energy is very sensitive to low-density regions. For the case of 3.0 Å, removing 2% of the density (as indicated in the outer contour in Fig. 7) removes 15% of the interaction energy; removing 8 % of the density (indicated in the inner counter of Fig. 7) removes 65% of the interaction energy. If we remove 50 % of the density, starting from the smallest density, 97% of the interaction energy is removed. For a vdW-DF like functional that weighted density equally, the non-local correlation would be quadratic in the density and the corresponding numbers would be 4%, 15%, and 75%. These widely different trends reflects how the kernel itself is a functional of the density .
Low density regions are even more important at the larger separation Å. To crudely explain this effect, we can consider that when the molecules move further apart, the density between the molecules is lowered, because the density overlap of the molecules is reduced.At the larger separations, we also find that vdW-DF1 is somewhat more sensitive to low density regions than vdW-DF2, while at 3.0 Å the curves are very similar. This trend agrees well with the separation analysis shown Fig. 10, since vdW-DF1 and vdW-DF2 differ more at larger separations.
Our analysis show that interactions between different low density regions (as connected by the two-point kernel ) dominate the contributions to the interaction energy. Clearly the response of low-density regions is much stronger that that of high density regions. Since part of the low-density regions are closer to each other (the region close to the interface) this result can be linked to the importance of shorter density separations (than atomic separations between fragments) established in the previous subsection.
Finally, we note that the importance of low-density regions has an important bearing on how we build our intuition of how strong a van der Waals interaction will be. According to the presented analysis, it is not the number of atoms that roughly scale with strength of van der Waals interactions but more the amount of low-density regions in a fragment. This observation can be for instance linked to the offset in the binding trends of PAHsChakarova and Schröder (2005); Björk et al. (2010) and alkanesLondero et al. (2012) on graphene and the big response of hydrogen atoms.Berland and Hyldgaard (2010); Björk et al. (2010)
V Corrugation results
The predicted corrugation is very sensitive to which version of vdW-DF and exchange variant is used. Corrugation is here understood as the potential energy variation of a molecule on different sites on the surface. We have recently shownLee et al. (2012) for H on Cu(111) that the corrugation is even more sensitive to switching between a vdW-DF to a DFT-D account. The adsorption energies and optimal separations, presented in section III, is also significantly affected by the vdW-DF variant used; for benzene on graphene, the optimal separations differ by and the energy by . Nevertheless, the vdW-DF sensitivity to corrugation stands out: for this system, the diffusion barriers can change by up to a factor of three as we switch the exchange variant and correlation version.
This huge sensitivity arise because the corrugation depends exponentially on the binding separation; slight changes in separation greatly affect the corrugation.
Figure 12 shows the how the optimal binding energies varies for the different high-symmetry sites on graphene and BN for benzene and C60 adsorbants, as calculated with the considerent variants of vdW-DF. It shows that the corrugation depends strongly on the correlation version used and in particular on the exchange choice. The upper panels show the results for benzene on graphene and on BN. The lower, show those for C60 with the carbon hexagon facing graphene and BN, which is preferred over the pentagon face. Graphene sites are identified in Fig. 3; the BN sites are similar, but with twice as many top sites: the hexagon can be centered on N or on B (top-N and top-B). The vertical scale in the lower panels are twice those of the upper; the corrugation of C60 adsorbants is roughly twice that of benzene.
Benzene can easily diffuse on graphene since the molecule can pass between top and bridge sites which have very similar energies. This open-path dynamics is similar to that of benzene on Cu(111).Berland et al. (2009); Lee et al. (2012) On BN this path is closed, since from a top-N site, the molecule must either move to a top-B site or a hollow (or between) to get to the next top-N site. The corrugation energies are also generally somewhat larger. Thus our calculations indicate that benzene will freeze at much higher temperatures on BN than on graphene. C60 has a similar, though less pronounced, low-energy path on graphene. Rotational barriers are in general small both for benzene and C60. We note that since C60 adsorbants induce a dipole in particular on graphene, in-plane dynamics must be accompanied by continuous charge redistribution. Since electron-hole excitations express such charge-relocations, this will open a damping mechanism causing frictionPersson (1978); Persson and Hellsing (1982); Persson and Persson (1980) for the dynamics for C60 on graphene.
There is a question of whether C60 might diffuse like a soccer ball on graphene. According to classical molecular dynamics calculations, this is not the case for thermally induced motion.Neek-Amal et al. (2010) However, such a classical analysis fail to account for dynamics of the friction that we argue will be present.
Comparing the functionals, we find that vdW-DF1 predicts the smallest corrugations followed by vdW-DF2, while vdW-DF1(C09) clearly predicts the largest. In between are vdW-DF2(C09), vdW-DF1(PW86r), and vdW-DF1(optPBE) which are more or less tied. The succession follows closely that of the binding separations from smallest to largest as listed in table 2. For benzene, vdW-DF1 clearly predicts a larger corrugation than vdW-DF2, but for C60 the difference is smaller. This trend is reflected in the more similar separations predicted for C60 adsorption.
The observation of a strong link between optimal separation and corrugation is somewhat clouded by fact that the separation is optimized separately for each site. To make a clearer analysis, we consider how the difference between the PECs of two specific sites (PEC) changes if we change the exchange variant or correlation version.
Fig. 13 shows, for benzene on graphene, that PEC barely changes with the exchange variant and correlation version used. In contrast, the full PEC might change much. Thus, the corrugation of graphene and BN is predominantly set by the binding separation (given by the minimum of the full PEC), making the huge sensitivity to functional version arise indirectly as consequence of the exponential sensitivity to binding separation.
In Fig. 13, the indirect sensitivity to corrugation is documented separately for the mostly-repulsive part of the PEC (upper panel) and for the non-local correlation part (lower panel). In the upper panel, the difference between of the top site and that of the bridge and hollow, for the exchange choice of revPBE, PW86r, and C09. And in the same manner in the lower panel for for vdW-DF1 and vdW-DF2 correlation. The curves conform closely to each other (the same goes for optPBE exchange which is not shown). The small difference in the curve is much bigger that the difference between the non-local correlation contribution curves . The vertical axis of the lower panel spans one-tenth of that of the upper panel.
The similarity of the tiny contributions to corrugation arising directly from vdW-DF1 and vdW-DF2 correlation can be explained by the similarity of their SNLI curves shown in Fig. 10 for short separations. The curves in the upper panel depend essentially exponentially on separation — for the largest separations noise contributes to the tiny difference between the top and hollow PEC. The curves document how corrugation is set by the optimal separation. It is the effect the exchange-correlation variant has on the binding separation that causes the corrugation to be strongly affected by the variant.
Vi Towards a surface-based experimental reference database
Adsorption on clean surfaces provide outstanding opportunities to compare theory and experiments. Specific molecules can be carefully deposited on surfacesBartels (2010) and imaged with scanning-tunneling (STM) or atomic-force microscopy (AFM). Several experimental probes can be used to gain detailed quantitative insight that can be compared with calculated numbers. One example is the comparison of theoretical and experimental vacuum-level shifts that can be used to deduce binding separations and charge transfer of weakly chemisorbed molecules.Toyoda et al. (2009, 2010, 2011) Another is the measurements of resonance levels arising in back-scattering of light molecules, which provide a powerful connection between experiments and theory.Andersson and Persson (1993a); Andersson et al. (1989); Andersson and Persson (1993b); Andersson et al. (2002); Lee et al. (2011, 2012) These experiments can be used to deduce the shape of the PECs as well as corrugation and could give detailed information on interactions at many different length scales.
vi.1 Proposed concerted theory-experiment approach
Based on our results and analysis, we propose a concerted theory-experiment approach to probe the magnitude of the various exchange-correlation contributions to physisorption on surfaces. The proposed approach is complimentary to others that emphasize a back-scattering based comparisonAndersson and Persson (1993a); Andersson et al. (1989); Andersson and Persson (1993b); Andersson et al. (2002); Lee et al. (2011, 2012) (limited to light molecules) or vacuum-level shift comparisonToyoda et al. (2009, 2010, 2011) (limited to systems with a sizeable charge transfer at the interface).
The suggested approach takes three steps: First, one compares experimental and calculated corrugation and thus determine the binding separation. Second, based on the separation, one estimates how much of the adsorption energy arise directly from the exchange-correlation energy. Third, comparing the adsorption of molecules with the same interface geometry but different shapes and volume, like benzene and C60, one probes how much of the adsorption energy arise from the van der Waals interaction.
The first step is motivated by the fact that the predicted corrugations are strongly sensitive to the vdW-DF variant —in particular to the exchange (Fig. 12)— while at the same time the difference between PECs of two sites is very insensitive to the variant used (Fig. 13). This result indicates that the strong corrugation sensitivity, at least for benzene and C60 adsorption on BN and graphene, arise almost exclusively from the shift in separation induced by the replacement of exchange-correlation functional. Turning this around, we can can conclude that the corrugation essentially determines the separation. If the corrugation can be measured, DFT calculations can determine the binding separation, since the specific exchange-correlation used for this comparison is unimportant.
Many other systems might also exhibit this kind of direct exchange-correlation insensitivity. The identification of additional systems of this kind merits further study, in particular if the corrugation has been accurately measured. The validity of the argument presented here is also limited by the fact that the exchange-correlation account affects the Kohn-Sham (KS) orbitals and thereby kinetic energy contributions. However, for purely van der Waals bonded systems, the direct part of this effect is likely small.
The second step also builds on the assumption that the Kohn-Sham orbitals are not much affected by the specific exchange-correlation functional used. This assumption is appropriate for weakly-bonded systems.Berland et al. (2013) If so, the total exchange-correlation contribution to the adsorption energy can be assessed if both the binding separation and energy is known. The separation also tells us how the exchange-correlation functional affects the steepness of the PEC, since the van der Waals forces must counter the kinetic-energy repulsion to make the molecule physisorb at a given separation.
The third step builds on the result that for molecules with similar interfaces to the surface (or sheet) like benzene and C60, but with different geometries and size, we can single out the effect of van der Waals forces (Fig. 10). Such results give us insight on the delicate balance between the exchange and correlation contributions to the binding. These results can be combined with results for the asymptotic interactions, since advanced methods can compute these to high accuracy.Kumar and Meath (1992); Chu and Dalgarno (2004); Hult et al. (1999); Spackman (1991); Ruzsinszky, Adrienn and Perdew, John P. and Tao, Jianmin and Csonka, Gábor I. and Pitarke, J. M. (2012); Perdew et al. (2012a); Lebègue, S. and Harl, J. and Gould, Tim and Ángyán, J. G. and Kresse, G. and Dobson, J. F. (2010) Such analysis tells us how the van der Waals forces are modulated at shorter separations compared to the asymptotic form.Berland and Hyldgaard (2010); Berland et al. (2011a); Perdew et al. (2012b); Dobson and Gould (2012); Dobson:_asymptotic
How to accurately determine the corrugation experimentally is a question beyond the scope of this paper, but we note that the advancements in STM and AFM is astonishing. One might drag molecules across the surface and measure the potential energy landscape,Ternes et al. (2008); Langewisch et al. (2013) or observe the diffusion of molecules on the surface.Wong et al. (2006); MartÃnez-Galera and GÃ³mez-RodrÃguez (2011) Based on diffusion barriers, corrugation can be estimated.
vi.2 Surface-based database
Many adsorption systems may have the structural relation as benzene and C60 on graphene and BN and that could make them good discriminators of exchange-correlation accounts. Other experimental-systems, like scattering of light molecules on noble metal surfaces or systems with vacuum level shifts also provide excellent links between theory and experiment. As experimental methods continue to develop and new experimental data becomes available, these could be linked to DFT calculations or parameterized in such a manner that it becomes straightforward to assess the performance of sparse matter methods. The construction of PEC based on resonance levels of H backscatteringAndersson and Persson (1993a, b) makes it simple to assess the theoretical methods.Lee et al. (2011, 2012) The relation between corrugation and binding separation discussed in section V provides a different example.
If we identify and organize several of these theory-experiment links for surfaces we can build up a database that could provide invaluable for the continued advancement of sparse-matter methods. Today datasets based on high-accuracy quantum-chemistry calculations play such a role,Björkman et al. (2012a); Jurečka et al. (2006); GraÌfovaÌ et al. (2010); ŘezaÌč et al. (2011) but surface electrons have a different collective nature than that of small molecular dimers.Rydberg (2001); Langreth et al. (2005) Surfaces are also especially suited to extract a range of experimental properties. The corrugation sensitivity to binding separation is central to the proposed approach. Other quantities sensitive to geometry, like the HOMO-LUMO gap of adsorbed molecules,Kubatkin et al. (2003); Neaton et al. (2006) may help uncover other links between theory and experiment useful for assessing sparse-matter methods.
In this paper, we identify a comparison of the adsorption of benzene and C60 on inert surfaces like graphene and BN as a good arena to evaluate exchange-correlation accounts of van der Waals forces and to gain insight into the nature of physisorption. Our study of benzene and C60 on graphene have shown that vdW-DF1 and vdW-DF2 accurately capture the increased van der Waals binding of the larger molecule.
A detailed analysis of adsorption results, exchange-correlation components to the PECs, and corrugation have been performed. Based on this analysis, we have made a number of observations. We have shown that C60 induces a moderate dipole on graphene, while for the other systems dipoles are tiny. We have shown that the binding separation is more sensitive to the exchange variant than the correlation version of vdW-DF. This result making the separation a better target for assessing or optimizing an exchange account. We have shown that the difference between benzene and C60 binding arise predominantly from the increased non-local correlation, displaying fundamental shortcomings of relying on unphysical exchange binding. With our density separation decomposition of the non-local correlation interaction, we documented that vdW-DF1 and vdW-DF2 are similar for short density separations, but differ much at large. The similar density-cutoff analysis reveals how the non-local correlation of vdW-DF is very sensitive to low-density regions.
Further, we have demonstrated that the corrugation is very sensitive to the vdW-DF variant used, but also that a particular low-energy path exists for benzene on graphene, making it easy for this molecule to diffuse. The corrugation sensitivity have been shown to arise almost exclusively as an indirect effect arising from the shift in binding separation caused by the exchange choice and correlation versions. This result has lead us to suggest a experiment-theory approach to assess variants of, and gain insight into, the exchange-correlation functional and its account of van der Waals forces.
The novel analysis approaches presented here, like density-separation and density-regime decomposition, can also be applied to other material systems and functionals. Such approaches can allow us to view materials and methods in new ways helping the development of non-local correlation functionals
Acknowledgements.We thank T. Löfwander and B. I. Lundqvist for discussions. The work was supported by the Swedish Research Council (VR) and the Chalmers Area of Advance. The Swedish National Infrastructure for Computing (SNIC) at the C3SE and HPC2N is acknowledged for computer time.
- Kelkkanen et al. (2009) A. K. Kelkkanen, B. I. Lundqvist, and J. K. Nørskov, J. Chem. Phys. 131, 046102 (2009).
- Møgelhøj et al. (2011) A. Møgelhøj, A. K. Kelkkanen, K. T. Wikfeldt, J. Schiotz, J. J. Mortensen, L. G. M. Pettersson, B. I. Lundqvist, K. W. Jacobsen, A. Nilsson, and J. K. Nørskov, J. Phys. Chem. B 115, 14149 (2011).
- Cooper et al. (2008) V. R. Cooper, T. Thonhauser, and D. C. Langreth, J. Chem. Phys. 128, 204102 (2008).
- Nabok et al. (2008a) D. Nabok, P. Puschnig, and C. Ambrosch-Draxl, Phys. Rev. B. 77, 245316 (2008a).
- Berland and Hyldgaard (2010) K. Berland and P. Hyldgaard, J. Chem. Phys. 132, 134705 (2010).
- Nabok et al. (2008b) D. Nabok, P. Puschnig, and C. Ambrosch-Draxl, Phys. Rev. B 77, 245316 (2008b).
- Berland et al. (2011a) K. Berland, Ø. Borck, and P. Hyldgaard, Comp. Phys. Comm. 182, 1800 (2011a).
- Kong et al. (2009) L. Kong, V. R. Cooper, N. Nijem, K. Li, J. Li, Y. J. Chabal, and D. C. Langreth, Phys. Rev. B 79, 081407 (2009).
- Li and Thonhauser (2012) Q. Li and T. Thonhauser, J. Phys.: Condens. Matter 24, 424204 (2012).
- Langreth and Perdew (1977) D. C. Langreth and J. P. Perdew, Phys. Rev. B 15, 2884 (1977).
- Gunnarsson and Lundqvist (1976) O. Gunnarsson and B. I. Lundqvist, Phys. Rev. B 13, 4274 (1976).
- Perdew and Yue (1986) J. P. Perdew and W. Yue, Phys. Rev. B 33, 8800 (1986).
- Becke (1986) A. D. Becke, J. Chem. Phys. 85, 7184 (1986).
- Perdew and Wang (1992) J. P. Perdew and Y. Wang, Phys. Rev. B 45, 13244 (1992).
- Perdew et al. (1996) J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996).
- Armiento and Mattsson (2005) R. Armiento and A. E. Mattsson, Phys. Rev. B 72, 085108 (2005).
- Wyrick et al. (2011) J. Wyrick, D.-H. Kim, D. Sun, Z. Cheng, W. Lu, Y. Zhu, K. Berland, Y. S. Kim, E. Rotenberg, M. Luo, et al., Nano Lett. 11, 2944 (2011).
- Berland et al. (2009) K. Berland, T. L. Einstein, and P. Hyldgaard, Phys. Rev. B 80, 155431 (2009).
- Chan et al. (2001) H. B. Chan, V. A. Aksyuk, R. N. Kleiman, D. J. Bishop, and F. Capasso, Science 291, 1941 (2001).
- Klimeš and Michaelides (2012) J. Klimeš and A. Michaelides, J. Chem. Phys. 137, 120901 (2012).
- Langreth et al. (2009) D. C. Langreth, B. I. Lundqvist, S. D. Chakarova-Käck, V. R. Cooper, M. Dion, P. Hyldgaard, A. Kelkkanen, J. Kleis, L. Kong, S. Li, et al., J. Phys.:Condens. Matt. 21, 084203+ (2009).
- Wu and Yang (2002) Q. Wu and W. Yang, J. Chem. Phys. 116, 515 (2002).
- Grimme (2004) S. Grimme, J. Comput. Chem. 25, 1463 (2004).
- Grimme (2006) S. Grimme, J. Comput. Chem. 27, 1787 (2006).
- Grimme et al. (2010) S. Grimme, J. Antony, S. Ehrlich, and H. Krieg, J. Chem. Phys. 132, 154104 (2010).
- Tkatchenko and Scheffler (2009) A. Tkatchenko and M. Scheffler, Phys. Rev. Lett. 102, 073005 (2009).
- Rydberg et al. (2000) H. Rydberg, B. I. Lundqvist, D. C. Langreth, and M. Dion, Phys. Rev. B 62, 6997 (2000).
- Rydberg et al. (2003) H. Rydberg, M. Dion, N. Jacobson, E. Schröder, P. Hyldgaard, S. I. Simak, D. C. Langreth, and B. I. Lundqvist, Phys. Rev. Lett. 91, 126402 (2003).
- Dion et al. (2004) M. Dion, H. Rydberg, E. Schröder, D. C. Langreth, and B. I. Lundqvist, Phys. Rev. Lett. 92, 246401 (2004).
- Langreth et al. (2005) D. C. Langreth, M. Dion, H. Rydberg, E. Schröder, P. Hyldgaard, and B. I. Lundqvist, Int. J. Quan. Chem. 101, 599 (2005).
- Thonhauser et al. (2007) T. Thonhauser, V. R. Cooper, S. Li, A. Puzder, P. Hyldgaard, and D. C. Langreth, Phys. Rev. B. 76, 125112 (2007).
- Lee et al. (2010a) K. Lee, E. D. Murray, L. Kong, B. I. Lundqvist, and D. C. Langreth, Phys. Rev. B 82, 081101 (2010a).
- Lundqvist (1967a) B. Lundqvist, Phys. kondens. Materie 6, 193 (1967a).
- Lundqvist (1967b) B. Lundqvist, Phys. kondens. Materie 6, 206 (1967b).
- Lundqvist (1967c) B. Lundqvist, Phys. kondens. Materie 7, 117 (1967c).
- Sernelius and Björk (1998) B. E. Sernelius and P. Björk, Phys. Rev. B 57, 6592 (1998).
- Perdew et al. (2012a) J. P. Perdew, J. Tao, P. Hao, A. Ruzsinszky, G. I. Csonka, and J. M. Pitarke, J. Phys: Condens. Matter 24, 424207 (2012a).
- Johnson and Becke (2006) E. R. Johnson and A. D. Becke, J. Chem. Phys. 124, 174104 (2006).
- Barask (1988) Y. U. Barask, Fiz. Tverd. Tela. 30, 1578 (1988).
- Perdew et al. (2012b) J. P. Perdew, A. Ruzsinszky, J. Sun, S. Glindmeyer, and G. I. Csonka, Phys. Rev. A 86, 062714 (2012b).
- Björkman et al. (2012a) T. Björkman, A. Gulans, A. V. Krasheninnikov, and R. M. Nieminen, J. Phys.: Conden. Mat. 24, 424218 (2012a).
- Jurečka et al. (2006) P. Jurečka, J. Šponer, J. Cerny, and P. Hobza, Phys. Chem. Chem. Phys. 8, (2006).
- GraÌfovaÌ et al. (2010) L. GraÌfovaÌ, M. PitoňaÌk, J. ŘezaÌč, and P. Hobza, J. Chem. Theory Comput. 6, 2365 (2010).
- ŘezaÌč et al. (2011) J. ŘezaÌč, K. E. Riley, and P. Hobza, J. Chem. Theory Comput. 7, 2427 (2011).
- Dobson and Gould (2012) J. F. Dobson and T. Gould, J. Phys.: Condens. Matter 24, 073201 (2012).
- Lee et al. (2011) K. Lee, A. K. Kelkkanen, K. Berland, S. Andersson, D. C. Langreth, E. Schröder, B. I. Lundqvist, and P. Hyldgaard, Phys. Rev. B 84, 193408 (2011).
- Lee et al. (2012) K. Lee, K. Berland, M. Yoon, S. Andersson, E. Schröder, P. Hyldgaard, and B. I. Lundqvist, J. Phys: Condensed Mat. 24, 424213 (2012).
- Cooper (2010) V. R. Cooper, Phys. Rev. B 81, 161104 (2010).
- Murray et al. (2009) É. D. Murray, K. Lee, and D. C. Langreth, J. Chem. Theory Comput. 5, 2754 (2009).
- Klimeš et al. (2010) J. Klimeš, D. R. Bowler, and A. Michaelides, J. Phys.:Condens. Matter 22, 022201 (2010).
- Zhang and Yang (1998) Y. Zhang and W. Yang, Phys. Rev. Lett. 80, 890 (1998).
- Puzder et al. (2006) A. Puzder, M. Dion, and D. C. Langreth, J. Chem. Phys. 124, 164105 (2006).
- Chakarova-Käck et al. (2010) S. D. Chakarova-Käck, A. Vojvodic, J. Kleis, P. Hyldgaard, and E. Schröder, New J. Phys. 12, 013017 (2010).
- Hamada (2010) I. Hamada, J. Chem. Phys. 133, 214503 (2010).
- Viitala et al. (2012) M. Viitala, M. Kuisma, and T. T. Rantala, Phys. Rev. B 85, 085412 (2012).
- Moses et al. (2009) P. G. Moses, J. J. Mortensen, B. I. Lundqvist, and J. K. Nørskov, J. Chem. Phys. 130, 104709 (2009).
- Lee et al. (2010b) K. Lee, Y. Morikawa, and D. C. Langreth, Phys. Rev. B 82, 155461 (2010b).
- Johnston et al. (2008) K. Johnston, J. Kleis, B. I. Lundqvist, and R. M. Nieminen, Phys. Rev. B. 77, 121404 (2008).
- Chakarova-Käck et al. (2006) S. D. Chakarova-Käck, E. Schröder, B. I. Lundqvist, and D. C. Langreth, Phys. Rev. Lett. 96, 146107 (2006).
- Berland et al. (2011b) K. Berland, S. D. Chakarova-Käck, V. R. Cooper, D. C. Langreth, and E. Schröder, J. Phys: Condens. Matter 23, 135001 (2011b).
- Londero et al. (2012) E. Londero, E. K. Karlson, M. Landahl, D. Ostrovskii, J. D. Rydberg, and E. Schröder, J. Phys.: Condens. Matter 24, 424212 (2012).
- Sony et al. (2007) P. Sony, P. Puschnig, D. Nabok, and C. Ambrosch-Draxl, Phys. Rev. Lett. 99, 176401 (2007).
- Wellendorff et al. (2010) J. Wellendorff, A. Kelkkanen, J. Mortensen, B. Lundqvist, and T. Bligaard, Top. Catal. 53, 378 (2010).
- Sun et al. (2010) D. Sun, D.-H. Kim, D. Le, Ø. Borck, K. Berland, K. Kim, W. Lu, Y. Zhu, M. Luo, J. Wyrick, et al., Phys. Rev. B 82, 201410 (2010).
- Li et al. (2012) G. Li, I. Tamblyn, V. R. Cooper, H.-J. Gao, and J. B. Neaton, Phys. Rev. B 85, 121409 (2012).
- Björkman et al. (2012b) T. Björkman, A. Gulans, A. V. Krasheninnikov, and R. M. Nieminen, Phys. Rev. Lett. 108, 235502 (2012b).
- Londero and Schröder (2010) E. Londero and E. Schröder, Phys. Rev. B 82, 054116 (2010).
- Wu et al. (2001) X. Wu, M. C. Vargas, S. Nayak, V. Lotrich, and G. Scoles, J. Chem. Phys. 115, 8748 (2001).
- Kannemann and Becke (2009) F. O. Kannemann and A. D. Becke, J. Chem. Theory Comput. 5, 719 (2009).
- Kelkkanen et al. (2011) A. K. Kelkkanen, B. I. Lundqvist, and J. K. Nørskov, Phys. Rev. B 83, 113401 (2011).
- Hamada and Tsukada (2011) I. Hamada and M. Tsukada, Phys. Rev. B 83, 245437 (2011).
- Hamada and Otani (2010) I. Hamada and M. Otani, Phys. Rev. B 82, 153412 (2010).
- Peköz et al. (2012) R. Peköz, K. Johnston, and D. Donadio, J. Phys. Chem. C 116, 20409 (2012).
- Atodiresei et al. (2009) N. Atodiresei, V. Caciuc, P. Lazić, and S. Blügel, Phys. Rev. Lett. 102, 136809 (2009).
- Busse et al. (2011) C. Busse, P. Lazić, R. Djemour, J. Coraux, T. Gerber, N. Atodiresei, V. Caciuc, R. Brako, A. T. N’Diaye, S. Blügel, et al., Phys. Rev. Lett. 107, 036101 (2011).
- Brako et al. (2010) R. Brako, D. Šokčević, P. Lazić, and N. Atodiresei, New J. Phys 12 (2010).
- Perdew et al. (2008) J. P. Perdew, A. Ruzsinszky, G. I. Csonka, O. A. Vydrov, G. E. Scuseria, L. A. Constantin, X. Zhou, and K. Burke, Phys. Rev. Lett. 100, 136406 (2008).
- Schwinger (1980) J. Schwinger, Phys. Rev. A 22, 1827 (1980).
- Schwinger (1981) J. Schwinger, Phys. Rev. A 24, 2353 (1981).
- Elliott and Burke (2009) P. Elliott and K. Burke, Can. J. Chem. 87, 1485 (2009).
- Harris (1985) J. Harris, Phys. Rev. B 31, 1770 (1985).
- Vydrov and Voorhis (2009a) O. A. Vydrov and T. V. Voorhis, J. Chem. Phys. 130, 104105 (2009a).
- Vydrov and Voorhis (2009b) O. A. Vydrov and T. V. Voorhis, Phys. Rev. Lett. 103, 063004 (2009b).
- Vydrov and Voorhis (2010) O. A. Vydrov and T. V. Voorhis, J. Chem. Phys. 133, 244103 (2010).
- (85) Open source code DACAPO, http://www.fysik.dtu.dk/CAMPOS/.
- Rohrer and Hyldgaard (2011) J. Rohrer and P. Hyldgaard, Phys. Rev. B 83, 165423 (2011).
- Kleis et al. (2008) J. Kleis, E. Schröder, and P. Hyldgaard, Phys. Rev. B. 77, 205422 (2008).
- Zacharia et al. (2004) R. Zacharia, H. Ulbricht, and T. Hertel, Phys. Rev. B 69, 155406 (2004).
- Ulbricht et al. (2003) H. Ulbricht, G. Moos, and T. Hertel, Phys. Rev. Lett. 90, 095501 (2003).
- Reinke et al. (2003) P. Reinke, H. Feldermann, and P. Oelhafen, J. Chem. Phys. 119, 12547 (2003).
- Caciuc et al. (2012) V. Caciuc, N. Atodiresei, M. Callsen, P. Lazić, and S. Blügel, J. Phys.: Conden. Mat. 24, 424214 (2012).
- Kresse and Joubert (1999) G. Kresse and D. Joubert, Phys. Rev. B 59, 1758 (1999).
- Pinto et al. (2009) H. Pinto, R. Jones, J. P. Goss, and P. R. Briddon, J. Phys.:Condens. Matter 21, 402001 (2009).
- Humphrey et al. (1996) W. Humphrey, A. Dalke, and K. Schulten, J. Mol. Graphics 14, 33 (1996).
- Berland et al. (2013) K. Berland, E. Londero, E. Schröder, and P. Hyldgaard (2013), unpublished.
Quantum Theory of Molecules and Solids:
Electronic structure of molecules, International series in pure and applied physics (McGraw-Hill, 1963).
- Lazić et al. (2010) P. Lazić, M. Alaei, N. Atodiresei, V. Caciuc, R. Brako, and S. Blügel, Phys. Rev. B 81, 045401 (2010).
- Vydrov and Van Voorhis (2010) O. A. Vydrov and T. Van Voorhis, Phys. Rev. A 81, 062708 (2010).
- Chakarova and Schröder (2005) S. D. Chakarova and E. Schröder, J. Chem. Phys. 122, 054102 (2005).
- Björk et al. (2010) J. Björk, F. Hanke, C.-A. Palma, P. Samori, M. Cecchini, and M. Persson, J Phys. Chem. Lett. 1, 3407 (2010).
- Persson (1978) B. N. J. Persson, J. Phys. C: Sol. Stat. Phys. 11, 4251 (1978).
- Persson and Hellsing (1982) M. Persson and B. Hellsing, Phys. Rev. Lett. 49, 662 (1982).
- Persson and Persson (1980) B. Persson and M. Persson, Surf. Sci. 97, 609 (1980).
- Neek-Amal et al. (2010) M. Neek-Amal, N. Abedpour, S. N. Rasuli, A. Naji, and M. R. Ejtehadi, Phys. Rev. E 82, 051605 (2010).
- Bartels (2010) L. Bartels, Nature Chem. 2, 87 (2010).
- Toyoda et al. (2009) K. Toyoda, Y. Nakano, I. Hamada, K. Lee, S. Yanagisawa, and Y. Morikawa, Surf. Sci. 603, 2912 (2009).
- Toyoda et al. (2010) K. Toyoda, I. Hamada, K. Lee, S. Yanagisawa, and Y. Morikawa, J. Chem. Phys. 132, 134703 (2010).
- Toyoda et al. (2011) K. Toyoda, I. Hamada, K. Lee, S. Yanagisawa, and Y. Morikawa, J. Phys. Chem. C 115, 5767 (2011).
- Andersson and Persson (1993a) S. Andersson and M. Persson, Phys. Rev. Lett. 70, 202 (1993a).
- Andersson et al. (1989) S. Andersson, L. Wilzén, M. Persson, and J. Harris, Phys. Rev. B 40, 8146 (1989).
- Andersson and Persson (1993b) S. Andersson and M. Persson, Phys. Rev. B 48, 5685 (1993b).
- Andersson et al. (2002) T. Andersson, F. Althoff, P. Linde, S. Andersson, and K. Burke, Phys. Rev. B 65, 045409 (2002).
- Kumar and Meath (1992) A. Kumar and W. J. Meath, Molecular Physics 75, 311 (1992).
- Chu and Dalgarno (2004) X. Chu and A. Dalgarno, J. Chem. Phys. 121, 4083 (2004).
- Hult et al. (1999) E. Hult, H. Rydberg, B. I. Lundqvist, and D. C. Langreth, Phys. Rev. B 59, 4708 (1999).
- Spackman (1991) M. A. Spackman, J. Chem. Phys. 94, 1288 (1991).
- Ruzsinszky, Adrienn and Perdew, John P. and Tao, Jianmin and Csonka, Gábor I. and Pitarke, J. M. (2012) Ruzsinszky, Adrienn and Perdew, John P. and Tao, Jianmin and Csonka, Gábor I. and Pitarke, J. M., Phys. Rev. Lett. 109, 233203 (2012).
- Lebègue, S. and Harl, J. and Gould, Tim and Ángyán, J. G. and Kresse, G. and Dobson, J. F. (2010) Lebègue, S. and Harl, J. and Gould, Tim and Ángyán, J. G. and Kresse, G. and Dobson, J. F., Phys. Rev. Lett. 105, 196401 (2010).
- Ternes et al. (2008) M. Ternes, C. P. Lutz, C. F. Hirjibehedin, F. J. Giessibl, and A. J. Heinrich, Science 319, 1066 (2008).
- Langewisch et al. (2013) G. Langewisch, J. Falter, H. Fuchs, and A. Schirmeisen, Phys. Rev. Lett. 110, 036101 (2013).
- Wong et al. (2006) K. L. Wong, K.-Y. Kwon, and L. Bartels, Appl. Phys. Lett. 88, 183106 (2006).
- MartÃnez-Galera and GÃ³mez-RodrÃguez (2011) A. J. MartÃnez-Galera and J. M. GÃ³mez-RodrÃguez, J. Phys. Chem. C 115, 23036 (2011).
- Rydberg (2001) H. Rydberg, Ph.D. thesis, Chalmers University of Technology (2001).
- Kubatkin et al. (2003) S. Kubatkin, A. Danilov, M. Hjort, J. Cornil, J.-L. Brédas, N. Stuhr-Hansen, P. Hedegård, and T. Bjørnholm, Nature 425, 698 (2003).
- Neaton et al. (2006) J. B. Neaton, M. S. Hybertsen, and S. G. Louie, Phys. Rev. Lett. 97, 216405 (2006).