Application of van der Waals Density Functional to an Extended System: Adsorption of Benzene and Naphthalene on Graphite

0.15 MB·English

by  Svetla D. Chakarova Kack

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).