Potassium intercalation in graphite: A van der Waals density-functional
  study
Potassium intercalation in graphite: A van der Waals density-functional study
Eleni Ziambaras,1 Jesper Kleis,1 Elsebeth Schröder,1 and Per Hyldgaard1, 2, ∗
Department of Applied Physics, Chalmers University of Technology, SE–412 96 Göteborg, Sweden
Microtechnology and Nanoscience, MC2, Chalmers University of Technology, SE–412 96 Göteborg, Sweden
(Dated: April 1, 2007)
Potassium intercalation in graphite is investigated by first-principles theory. The bonding in
the potassium-graphite compound is reasonably well accounted for by traditional semilocal density
functional theory (DFT) calculations. However, to investigate the intercalate formation energy
from pure potassium atoms and graphite requires use of a description of the graphite interlayer
binding and thus a consistent account of the nonlocal dispersive interactions. This is included
seamlessly with ordinary DFT by a van der Waals density functional (vdW-DF) approach [Phys.
Rev. Lett. 92, 246401 (2004)]. The use of the vdW-DF is found to stabilize the graphite crystal,
with crystal parameters in fair agreement with experiments. For graphite and potassium-intercalated
graphite structural parameters such as binding separation, layer binding energy, formation energy,
and bulk modulus are reported. Also the adsorption and sub-surface potassium absorption energies
are reported. The vdW-DF description, compared with the traditional semilocal approach, is found
to weakly soften the elastic response.
I. INTRODUCTION
Graphite with its layered structure is easily interca-
lated by alkali metals (AM) already at room tempera-
ture. The intercalated compound has two-dimensional
layers of AM between graphite layers,1,2,3,4,5 giving rise
to interesting properties, such as superconductivity.6,7
The formation of an AM-graphite intercalate proceeds
with adsorption of AM atoms on graphite and absorption
of AM atoms below the top graphite layer, after which
further exposure to AM atoms leads the AM intercalate
compound.
Recent experiments8,9 on the structure and elec-
tronic properties of AM/graphite systems use samples
of graphite that are prepared by heating SiC crystals
to temperatures around ∼ 1400◦ C.10 This heat-induced
graphitization is of great value for spectroscopic studies
of graphitic systems, since the resulting graphite overlay-
ers are of excellent quality.11 The nature of the bonding
between the SiC surfaces and graphite has been explored
experimentally with photoemission spectroscopy12 and
theoretically13 with a van der Waals density functional
(vdW-DF) theory approach that accounts for the van der
Waals (vdW) forces.14,15,16,17
Here we investigate with density functional theory
(DFT) the effects on the graphite structure and the
energetics and the elastic response when potassium
is intercalated. The final intercalate compound is
C8K. The AM intercalate system is interesting in it-
self and has been the focus of numerous experimen-
tal investigations.18,19,20,21,22 Graphitic systems are also
ideal test materials in ongoing theory development that
aims at improving the description of the nonlocal inter-
layer bonds in sparse systems.14,23,24 Standard DFT ap-
proaches are based on local (local density approximation,
LDA) and semilocal approximations (generalized gradi-
ent approximation, GGA)25,26,27,28 for the electron ex-
change and correlation. Such regular DFT tools do not
treat correctly the weak vdW binding, e.g., the cohe-
sion between (adjacent) graphite layers. The failure of
traditional DFT for graphite makes it impossible to ob-
tain a meaningful comparison of the energetics in on-
surface AM adsorption and subsurface AM absorption.
Conversely, investigations of graphitic systems like C8K
permit us to test the accuracy of our vdW-DF develop-
ment work.
We explore the nature of the bonding of graphite,
the process leading to intercalation via adsorption and
absorption of potassium, and the nature of potassium-
intercalated graphite C8K using a recently developed
vdW-DF density functional.16 This choice of functional is
essential for a comparison of graphite and C8K properties
because of the inability of traditional GGA-based DFT
to describe graphite. We calculate the structure and elas-
tic response (bulk modulus B0) of pristine graphite and
potassium intercalated graphite and we present results
for the formation energies of the C8K system.
The intercalation of potassium in graphite is preceded
by the adsorption of potassium on top of a graphite
surface and potassium absorption underneath the top
graphite layer of the surface. In this work we study
how potassium bonds to graphite in these two parts of
the process towards intercalation. Our vdW-DF inves-
tigations of the binding of potassium in or on graphite
supplements corresponding vdW-DF studies of the bind-
ing of polycyclic aromatic hydrocarbon dimers, of the
polyethylene crystal, of benzene dimers, and of poly-
cyclic aromatic hydrocarbon and phenol molecules on
graphite.29,30,31,32,33,34
The outline of the paper is as follows. Section II con-
tains a short description of the materials of interest here:
graphite, C8K, and graphite with an adsorbed or ab-
sorbed K atom layer. The vdW-DF scheme is described
in Sec. III. Section IV presents our results, Sec. V the
discussion, and conclusions are drawn in Sec. VI.
http://arxiv.org/abs/0704.0055v1
FIG. 1: (Color online) Simple hexagonal graphite (AA stack-
ing) and natural hexagonal graphite (AB stacking). The two
structures differ by that each second carbon layer in AB-
stacked graphite is shifted, whereas in AA-stacked graphite
all planes are directly above each other. The experimentally
obtained in-plane lattice constant and sheet separation of nat-
ural graphite is (Ref. 40) a = 2.459 Å and dC-C = 3.336 Å,
respectively.
II. MATERIAL STRUCTURE
Graphite is a semimetallic solid with strong intra-plane
bonds and weakly coupled layers. The presence of these
two types of bonding results in a material with different
properties along the various crystallographic directions.35
For example, the thermal and electrical conductivity
along the carbon sheets is two orders of magnitude higher
than that perpendicular to the sheets. This specific prop-
erty allows heat to move directionally, which makes it
possible to control the heat transfer. The relatively weak
vdW forces between the sheets contribute to another in-
dustrially important property: graphite is an ideal lubri-
cant. In addition, the anisotropic properties of graphite
make the material suitable as a substrate in electronic
studies of ultrathin metal films.36,37,38,39
The natural structure of graphite is an AB stacking,
with the graphite layers shifted relative to each other,
as illustrated in Fig. 1. The figure also shows hexagonal
graphite, consisting of AA-stacked graphite layers. The
in-plane lattice constant a and the layer separation dC-C
is also illustrated. In natural graphite the primitive unit
cell is hexagonal, includes four carbon atoms in two lay-
ers, and has unit cell side lengths a and height c = 2dC-C.
The physical properties of graphite have been studied
in a variety of experimental40,41,42 and theoretical43,44
work. Some of the DFT work has been performed in
LDA, which does not provide a physically meaningful
account of binding in layered systems.15,45 At the same
time, using GGA is not an option because it does not
bind the graphite layers. For a good description of the
FIG. 2: (Color online) Crystalline structure of C8K show-
ing the AA-stacking of the carbon layers (small balls) and
the αβγδ-stacking of the potassium layers (large balls) per-
pendicular to the graphene sheets. The potassium layers are
arranged in a p(2× 2) structure, with the K atoms occupying
the sites over the hollows of every fourth carbon hexagon.
graphite structure and nature the vdW interactions must
be included.45
Alkali metals (AM), except Na, easily penetrate the
gallery of the graphite forming alkali metal graphite in-
tercalation compounds. These intercalation compounds
are formed through electron exchange between the inter-
calated layer and the host carbon layers, resulting in a
different nature of the interlayer bonding type than that
of pristine graphite. The intercalate also affects the con-
ductive properties of graphite, which becomes supercon-
ductive in the direction parallel to the planes at critical
temperatures below 1 K.6,7
The structure of AM graphite intercalation compounds
is characterized by its stage n, where n is the number of
graphite sheets located between the AM layers. In this
work we consider only stage-1 intercalated graphite C8K,
in which the layers of graphite and potassium alternate
throughout the crystal. The primitive unit cell of C8K
is orthorhombic and contains sixteen C atoms and two
K atoms. In the C8K crystal the K atoms are ordered
in a p(2× 2) registry with K-K separation 2a, where a is
the in-plane lattice constant of graphite. This separation
of the potassium atoms is about 8% larger than that in
the natural K bcc crystal (based on experimental values).
The carbon sheet stacking in C8K is of AA type, with the
K atoms occupying the sites over the hollows of every
fourth carbon hexagon, each position denoted by α, β,
γ, or δ, and the stacking of the K atoms perpendicular
to the planes being described by the αβγδ-sequence as
illustrated in Fig. 2.
III. COMPUTATIONAL METHODS
The first-principle total-energy and electronic struc-
ture calculations are performed within the framework
of DFT. The semilocal Perdew-Burke-Ernzerhof (PBE)
flavor26 of GGA is chosen for the exchange-correlation
functional for the traditional self-consistent calculations
underlying the vdW-DF calculations. For all GGA cal-
culations we use the open source DFT code Dacapo,46
which employs Vanderbilt ultrasoft pseudopotentials,47
periodic boundary conditions, and a plane-wave basis set.
An energy cut-off of 500 eV is used for the expansion of
the wave functions and the Brillouin zone (BZ) of the
unit cells is sampled according to the Monkhorst-Pack
scheme.48 The self-consistently determined GGA valence
electron density n(r) as well as components of the energy
from these calculations are passed on to the subsequent
vdW-DF calculation of the total energy.
For the adsorption and absorption studies a graphite
surface slab consisting of 4 layers is used, with a surface
unit cell of side lengths twice those in the graphite bulk
unit cell (i.e., side lengths 2a). The surface calculations
are performed with a 4×4×1 k-point sampling of the BZ.
The (pure) graphite bulk GGA calculations are per-
formed with a 8×8×4 k-point sampling of the BZ,
whereas for the C8K bulk structure, in a unit cell at least
double the size in any direction, 4×4×2 k-points are used,
consistent also with the choice of k-point sampling of the
surface slabs.
We choose to describe C8K by using a hexagonal unit
cell with four formula units, lateral side lengths approxi-
mately twice those of graphite and with four graphite and
four K-layers in the direction perpendicular to the layers.
C8K can also be described by the previously mentioned
primitive orthorhombic unit cell containing two formula
units of atoms but we retain the orthorhombic cell for
ease of description and for simple implementation of nu-
merically robust vdW-DF calculations.
In all our studies, except test cases, the Fast Fourier
Transform (FFT) grids are chosen such that the separa-
tion of neighboring points is maximum ∼0.13 Å in any
direction in any calculation.
A. vdW density function calculations
In graphite, the carbon layers bind by vdW interac-
tions only. In the intercalated compound a major part of
the attraction is ionic, but also here the vdW interactions
cannot be ignored. In order to include the vdW interac-
tions systematically in all of our calculations we use the
vdW-DF of Ref. 16. There, the correlation energy func-
tional is divided into a local and a nonlocal part,
Ec ≈ E
c + E
c , (1)
where the local part is approximated in the LDA and the
nonlocal part Enlc is consistently constructed to vanish
for a homogeneous system. The nonlocal correlation Enlc
is calculated from the GGA-based n(r) and its gradients
by using information about the many-body response of
the weakly inhomogeneous electron gas:
Enlc =
dr′n(r)φ(r, r′)n(r′). (2)
The nonlocal kernel φ(r, r′) can be tabulated in terms
of the separation |r − r′| between the two fragments at
positions r and r′ through the parameters D = (q0 +
q′0)|r − r
′|/2 and δ = (q0 − q
0)/(q0 + q
0). Here q0 is a
local parameter that depends on the electron density and
its gradient at position r. The analytic expression for the
kernel φ in terms of D and δ can be found in Ref. 16.
For periodic systems, such as bulk graphite, C8K,
and the graphite surface (with adsorbed or absorbed K-
atoms), the nonlocal correlation per unit cell is simply
evaluated from the interaction of the points in the unit
cell V0 with points everywhere in space (V ) in the three
(for bulk graphite and C8K) or two (for the graphite
surface) dimensions of periodicity. Thus, the V -integral
in Eq. (2) in principle requires a representation of the
electron density infinitely repeated in space. In prac-
tice, the nonlocal correlation rapidly converges31 and it
suffices with repetitions of the unit cell a few times in
each spatial direction. For graphite bulk the V -integral
is converged when we use a V that extends 9 (7) times
the original unit cell in directions parallel (perpendicular)
to the sheets. For the potassium investigation a signif-
icantly larger original unit cell is adopted (see Fig. 2);
here a fully converged V corresponds to a cell extending
five (three) times the original cell in the direction parallel
(perpendicular) to the sheets for C8K bulk. To describe
the nonlocal correlations (2) for the graphite surface a
sufficient V extends five times the original unit cell along
the carbon sheets.
For the exchange energy Ex we follow the choice of
Ref. 16 of using revPBE27 exchange. Among the func-
tionals that we have easy access to, the revPBE has
proved to be the best candidate for minimizing the ten-
dency of artificial exchange binding in graphite.15
Using the scheme described above to evaluate Enlc , the
total energy finally reads:
EvdW−DF = EGGA − EGGAc + E
c + E
c , (3)
where EGGA is the GGA total energy with the revPBE
choice for the exchange description and EGGAc (E
the GGA (LDA) correlation energy. As our GGA calcu-
lations in this specific application of vdW-DF are carried
out in PBE, not revPBE, we further need to explicitly
replace the PBE exchange in EGGA by that of revPBE
for the same electron charge density distribution.
B. Convergence of the local and nonlocal energy
variation
DFT calculations provide physically meaningful results
for energy differences between total energies (3). To un-
derstand materials and processes we must compare total
energy differences between a system with all constituents
at relatively close distance and a system of two or more
fragments at “infinite” separation (the reference system).
Since the total energy (3) consists both of a long-range
term and shorter-ranged GGA and LDA terms it is nat-
ural to choose different ways to represent the separated
fragments for these different long- or short-range energy
terms.
For the shorter-range energy parts (LDA and GGA
terms) the reference system is a full system with vacuum
between the fragments. For LDA and GGA calculations
it normally suffices to make sure that the charge den-
sity tails of the fragments do not overlap, but here we
find that the surface dipoles cause a slower convergence
with layer separation. We use a system with the layer
separation between the potassium layer and the nearest
graphite layer(s) dC-K = 12 Å (8 Å) as reference for the
adsorption (absorption) study.
The evaluation of the nonlocal correlationsEnlc requires
additional care. This is due to technical reasons per-
taining to numerically stability in basing the Enlc eval-
uation on the FFT grid used to converge the underly-
ing traditional-DFT calculations. The evaluation of the
nonlocal correlation energy, Eq. (2), involves a weighted
double integral of a kernel with a significant short-range
variation16. The shape of the kernel makes the Enlc eval-
uation sensitive to the particulars of FFT-type griding,49
for example, to the relative position of FFT grid points
relative to the nuclei position (for a finite grid-point spac-
ing).
However, robust evaluation of binding- or cohesive-
energy contributions by nonlocal correlations can gen-
erally be secured by a further splitting of energy differ-
ences into steps that minimize the above-mentioned grid
sensitivity. The problem of FFT sensitivity of the Enlc
evaluation is accentuated because the binding in the Enlc
channel arises as a smaller energy difference between siz-
ableEnlc contributions of the system and of the fragments.
Conversely, convergence in vdW-DF calculations of bind-
ing and cohesive energies can be obtained even at a mod-
erate FFT grid accuracy (0.13 Å used here) by devising a
calculational scheme that always maintains identical po-
sition of the nuclei relative to grid points in the combined
systems as well as in the fragment reference system.
Thus we obtain a numerically robust evaluation of the
Enlc energy differences by choosing steps for which we
can explicitly control the FFT griding. For adsorption
and absorption cases we calculate the reference systems
as a sum of Enlc -contributions for each fragment and we
make sure to always position the fragment at the exact
same position in the system as in the interacting system.
For bulk systems we choose steps in which we exclusively
adjust the inter-plane or in-plane lattice constant. Here
the reference system is then simply defined as a system
with either double (or in some cases quadruple) lattice
constant and with a corresponding doubling of the FFT
griding along the relevant unit-cell vector.
The cost of full convergence is that, in practice, we of-
ten do three or more GGA calculations and subsequent
Enlc calculations for each point on the absorption, absorp-
tion, or formation-energy curve. In addition to the cal-
culations for the full system we have to do one for each of
the isolated fragments at identical position in the adsorp-
tion/absorption cases and one or more for fragments in
the doubled unit-cell and doubled griding reference. We
have explicitly tested that using a FFT grid spacing of
< 0.13 Å (but not larger) for such reference calculations
is sufficient to ensure full convergence in the reported
Enlc (and E
vdW−DF total) energy variation for graphitic
systems.
C. Material formation and sorption energies
The cohesive energy of graphite (G) is the energy gain,
per carbon atom, of creating graphite at in-plane lattice
constant a and layer separation dC-C from isolated (spin-
polarized) carbon atoms.
EG,coh(a, dC-C) = EG,tot(a, dC-C)− EC-atom,tot (4)
where EG,tot and EC-atom,tot are total energies per carbon
atom. The graphite structure is stable at the minimum
of the cohesive energy, at lattice constants a = aG and
2dC-C = cG.
The adsorption (absorption) energy for a p(2 × 2) K-
layer over (under) the top layer of a graphite surface is
the difference in total energy [from Eq. (3)] for the system
at hand minus the total energy of the initial system, i.e.,
a clean graphite surface and isolated gas-phase potassium
atoms. However, due to the above mentioned technical
issues in using the vdW-DF we calculate the adsorption
and absorption energy as a sum of (artificial) stages lead-
ing to the desired system: First the initially isolated,
spin-polarized potassium atoms are gathered into a free
floating potassium layer with the structure correspond-
ing to a full cover of potassium atoms. By this the total
system gains the energy ∆EK-layer(aG), with
∆EK-layer(a) = EK,tot(a)− EK-atom,tot . (5)
In adsorption the potassium layer is then simply placed
on top of the four-layer (2 × 2) graphite surface (with
the K atoms above graphite hollows) at distance dC-K.
The system thereby gains a further energy contribution
∆EK-G(dC-K). This leads to an adsorption energy per
K-atom
Eads(dC-K) = ∆EK-layer(aG) + ∆EK-G(dC-K) . (6)
In absorption the top graphite layer is peeled off the
(2 × 2) graphite surface and moved to a distance far
from the remains of the graphite surface. This process
costs the system an (“exfoliation”) energy −∆EC-G =
−[Etot,C-G(dC-C = cG/2)− Etot,C-G(dC-C → ∞)]. At the
far distance the isolated graphite layer is moved into AA
stacking with the surface, at no extra energy cost. Then,
the potassium layer is placed midway between the far-
away graphite layer and the remains of the graphite sur-
face. Finally the two layers are gradually moved towards
the surface. At distance 2dC-K between the two topmost
graphite layers (sandwiching the K-layer) the system has
further gained an energy ∆EC-K-G(dC-K). The absorp-
tion energy per K-atom is thus
Eabs(dC-K) = −∆EC-G+∆EK-layer(aG)+∆EC-K-G(dC-K) .
Similarly, the C8K intercalate compound is formed
from graphite by first moving the graphite layers far
apart accordion-like (and there shift the graphite stack-
ing from ABA . . . to AAA . . . at no energy cost), then
changing the in-plane lattice constant of the isolated
graphene layers from aG to a, then intercalating K-layers
(in stacking αβγδ) between the graphite layers, and fi-
nally moving all the K- and graphite layers back like an
accordion, with in-plane lattice constant a (which has the
value aC8K at equilibrium).
In practice, a unit cell of four periodically repeated
graphite layers is used in order to accommodate the
potassium αβγδ-stacking. The energy gain of creating
a (2× 2) graphene sheet from 8 isolated carbon atoms is
defined similarly to that of the K-layer:
∆EC-layer(a) = EC-layer,tot(a)− 8EC-atom,tot . (8)
The formation energy for the C8K intercalate com-
pound per K atom or formula unit, Eform, is thus found
from the energy cost of moving four graphite layers
apart by expanding the (2 × 2) unit cell to large height,
−∆EG-acc, the cost of changing the in-plane lattice con-
stant from aG to a in each of the four isolated graphene
layers, 4(∆EC-layer(a)−∆EC-layer(aG)), the gain of creat-
ing four K-layers from isolated K-atoms, 4∆EK-layer(a),
plus the gain of bringing four K-layers and four graphite
layers together in the C8K structure, ∆EC8K-acc(a, dC-K),
yielding
Eform(a, dC-K)
−∆EG-acc + 4∆EC-layer(a)− 4∆EC-layer(aG)
+ 4∆EK-layer(a) + ∆EC8K-acc(a, dC-K)
. (9)
The relevant energies to use for comparing the three
different mechanisms of including potassium (adsorp-
tion, absorption and intercalation) are thus Eads(dC-K),
Eabs(dC-K) and Eform(a, dC-K) at their respective mini-
mum values.
IV. RESULTS
Experimental observations indicate that the intercala-
tion of potassium into graphite starts with the absorption
of evaporated potassium into an initially clean graphite
surface.50 This subsurface absorption is preceded by ini-
tial, sparse potassium adsorption onto the surface, and
proceeds with further absorption into deeper graphite
voids. The general view is that the K atoms enter
graphite at the graphite step edges.20 The amount and
position of intercalated K atoms is controlled by the tem-
perature and time of evaporation.
Below, we first describe the initial clean graphite sys-
tem, and the energy gain in (artificially) creating free-
floating K-layers from isolated K-atoms. Then we present
and discuss our results on potassium adsorption and sub-
surface absorption, followed by a characterization of bulk
For the adsorption (absorption) system we calculate
the adsorption (absorption) energy curve, including the
equilibrium structure. As a demonstration of the need
for a relatively fine FFT griding in the vdW-DF cal-
culations we also calculate and compare the absorption
curve for a more sparse FFT grid. For the bulk sys-
tems (graphite and C8K) we determine the lattice pa-
rameters and the bulk modulus. We also calculate the
formation energy of C8K and the energy needed to peel
off one graphite layer from the graphite surface and com-
pare with experiment.51
A. Graphite bulk structure
The present calculations on pure graphite are for the
natural, AB-stacked graphite (lower panel of Fig. 1). The
cohesive energy is calculated at a total of 232 structure
values (a, dC-C) and the equilibrium structure and bulk
modulus B0 are then evaluated using the method de-
scribed in Ref. 52.
Figure 3 shows a contour plot of the graphite cohesive
energy variation EG,coh as a function of the layer sep-
aration dC-C and the in-plane lattice constant a, calcu-
lated within the vdW-DF scheme. The contour spacing
is 5meV per carbon atom, shown relative to the energy
minimum located at (a, dC-C) = (aG, cG/2) =(2.476 Å,
3.59 Å). These values are summarized in Table I together
with the results obtained from a semilocal PBE cal-
culation. As expected, and discussed in Ref. 14, the
semilocal PBE calculation yields unrealistic results for
the layer separation. The table also presents the cor-
responding experimental values. Our calculated lattice
values obtained using vdW-DFT are in good agreement
with experiment,40 and close to those found from the
older vdW-DF of Refs. 14 and 15, (in which we for Enlc as-
sume translational invariance of n(r) along the graphite
planes,) at (2.47 Å, 3.76 Å).
Consistent with experimental reports18 and our previ-
ous calculations14,15,45 we find graphite to be rather soft,
indicated by the bulk modulus B0 value. Since in-plane
compression is very hard in graphite most of the softness
suggested by (the isotropic) B0 comes from compression
perpendicular to the graphite layers, and the value of
B0 is expected to be almost identical to the C33 elastic
3 3.5 4 4.5 5 5.5 6
dC−C [Å]
FIG. 3: Graphite cohesive energy EG,coh (AB-stacked), based
on vdW-DF, as a function of the carbon layer separation dC-C
and the in-plane lattice constant a. The energy contours are
spaced by 5meV per carbon atom.
TABLE I: Optimized structure parameters and elastic
properties for natural hexagonal graphite (AB-stacking)
and the potassium-intercalated graphite structure C8K in
AαAβAγAδAα . . . stacking. The table shows the calcu-
lated optimal values of the in-plane lattice constant a, the
(graphite-)layer-layer separation dC-C, and the bulk modulus
B0. In C8K the value if dC-C is twice the graphite-potassium
distance dC-K.
Graphite C8K
PBE vdW-DF Exp. PBE vdW-DF Exp.
a (Å) 2.473 2.476 2.459a 2.494 2.494 2.480b
dC-C (Å) ≫ 4 3.59 3.336
a 5.39 5.53 5.35c
B0 (GPa) 27 37
de 37 26 47de
aRef. 40. bRef. 53. cRef. 4. dRef. 18.
eValue presented is for C33; for laterally rigid materials, like
graphite and C8K, C33 is a good approximation of B0.
coefficient.14,18
We find the energy cost of peeling off a graphite layer
from the graphite surface (the exfoliation energy) to be
∆EC-G = −435 meV per (2 × 2) unit cell, i.e., −55
meV per surface carbon atom (Table II). A recent
experiment51 measured the desorption energy of poly-
cyclic aromatic hydrocarbons (basically flakes of graphite
sheets) off a graphite surface. From this experiment
the energy cost of peeling off a graphite layer from the
graphite surface was deduced to −52± 5meV/atom.
Our value −55 meV/C-atom is also consistent with
a separate vdW-DF determination29 of the binding
(−47meV per in-plane atom) between two (otherwise)
isolated graphene sheets.
For the energies of the absorbate system and of the
C8K intercalate a few other graphite-related energy con-
tributions are needed. The energy of collecting C atoms
to form a graphene sheet at lattice constant a from iso-
lated (spinpolarized) atoms is given by ∆EC-layer(a); we
find that changing the lattice constant a from aG to
the equilibrium value aC8K of C8K causes this energy
to change a mere 30 meV per (2 × 2) sheet. The contri-
bution ∆EG-acc is the energy of moving bulk graphite
layers (in this case four periodically repeated layers)
far away from each other, by expanding the unit cell
along the direction perpendicular to the layers. Thus,
∆EG-acc = 32∆EG,coh(aG, cG/2) − 4∆EC-layer(aG) tak-
ing the number of atoms and layers per unit cell into
account. We find the value ∆EG-acc = −1600 meV per
(2×2) four-layer unit cell. This corresponds to −50 meV
per C atom, again consistent with our result for the ex-
foliation energy, ∆EC-G/8 = −55 meV.
B. Creating a layer of K-atoms
The (artificial) step of creating a layer of potassium
atoms from isolated atoms releases a significant energy
∆EK-layer. This energy contains the energy variation
with in-plane lattice constant and the energy cost of
changing from a spin-polarized to a spin-balanced elec-
tron configuration for the isolated atom.54
The creation of the K-layer provides an energy gain
which is about half an eV per potassium atom, depending
on the final lattice constant. With the graphite lattice
constant aG the energy change, including the spin-change
cost, is ∆EK-layer(aG) = −476meV per K atom in vdW-
DF (−624meV when calculated within PBE), whereas
∆EK-layer(aC8K) = −473meV in vdW-DF.
C. Graphite-on-surface adsorption of potassium
The potassium atoms are adsorbed on a usual
ABA . . .-stacked graphite surface. We consider here full
(one monolayer) coverage, which is one potassium atom
per (2 × 2) graphite surface unit cell. This orders the
potassium atoms in a honeycomb structure with lattice
constant 2aG, and a nearest-neighbor distance within the
K-layer of aG.
The unit cell used in the standard DFT calculations
for adsorption and absorption has a height of 40 Å and
includes a vacuum region sufficiently big that no interac-
tions (within GGA) can occur between the top graphene
sheet and the slab bottom in the periodically repeated
image of the slab. The vacuum region is also large in
order to guarantee that the separation from any atom to
the dipole layer55 always remains larger than 4 Å.
In the top panel of Fig. 4 we show the adsorption en-
ergy per potassium atom. The adsorption energy at equi-
librium is −937 meV per K atom at distance dC-K = 3.02
Å from the graphite surface.
For comparison we also show the adsorption curve cal-
culated in a PBE-only traditional DFT calculation. Since
vdW−DF
2 3 4 5 6 7
dC−K [Å]
vdW−DF
43.532.5
dC−K [Å]
sparse
FIG. 4: Potassium adsorption and absorption energy at the
graphite surface as a function of the separation dC-K of the
K-atom layer and the nearest graphite layer(s) (at in-plane
lattice constant corresponding to that of the surface, aG).
Top panel: Adsorption curve based on vdW-DF calculations
(solid line with black circles) and PBE GGA calculations
(dashed line). The horizontal lines to the left show the en-
ergy gain in creating the isolated K layer from isolated atoms,
∆EK-layer(aG), the asymptote of Eads(dC-K) in this plot.
Bottom panel: Absorption curve based on vdW-DF calcu-
lations. The asymptote is here the sum ∆EK-layer(aG) −
∆EC-G. Inset: Binding energy of the K-layer and the
top graphite layer (“C-layer”) on top of the graphite slab,
∆EC-K-G. The dashed curve shows our results when in E
ignoring every second FFT grid point (in each direction) of
the charge density from the underlying GGA calculations, the
solid curve with black circles shows the result of using every
available FFT grid point.
the interaction between the K-layer and the graphite sur-
face has a short-range component to it, even GGA calcu-
lations, such as the PBE curve, show significant binding
(−900 meV/K-atom at dC-K = 2.96 Å). This is in con-
trast to the pure vdW binding between the layers in clean
graphite.14,15 Note that the asymptote of the PBE curve
is different from that of the vdW-DF curve, this is due
to the different energy gains (∆EK-layer) in collecting a
potassium layer from isolated atoms when calculated in
PBE or in vdW-DF.
For K-adsorption the vdW-DF and PBE curves agree
reasonably well, and the use of vdW-DF for this spe-
cific calculation is not urgently necessary. However, in
order to compare the adsorption results consistently to
absorption, intercalation and clean graphite, it is neces-
sary to include the long-range interactions through vdW-
DF. As shown for the graphite bulk results above, PBE
yields quantitatively and qualitatively wrong results for
the layer separation.
D. Graphite-subsurface absorption of potassium
The first subsurface adsorption of K takes place in the
void under the top-most graphite layer. The surface ab-
sorption of the first K-layer causes a lateral shift of the
top graphite sheet, resulting in a A/K/ABAB . . . stack-
ing of the graphite. We have studied the bonding nature
of this absorption process by considering a full p(2× 2)-
intercalated potassium layer in the subsurface of a four
layer thick graphite slab.
Following the receipt of Section III for the absorption
energy (7) the energies ∆EC-K-G are approximated by
those from a four-layer intercalated graphite slab with
the stacking A/K/ABA, and the values are shown in the
inset of Fig. 4. The absorption energy Eabs is given by
the curve in the bottom panel of Fig. 4, and its minimum
is −952meV per K atom at dC-K = 2.90 Å.
To investigate what grid spacing is sufficiently dense
to obtain converged total-energy values in vdW-DF we
do additional calculations in the binding distance region
with a more sparse grid. Specifically, the inset of Fig. 4
compares the vdW-DF calculated at full griding with one
that uses only every other FFT grid point in each direc-
tion, implying a grid spacing for Enlc (but not for the lo-
cal terms) which is maximum 0.26 Å. We note that using
the full grid yields smaller absolute values of the absorp-
tion energy. We also notice that the effect is more pro-
nounced for small separations than for larger distances.
Thus given resources, the dense FFT grid calculations
are preferred, but even the less dense FFT grid calcu-
lations yield reasonably well-converged results. In all
calculations (except tests of our graphitic systems) we
use a spacing with maximum 0.13 Å between grid points.
This is a grid spacing for which we have explicitly tested
convergence of the vdW-DF for graphitic systems given
the computational strategy described and discussed in
Sec. III.
E. Potassium-intercalated graphite
When potassium atoms penetrate the gallery of the
graphite, they form planes that are ordered in a p(2× 2)
fashion along the planes. The K intercalation causes a
shift of every second carbon layer resulting into an AA
stacking of the graphite sheets. The K atoms then simply
occupy the sites over the hollows of every fourth carbon
hexagon. The order of the K atoms perpendicular to the
planes is described by the αβγδ stacking, illustrated in
Fig. 2.
For the potassium intercalated compound C8K we cal-
culate in standard DFT using PBE the total energy at
132 different combinations of the structural parameters a
TABLE II: Comparison of the graphite exfoliation energy per
surface atom, EC-G/8, graphite layer binding energy per car-
bon atom, ∆EC-acc/32, the energy gain per K atom of col-
lecting K- and graphite-layers at equilibrium to form C8K,
∆EC8K-acc/4, and the equilibrium formation energy of C8K,
Eform.
∆EC-G/8 ∆EC-acc/32 ∆EC8K-acc/4 Eform
[meV/atom] [meV/atom] [meV/C8K] [meV/C8K]
vdW-DF −55 −50 −818 −861
PBE − − −511 −
Exp. −52± 5a −1236b
aRef. 51. bRef. 1.
 5  5.2  5.4  5.6  5.8  6
dC−C [Å]
FIG. 5: Formation energy of C8K, Eform, as a function of
the carbon-to-carbon layer separation dC-C and half the in-
plane lattice constant, a. The energy contours are spaced by
20meV per formula unit.
and dC-C. The charge densities and energy terms of these
calculations are then used as input to vdW-DF. The equi-
librium structure and elastic properties (B0) both for the
vdW-DF results and for the PBE results are then evalu-
ated with the same method as in the graphite case.52
Figure 5 shows a contour plot of the C8K formation
energy, calculated in vdW-DF, as a function of the C-C
layer separation (dC-C) and the in-plane periodicity (a)
of the graphite-layer structure. The contour spacing is
20 meV per formula unit and are shown relative to the
energy minimum at (a, dC-C) = (2.494 Å, 5.53 Å).
V. DISCUSSION
Table I presents an overview of our structural results
obtained with the vdW-DF for graphite and C8K. The
table also contrasts the results with the corresponding
values calculated with PBE where available. The vdW-
DF value dC-C = 5.53 Å for the C8K C − C layer sep-
aration is 3% larger than the experimentally observed
value whereas the PBE value corresponds to less than a
1% expansion. Our vdW-DF result for the C8K bulk
modulus (26 GPa) is also softer than the PBE result
(37 GPa) and further away from the experimental esti-
mates (47 GPa) based on measurements of the C33 elastic
response.18 A small overestimation of atomic separation
is consistent with the vdW-DF behavior that has been
documented in a wide range of both finite and extended
systems.14,15,16,17,29,30,33,34 This overestimation results,
at least in part, from our choice of parametrization of
the exchange behavior — an aspect that lies beyond
the present vdW-DF implementation which focuses on
improving the account of the nonlocal correlations, per
se. It is likely that systematic investigations of the ex-
change effects can further refine the accuracy of vdW-DF
implementations.56 In any case, vdW-DF theory calcu-
lations represent, in contrast to PBE, the only approach
to obtain a full ab initio characterization of the AM in-
tercalation process.
The C8K system is more compact than graphite and
this explains why PBE alone can here provide a good de-
scription of the materials structure and at least some ma-
terials properties, whereas it fails completely for graphite.
The distance between the graphene sheets upon interca-
lation of potassium atoms is stretched compared to that
of pure graphite, but the (K-)layer to (graphite-)layer
separation, dC-K = dC-C/2 = 2.77 Å, is significantly less
than the layer-layer separation in pure graphite. This in-
dicates that C8K is likely held together, at least in part,
by shorter-ranged interactions.
Table II documents that the vdW binding neverthe-
less plays an important role in the binding and forma-
tion of C8K. The table summarizes and contrasts our
vdW-DF and PBE results for graphite exfoliation and
layer binding energies as well as C8K interlayer binding
and formation energies. The vdW-DF result for the C8K
formation energy is smaller than experimental measure-
ments by 31% but it nevertheless represent a physically
motivated ab initio calculation. In contrast, the C8K
formation energy is simply unavailable in PBE because
PBE, as indicated, fails to describe the layer binding in
graphite. Moreover, for the vdW-DF/PBE comparisons
that we can make — for example, of the C8K layer inter-
action ∆EC8K-acc — the vdW-DF is found to significantly
strengthen the bonding compared with PBE.
It is also interesting to note that the combination of
shorter-ranged and vdW bonding components in C8K
yields a layer binding energy that is close to that of
the graphite case. In spite of the difference in nature
of interactions, we find almost identical binding energies
per layer for the case of the exfoliation and accordion in
graphite and for the accordion in C8K. This observation
testifies to a perhaps surprising strength of the so-called
soft-matter vdW interactions.
In a wider perspective our vdW-DF permits a first
comparison of the range of AM-graphite systems from
adsorption over absorption to full intercalation and thus
insight on the intercalation progress. Assuming a dense
2×2 configuration, we find that the energy for potassium
adsorption and absorption is nearly degenerate with an
indication that absorption is slightly preferred, consis-
tent with experimental behavior. We also find that the
potassium absorption may eventually proceeds towards
full intercalation thanks to a significant release of forma-
tion energy.
VI. CONCLUSIONS
The potassium intercalation process in graphite has
been investigated by means of the vdW-DF density func-
tional method. This method includes the dispersive in-
teractions needed for a consistent investigation of the
intercalation process. For clean graphite the vdW-DF
predicts — contrary to standard semilocal DFT imple-
mentations — a stabilized bulk system with equilibrium
crystal parameters in close agreement with experiments.
Two limits of the absorption process have been inves-
tigated by the vdW-DF, namely single layer subsurface
absorption and the fully potassium intercalated stage-1
crystal C8K. Here the vdW-DF is shown to enhance the
(semi-)local type of bonding described by traditional ap-
proaches. The significant impact on the materials behav-
ior indicates that the vdW-DF is needed not only for a
consistent description of sparse matter systems that are
solely stabilized by dispersion forces, but also for their
intercalates.
We thank D.C. Langreth and B.I. Lundqvist for stim-
ulating discussions. Partial support from the Swedish
Research Council (VR), the Swedish National Graduate
School in Materials Science (NFSM), and the Swedish
Foundation for Strategic Research (SSF) through the
consortium ATOMICS is gratefully acknowledged, as
well as allocation of computer time at UNICC/C3SE
(Chalmers) and SNIC (Swedish National Infrastructure
for Computing).
∗ Electronic address: hyldgaar@chalmers.se
1 S. Aronson, F.J. Salzano, and D. Ballafiore, J. Chem. Phys.
49, 434 (1968).
2 D.E. Nixon and G.S. Parry, J. Phys. D 1, 291 (1968).
3 R. Clarke, N. Wada, and S.A. Solin, Phys. Rev. Lett. 44,
1616 (1980).
4 M.S. Dresselhaus and G. Dresselhaus, Adv. Phys. 30, 139
(1981).
5 D.P. DiVincenzo and E.J. Mele, Phys. Rev. B 32, 2538
(1985).
6 N.B. Hannay, T.H. Geballe, B.T. Matthias, K. Andreas, P.
Schmidt, and D. MacNair, Phys. Rev. Lett. 14, 225 (1965).
7 R.A. Jishi and M.S. Dresselhaus, Phys. Rev. B 45, 12465
(1992).
8 T. Kihlgren, T. Balasubramanian, L. Walldén, and R.
Yakimova, Surf. Sci. 600, 1160 (2006).
9 M. Breitholtz, T. Kihlgren, S.-Å. Lindgren, and L.
Walldén, Phys. Rev. B 66, 153401 (2002).
10 I. Forbeaux, J.-M. Themlin, and J.-M. Debever, Phys. Rev.
B 58, 16396 (1998).
11 T. Kihlgren, T. Balasubramanian, L. Walldén, and R.
Yakimova, Phys. Rev. B 66, 235422 (2002).
12 I. Forbeaux, J.-M. Themlin, A. Charrier, F. Thibaudau,
and J.-M. Debever, Appl. Surf. Sci. 162–163, 406 (2000).
13 E. Ziambaras, Ph.D. thesis, Chalmers (2006).
14 H. Rydberg, M. Dion, N. Jacobson, E. Schröder, P.
Hyldgaard, S.I. Simak, D.C. Langreth, and B.I. Lundqvist,
Phys. Rev. Lett. 91, 126402 (2003).
15 D.C. Langreth, M. Dion, H. Rydberg, E. Schröder, P.
Hyldgaard, and B.I. Lundqvist, Int. J. Quantum Chem.
101, 599 (2005).
16 M. Dion, H. Rydberg, E. Schröder, D.C. Langreth, and
B.I. Lundqvist, Phys. Rev. Lett. 92, 246401 (2004); 95,
109902(E) (2005).
17 T. Thonhauser, V.R. Cooper, S. Li, A. Puzder,
P. Hyldgaard, and D.C. Langreth, Van der
Waals density functional: Self-consistent poten-
tial and the nature of the van der Waals bond,
http://arxiv.org/abs/cond-mat/0703442
18 N. Wada, R. Clarke, and S.A. Solin, Solid State Comm.
35, 675 (1980).
19 H. Zabel and A. Magerl, Phys. Rev. B 25, 2463 (1982).
20 J.C. Barnard, K.M. Hock and R.E. Palmer, Surf. Science
287–288, 178 (1993).
21 K. M. Hock and R. E. Palmer, Surf. Science 284, 349
(1993).
22 Z.Y. Li, K.M. Hoch, and R.E. Palmer, Phys. Rev. Lett.
67, 1562 (1991).
23 S.D. Chakarova and E. Schröder, Materials Science and
Engineering C 25, 787 (2005).
24 L.A. Girifalco and M. Hodak, Phys. Rev. B 65, 125404
(2002).
25 J.P. Perdew, J.A. Chevary, S.H. Vosko, K.A. Jackson,
M.R. Pederson, D.J. Singh, and C. Fiolhais, Phys. Rev.
B 48, 6671 (1992).
26 J.P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett.
77, 3865 (1996).
27 Y. Zhang and W. Yang, Phys. Rev. Lett. 80, 890 (1998).
28 B. Hammer, L.B. Hansen, and J.K. Nørskov, Phys. Rev.
B 59, 7413 (1999).
29 S.D. Chakarova-Käck, J. Kleis, and E. Schröder, Appl.
Phys. Rep. 2005-16 (2005).
30 J. Kleis, B.I. Lundqvist, D.C. Langreth, and E. Schröder,
Towards a working density-functional theory for polymers:
First-principles determination of the polyethylene crystal
structure, http://arxiv.org/abs/cond-mat/0611498
31 S.D. Chakarova-Käck, E. Schröder, B.I. Lundqvist, and
D.C. Langreth, Phys. Rev. Lett. 96, 146107 (2006).
32 S.D. Chakarova-Käck, Ø. Borck, E. Schröder, and B.I.
Lundqvist, Phys. Rev. B 74, 155402 (2006).
33 A. Puzder, M. Dion, and D.C. Langreth, J. Chem. Phys.
124, 164105 (2006).
34 T. Thonhauser, A. Puzder, and D.C. Langreth, J. Chem.
Phys. 124, 164106 (2006).
35 D.D.L. Chung, J. Mat. Sci. 37, 1475 (2002).
36 M. Breitholtz, T. Kihlgren, S.-Å. Lindgren, H. Olin, E.
mailto:hyldgaar@chalmers.se
http://arxiv.org/abs/cond-mat/0703442
http://arxiv.org/abs/cond-mat/0611498
Wahlström, and L. Walldén, Phys. Rev. B 64, 073301
(2001).
37 Z.P. Hu, N.J. Wu, and A. Ignatiev, Phys. Rev. B 33, 7683
(1986).
38 J. Cui, J.D. White, R.D. Diehl, J.F. Annett, and M.W.
Cole, Surf. Sci. 279, 149 (1992).
39 L. Österlund, D.V. Chakarov, and B. Kasemo, Surf. Sci.
420, L437 (1991).
40 Y. Baskin and L. Meyer, Phys. Rev. 100, 544 (1955).
41 W. Eberhardt, I.T. McGovern, E.W. Plummer, and J.E.
Fisher, Phys. Rev. Lett. 44, 200 (1980).
42 A.R. Law, J.J. Barry, and H.P. Hughes, Phys. Rev. B 28,
5332 (1983).
43 R. Ahuja, S. Auluck, J. Trygg, J.M. Wills, O. Eriksson,
and B. Johansson, Phys. Rev. B 51, 4813 (1995).
44 N.A.W. Holzwarth, S.G. Louie, and S. Rabii, Phys. Rev.
B 26, 5382 (1982).
45 H. Rydberg, N. Jacobson, P. Hyldgaard, S.I. Simak, B.I.
Lundqvist, and D.C. Langreth, Surf. Sci. 532-535, 606
(2003).
46 Open-source plane-wave DFT computer code Dacapo,
http://www.fysik.dtu.dk/CAMPOS/
47 D. Vanderbilt, Phys. Rev. B 41, 7892 (1990).
48 H.J. Monkhorst and J.D. Pack, Phys. Rev. B 13, 5188
(1976).
49 D.C. Langreth, private communication; J. Kleis and P.
Hyldgaard, unpublished.
50 The transition from on-surface adsorption to subsurface
absorption is identified in experiment by a work function
change, Refs. 20 and 21.
51 R. Zacharia, H. Ulbricht, and T. Hertel, Phys. Rev. B 69,
155406 (2004).
52 E. Ziambaras and E. Schröder, Phys. Rev. B 68, 064112
(2003).
53 D.E. Nixon and G.S. Parry, J. Phys. C 2, 1732 (1969).
54 O. Gunnarsson, B.I. Lundqvist, and J.W. Wilkins, Phys.
Rev. B 10, 1319 (1974). Since no spin-polarized version
of vdW-DF exists at present, we calculate the the energy
cost for changing the spin of isolated potassium atoms in
PBE. The spin-change cost is thus determined to be 26
meV/K-atom.
55 L. Bengtsson, Phys. Rev. B 59, 12301 (1999), and refer-
ences therein.
56 The choice of exchange flavor in vdW-DF was set in Ref. 15
to avoid artificial bonding in noble-gas systems and to
better mimic exact exchange calculations for those sys-
tems. However, it is far from certain and even unlikely that
the conclusions drawn for noble-gas systems carry over to
bonding separations smaller than 3 Å.
http://www.fysik.dtu.dk/CAMPOS/
