Application of van der Waals Density Functional to an Extended System: Adsorption of Benzene and Naphthalene on Graphite Svetla D. Chakarova K¨ack, Elsebeth Schr¨oder, Bengt I. Lundqvist Department of Applied Physics, Chalmers University of Technology, SE-412 96 G¨oteborg, Sweden David C. Langreth 6 Department of Physics and Astronomy, Rutgers University, Piscataway, N. J. 08854-8019, USA 0 (Dated: February 4, 2008) 0 2 It is shown that it is now possible to include van der Waals interactions via a nonempirical n implementationofdensityfunctionaltheorytodescribethecorrelationenergyinelectronicstructure a calculationsoninfinitesystemsofnoparticularsymmetry. ThevdW-DFfunctional[Phys.Rev.Lett. J 92, 246401 (2004)] is applied to the adsorption of benzene and naphthalene on an infinite sheet of 9 graphite, as well as the binding between two graphite sheets. Comparison with recent thermal desorption data [Phys.Rev. B 69, 535406 (2004)] shows great promise for thevdW-DFmethod. ] i c s Arecentstudyoftheinteractionofpolycyclicaromatic Condensed matter is held together by several kinds of - hydrocarbon molecules (PAH’s) with the basal plane of interatomic forces, including the ubiquitous vdW forces, l r graphite [1] provides experimental benchmark data that which are particularly significant in sparse matter. For t m constitute an ideal challenge for our recently proposed densematterDFThaswellprovenitsvalue,state-of-the- densityfunctional(vdW-DF)[2]whichbothincludesvan artversionsofitgivingvaluesforground-stateproperties . t a der Waals (vdW) interactions and all the benefits of the of covalent molecules and hard materials close to exper- m earlier state-of-the-art versions of the density-functional imental data. The key to success in DFT is the func- theory (DFT). Aiming at a better experimental charac- tional for exchange and correlation, and DFT calcula- - d terizationof the weak interlayerinteractionsingraphite, tionstodaytypicallyapplysomeflavorofthegeneralized- n carefulanalysisofthermal-desorptionkinetics yieldacti- gradient approximation (GGA) [10, 11]. For sparse o vation energies for benzene and PAH’s at submonolayer matter, including soft matter, vdW complexes, poly- c coverageswithexpliciterrorbars[1]. Ourcalculatedval- mers, biomolecules, and other abundant systems, how- [ ues for the adsorption energy of benzene and naphtha- ever,DFTinGGAperformsbadly. Forinstance,itgives 1 lene on graphene and for the weak interlayer interaction unphysicalresultsfortheinterlayerbondofgraphite[12], v energy of graphene agree with the values deduced from a canonical vdW case. The vdW interaction stems from 1 experiment. From this we conclude that the vdW-DF truly nonlocal correlation and is contained in neither 6 1 is, indeed, very promising, and that it can be applied to thelocal-densityapproximation(LDA)northesemilocal 1 systems that are neither periodic nor finite. This distin- GGA. We have subsequently proposed a density func- 0 guishes it from the various wave-function methods that tional(vdW-DF) for generalgeometries[2] thatincludes 6 are often applied to vdW complexes. both vdW interactions and all GGA achievements. It 0 has been applied to the benzene dimer [2, 13] and to t/ Our method differs also from a newly published study complexes of benzene with phenol, toluene, fluoroben- a of adenine on graphite [3], which treats the vdW part of zene,orbenzonitrile[14]. Allthesepreviousapplications m thecorrelationenergybyafrequentlyused[4,5,6,7,8,9] have been to finite systems. We now apply it to single - semi-empirical method. This method introduces empir- benzene and naphthalene molecules in interaction with d n ical damping functions applied to an asymptotic attrac- an extended graphene sheet. This establishes the appli- o tive 1/R6 interaction assumed to occur between each cabilityofthevdW-DFfunctionaltoamoregeneralclass c pair of nuclei. At shorter distances this interaction is of systems. : damped by a physically motivated, but arbitrary and v varying functional form, which introduces one empirical The derivation of the vdW-DF starts from the adia- i X parameterfor everypair of atomic types in the complex. batic connection formula and benefits from insights into r On the other hand our method (vdW-DF) for the cor- the polarization properties of the inhomogeneous elec- a relation energy is completely free from empiricism, and tron gas. It involves several approximations, which sat- although containing approximations, represents a first- isfy known limits, sum rules, and invariances [2]. To principlesdensityfunctional,whichsincetheappearance assess the validity of vdW-DF, internal criteria must ofRef.2isapplicabletoarbitrarygeometries,andwhich be supplemented with external ones, that is, we must is seamless as two fragments merge into a single one. As evaluate its performance in actual applications by com- discussed later,it has been applied to a number of phys- parison with experimental findings. Systems where the ical systems with promising results. The present appli- vdW forces play a prominent role or dominate the in- cation is particularly pertinent, however, as alternative teractions entirely are to be preferred. Being based on first-principles methods for including vdW interactions inhomogeneous-electron-gas polarization properties, the are lacking for extended systems. vdW-DF should work best on systems with similar po- 2 larizationproperties,thatis,ondelocalizedsystemswith dense excitation spectra. Dimers of benzene and PAH’s aresuchsystems. Experimentaldataforthesearescarce, and in practice the assessment has to be made by com- paring to accurate wave-function calculations [2, 13], for which, however, the computational effort grows strongly with size, thus making benchmark results unavailable for really large systems. Further, to properly assess the vdW-DFmethod,extendedsystemsshouldbeaddressed. The adsorption problem, with its extended substrate, is far beyond the possibilities of such wave-function meth- ods. The agreement between values calculated with FIG. 1: The lateral configuration of the AB adsorption vdW-DFandmeasuredvaluesofkeyquantitiesisdemon- positions for the molecules benzene and naphthalene on a graphene sheet. Carbon atoms are shown in black, when strated below. belonging to the substrate and gray, when belonging to the Recentthermal-desorptionstudiesoftheinteractionof adsorbate. Hydrogen atoms are shown in white. PAH’swiththe basalplaneofgraphite[1]yield,through the analysis of desorption kinetics, activation energies of 0.50, 0.85, 1.40, and 2.1 eV for benzene, naphtha- SC calculations are performed with the pseudo- lene, coronene, and ovalene, respectively, at submono- potential-plane-wave-based DFT code dacapo [18], us- layer coverages. Using a force-field model [15], the au- ingthe revPBEfunctionaltogetGGAelectrondensities thorsofRef.1predictanexfoliation energy of52±5meV and total energies. The hexagonalcell used for the stan- per/atom to peel a single graphene layer off a graphite dard part of the DFT calculations has the size (12.32 ˚A, surface. Another determination [16] involving collapsed 12.32 ˚A, 26 ˚A) for the benzene calculations and (14.783 nanotubes,whichobtained35+15meV/atomforarelated −10 ˚A, 12.32 ˚A, 26 ˚A) for the naphthalene calculations. We quantity(theenergytoseparatetwographenesheets),re- use a Monkhorst-Pack grid with (2, 2, 1) k-points and a quired an equally tenuous theoretical model to extract a plane-wave cutoff energy of 450 eV. The in-plane struc- prediction. HereweapplyvdW-DFtocalculatethebind- ture of the adsorbate is found through a monomer cal- ing energyofa benzene anda naphthalene molecule to a culation, where the molecule is fully relaxed within a singlegraphenesheet,aswellasthebindingenergyoftwo GGA-revPBE SC calculation. This results in an opti- graphene sheets to each other. We also present evidence mumstructurewithcarbon-carbonandcarbon-hydrogen thatsecond-layerinteractionsaresmall,whichimplythat bondlengths in agreement with experimental data [19]. our results also apply approximately to molecular des- In the same way, the structure of the graphene sheet is orption or to graphene exfoliation from a basal graphite obtainedfromanisolatedcalculation,resultinginanop- surface, as relevant for the experiment mentioned first timumintralayerlatticeconstantof2.46˚A,inagreement above. withexperiment[20]. Afterwardsthe in-planestructures The generalgeometry(gg)vdW-DF [2]basically“cor- arekeptfixedwhilewevarythedistancefromthe adsor- rects” the correlation part of the energy of a standard bate to the surface and map out an energy profile as a self-consistent(SC) DFT calculation,using the standard function of the separation. We place the molecule in an exchange from the GGA of the revPBE-flavor [17], cho- ABconfiguration(Fig.1),whichweexpecttobeenerget- sen for its being closest to the Hartree-Fock exchange icallymostfavorable,asisthecaseforbothbenzeneand in some key applications [2, 13]. The vdW-DF energy naphthalenedimers,andforgraphenelayers. Thatis,we functional is written place the center of the benzene ring exactly above a car- E =E −E +E +Enl, (1) bonatominthegraphenesheet,andfornaphthalene,the vdW-DF GGA GGA,c LDA,c c moleculecenter-of-massispositionedabovethecenter-of- mass of the two graphite atoms below the naphthalene where the LDA correlation energy E is substituted LDA,c rings, as seen in Fig. 1 [21]. We know that for naphtha- for the GGA correlation energy, E , and where the GGA,c nonlocal correlation energy Enl is added. The last three lene dimers in the AB stacking, small shifts in the exact c lateral position yield minor changes in the total energy, terms in (1) are treated as a post-GGA perturbation, while the separation in the perpendicular direction is of utilizing their low sensitivity to the choice of GGA SC much greater importance [22], and we believe the same electronicdensity. The nonlocalcorrelationenergyis ex- should hold also for small lateral shifts with respect to pressed as the graphene layer in Fig. 1. 1 Fortheadsorptionpositiongivenabove,fullGGAcal- Enl = d3rd3r′n(~r)φ(~r,~r′)n(~r′), (2) c 2Z culations are performed, thus obtaining the first three terms in Eq. (1). The last term, given by Eq. (2), re- ′ ′ where φ(~r,~r ) is a function depending on |~r −~r |, the quires some extra considerations, as (lateral) cell size charge density n and its gradient at the positions ~r and becomes more important when evaluating the nonlocal ′ ~r , respectively, given in detail by Eq. (14) of Ref. 2. vdW energy correction than for the standard calcula- 3 tions. Forthispurposeweenlargethesystemsizeforthe Enlevaluationasfollows. TheSCelectrondensityforthe 20 c adsorbate-substratesystemobtainedwith GGA-revPBE 10 Exp. (Ref. 1) isusedwithintheunitcell,whileoutsidethisunitcellthe Exp. (Ref. 16) graphene electronic density is simulated by that of pure m] 0 graphene, enabling the extension of the substrate to, in o at−10 principle, arbitrary size [23]. Thus we can obtain the V/ e binding energy for benzene and naphthalene interacting m−20 with an increasingly larger circle of the graphene layer y [ g−30 belowthecenteroftheadsorbedmolecule. While91%of er n E the interactionis obtained for the originalunit cell, 99% −40 isobtainedfora10˚Aradiuscircle,andfullconvergenceis −50 reachedfora14˚Acircle. Thesamemethodisusedforthe calculationontheinteractionsofthetwographenelayers. −60 3 4 5 6 7 8 Thus calculationof Ecnl allows anincrease of system size separation [Å] withoutmuchincreasedcomputationalcostcomparedto standard DFT calculations. FIG.2: (Coloronline)Binding-energycurvefortwographene Binding-energy results for two interacting graphene sheets in AB stacking. The calculated general geometry sheets in the AB stacking, which is the structure occur- vdW-DFcurve(solid line) indicates a binding energy of 45.5 ringinbulkgraphite,areshowninFig.2. Thegraphene– meV/atomattheequilibriumseparation3.6˚A,whichisclose to the experimental energy estimates (diamonds with energy graphene binding-energy curve is somewhat deeper than error bars; see text). The two vertical dashed lines mark the that obtained with our earlier layered-functional version positionsofthepotential-energyminimum(3.6˚A)andtheap- [12]. Thegeneralgeometryversionemploysimprovement proximate point where a next-neighbor graphite layer would made possible by an approximate expansion not made havebeen placed (7.2 ˚A). in the layered-functional version. Further discussion of the similarities anddifferences between the two methods is given separately [24]. The gg-vdW-DF curve deviates tions as shown in Fig. 1. We judge these calculations onlyslightlyfromthecorrespondingoneinRef.25,where equivalent to those of adsorption on the basal plane of an averaging procedure was used. graphite, when a small correction for a second graphite OurresultsfortwographenesheetsintheABstacking layer (3.0% increase for benzene and 3.3% for naph- may be compared with the results of the experimental thalene) is taken into account, as discussed above for studiesdiscussedabove[1,16]. ThesearemarkedinFig. graphene-graphene. Figure 3 shows the calculated gg- 2withdiamondsatthelayerseparationinbulkgraphite, vdW-DF potential-energy curves of adsorbed benzene whichis3.36˚A[20]. WeexpecttheseparationoftheAB and naphthalene, respectively, with their minima well graphenesheetstobesimilar,oronlyslightlylarger. The within the ranges of measured binding-energy values energy error bars from the two experiments merge into [1, 28]. For benzene, we find a binding energy of 495 one. meV (to be compared to the values 500± 80 meV [1] While the exfoliation energy from a graphite basal and 480 meV [28]) at the equilibrium separation 3.6 ˚A. plane is not the same as the energy to separate two The same separation is found optimal for naphthalene, graphene sheets, we may use the value of our curve at a withtheadsorptionenergy763meV,which(particularly, second layer distance to estimate the change in our pre- when including the second-layer 3.3% contribution) can diction if an exfoliation calculation were done instead. be compared to the experimental values 800±100 meV From the value of our curve of −2.5 meV at 7.2 ˚A in and900±70meVresultingfromslightlydifferentanaly- Fig. 2, we can estimate that our method will give a sisofthesameexperiment[1],summedupbytheauthors 48 meV/atom exfoliation energy, and a 50.5 meV/atom themselvesbyreportingthenumber850meV.Theexper- graphitecleavageenergy,a5.5%increaseeach [26]. This imentalvalueswitherrorbarsontheenergyareshownas conclusion requires the forces from subsequent layers to shaded regions in the figures. The separations were not be negligible, which is certainly true for our functional, measured in the experiment. However, it is reasonable andprobablyalsowhenthe anomalousasymptoticeffect toexpecttheseparationtobe similartothatingraphite of graphite’s semimetallic nature [27] is accounted for. also in these two cases, which is what we find. Forgraphite,thermalvibrationsalsogivearelevantcon- The analysis of the desorption experiments on ben- tributiontotheenergyfromthemotionperpendicularto zene and naphthalene and the close agreement between the sheets [7]. In any case the comparisonof our curves experimentaland theoretical results let us conclude that with the points given by the two experimental groups is the vdW-DF is a very promising functional to account satisfying,despite thedependence ofthe latteronearlier for the nonlocal vdW forces in these very typical cases. theoretical methods that cannot be fully justified. While the experimental evidence is not so direct in the Binding-energy curves are also calculated for benzene graphene–graphene binding case, the values found are and naphthalene adsorbed on graphene, with configura- probably representative of what more direct experimen- 4 0.4 0.4 0.2 0.2 0 0 V] V] e−0.2 e−0.2 y [ y [ g g er−0.4 er−0.4 n n E E −0.6 −0.6 −0.8 −0.8 −1 −1 3 4 5 6 7 8 3 4 5 6 7 8 separation [Å] separation [Å] FIG.3: (Coloronline)Benzene-graphene(leftpanel)andnaphthalene-graphene(rightpanel)binding-energycurves. Likethat forgraphene-graphene,thegeneralgeometry vdW-DFcalculations (solid lines) giveresultscomparabletoexperiment(boxes); see text. The two vertical dashed lines mark theminima and thenext-neighborgraphite layer, as in Figure 2. tal methods would find, and give additional support to ATOMICS consortium and the Swedish Research Coun- the vdW-DF method. cil, as well as allocation of computer time at UNICC Valuable exchanges with P. Hyldgaard and T. Her- (Chalmers) and SNIC (Swedish National Infrastructure tel are gratefully acknowledged, as is support from for Computing). Work by D.C.L. was supported in part the Swedish Foundation for Strategic Research via the by NSF Grant DMR-0456937. [1] R. Zacharia, H. Ulbricht, and T. Hertel, Phys. Rev. B [18] Computer code dacapo, http://www.fysik.dtu.dk 69, 155406 (2004). /CAMPOS/ [2] M. Dion, H. Rydberg,E. Schr¨oder, D. C. Langreth, and [19] S. Chakarova, Licenciate Thesis, Chalmers University of B.I.Lundqvist,Phys.Rev.Lett. 92,246401 (2004); 95, Technology, G¨oteborg, Sweden (2004). 109902 (2005). [20] Landolt-B¨ornstein search (Springer, Berlin, 2003), [3] F. Ortmann, W.G. Schmidt, and F. Bechstedt, Phys. http://link.springer.de. Rev.Lett. 95, 186101 (2005). [21] The C–C bondlengths in naphthalene differ from those [4] X.Wu et al., J. Chem. Phys. 115, 8748 (2001). in graphite by up to 3%. [5] M. Elstner et al., J. Chem. Phys. 114, 5149 (2001). [22] S. Tsuzuki, K. Honda, T. Uchimaru, and M. Mikami, J. [6] Q.Wu and W. Yang, J. Chem. Phys. 116, 515 (2002). Chem. Phys.120, 647 (2004). [7] M.HasegawaandK.Nishidate,Phys.Rev.B70,205431 [23] For Ecnl we sampled the electron density only on ev- (2004). erysecondstandard-DFTFast-Fourier-Transform(FFT) [8] U. Zimmerli, M. Parrinello, and P. Koumoutsakos, J. gridpoint in each direction. (Ourstandard-DFTcalcula- Chem. Phys.120, 2693 (2004). tions have approx. 0.14 ˚A/FFT-gridpoint). A finer sam- [9] S.Grimme, J. Comp. Chem. 25, 1463 (2004). pling scheme using every FFT gridpoint was tested and [10] V. N. Staroverov, G. E. Scuseria, J. Tao, and J. P. found unnecessary. Perdew, J. Chem. Phys. 119, 12129 (2003). [24] S.D.ChakarovaK¨ack,J.Kleis,andE.Schr¨oder,Dimers [11] V. N. Staroverov, G. E. Scuseria, J. Tao, and J. P. of polycyclic aromatic hydrocarbons indensity functional Perdew, Phys. Rev.B 69, 075102 (2004). theory, to be published. [12] H. Rydberg, M. Dion, N. Jacobson, E. Schr¨oder, P. [25] M. Dion, Ph.D. Thesis, Rutgers University, New Hyldgaard, S. I. Simak, D. C. Langreth, and B. I. Brunswick, New Jersey, (2004). Lundqvist,Phys. Rev.Lett. 91, 126402 (2003). [26] Earlier theory [L. A. Girifalco and R. A. Lad, J. Chem. [13] A. Puzder, M. Dion, and D. C. Langreth, Phys. 25, 693 (1956)] implied a larger 18% difference, cond-mat/0509421 (submitted to J. Chem. Phys.). later used in [1] to predict a cleavage energy of 61 [14] T. Thonhauser, A. Puzder, and D. C. Langreth, meV/atom. cond-mat/0509426 (submitted to J. Chem. Phys.). [27] J. F. Dobson, A. White, and A. Rubio, [15] N. L. Allinger, Y. H. Yuh, and J.-H. Lii, J. Am. Chem. cond-mat/0507156. Soc. 111, 8551 (1989). [28] C. Pierce, Journ. Phys. Chem. 73, 813 (1969). [16] L.X.Benedictetal.,Chem.Phys.Lett.286,490(1998). [17] Y.ZhangandW.Yang,Phys.Rev.Lett.80,890(1998).