Computational Materials Science 47 (2009) 237–253

Contents lists available at ScienceDirect

Computational Materials Science

journal homepage: www.elsevier.com/locate/commatsci

Density-functional tight-binding for beginners

Pekka Koskinen *, Ville Mäkinen
NanoScience Center, Department of Physics, University of Jyväskylä, 40014 Jyväskylä, Finland

a r t i c l e i n f o a b s t r a c t

Article history: This article is a pedagogical introduction to density-functional tight-binding (DFTB) method. We derive it
Received 1 June 2009 from the density-functional theory, give the details behind the tight-binding formalism, and give practi-
Accepted 29 July 2009 cal recipes for parametrization: how to calculate pseudo-atomic orbitals and matrix elements, and espe-
Available online 22 August 2009
cially how to systematically fit the short-range repulsions. Our scope is neither to provide a historical
review nor to make performance comparisons, but to give beginner’s guide for this approximate, but
PACS: in many ways invaluable, electronic structure simulation method—now freely available as an open-
31.10.+z
source software package, hotbit.
71.15.Dx
71.15.Nc
Ó 2009 Elsevier B.V. All rights reserved.
31.15.E

Keywords:
DFTB
Tight-binding
Density-functional theory
Electronic structure simulations

1. Introduction sizes were already possible. With DFT as a competitor, the DFTB
community never grew large.
If you were given only one method for electronic structure calcu- Despite being superseded by DFT, DFTB is still useful in many
lations, which method would you choose? It certainly depends on ways: (i) In calculations of large systems [4,5]. Computational scal-
your field of research, but, on average, the usual choice is probably ing of DFT limits system sizes, while better scaling is easier to achieve
density-functional theory (DFT). It is not the best method for every- with DFTB. (ii) Accessing longer time scales. Systems that are limited
thing, but its efficiency and accuracy are suitable for most purposes. for optimization in DFT can be used for extensive studies on dynam-
You may disagree with this argument, but already DFT commu- ical properties in DFTB [6]. (iii) Structure search and general trends
nity’s size is a convincing evidence of DFT’s importance—it is even [7]. Where DFT is limited to only a few systems, DFTB can be used
among the few quantum mechanical methods used in industry. to gather statistics and trends from structural families. It can be used
Things were different before. When computational resources also for pre-screening of systems for subsequent DFT calculations
were modest and DFT functionals inaccurate, classical force fields, [8,9]. (iv) Method development. The formalism is akin to that of
semiempirical tight-binding, and jellium DFT calculations were DFT, so methodology improvements, quick to test in DFTB, can be
used. Today tight-binding is mostly familiar from solid-state text- easily exported and extended into DFT[10]. (v) Testing, playing
books as a method for modeling band-structures, with one to sev- around, learning and teaching. DFTB can be used simply for playing
eral fitted hopping parameters [1]. However, tight-binding could around, getting the feeling of the motion of atoms at a given temper-
be used better than this more often even today—especially as a ature, and looking at the chemical bonding or realistic molecular
method to calculate total energies. Particularly useful for total en- wavefunctions, even with real-time simulations in a classroom—
ergy calculations is density-functional tight-binding, which is DFTB runs easily on a laptop computer.
parametrized directly using DFT, and is hence rooted in first prin- DFTB is evidently not an ab initio method since it contains
ciples deeper than other tight-binding flavors. It just happened that parameters, even though most of them have a theoretically solid
density-functional tight-binding came into existence some ten basis. With parameters in the right place, however, computational
years ago [2,3] when atomistic DFT calculations for realistic system effort can be reduced enormously while maintaining a reasonable
accuracy. This is why DFTB compares well with full DFT with
* Corresponding author. Tel.: +358 14 260 4717.
minimal basis, for instance. Semiempirical tight-binding can be
E-mail address: pekka.koskinen@iki.fi (P. Koskinen). accurately fitted for a given test set, but transferability is usually

0927-0256/$ - see front matter Ó 2009 Elsevier B.V. All rights reserved.
doi:10.1016/j.commatsci.2009.07.013
238 P. Koskinen, V. Mäkinen / Computational Materials Science 47 (2009) 237–253

worse; for general purposes DFTB is a good choice among the dif- length (aB ¼ 0:5292 Å) and Hartree as the unit of energy (Ha ¼
ferent tight-binding flavors. 27:2114 eV). Selecting atomic mass unit (u ¼ 1:6605  1027 kg)
Despite having its origin in DFT, one has to bear in mind that DFTB the unit of mass, the unit of time becomes 1:0327 fs, appropriate
is still a tight-binding method, and should not generally be consid- for molecular dynamicspsimulations.
ffiffiffiffiffiffiffiffiffiffiffiffi Some useful fundamental
ered to have the accuracy of full DFT. Absolute transferability can constants are then h  ¼ me =u ¼ 0:0234, kB ¼ 3:1668  106 , and
never be achieved, as the fundamental starting point is tightly bound e0 ¼ 1=ð4pÞ, for instance.
electrons, with interactions ultimately treated perturbatively. DFTB Electronic eigenstates are denoted by wa , and (pseudo-atomic)
is hence ideally suited for covalent systems such as hydrocarbons basis states ul , occasionally adopting Dirac’s notation. Greek let-
[2,11]. Nevertheless, it does perform surprisingly well describing ters l, m are indices for basis states, while capital Roman letters I,
also metallic bonding with delocalized valence electrons [12,8]. J are indices for atoms; notation l 2 I stands for orbital l that be-
This article is neither a historical review of different flavors of longs to atom I. Capital R denotes nuclear positions, with the posi-
tight-binding, nor a review of even DFTB and its successes or failures. tion of atom I at RI , and displacements RIJ ¼ jRIJ j ¼ jRJ  RI j. Unit
For these purposes there are several good reviews, see, for example vectors are denoted by R b ¼ R=jRj.
Refs. [13–16]. Some ideas were around before [17–19], but the DFTB In other parts our notation remains conventional; deviations are
method in its present formulation, presented also in this article, was mentioned or made self-explanatory.
developed in the mid-1990s [2,3,20–23]. The main architects behind
the machinery were Eschrig, Seifert, Frauenheim, and their co-work- 2.2. Starting point: full DFT
ers. We apologize for omitting other contributors—consult Ref. [13]
for a more organized literature review on DFTB. The derivation of DFTB from DFT has been presented several
Instead of reviewing, we intend to present DFTB in a pedagogi- times; see, for example Refs. [18,13,3]. We do not want to be
cal fashion. By occasionally being more explicit than usual, we de- redundant, but for completeness we derive the equations briefly;
rive the approximate formalism from DFT in a systematic manner, our emphasis is on the final expressions.
with modest referencing to how the same approximations were We start from the total energy expression of interacting elec-
done before. The approach is practical: we present how to turn tron system
DFT into a working tight-binding scheme where parametrizations
E ¼ T þ Eext þ Eee þ EII ; ð1Þ
are obtained from well-defined procedures, to yield actual num-
bers—without omitting ugly details hiding behind the scenes. Only where T is the kinetic energy, Eext the external interaction (including
basic quantum mechanics and selected concepts from density- electron–ion interactions), Eee the electron–electron interaction, and
functional theory are required as pre-requisite. EII ion–ion interaction energy. Here EII contains terms like
The DFTB parametrization process is usually presented as Z vI Z vJ =jRI  RJ j, where Z vI is the valence of the atom I, and other contri-
superficially easy, while actually it is difficult, especially regarding butions from the core electrons. In density-functional theory the en-
the fitting of short-range repulsion (Section 4). In this article we ergy is a functional of the electron density nðrÞ, and for Kohn–Sham
want to present a systematic scheme for fitting the repulsion, in or- system of non-interacting electrons the energy can be written as
der to accurately document the way the parametrization is done.
We discuss the details behind tight-binding formalism, like Sla-
E½nðrÞ ¼ T s þ Eext þ EH þ Exc þ EII ; ð2Þ
ter–Koster integrals and transformations, largely unfamiliar for where T s is the non-interacting kinetic energy, EH is the Hartree en-
density-functional community; some readers may prefer to skip ergy, and Exc ¼ ðT  T s Þ þ ðEee  EH Þ is the exchange-correlation (xc)
these detailed appendices. Because one of DFTB’s strengths is the energy, hiding all the difficult many-body effects. More explicitly,
transparent electronic structure, in the end we also present se- Z !
X 3
lected analysis tools. 1 1 nðr 0 Þd r 0
E½n ¼ fa hwa j  r2 þ V ext þ jwa i þ Exc ½n þ EII ;
We concentrate on ground-state DFTB, leaving time-dependent a
2 2 jr0  rj
[24–27] or linear response [28] formalisms outside the discussion.
ð3Þ
We do not include spin in the formalism. Our philosophy lies in the
limited benefits of improving upon spin-paired self-consistent - where fa 2 ½0; 2 is the occupation of a single-particle state wa with
charge DFTB. Admittedly, one can adjust parametrizations for cer- energy ea , usually taken from the Fermi-function (with factor 2
tain systems, but the tight-binding formalism, especially the pres- for spin)
ence of the repulsive potential, contains so many approximations 1
that the next level in accuracy, in our opinion, is full DFT. fa ¼ f ðea Þ ¼ 2  ½expðea  lÞ=kB T þ 1 ð4Þ
This philosophy underlies hotbit [29] software. It is an open- P
with chemical potential l chosen such that a fa ¼ number of elec-
source DFTB package, released under the terms of GNU general trons. The Hartree potential
public license [30]. It has an interface with the atomic simulation Z 0
environment (ASE) [31], a python module for multi-purpose atom- nðr0 Þ
V H ½nðrÞ ¼ ð5Þ
istic simulations. The ASE interface enables simulations with dif- jr0  rj
ferent levels of theory, including many DFT codes or classical
is a classical electrostatic potential from given nðrÞ; for brevity we
potentials, with DFTB being the lowest-level quantum-mechanical R 3 R R 3 0 R0
will use the notation d r ! , d r ! , nðrÞ ! n, and
method. hotbit is built upon the theoretical basis described 0 0
nðr Þ ! n . With this notation the Kohn–Sham DFT energy is, once
here—but we avoid technical issues related practical implementa-
more,
tions having no scientific relevance.  
X Z
1
E½n ¼ fa hwa j  r2 þ V ext ðrÞnðrÞ jwa i
a
2
2. The origins of DFTB Z Z 0
1 nn0
þ þ Exc ½n þ EII : ð6Þ
2 jr  r 0 j
2.1. Warm-up
So far everything is exact, but now we start approximating.
We begin by commenting on practical matters. The equations Consider a system with density n0 ðrÞ that is composed of atomic
2
h =me ¼ 4pe0 ¼ e ¼ 1. This gives Bohr radius as the unit of
use  densities, as if atoms in the system were free and neutral. Hence
 P. Koskinen, V. Mäkinen / Computational Materials Science 47 (2009) 237–253 239

X
n0 ðrÞ contains (artificially) no charge transfer. The density n0 ðrÞ Erep ¼ V IJrep ðRIJ Þ: ð14Þ
does not minimize the functional E½nðrÞ, but neighbors the true I<J

minimizing density nmin ðrÞ ¼ n0 ðrÞ þ dn0 ðrÞ, where dn0 ðrÞ is sup-
For each pair of atoms IJ we have a repulsive function V IJrep ðRÞ
posed to be small. Expanding E½n at n0 ðrÞ to second order in fluc-
depending only on atomic numbers. Note that Erep contains also
tuation dnðrÞ the energy reads
on-site contributions, not only the atoms’ pair-wise interactions,
X 1 but these depend only on n0 ðrÞ and shift the total energy by a
E½dn  fa hwa j  r2 þ V ext þ V H ½n0  þ V xc ½n0 jwa i
a
2 constant.
Z Z 0 2 ! Z The pair-wise repulsive functions V IJrep ðRÞ are obtained by fitting
1 d Exc ½n0  1 1
þ þ 0
dndn0  V H ½n0 ðrÞn0 ðrÞ to high-level theoretical calculations; detailed description of the
2 dndn 0 jr  r j 2 fitting process is discussed in Section 4.
Z
þ Exc ½n0  þ EII  V xc ½n0 ðrÞn0 ðrÞ; ð7Þ
2.4. Charge fluctuation term
while linear terms in dn vanish. The first line in Eq. (7) is the band-
structure energy Let us first make a side road to recall some general concepts
X from atomic physics. Generally, the atom energy can be expressed
EBS ½dn ¼ fa hwa jH½n0 jwa i; ð8Þ as a function of Dq extra electrons as [32]
a
  !
where the Hamiltonian H0 ¼ H½n0  itself contains no charge transfer. @E 1 @2E
EðDqÞ  E0 þ Dq þ Dq 2
The second line in Eq. (7) is the energy from charge fluctuations, @ Dq 2 @ Dq 2
being mainly Coulomb interaction but containing also xc- 1
contributions ¼ E0  vDq þ U Dq2 : ð15Þ
2
Z Z !
1 0
d2 Exc ½n0  1 The (negative) slope of EðDqÞ at Dq ¼ 0 is given by the (positive)
Ecoul ½dn ¼ þ dndn0 : ð9Þ
2 dndn0 jr  r 0 j electronegativity, which is usually approximated as

The last four terms in Eq. (7) are collectively called the repulsive v  ðIE þ EAÞ=2; ð16Þ
energy where IE the ionization energy and EA the electron affinity. The (up-
Z Z ward) curvature of EðDqÞ is given by the Hubbard U, which is
1
Erep ¼ V H ½n0 ðrÞn0 ðrÞ þ Exc ½n0  þ EII  V xc ½n0 ðrÞn0 ðrÞ;
2 U  IE  EA ð17Þ
ð10Þ
and is twice the atom absolute hardness g (U ¼ 2g) [32]. Electro-
because of the ion–ion repulsion term. Using this terminology the negativity comes mainly from orbital energies relative to the vac-
energy is uum level, while curvature effects come mainly from Coulomb
E½dn ¼ EBS ½dn þ Ecoul ½dn þ Erep : ð11Þ interactions.
Let us now return from our side road. The energy in Eq. (9)
Before switching into tight-binding description, we discuss Erep and comes from Coulomb and xc-interactions due to fluctuations
Ecoul separately and introduce the main approximations. dnðrÞ, and involves a double integrals over all space. Consider the
space V divided into volumes VI related to atoms I, such that
2.3. Repulsive energy term X Z XZ
VI ¼ V and ¼ : ð18Þ
In Eq. (10) we lumped four terms together and referred them to I V I VI

as repulsive interaction. It contains the ion–ion interaction so it is

We never precisely define what these volumes VI exactly are—they
repulsive (at least at small atomic distances), but it contains also
are always used qualitatively, and the usage is case-specific. For
xc-interactions, so it is a complicated object. At this point we adopt
example, volumes can be used to calculate the extra electron popu-
manners from DFT: we sweep the most difficult physics under the
lation on atom I as
carpet. You may consider Erep as practical equivalent to an xc-func- Z
tional in DFT because it hides the cumbersome physics, while we 3
DqI  dnðrÞd r: ð19Þ
approximate it with simple functions. VI
For example, consider the total volumes in the first term, the
Hartree term By using these populations we can decompose dn into atomic
Z Z contributions
1 n0 ðrÞn0 ðr 0 Þ 3 3 0 X
 d rd r ; ð12Þ dnðrÞ ¼ DqI dnI ðrÞ; ð20Þ
2 jr  r 0 j
I
divided into atomic volumes; the integral becomes a sum over atom R 3
such that each dnI ðrÞ is normalized, VI dnI ðrÞd r ¼ 1. Note that Eqs.
pairs with terms depending on atomic numbers alone, since n0 ðrÞ
(19) and (20) are internally consistent. Ultimately, this division is
depends on them. We can hence approximate it as a sum of terms
used to convert the double integral in Eq. (9) into sum over atom
over atom pairs, where each term depends only on elements and R R
pairs IJ, and integrations over volumes VI VJ .
their distance, because n0 ðrÞ is spherically symmetric for free
First, terms with I ¼ J are
atoms. Similarly ion–ion repulsions
Z Z !
Z vI Z vJ Z vI Z vJ 1 2 0
d2 Exc ½n0  1
¼ ; ð13Þ Dq þ dnI dn0I : ð21Þ
jRJ  RI j RIJ 2 I VI VI dndn0 jr  r0 j

depend only on atomic numbers via their valence numbers Z vI . Term depends quadratically on DqI and by comparing to Eq. (15) we
Using similar reasoning for the remaining terms the repulsive en- can see that integral can be approximated by U. Hence, terms with
ergy can be approximated as I ¼ J become 12 U I Dq2I .
240 P. Koskinen, V. Mäkinen / Computational Materials Science 47 (2009) 237–253

Second, when I – J xc-contributions will vanish for local xc- and FWHMI is the full width at half maximum for the profile. This
functionals for which choice of profile is justified for a carbon atom in Fig. 1a. With these
assumptions, the Coulomb energy of two spherically symmetric
d2 Exc Gaussian charge distributions in Eq. (23) can be calculated analyti-
/ dðr  r0 Þ ð22Þ
dnðrÞdnðr0 Þ cally to yield
Z Z
and the interaction is only electrostatic, 0 dnI dn0J erfðC IJ RIJ Þ
Z Z ¼  cIJ ðRIJ Þ; ð26Þ
1 0 dn dn0 V V jr  r 0 j RIJ
I J
Dq Dq ð23Þ where
2 I J VI
0
VJ jr  r j sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
4 ln 2
between extra atomic populations DqI and DqJ . Strictly speaking, we C IJ ¼ : ð27Þ
do not know what the functions dnI ðrÞ are. However, assuming FWHM2I þ FWHM2J
spherical symmetry, they tell how the density profile of a given
In definition (26) we integrate over whole spaces because dnI ’s are
atom changes upon charging. By assuming functional form for the
strongly localized. The function (26) is plotted in Fig. 1b, where we
profiles dnI ðrÞ, the integrals can be evaluated. We choose a Gaussian
see that when R  FWHM we get point-like 1=R-interaction. Fur-
profile [33], pffiffiffiffi
thermore, as R ! 0, c ! C  2= p, which gives us a connection to
 
1 r2 on-site interactions: if I ¼ J
dnI ðrÞ ¼  3=2 exp  2 ; ð24Þ rffiffiffiffiffiffiffiffiffiffiffiffi
2pr2I 2rI
8 ln 2 1
cII ðRII ¼ 0Þ ¼ : ð28Þ
where p FWHMI
FWHM This is the on-site Coulomb energy of extra population on atom I.
rI ¼ pffiffiffiffiffiffiffiffiffiffiffiffiI ð25Þ Comparing to the I ¼ J case above, we can approximate
8 ln 2
rffiffiffiffiffiffiffiffiffiffiffiffi
8 ln 2 1 1:329
FWHMI ¼ ¼ : ð29Þ
p UI UI
We interpret Eq. (28) as: narrower atomic charge distributions
causes larger costs to add or remove electrons; for a point charge
charging energy diverges, as it should.
Hence, from absolute hardness, by assuming only Coulombic
origin, we can estimate the sizes of the charge distributions, and
these sizes can be used to estimate Coulomb interactions also be-
tween the atoms. U and FWHM are coupled by Eq. (29), and hence
for each element a single parameter U I , which can be found from
standard tables, determines all charge transfer energetics.
To conclude this subsection, the charge fluctuation interactions
can be written as
1X
Ecoul ¼ c ðRIJ ÞDqI DqJ ; ð30Þ
2 IJ IJ

where
(
UI ; I¼J
cIJ ðRIJ Þ ¼ erfðC IJ RIJ Þ ð31Þ
RIJ
; I – J:

2.5. TB formalism

So far the discussion has been without references to tight-bind-

ing description. Things like eigenstates jwa i or populations DqI can
be understood, but what are they exactly?
As mentioned above, we consider only valence electrons; the
repulsive energy contains all the core electron effects. Since
tight-binding assumes tightly bound electrons, we use minimal lo-
cal basis to expand
X
wa ðrÞ ¼ cal ul ðrÞ: ð32Þ
l

Minimality means having only one radial function for each angular
Fig. 1. (a) The change in density profile for C upon charging. Atom is slightly
charged (2=3 e and þ1 e), and we plot the averaged dnðrÞ ¼ jn
ðrÞ  n0 ðrÞj momentum state: one for s-states, three for p-states, five for
(shadowed), where n
ðrÞ is the radial electron density for charged atom, and n0 ðrÞ is d-states, and so on. We use real spherical harmonics familiar
the radial electron density for neutral atom. This is compared to the Gaussian from chemistry—for completeness they are listed in Table 2 in
profile of Eq. (24) with FWHM ¼ 1:329=U where U is given by Eq. (17). The change Appendix A.
in density near the core is irregular, but the behavior is smooth for up to  90 % of
the density change. (b) The interaction energy of two spherically symmetric
With this expansion the band-structure energy becomes
Gaussian charge distributions with equal FWHMI ¼ FWHMJ ¼ 1:329=U with U ¼ 1,
X X
a 0
EBS ¼ fa ca
l cm Hlm ; ð33Þ
as given by Eq. (26). With RIJ  FWHMI interaction is Coulomb-like, and approaches
a lm
U as RIJ ! 0.
 P. Koskinen, V. Mäkinen / Computational Materials Science 47 (2009) 237–253 241

where where
H0lm ¼ hul jH0 jum i: ð34Þ 1 1
hlm ¼ ðI þ J Þ; l 2 I m 2 J: ð42Þ
2
The tight-binding formalism is adopted by accepting the matrix ele-
ments H0lm themselves as the principal parameters of the method. This expression suggests a reasonable interpretation: charge fluctu-
This means that in tight-binding spirit the matrix elements H0lm ations shift the matrix element Hlm according to the averaged elec-
are just numbers. Calculation of these matrix elements is discussed trostatic potentials around orbitals l and m. As in Kohn–Sham
in Section 3, with details left in Appendix C. equations in DFT, also Eqs. (39) and (40) have to be solved self-con-
1
How about the atomic populations DqI ? Using the localized sistently: from a given initial guess for fDqI g one obtains hlm and
a
basis, the total number of electrons on atom I is Hlm , then by solving Eq. (39) one obtains new fcl g, and, finally,
X Z X X Z new fDqI g, iterating until self-consistency is achieved. The number
3
qI ¼ fa jwa ðrÞj2 d r ¼ fa ca a
l cm u l ðrÞum ðrÞd3 r: ð35Þ of iterations required for convergence is usually markedly less than
a VI a lm VI in DFT, albeit similar convergence problems are shared.
Atomic forces can be obtained directly by taking gradients of Eq.
If neither l nor m belong to I, the integral is roughly zero, and if both
(38) with respect to coordinates (parameters) RI . We get (with
l and m belong to I, the integral is approximately dlm since orbitals
rJ ¼ @=@RJ )
on the same atom are orthonormal. If l belongs to I and m to some
X X h i
other atom J, the integral becomes ca a 0 1
FI ¼  fa l cm rI Hlm  ðea  hlm ÞrI Slm
Z Z a lm
1 1 X
u l ðrÞum ðrÞ  u l ðrÞum ðrÞ ¼ Slm ; ð36Þ
 DqI ðrI cIJ ÞDqJ  rI Erep ; ð43Þ
VI 2 V 2
J
as suggested by Fig. 2, where Slm ¼ hul jum i is the overlap of orbitals
l and m. Charge on atom I is where the gradients of cIJ are obtained analytically from Eq. (26),
and the gradients of H0lm and Slm are obtained numerically from an
X X X 1 
qI ¼ fa ca a interpolation, as discussed in Appendix C.
l cm þ c:c: Slm ; ð37Þ
a l2I m
2
3. Matrix elements
where c.c. stands for complex conjugate. Hence DqI ¼ qI  q0I , where
q0I is the number of valence electrons for a neutral atom. This ap-
Now we discuss how to calculate the matrix elements H0lm and
proach is called the Mulliken population analysis [34].
Slm . In the main text we describe only main ideas; more detailed
Now were are ready for the final energy expression,
issues are left to Appendices Appendices A–C
X X 1X X IJ
a 0
E¼ fa ca
l cm Hlm þ c ðRIJ ÞDqI DqJ þ V rep ðRIJ Þ; ð38Þ
a lm 2 IJ IJ I<J 3.1. The pseudo-atom

where everything is, in principle, defined. We find the minimum of

P The minimal basis functions ul in the expansion (32),
this expression by variation of dðE  a ea hwa jwa iÞ, where ea are
undetermined Lagrange multipliers, constraining the wave function ul ðr0 Þ ¼ Rl ðrÞ Ye l ðh; uÞ ðr0 ¼ RI þ r; l 2 IÞ ð44Þ
norms, and obtain
X with real spherical functions Y e l ðh; uÞ as defined in Appendix A,
cam ðHlm  ea Slm Þ ¼ 0 ð39Þ should robustly represent bound electrons in a solid or molecule,
m which is what we ultimately want to simulate. Therefore orbitals
for all a and l. This equation for the coefficients cal is the Kohn– should not come from free atoms, as they would be too diffuse.
Sham equation-equivalent in DFTB. Here To this end, we use the orbitals from a pseudo-atom, where an addi-
X tional confinement potential V conf ðrÞ is added to the Hamiltonian
1
Hlm ¼ H0lm þ Slm ðcIK þ cJK ÞDqK ; l 2 I m 2 J: ð40Þ 1 Z
2 K  r2  þ V H ðrÞ þ V xc ðrÞ þ V conf ðrÞ: ð45Þ
2 r
By noting that the electrostatic potential on atom I due to charge
P This additional, spherically symmetric confinement cuts the orbi-
fluctuations is I ¼ K cIK DqK , the equation above can be written as
tals’ diffuse tails off and makes a compact basis—and ultimately bet-
1
Hlm ¼ H0lm þ hlm Slm ; ð41Þ ter basis[17]—for the wave function expansion.
A general, spherically symmetric environment can be repre-
sented by a potential
X
1
V conf ðrÞ ¼ v2i r2i ; ð46Þ
i¼0

where the odd terms disappear because the potential has to be

smooth at r ¼ 0. Since the first v0 term is just a constant shift, the
first non-trivial term is v2 r 2 . To first approximation we hence
choose the confining potential to be of the form
 2
r
V conf ðrÞ ¼ ; ð47Þ
r0
where r 0 is a parameter. The quadratic form for the confinement has
Fig. 2. Integrating the overlap of local orbitals. The large shaded area denotes the appeared before [2], but also other forms have been analyzed [35].
volume VI of atom I, spheres represent schematically the spatial extent of orbitals,
and the small hatched area denotes the overlap region. With m 2 J and l 2 I, the
While different forms can be considered for practical reasons, they
integration of ul ðrÞ um ðrÞ over atom I’s volume VI misses half of the overlap, since have only little effect on DFTB performance. The adjustment of the
the other half is left approximately to VJ . parameter r0 is discussed in Section 5.
242 P. Koskinen, V. Mäkinen / Computational Materials Science 47 (2009) 237–253

The pseudo-atom is calculated with DFT only once for a given is the effective potential evaluated at the (artificial) neutral density
confining potential. This way we get ul ’s (more precisely, Rl ’s), n0 ðrÞ of the system. The density n0 is determined by the atoms in
the localized basis functions, for later use in matrix element the system, and the above matrix element between basis states l
calculations. and m, in principle, depends on the positions of all atoms. However,
One technical detail we want to point out here concerns orbital since the integrand is a product of factors with three centers, two
conventions. Namely, once the orbitals ul are calculated, their sign wave functions and one potential (and kinetic), all of which are
and other conventions should never change. Fixed convention non-zero in small spatial regions only, reasonable approximations
should be used in all simultaneously used Slater–Koster tables; can be made.
using different conventions for same elements gives inconsistent First, for diagonal elements Hll one can make a one-center
tables that are plain nonsense. The details of our conventions, approximation where the effective potential within volume VI is
along with other technical details of the pseudo-atom calculations,
V s ½n0 ðrÞ  V s;I ½n0;I ðrÞ; ð51Þ
are discussed in Appendix A.
where l 2 I. This integral is approximately equal to the eigenener-
3.2. Overlap matrix elements gies el of free atom orbitals. This is only approximately correct since
the orbitals ul are from the confined atom, but is a reasonable
Using the orbitals from pseudo-atom calculations, we need to approximation that ensures the correct limit for free atoms.
calculate the overlap matrix elements Second, for off-diagonal elements we make the two-center
Z approximation: if l is localized around atom I and m is localized
Slm ¼ ul ðrÞ um ðrÞd3 r: ð48Þ around atom J, the integrand is large when the potential is local-
ized either around I or J as well; we assume that the crystal field
contribution from other atoms, when the integrand has three dif-
Since orbitals are chosen real, the overlap matrix is real and ferent localized centers, is small. Using this approximation the
symmetric. effective potential within volume VI þ VJ becomes
The integral with ul at RI and um at RJ can be calculated also
with ul at the origin and um at RIJ . Overlap will hence depend on V s ½n0 ðrÞ  V s;I ½n0;I ðrÞ þ V s;J ½n0;J ðrÞ; ð52Þ
RIJ , or equivalently, on RIJ and R b IJ separately. Fortunately, the
b IJ is fully governed by Slater–Koster transformation where V s;I ½n0;I ðrÞ is the Kohn-Sham potential with the density of a
dependence on R
neutral atom. The Hamiltonian matrix element is
rules [36]. Only one to three Slater–Koster integrals, depending on
Z  
the angular momenta of ul and um , are needed to calculate the 1
b IJ for fixed RIJ . These rules originate from the H0lm ¼ ul ðrÞ  r2 þ V s;I ðrÞ þ V s;J ðrÞ um ðrÞ; ð53Þ
integral with any R 2
properties of spherical harmonics.
The procedure is hence the following: we integrate numerically where l 2 I and m 2 J. Prior to calculating the integral, we have to
the required Slater–Koster integrals for a set of RIJ , and store them apply the Hamiltonian to um . But in other respects the calculation
in a table. This is done once for all orbital pairs. Then, for a given is similar to overlap matrix elements: Slater–Koster transforma-
orbital pair, we interpolate this table for RIJ , and use the Slater– tions apply, and only a few integrals have to be calculated numeri-
Koster rules to get the overlap with any geometry—fast and cally for each pair of orbitals, and stored in tables for future
accurately. reference. See Appendix C.2 for details of numerical integration of
Readers unfamiliar with the Slater–Koster transformations can the Hamiltonian matrix elements.
read the detailed discussion in Appendix B. The numerical integra-
tion of the integrals is discussed in Appendix C. 4. Fitting the repulsive potential
Before concluding this subsection, we make few remarks about
non-orthogonality. In DFTB it originates naturally and inevitably In this section we present a systematic approach to fit the repul-
from Eq. (48), because non-overlapping orbitals with diagonal sive functions V IJrep ðRIJ Þ that appear in Eq. (38)—and a systematic
overlap matrix would yield also diagonal Hamiltonian matrix, way to describe the fitting. But first we discuss some difficulties re-
which would mean chemically non-interacting system. The lated to the fitting process.
transferability of a tight-binding model is often attributed to The first and straightforward way of fitting is simple: calculate
non-orthogonality, because it accounts for the spatial nature of dimer curve EDFT ðRÞ for the element pair with DFT, require
the orbitals more realistically. EDFT ðRÞ ¼ EDFTB ðRÞ, and solve
Non-orthogonality requires solving a generalized eigenvalue
V rep ðRÞ ¼ EDFT ðRÞ  ½EBS ðRÞ þ Ecoul ðRÞ: ð54Þ
problem, which is more demanding than normal eigenvalue prob-
lem. Non-orthogonality complicates, for instance, also gauge trans- We could use also other symmetric structures with N bonds having
formations, because the phase from the transformations is not well equal RIJ ¼ R, and require
defined for the orbitals due to overlap. The Peierls substitution
[37], while gauge invariant in orthogonal tight-binding [38,39], is N  V rep ðRÞ ¼ EDFT ðRÞ  ½EBS ðRÞ þ Ecoul ðRÞ: ð55Þ
not gauge invariant in non-orthogonal tight-binding (but affects In practice, unfortunately, it does not work out. The approximations
only time-dependent formulation). made in DFTB are too crude, and hence a single system is insuffi-
cient to provide a robust repulsion. As a result, fitting repulsive
3.3. Hamiltonian matrix elements potentials is difficult task, and forms the most laborous part of
parametrizing in DFTB.
From Eq. (34) the Hamiltonian matrix elements are Let us compare things with DFT. As mentioned earlier, we tried
Z  
1 to dump most of the difficult physics into the repulsive potential,
H0lm ¼ ul ðrÞ  r2 þ V s ½n0 ðrÞ um ðrÞ; ð49Þ and hence V rep in DFTB has practical similarity to Exc in DFT. In
2
DFTB, however, we have to make a new repulsion for each pair
where of atoms, so the testing and fitting labor compared to DFT function-
als is multifold. Because xc-functionals in DFT are well docu-
V s ½n0 ðrÞ ¼ V ext ðrÞ þ V H ½n0 ðrÞ þ V xc ½n0 ðrÞ ð50Þ mented, DFT calculations of a reasonably documented article can
 P. Koskinen, V. Mäkinen / Computational Materials Science 47 (2009) 237–253 243


be reproduced, whereas reproducing DFTB calculations is usually E0 ðRAB Þ  E0BS ðRAB Þ þ E0coul ðRAB Þ
V 0rep ðRAB Þ ¼ DFT ; ð57Þ
harder. Even if the repulsive functions are published, it would be N
a great advantage to be able to precisely describe the fitting pro-
where the prime stands for a derivative with respect to RAB . The eas-
cess; repulsions could be more easily improved upon.
iest way is first to calculate the energy curve and use finite differ-
Our starting point is a set of DFT structures, with geometries R,
ences for derivatives. In fact, systems treated this way can have
energies EDFT ðRÞ, and forces F DFT (zero for optimized structures). A
even different RAB ’s if only the ones that are equal are chosen to vary
natural approach would be to fit V rep so that energies EDFTB ðRÞ and
(e.g. a complex system with one appropriate AB bond on its surface).
forces F DFTB are as close to DFT ones as
possible. In other words, we
For each system, this gives a family of data points for the fitting; the
want to minimize force differences F DFT  F DFTB and energy dif-
number of points in the family does not affect fitting, as explained
ferences jEDFT  EDFTB j on average for the set of structures. There
later. The dimer curve, with N ¼ 1, is clearly one system where this
are also other properties such as basis set quality (large overlap
method can be applied. For any equilibrium DFT structure things
with DFT and DFTB wave functions), energy spectrum (similarity
simplify into
of DFT and DFTB density of states), or charge transfer to be com-
pared with DFT, but these originate already from the electronic
   
  E0BS R0AB  E0coul R0AB
part, and should be modified by adjusting V conf ðrÞ and Hubbard V 0rep R0AB ¼ ; ð58Þ
U. Repulsion fitting is always the last step in the parametrizing, N
 
and affects only energies and forces. where R0AB is the distance for which E0DFT R0AB ¼ 0.
In practice we shall minimize, however, only force differ-
ences—we fit repulsion derivative, not repulsion directly. The fit-
4.1.2. Homonuclear systems
ting parameters, introduced shortly, can be adjusted to get energy
differences qualitatively right, but only forces are used in the If a cluster or a solid has different bond lengths, the energy
curve method above cannot be applied (unless a subset of bonds
practical fitting algorithm. There are several reasons for this. First,
forces are absolute, energies only relative. For instance, since we are selected). But if the system is homonuclear, the data points
can be obtained the following way. The force on atom I is
do not consider spin, it is ambiguous whether to fit to DFT dimer
curve with spin-polarized or spin-paired free atom energies. We X
F I ¼ F 0I þ b IJ
V 0rep ðRIJ Þ R ð59Þ
could think that lower-level spin-paired DFT is the best DFTB J–I
can do, so we compare to spin-paired dimer curve—but we should X
fit to energetics of nature, not energetics in some flavors of DFT. ¼ F 0I þ IJ Rb IJ ; ð60Þ
J–I
Second, for faithful dynamics it is necessary to have right forces
and right geometries of local energy minima; it is more important where F 0I is the force without repulsions. Then we minimize the
first to get local properties right, and afterwards look how the sum
global properties, such as energy ordering of different structural X
motifs, come out. Third, the energy in DFTB comes mostly from F DFT;I  F I 2 ð61Þ
the band-structure part, not repulsion. This means that if already I

the band-structure part describes energy wrong, the short-ranged

with respect to IJ , with IJ ¼ 0 for pair distances larger than the cut-
repulsions cannot make things right. For instance, if EDFT ðRÞ and
off. The minimization gives optimum IJ , which can be used directly,
EDFTB ðRÞ for dimer deviates already with large R, short-range
together with their RIJ ’s, as another family of data points in the
repulsion cannot cure the energetics anymore. For transferability
fitting.
repulsion has to be monotonic and smooth, and if repulsion is ad-
justed too rapidly catch up with DFT energetics, the forces will go
wrong. 4.1.3. Other algorithms
For the set of DFT structures, we will hence minimize DFT and Fitting algorithms like the ones above are easy to construct, but
DFTB force differences, using the recipes below. a few general guidelines are good to keep in mind.
While pseudo-atomic orbitals are calculated with LDA-DFT, the
systems to fit the repulsive potential should be state-of-the-art cal-
4.1. Collecting data
culations; all structural tendencies—whether right or wrong—are
directly inherited by DFTB. Even reliable experimental structures
To fit the derivative of the repulsion for element pair AB, we
can be used as fitting structures; there is no need to think DFTB
need a set of data points fRi ; V 0rep ðRi Þg. As mentioned before, fitting
should be parametrized only from theory. DFTB will not become
to dimer curve alone does not give a robust repulsion, because the
any less density-functional by doing so.
same curve is supposed to work in different chemical environ-
As data points are calculated by stretching selected bonds (or
ments. Therefore it is necessary to collect the data points from sev-
calculating static forces), also other bonds may stretch (dimer is
eral structures, to get a representative average over different types
one exception). These other bonds should be large enough to ex-
of chemical bonds. Here we present examples on how to acquire
clude repulsive interactions; otherwise fitting a repulsion between
data points.
two elements may depend on repulsion between some other ele-
ment pairs. While this is not illegal, the fitting process easily be-
4.1.1. Force curves and equilibrium systems comes complicated. Sometimes the stretching can affect chemical
This method can be applied to any system where all the bond interactions between elements not involved in the fitting; this is
lengths between the elements equal RAB or otherwise are beyond worth avoiding, but sometimes it may be inevitable.
the selected cutoff radius Rcut . In other words, the only energy com-
ponent missing from these systems is the repulsion from N bonds
4.2. Fitting the repulsive potential
between elements A and B with matching bond lengths. Hence,

EDFTB ðRAB Þ ¼ EBS ðRAB Þ þ Ecoul ðRAB Þ þ e

E rep þ N  V rep ðRAB Þ; ð56Þ Transferability requires the repulsion to be short-ranged, and
we choose a cutoff radius Rcut for which V rep ðRcut Þ ¼ 0, and also
where eE rep is the repulsive energy independent of RAB . This setup V 0rep ðRcut Þ ¼ 0 for continuous forces. Rcut is one of the main
allows us to change RAB , and we will require parameters in the fitting process. Then, with given Rcut , after having
244 P. Koskinen, V. Mäkinen / Computational Materials Science 47 (2009) 237–253

collected enough data points fRi ; V 0rep;i g, we can fit the function to be independent of the number of points in each family; a fit with
V 0rep ðRÞ. The repulsion itself is dimer force curve should yield the same result with 10 or 100
Z points in the curve. Hence for each system
Rcut
V rep ðRÞ ¼  V 0rep ðrÞdr: ð62Þ pffiffiffiffiffiffi
R ri ¼ rs Ns ; ð64Þ
Fitting of V 0rep
using the recipe below provides a robust and
where rs is the uncertainty given for system s, with N s points in the
unbiased fit to the given set of points, and the process is easy to
family. This means that systems with the same rs ’s have the same
control. We choose a standard smoothing spline [40] for
significance in the process, irregardless of the number of data points
V 0rep ðRÞ  UðRÞ, i.e. we minimize the functional
in each system. The effect of Eq. (64) is the same as putting a weight
!2 Z 1=Ns for each data point in the family. Note that k has nothing to do
X
M
V 0rep;i  UðRi Þ Rcut
with the number of data points, and is more universal parameter.
S½UðRÞ ¼ þk U 00 ðRÞ2 dR ð63Þ
i¼1
ri The cutoff is set by adding a data point at UðRcut Þ ¼ 0 with a tiny r.
Fig. 3 shows an example of fitting carbon–hydrogen repulsion.
for total M data points fRi ; V 0rep;i g, where UðRÞ is given by a cubic The parameters Rcut and k, as well as parameters ri , are in practice
spline. Spline gives an unbiased representation for UðRÞ, and the chosen to yield visually satisfying fit; the way of fitting should not
smoothness can be directly controlled by the parameter k. Large k affect the final result, and in this sense it is just a technical neces-
means expensive curvature and results in linear UðRÞ (quadratic sity—the simple objective is to get a smooth curve going nicely
V rep ) going through the data points only approximately, while small through the data points. Visualization of the data points can be also
k considers curvature cheap and may result in a wiggled UðRÞ pass- generally invaluable: deviation from a smooth behavior tells how
ing through the data points exactly. The parameter k is the second well you should expect DFTB to perform. If data points lay nicely
parameter in the fitting process. Other choices for UðRÞ can be used, along one curve, DFTB performs probably well, but scattered data
such as low-order polynomials [2], but they sometimes behave sur- points suggest, for instance, the need for improvements in the elec-
prisingly while continuously tuning Rcut . For transferability the tronic part. In the next section we discuss parameter adjustment
behavior of the derivative should be as smooth as possible, prefer- further.
ably also monotonous (the example in Fig. 3a is slightly non-monot-
onous and should be improved upon).
The parameters ri are the data point uncertainties, and can be 5. Adjusting parameters
used to weight systems differently. With the dimension of force,
ri ’s have also an intuitive meaning as force uncertainties, the In this section we summarize the parameters, give practical
lengths of force error bars. As described above, each system may instructions for their adjustment, and give a demonstration of their
produce a family of data points. We would like, however, the fitting usage. The purpose is to give an overall picture of the selection of
knobs to turn while adjusting parametrization.
Occasionally one finds published comments about the perfor-
mance of DFTB. While DFTB shares flaws and failures characteristic
to the method itself, it should be noted that DFTB parametrizations
are even more diverse than DFT functionals. A website in Ref. [41],
maintained by the original developers of the method, contains var-
ious sets of parametrizations. While these parametrizations are of
good quality, they are not unique. Namely, there exists no auto-
mated way of parametrizing, so that a straightforward process
would give all parameters definite values. This is not necessarily
a bad thing, since some handwork in parametrizing also gives feel-
ing what is to be expected in the future, how well parametrizations
are expected to perform.

5.1. Pseudo-atoms

The basis functions, and, consequently, the matrix elements are

determined by the confinement potentials V conf ðrÞ containing the
parameters r 0 for each element. The value r0 ¼ 2  r cov , where r cov
is the covalent radius, can be used as a rule of thumb [2]. Since
the covalent radius is a measure for binding range, it is plausible
that the range for environmental confining potential should de-
pend on this scale—the number 2 in this rule is empirical.
With this rule of thumb as a starting point, the quality of the ba-
sis functions can be inspected and r0 adjusted by looking at (i)
band-structure (for solids), (ii) densities of states (dimer or other
simple molecules), or (iii) amount of data point scatter in repulsion
fit (see Section 4.2, especially Fig. 3). Systems with charge transfer
should be avoided herein, since the properties listed above would
depend on electrostatics as well, which complicates the process;
adjustment of r0 should be independent of electrostatics.
Fig. 3. (a) Fitting the derivative of repulsive potential. Families of points from
various structures, obtained by stretching C–H bonds and using Eq. (57). Here
The inspections above are easiest to make with homonuclear
Rcut ¼ 1:8 Å, for other details see Section 5. (b) The repulsive potential V rep ðRÞ, which interactions, even though heteronuclear interactions are more
is obtained from by integration of the curve in (a). important for some elements, such as for hydrogen. Different
 P. Koskinen, V. Mäkinen / Computational Materials Science 47 (2009) 237–253 245

chemical environments can affect the optimum value of r 0 , but this is why the parameter k is not given any guidelines here. Large
usually it is fixed for all interactions of a given element. r’s can be used to give less weight for systems with (i) marginal
importance for systems of interest, (ii) inter-dependence on other
5.2. Electrostatics parameters (dependence on other repulsions, on electrostatics,
or on other chemical interactions), (iii) statistically peculiar stick-
Electrostatic energetics, as described in Section 2.4, are deter- ing out from the other systems (reflecting situation that cannot
mined by the Hubbard U parameter, having the default value be described by tight-binding or the urge to improve electronic
U ¼ IE  EA. Since U is a value for a free atom, while atoms in mol- part).
ecules are not free, it is permissible to adjust U in order to improve
(i) charge transfer, (ii) density of states, (iii) molecular ionization 5.4. Hydrocarbon parametrization
energies and electron affinities, or (iv) excitation spectra for se-
lected systems. Since U is an atomic property, one should beware To demonstrate the usage of the parameters, we present the
when using several different elements in fitting—charge transfer hydrocarbon parametrization used in this article (this was first
depends on U’s of other elements as well. shot fitting without extensive adjustment, but works reasonably
Eq. (29) relates U and FWHM of given atom together. But since well). The Hubbard U is given by Eq. (17) and FWHM by Eq. (29),
Eq. (21) contains, in principle, also xc-contributions, the relation both for hydrogen and carbon. The force curves have been calcu-
can be relaxed, if necessary. Since FWHM affects only pair-interac- lated with GPAW [42,43] using the PBE xc-functional [44].
tions, it is better to adjust the interactions directly like Hydrogen: r 0 ¼ 1:08, U ¼ 0:395. Carbon: r 0 ¼ 2:67, U ¼ 0:376.
Hydrogen–carbon repulsion: Rcut ¼ 3:40, k ¼ 35, and systems with
C IJ ! C IJ =xIJ ; ð65Þ
force curve: CH , ethyne, methane, benzene; all with ri ¼ 1. Car-
where xIJ (being close to one) effectively scales both FWHMI and bon–carbon repulsion: Rcut ¼ 3:80, k ¼ 200, and systems with
FWHMJ . If atomic FWHMs would be changed directly, it would af- force curve: C2 , CC2 , C3 , C2 2
4 with ri ¼ 1, and C3 with ri ¼ 0:3.
fect all interactions and complicate the adjusting process. Note that
FWHM affects only nearest neighbor interactions (see Fig. 1) and al- 6. External fields
ready next-nearest neighbors have (very closely) the pure 1=R
interaction. Including external potentials to the formalism is straightfor-
An important principle, general for all parameters but particu- ward: just add one more term in the Hamiltonian of Eq. (53). Ma-
larly for electrostatics, is this: all parameters should be adjusted trix element with external (scalar) potential becomes, with
within reasonable limits. This means that, since all parameters plausible approximations,
have a physical meaning, if a parameter is adjusted beyond a rea-
Z
sonable and physically motivated limit, the parametrization will in
general not be transferable. If a good fit should require overly large
Hlm ! Hlm þ u l ðrÞV ext ðrÞum ðrÞd3 r
parameter adjustments (precise ranges are hard to give), the origi-
Z Z
nal problem probably lies in the foundations of tight-binding.  Hlm þ V ext ðRI Þ u l um þ V ext ðRJ Þ u l um
VI VJ

1 I 
5.3. Repulsive potentials  H lm þ V þ V Jext Slm ð66Þ
2 ext
The last step in the parametrization is to fit the repulsive poten-
for a smoothly varying V ext ðrÞ. The electrostatic part in the Hamilto-
tial. Any post-adjustment of other parameters calls for re-fitting of
nian is
the repulsive potential.
The most decisive part in the fitting is choosing the set of struc- 1 1 
hlm ¼ I þ J þ V Iext þ V Jext ð67Þ
tures and bonds to fit. Parameters Rcut , r and k are necessary, but 2
they play only a limited part in the quality of the fit—the quality and naturally extends Eq. (42).
and transferability is determined by band structure energy, elec-
trostatic energy, and the chosen structures. In fact, the repulsion
7. van der Waals forces
fitting was designed such that the user has as little space for adjust-
ment as possible.
Accurate DFT xc-functionals, which automatically yield the R6
The set of structures should contain the fitted interaction in dif-
long range attractive van der Waals interactions, are notoriously
ferent circumstances, with (i) different bond lengths, (ii) varying
hard to make [45]. Since DFT in other respects is accurate with
coordination, and (iii) varying charge transfer. In particular, if
short-range interactions, it would be wrong to add van der Waals
charge transfer is important for the systems of interest, calculation
interactions by hand—addition inevitably modifies short-range
of charged molecules is recommended.
parts as well.
A reasonable initial guess for the cutoff radius is Rcut ¼ 1:5
DFTB, on the other hand, is more approximate, and adding
Rdimer , being half-way between nearest and next-nearest neighbors
physically motivated terms by hand is easier. In fact, van der Waals
for homonuclear systems. It is then adjusted to yield a satisfying fit
forces in DFTB can conceptually be thought of as modifications of
for the derivative of the repulsion, as in Fig. 3, while remembering
the repulsive potential. Since dispersion forces are due to xc-con-
that it has to be short ranged (Rcut ¼ 2  Rdimer , for instance, is too
tributions, one can see that for neutral systems, where dnðrÞ  0,
large, lacks physical motivation, and makes fitting hard). Increasing
dispersion has to come from Eq. (10). However, in practice it is bet-
Rcut will increase V rep ðRÞ at given R < Rcut , which is an aspect that
ter to leave V IJrep ’s short-ranged and add the dispersive forces as
can be used to adjust energies (but not much forces). Usually the
additional terms
parameters r are used in the sense of relative weights between
systems, as often they cannot be determined in the sense of abso- X C IJ6
lute force uncertainties. The absolute values do not even matter,
EvdW ¼  fIJ ðRIJ Þ ð68Þ
I<J R6IJ
since the scale of r’s merely sets the scale for the smoothness
parameter k (you can start with r ¼ 1 for the first system; if you in the total energy expression. Here f ðRÞ is a damping function with
multiply r’s by x, the same fit is obtained with k multiplied by x2 )— the properties
246 P. Koskinen, V. Mäkinen / Computational Materials Science 47 (2009) 237–253


 1; R J R0 for each k-point from a chosen set, such as Monkhorst–Pack sam-
f ðRÞ ¼ ð69Þ
 0; R K R0 ; pled [50]. Above we have

because the idea is to switch off van der Waals interactions for dis- 1 1
hlm ¼ ðI þ J Þ l 2 I; m 2 J; ð78Þ
tances smaller than R0 , a characteristic distance where chemical 2
interactions begin to emerge.
as in Eq. (42), and Mulliken charges are extensions of Eq. (37),
The C 6 -parameters depend mainly on atomic polarizabilities
and have nothing to do with DFTB formalism. Care is required to XX X 1h i
qI ¼ fa ðkÞ ca a
l ðkÞcm ðkÞSlm ðkÞ þ c:c: : ð79Þ
avoid large repulsive forces, coming from abrupt behavior in f ðRÞ 2
a k l2I;m
near R  R0 , which could result in local energy minima. For a de-
tailed descriptions about the C 6 -parameters and the form of f ðRÞ The sum for the electrostatic energy per unit cell,
we refer to original Refs. [46,47]; in this section we merely demon-
strate how straightforward it is, in principle, to include van der X
1 unit cell X

Waals forces in DFTB.

Ecoul ¼ cIJ ðRIJ  TÞDqI DqJ ; ð80Þ
2 IJ T

8. Periodic boundary conditions can be calculated with standard methods, such as Ewald summation
[51], and the repulsive part,
8.1. Bravais lattices
X X
1 unit cell X
V IJrep ðRIJ Þ ¼ V IJrep ðRIJ  TÞ ð81Þ
Calculation of isolated molecules with DFTB is straightforward, I<J
2 IJ T
but implementation of periodic boundary conditions and calcula-
tion of electronic band-structures is also easy[48]. As mentioned is easy because repulsions are short-ranged (and V rep ð0Þ ¼ 0 is
in the introduction, this is usually the first encounter with tight- understood).
binding models for most physicists; our choice was to discuss peri-
odic systems at later stage. 8.2. General symmetries
In a crystal periodic in translations T, the wave functions have
the Bloch form Thanks to the transparent formalism of DFTB, it is easy to con-
ikr struct more flexible boundary conditions, such as the ‘‘wedge
wa ðk; rÞ ¼ e ua ðk; rÞ; ð70Þ
boundary condition” introduced in Ref. [52]. This is one example
where ua ðk; rÞ is function with crystal periodicity [49]. This means of DFTB in method development.
that a wave function wa ðk; rÞ changes by a phase eikT in translation General triclinic unit cells are copied by translations, and DFT
T. We define new basis functions, not as localized orbitals anymore, implementation is easy with plane waves, real-space grids or local-
but as Bloch waves extended throughout the whole crystal ized orbitals with fixed quantization axis. But if we allow the quan-
tization axis of localized orbitals to be position-dependent, we can
1 X
ul ðk; rÞ ¼ pffiffiffiffi eikT ul ðr  TÞ; ð71Þ treat more general symmetries which have rotational symmetry
N T [53] or even combined rotational and translational (chiral) symme-
where N is the (infinite) number of unit cells in the crystal. The tries [54].
eigenfunction ansatz The basic idea is to enforce the orbitals to have the the same
X symmetry as the system. This requires that basis functions not only
wa ðk; rÞ ¼ cal ðkÞul ðk; rÞ ð72Þ depend on atom positions like
l
ul ðr  Rl Þ; ð82Þ
is then also an extended Bloch wave, as required by Bloch theorem,
because k is the same for all basis states. Matrix elements in this as usual, but more generally like
new basis are
DðRl Þul ðrÞ; ð83Þ
0 0
Slm ðk; k Þ ¼ dðk  k ÞSlm ðkÞ ð73Þ
where DðRl Þ is an operator transforming the orbitals in any posi-
and tion-dependent manner, including both translations and rotations.
The only requirement is that the orbitals are complete and ortho-
0 0
Hlm ðk; k Þ ¼ dðk  k ÞHlm ðkÞ; ð74Þ normal for a given angular momentum. If the quantization axes
change, things become unfortunately messy. However, suitably de-
where fined basis orbitals yield well-defined Hamiltonian and overlap
X Z  X matrices, and enable simulations of systems like bent tubes or slabs,
Slm ðkÞ ¼ eikT u l ðrÞum ðr  TÞ  eikT Slm ðTÞ ð75Þ helical structures such as DNA, or a piece of spherical surface—with
T T
a greatly reduced number of atoms. Similar concepts are familiar
and similarly for H. Obviously the Hamiltonian conserves k—that is from chemistry, where symmetry-adapted molecular orbitals are
why k labels the eigenstates in the first place. Note that the new ba- constructed from the atomic orbitals, and computational effort is
sis functions are usually not normalized. hereby reduced [55]. Detailed treatment of these general symme-
Inserting the trial wave function (72) into Eq. (38) and by using tries is a work in progress [56].
the variational principle we obtain the secular equation
X

cam ðkÞ Hlm ðkÞ  ea ðkÞSlm ðkÞ ¼ 0; ð76Þ 9. Density-matrix formulation
m
In this section we introduce DFTB using density-matrix formu-
where
lation. We do this because not only does the formulation simplify
Hlm ðkÞ ¼ H0lm ðkÞ þ hlm Slm ðkÞ
1
ð77Þ expressions, but it also makes calculations faster in practice. This
practical advantage comes because the density-matrix,
 P. Koskinen, V. Mäkinen / Computational Materials Science 47 (2009) 237–253 247

X X
qlm ¼ fa ðcal ca
m Þ¼ fa ðqalm Þ; ð84Þ is easy to partition into smaller pieces. Population of a single orbital
a a l is
X
contains a loop over eigenstates; quantities calculated with qlm sim- qðlÞ ¼ q~ lm Sml ¼ ðq~ SÞll ; ð98Þ
ply avoid this extra loop. It has the properties m

qSq ¼ q ð idempotencyÞ; ð85Þ while population on atom I due to eigenstate wa alone is

XX
qlm ¼ q ml ðq ¼ qy Þ; ð86Þ ~ a SÞ ¼
qI;a ¼ TrI ðq q~ alm Sml ; ð99Þ
Trðqa SÞ ¼ 1 ðeigenfunction normalizationÞ: ð87Þ l2I m

We define also the energy-weighted density-matrix so that

X !
e a a
qlm ¼ ea fa cl cm ; ð88Þ X X X
a fa qI;a ¼ ðqI Þ ¼ Nel : ð100Þ
I a I
and symmetrized density matrix
1 Population on orbitals of atom I with angular momentum l is,
q~ lm ¼ ðqlm þ q lm Þ; ð89Þ similarly,
2
which is symmetric and real. Using qlm we obtain simple expres- X
qlI ¼ ðq
~ SÞll : ð101Þ
sions, for example, for l2Iðll ¼lÞ
EBS ¼ TrðqH0 Þ ð90Þ
XX The partial Mulliken populations introduced above are simple, but
Nel ¼ Trðq
~ SÞ ¼ q~ lm Sml ; ð91Þ enable surprisingly rich analysis of the electronic structure, as dem-
l m
XX onstrated below.
qI ¼ TrI ðq
~ SÞ ¼ q~ lm Sml ; ð92Þ
l2I m 10.2. Analysis beyond Mulliken charges
where EBS is the band-structure energy, Nel is the total number of
electrons, and qI is the Mulliken population on atom I. TrI is partial At this point, after discussing Mulliken population analysis, we
trace over orbitals of atom I alone. comment on the role of wave functions in DFTB. Namely, internally
It is practical to define also matrices’ gradients. They do not di- DFTB formalism uses atom resolution for any quantity, and the
rectly relate to density matrix formulation, but equally simplify tight-binding spirit means that the matrix elements Hlm and Slm
notation, and are useful in practical implementations. We define are just parameters, nothing more. Nonetheless, the elements Hlm
(with rJ ¼ @=@RJ ) and Slm are obtained from genuine basis orbitals ul ðrÞ using
well-defined procedure—these basis orbitals remain constantly
dS lm ¼ rJ Slm ðl 2 I; m 2 JÞ available for deeper analysis. The wave functions are
Z
X
 u l ðr  RI ÞrJ um ðr  RJ Þ ¼ hul jrJ um i ð93Þ wa ðrÞ ¼ cal ul ðrÞ ð102Þ
l
and
and the total electron density is
dH lm ¼ rJ Hlm hul jHjrJ um i ðl 2 I; m 2 JÞ ð94Þ X X
nðrÞ ¼ fa jwðrÞj2 ¼ qlm u m ðrÞul ðrÞ; ð103Þ

with the properties dS lm ¼ dS ml and dH lm ¼ dH ml . From these a lm
definitions we can calculate analytically, for instance, the time
derivative of the overlap matrix for a system in motion awaiting for inspection with tools familiar from DFT. One should,
however, use the wave functions only for analysis[58]. The formal-
S_ ¼ ½dS; V ð95Þ ism itself is better off with Mulliken charges. But for visualization
with commutator ½A; B and matrix V lm ¼ dlm R_ I , l 2 I. Force from the and for gaining understanding this is a useful possibility. This dis-
band-energy part, the first line in Eq. (43), for atom I can be ex- tinguishes DFTB from semiempirical methods, which—in princi-
pressed as ple—do not possess wave functions but only matrix elements
(unless made ad hoc by hand).
F I ¼ TrI ðqdH  qe dSÞ þ c:c:; ð96Þ
which is, besides compact, useful in implementation. The density- 10.3. Densities of states
matrix formulation introduced here is particularly useful in elec-
tronic structure analysis, discussed in the following section. Mulliken populations provide intuitive tools to inspect elec-
tronic structure. Let us first break down the energy spectrum into
10. Simplistic electronic structure analysis various components. The complete energy spectrum is given by the
density of states (DOS),
One great benefit of tight-binding is the ease in analyzing the X r
DOSðeÞ ¼ d ðe  ea Þ; ð104Þ
electronic structure. In this section we present selected analysis a
tools, some old and renowned, others casual but intuitive. Other
simple tools for chemical analysis of bonding can be found from where dr ðeÞ can be either the peaked Dirac delta-function, or some
Ref. [57]. function—such as a Gaussian or a Lorentzian—with broadening
parameter r. DOS carrying spatial information is the local density
10.1. Partial Mulliken populations of states,
X r
The Mulliken population on atom I, LDOSðe; rÞ ¼ d ðe  ea Þjwa ðrÞj2 ð105Þ
a
XX X R
qI ¼ TrI ðq
~ SÞ ¼ q~ lm Sml ¼ ðq
~ SÞll ð97Þ with integration over
3
d r yielding DOSðeÞ. Sometimes it is defined
l2I m l2I as
248 P. Koskinen, V. Mäkinen / Computational Materials Science 47 (2009) 237–253

X
LDOSðrÞ ¼ fa0 jwa ðrÞj2 ; ð106Þ Let us start by discussing promotion energy. When free atoms
a coalesce to form molecules, higher energy orbitals get occupied—
electrons get promoted to higher orbitals. This leads to the natural
where fa0 are weights chosen to select states with given energies, as
definition
in scanning tunneling microscopy simulations[59]. Mulliken
charges, pertinent to DFTB, yield LDOS with atom resolution, X 
Eprom ¼ qðlÞ  qfree
ðlÞ
atom
H0ll : ð111Þ
X r
l
LDOSðe; IÞ ¼ d ðe  ea ÞqI;a ; ð107Þ
a
Promotion energy is the price atoms have to pay to prepare them-
which can be used to project density for group of atoms R as selves for bonding. Noble gas atoms, for instance, cannot bind to
X other atoms, because the promotion energy is too high due to the
LDOSR ðeÞ ¼ LDOSðe; IÞ: ð108Þ large energy gap; any noble gas atom could in principle promote
I2R
electrons to closest s-state, but the gain from bonding compared
For instance, if systems consists of surface and adsorbed molecule, to the cost in promotion is too small.
we can plot LDOSmol ðeÞ and LDOSsurf ðeÞ to see how states are distrib- Covalent bond energy, on the other hand, is the energy reduc-
uted; naturally LDOSmol ðeÞ þ LDOSsurf ðeÞ ¼ LDOSðeÞ. tion from bonding. We define covalent bond energy as [61]
Similar recipes apply for projected density of states, where DOS
Ecov ¼ ðEBS  Efree atoms Þ  Eprom : ð112Þ
is broken into angular momentum components,
X r X This definition can be understood as follows. The term
PDOSðe; lÞ ¼ d ðe  ea Þ qlI;a ; ð109Þ
a I
ðEBS  Efree atoms Þ is the total gain in band-structure energy as atoms
P coalesce; but atoms themselves have to pay Eprom , an on-site price
such that, again l PDOSð e; lÞ ¼ DOSðeÞ. that does not enhance binding itself. Subtraction gives the gain
the system gets in bond energies as it binds together. More
10.4. Mayer bond-order explicitly,
X
Bond strengths between atoms are invaluable chemical infor- Ecov ¼ EBS  qðlÞ H0ll ð113Þ
mation. Bond order is a dimensionless number attached to the l
X
bond between two atoms, counting the differences of electron ¼ qlm ðH0lm  elm Slm Þ; ð114Þ
pairs on bonding and antibonding orbitals; ideally it is one for sin- lm
gle, two for double, and three for triple bonds. In principle, any
bond strength measure is equally arbitrary; in practice, some mea- where
sures are better than others. A measure suitable for many purposes
in DFTB is Mayer bond-order [57], defined for bond IJ as 1 0
X elm ¼ ðH þ H0mm Þ: ð115Þ
M IJ ¼ ðq
~ SÞlm ðq
~ SÞml : ð110Þ 2 ll
l2I;m2J
This can be resolved with respect to orbital pairs and energy as
The off-diagonal elements of q ~ S can be understood as Mulliken X r
overlap populations, counting the number of electrons in the over- Ecov;lm ðeÞ ¼ d ðe  ea Þqalm ðH0ml  eml Sml Þ: ð116Þ
a
lap region—the bonding region. It is straightforward, if necessary, to
partition Eq. (110) into pieces, for inspecting angular momenta or Ecov;l;m can be viewed as the bond strength between orbitals l and
eigenstate contributions in bonding. Look at Refs. [60,57] for further m—with strength directly measured in energy; negative energy
details, and Table 1 for examples of usage. means bonding and positive energy antibonding contributions.
Sum over atom orbitals yields bond strength information for atom
10.5. Covalent bond energy pairs
XX
Another useful bonding measure is the covalent bond energy, Ecov;IJ ðeÞ ¼ ðEcov;lm ðeÞ þ c:c:Þ ð117Þ
which is not just a dimensionless number but measures bonding l2I m2J
directly using energy [61].
and sum over angular momentum pairs
Table 1
ll ¼la lm ¼lb
Simplistic electronic structure and bonding analysis for selected systems: C2 H2 (triple X X  
CC bond), C2 H4 (double CC bond), C2 H6 (single CC bond), benzene, and graphene (C- Elcov
a lb
ðeÞ ¼ Ecov;lm ðeÞ þ ½c:c: ð118Þ
point calculation with 64 atoms). We list most energy and bonding measures l m
introduced in the main text (energies in eV).

Property C2 H2 C2 H 4 C2 H 6 Benzene Graphene

gives bonding between states with angular momenta la and lb (no
c.c. for la ¼ lb ). For illustration, we plot covalent bonding contribu-
qH 0.85 0.94 0.96 0.95
AH 1.23 0.45 0.26 0.32
tions in graphene in Fig. 4.
ABH 2.75 3.27 3.34 3.39
Eprom;H 0.97 0.41 0.25 0.30 10.6. Absolute atom and bond energies
qC 4.15 4.13 4.12 4.05 4.00
AC 5.86 5.86 5.78 6.48 6.94
ABC 9.16 9.74 10.45 9.74 9.78 While Mayer bond-order and Ecov are general tools for any tight-
Eprom;C 5.62 5.69 5.64 6.46 6.94 binding flavor, neither of them take electrostatics or repulsions be-
BCH 8.14 7.90 7.84 7.93 tween atoms into account. Hence, to conclude this section, we
BCC 22.01 15.82 9.39 13.01 12.16 introduce a new analysis tool that takes also these contributions
M CH 0.96 0.95 0.97 0.96
into account.
M CC 2.96 2.02 1.01 1.42 1.25
The DFTB energy with subtracted free-atom energies,
 P. Koskinen, V. Mäkinen / Computational Materials Science 47 (2009) 237–253 249

11. Conclusions

Here ends our journey with DFTB for now. The road up to this
point may have been long, but the contents have made it worth
the effort: we have given a detailed description of a method to
do realistic electronic structure calculations. Especially the trans-
parent chemistry and ease of analysis makes DFTB an appealing
method to support DFT simulations. With these features tight-
binding methods will certainly remain an invaluable supporting
method for years to come.

Acknowledgements

One of us (PK) is greatly indebted for Michael Moseler, for intro-

ducing molecular simulations in general, and DFTB in particular.
Academy of Finland is acknowledged for funding though Projects
121701 and 118054. Matti Manninen, Hannu Häkkinen, and Lars
Fig. 4. Covalent bonding energy contributions in graphene (C-point calculation Pastewka are greatly acknowledged for commenting and proof-
with 64 atoms in unit cell). Both bonding and antibonding ss bonds are occupied, reading the manuscript.
but sp and pp have only bonding contributions. At Fermi-level (zero-energy) the pp-
bonding (the p-cloud above and below graphene) transforms into antibonding
(p )—hence any addition or removal of electrons weakens the bonds. Appendix A. Calculating the DFT pseudo-atom

The pseudo-atom, and also the free atom without the confine-
ment, is calculated with LDA-DFT [62]. More recent xc-functionals
E0 ¼ E  Efree atoms ð119Þ could be used, but they do not improve DFTB parametrizations,
1 X X X whereas LDA provides a fixed level of theory to build foundation.
¼ TrðqH0 Þ þ c Dq Dq þ V IJrep  qfree
ðlÞ
atom
Hll ; ð120Þ
2 IJ IJ I J I<J l
However, better DFT functionals can—and should be—used in the
repulsive potential fitting; see Section 4.
can be rearranged as Elements with small atomic numbers are calculated using non-
relativistic radial Schrödinger equation. But for some elements,
X X such as gold, chemistry is greatly modified by relativistic effects,
E0 ¼ AI þ BIJ ; ð121Þ
I I<J
which have to be included in the atom calculation.
In the four-component Dirac equation with central potential
good quantum numbers are energy, total angular momentum j,
where
its z-component jz , and j which is the eigenvalue of the operator
 
1 X  RLþ1 0
AI ¼ cII Dq2I þ qðlÞ  qfree
ðlÞ
atom
H0ll ð122Þ K¼ ; ð126Þ
2 l2I
0 R  L  1
1 where L is the orbital angular momentum operator and the compo-
¼ c Dq2 þ Eprom;I ð123Þ
2 II I nents of 4  4 relativistic spin-matrix R are
 
rk 0
and Rk ¼ ð127Þ
0 rk
X  
BIJ ¼ V IJrep þ cIJ DqI DqJ þ qml H0lm  Slm elm þ c:c: ð124Þ with 2  2 Pauli spin-matrices rk . Remember that angular momen-
l2I;m2J tum l is not a good quantum number; the upper and lower two com-
ponents are separately eigenstates of L2 with different angular
We call AI the absolute atom energy of atom I, and BIJ the absolute momenta, and coupled by spin–orbit interaction. In other words,
bond energy of bond IJ. Measuring atom energies with AI and bond a given l (that we are interested in) appears in 4-component spinors
energies with BIJ explicitly takes into account all energetics. Eq. states with different j.
(121) can be further simplified into The intention is to include relativistic effects from the Dirac
equation, but still use the familiar non-relativistic machinery. This
!
X 1X X can be achieved by ignoring the spin–orbit interaction, decoupling
0
E ¼ AI þ BIJ ¼ ABI ; ð125Þ upper and lower components. By considering only the upper com-
I
2 J–I I
ponents as a non-relativistic limit, l becomes again a good quan-
tum number. The radial equation transforms into the scalar-
where ABI measures how much atom I contributes to total binding
relativistic equation [63,64]
energy—in electrostatic, repulsive, promotive, and bonding sense.
For homonuclear crystals the binding energy per atom is directly 2  
d Rnl ðrÞ lðl þ 1Þ
ABI , and for heteronuclear systems the binding energy per atom  þ 2MðrÞðV s ðrÞ  e nl Þ Rnl ðrÞ
dr2 r2
(the negative of cohesion energy) is averaged ABI ; positive ABI  
1 dMðrÞ dRnl ðrÞ Rnl ðrÞ
means atom I would rather be a free atom, even though for charged  þ hji ¼ 0: ð128Þ
MðrÞ dr dr r
systems the interpretation of these numbers is more complicated.
Visualizing AI , BIJ , and ABI gives a thorough and intuitive measure Here a ¼ 1=137:036 is the fine structure constant,
of energetics; see Table 1 for illustrative examples. Note that AI ,
Z
BIJ , and ABI come naturally from the exact DFTB energy expres- V s ðrÞ ¼  þ V H ðrÞ þ V LDA
xc ðrÞ þ V conf ðrÞ ð129Þ
r
sion—they are not arbitrary definitions.
250 P. Koskinen, V. Mäkinen / Computational Materials Science 47 (2009) 237–253

with or without the confinement, and we defined

a2
MðrÞ ¼ 1 þ ½enl  V s ðrÞ: ð130Þ
2
The reminiscent of the lower two components of Dirac equation is
hji, which is the quantum number j averaged over states with dif-
ferent j, using the degeneracy weights jðj þ 1Þ; a straightforward
calculation gives hji ¼ 1.
Summarizing shortly, for given l the potential in the radial
equation is the weighted average of the potentials in full Dirac the-
ory, with ignored spin–orbit interaction. This scalar-relativistic
treatment is a standard trick, and is, for instance, used routinely
for generating DFT pseudo-potentials[64].
The pseudo-atom calculation, as described, yields the localized
basis orbitals we use to calculate the matrix elements. Our conven-
tions for the real angular part Y e l ðh; uÞ of the orbitals u ðrÞ ¼
l
e
Rl ðrÞ Y l ðh; uÞ are shown in Table 2. The sign of Rl ðrÞ is chosen, as
usual, such that the first antinode is positive.

Appendix B. Slater–Koster transformations

Consider calculating the overlap integral

Z
3
Sxx ¼ px ðrÞpx ðr  RÞd r ð131Þ

for orbital px at origin and another px orbital at R, as shown in

Fig. 5a. We rotate the coordinate system passively clockwise, such
that the orbital previously at R shifts to R ¼ R^z in the new coordi-
nate system. Both orbitals become linear combinations of px and Fig. 5. Illustrating overlap integral calculation. (a) Originally one px orbital locates
pz in the new coordinate system, px ! px sin a þ pz cos a, and the at origin, another px orbital at R. (b) Coordinate system is rotated so that another
overlap integral becomes a sum of four terms orbital shifts to R^z; this causes the orbitals in the new coordinate system to become
Z linear combinations of px and pz . Hence the overlap can be calculated as the sum of
2 the so-called SðpppÞ-integral in (c), and SðpprÞ-integral in (d). The integrals in (e)
px ðrÞpx ðr  R^zÞ  sin a ðFig:5cÞ and (f) are zero by symmetry.
Z
þ pz ðrÞpz ðr  R^zÞ  cos2 a ðFig:5dÞ where x ¼ cos a is the direction cosine of R. For simplicity we
Z ð132Þ assumed y ¼ 0, but the equation above applies also for y – 0
þ px ðrÞpz ðr  R^zÞ  cos a sin a ðFig:5eÞ (using hindsight we wrote ð1  x2 Þ instead of z2 ). The integrals
Z Z
3
þ pz ðrÞpx ðr  R^zÞ  cos a sin a ðFig:5fÞ: SðpprÞ ¼ pz ðrÞpz ðr  R^zÞd r
Z ð134Þ
3
The last two integrals are zero by symmetry, but the first two terms SðpppÞ ¼ px ðrÞpx ðr  R^zÞd r
can be written as
are called Slater–Koster integrals and Eq. (133) is called the
Sxx ¼ x2 SðpprÞ þ ð1  x2 ÞSðpppÞ; ð133Þ
Slater–Koster transformation rule for the given orbital pair (orbitals

Table 2
R
Spherical functions, normalized to unity ( Y e / Y lm
Y . Note that angular momentum l remains a good
e ðh; uÞ sin hdhdu ¼ 1) and obtained from spherical harmonics as Y
lm
quantum number for all states, but magnetic quantum number m remains a good quantum number only for s-, pz -, and d3z2 r2 -states.

Visualization Symbol e ðh; uÞ

Y
sðh; uÞ p1
ffiffiffiffiffi
4p
qffiffiffiffiffi
px ðh; uÞ 3
u
4p sin h cos
qffiffiffiffiffi
py ðh; uÞ 3
u
4p sin h sin
qffiffiffiffiffi
pz ðh; uÞ 3
4p cos h
qffiffiffiffiffiffiffi
d3z2 r2 ðh; uÞ 5 2
16pð3 cos h  1Þ
qffiffiffiffiffiffiffi
dx2 y2 ðh; uÞ 15 2
h cos 2u
16p sin
qffiffiffiffiffiffiffi
dxy ðh; uÞ 15 2
h sin 2u
16p sin
qffiffiffiffiffiffiffi
dyz ðh; uÞ 15
u
16p sin 2h sin
qffiffiffiffiffiffiffi
dzx ðh; uÞ 15
u
16p sin 2h cos
 P. Koskinen, V. Mäkinen / Computational Materials Science 47 (2009) 237–253 251

may have different radial parts; the notation Slm ðpprÞ stands for ra- tions with an asterisk can be manipulated to yield another 16
dial functions Rl ðrÞ and Rm ðrÞ in the basis functions l and m). Similar transformations and the remaining ones can be obtained by inver-
reasoning can be applied for other combinations of p-orbitals as sion, which effectively changes the order of the orbitals. This inver-
well—they all reduce to Slater–Koster transformation rules involv- sion changes the sign of the integral according to orbitals’ angular
ing SðpprÞ and SðpppÞ integrals alone. This means that only two momenta ll and lm ,
integrals with a fixed R is needed for all overlaps between any
two p-orbitals from a given element pair. Slm ðsÞ ¼ Sml ðsÞ  ð1Þll þlm ; ð136Þ
Finally, it turns out that 10 Slater–Koster integrals, labeled ddr,
because orbital parity is ð1Þl .
ddp, ddd, pdr, pdp, ppr, ppp, sdr, spr, and ssr, are needed to trans-
The discussion here concentrated only on overlap, but the Sla-
form all s-, p-, and d- matrix elements. The last symbol in the nota-
ter–Koster transformations apply equally for Hamiltonian matrix
tion, r, p, or d, refers to the angular momentum around the
elements.
symmetry axis, and is generalized from the atomic notation s, p,
d for l ¼ 0; 1; 2.
Table 3 shows how to select the angular parts for calculating Appendix C. Calculating the Slater–Koster integrals
these 10 Slater–Koster integrals. The integrals can be obtained in
many ways; here we used our setup with the first orbital at the ori- C.1. Overlap integrals
gin and the second orbital at R^z. This means that we set the direc-
tion cosines z ¼ 1 and x ¼ y ¼ 0 in Table 4, and chose orbital pairs We want to calculate the Slater–Koster integral
accordingly.
Finally, Table 4 shows the rest of the Slater–Koster transforma- Slm ðsÞ ¼ huls1 jums2 i ð137Þ
tions. The overlap Slm between orbitals ul at Rl ¼ 0 and um at
Rm ¼ R is the sum with
X
Slm ¼ cs Slm ðsÞ; ð135Þ
s uls1 ðrÞ ¼ Rl ðrÞ Ye s1 ðh; uÞ ¼ Rl ðrÞHs1 ðhÞUs1 ðuÞ ð138Þ

for the pertinent Slater–Koster integrals s (at most three). The gra- and
dients of the matrix elements come directly from the above expres-
sion by chain rule: Slater–Koster integrals Slm ðsÞ depend only on R ums2 ðrÞ ¼ Rm ðrÞ Ye s2 ðh; uÞ ¼ Rm ðrÞHs2 ðhÞUs2 ðuÞ; ð139Þ
and the coefficients cs only on R.b
For 9 orbitals (one s-, three p-, and five d-orbitals) 81 transfor- e s ðh; uÞ are
where RðrÞ is the radial function, and the angular parts Y i

mations are required, whereas only 29 are in Table 4. Transforma- chosen from Table 3 and depend on the Slater–Koster integral s in

Table 3
Calculation of the 10 Slater–Koster integrals s. The first orbital (with angular part s1 ) is at origin and the second orbital (with angular part s2 ) at R^z. /s ðh1 ; h2 Þ is the function
R 2p e ðh1 ; uÞ Y
e s ðh2 ; uÞ, and is used in Eq. (142). The images in the first row
resulting from azimuthal integration of the real spherical harmonics (of Table 2) /s ðh1 ; h2 Þ ¼ u¼0 du Y s1 2

visualize the setup: orbital (o) is centered at origin, and orbital (x) is centered at R^z; shown are wave function isosurfaces where the sign is color-coded, red (light grey) is positive
and blue (dark grey) negative.

s e s ðh1 ; uÞ
Y 1
e s ðh2 ; uÞ
Y 2
/s ðh1 ; h2 Þ

ddr d3z2 r2 d3z2 r2 5

8 ð3 cos
2
h1  1Þð3 cos2 h2  1Þ

ddp dzx dzx 15

4 sin h1 cos h1 sin h2 cos h2

ddd dxy dxy 15 2 2

16 sin h1 sin h2

pffiffiffiffi
pdr pz d3z2 r2 15 2
4 cos h1 ð3 cos h2  1Þ

pffiffiffiffi
pdp px dzx 45
4 sin h1 sin h2 cos h2

ppr pz pz 3
2 cos h1 cos h2

ppp px px 3
4 sin h1 sin h2

pffiffi
sdr s d3z3 r2 5 2
4 ð3 cos h2  1Þ

pffiffi
spr s pz 3
2 cos h2

ssr s s 1
2
252 P. Koskinen, V. Mäkinen / Computational Materials Science 47 (2009) 237–253

Table 4 Hlm ðsÞ ¼ huls1 jH0 jums2 i ð143Þ

Slater–Koster transformations for s-, p-, and d-orbitals, as first published in Ref. [36].
To shorten the notation we used a ¼ x2 þ y2 and b ¼ x2  y2 . Here x, y, and z are the is mostly similar to overlap. The potentials V s;I ½n0;I ðrÞ in the
b with x2 þ y2 þ z2 ¼ 1. Missing transformations are obtained by
direction cosines of R,
Hamiltonian
manipulating transformations having an asterisk ( ) in the third column, or by
inversion. Here m is ul ’s angular momentum and n is um ’s angular momentum. 1
H0 ¼  r2 þ V s;I ½n0;I ðrÞ þ V s;J ½n0;J ðrÞ ð144Þ
l At 0 m At R cmnr cmnp cmnd 2
s s
*
1 with l 2 I and m 2 J, are approximated as
s px x
* pffiffiffi
s dxy 3xy V s;I ½n0;I ðrÞ  V conf
s;I ðrÞ  V conf;I ðrÞ; ð145Þ
pffiffiffi
s dx2 y2 1
2 3b
s d3z2 r2 z2  12 a where V confis the self-consistent effective potential from the con-
s;I ðrÞ
px px * fined pseudo-atom, but without the confining potential. The reason-
x2 1  x2
px py * xy xy ing behind this is that while V conf ðrÞ yields the pseudo-atom and the
*
px pz xz
pffiffiffi
xz pseudo-atomic orbitals, the potential V s ½n0 ðrÞ in H0 should be the
px dxy 3x2 y yð1  2x2 Þ approximation to the true crystal potential, and should not be aug-
pffiffiffi
px dyz 3xyz 2xyz mented by confinements anymore. The Hamiltonian becomes
pffiffiffi
px dzx 3x2 z zð1  2x2 Þ
pffiffiffi 1
px dx2 y2 1 xð1  bÞ
2 3xb
pffiffiffi H0 ¼  r2 þ V conf
s;I ðrÞ  V conf;I ðrÞ þ V conf
s;J ðrÞ  V conf;J ðrÞ ð146Þ
py dx2 y2 1
2 3yb
yð1 þ bÞ 2
pffiffiffi
pz dx2 y2 1 zb
2 3zb and the matrix element
pffiffiffi
px d3z2 r2 xðz2  12 aÞ  3xz2
pffiffiffi conf
py d3z2 r2 yðz2  12 aÞ  3yz2 Hlm ðsÞ ¼ econf
l Slm ðsÞ þ huls1 jV s;J  V conf;I  V conf;J jums2 i ð147Þ
pffiffiffi
pz d3z2 r2 zðz2  12 aÞ 3za conf conf
* ¼ em Slm ðsÞ þ hu ls1 jV s;I  V conf;I  V conf;J jums2 i; ð148Þ
dxy dxy 3x2 y2 a  4x2 y2 z2 þ x2 y2
*
dxy dyz 3xy2 z xzð1  4y2 Þ xzðy2  1Þ depending whether we operate left with r2 =2 þ V conf or right
* s;I
dxy dzx 3x2 yz yzð1  4x2 Þ yzðx2  1Þ
3 1
with r2 =2 þ V conf
s;J ; ul ’s are eigenstates of the confined atoms with
dxy dx2 y2 2 xyb
2xyb 2 xyb
eigenvalues econf
l (including the confinement energy contribution
dyz dx2 y2 3
2 yzb
yzð1 þ 2bÞ yzð1 þ 12 bÞ
which is then subtracted). The form used in numerical integration is
dzx dx2 y2 3
2 zxb
zxð1  2bÞ xzð1  12 bÞ
pffiffiffi pffiffiffi pffiffiffi Z Z
dxy d3z2 r2 3xyðz2  12 aÞ 2 3xyz2 1 2
2 3xyð1 þ z Þ
pffiffiffi pffiffiffi pffiffiffi Hlm ðsÞ ¼ econf
l Slm ðsÞ þ dzqdqRl ðr 1 ÞRm ðr 2 Þ/s ðh1 ; h2 Þ
dyz d3z2 r2 3yzðz2  12 aÞ 3yzða  z2 Þ  12 3yza
dzx d3z2 r2
pffiffiffi pffiffiffi pffiffiffi h i
3xzðz2  12 aÞ 3xzða  z2 Þ  12 3xza
3 2
 V conf
s;J ðr 2 Þ  V conf;I ðr 1 Þ  V conf;J ðr 2 Þ : ð149Þ
dx2 y2 dx2 y2
4b a  b2 z2 þ 14 b2
pffiffiffi pffiffiffi pffiffiffi
dx2 y2 d3z2 r2 1 2 1
aÞ  3z2 b 1 2
2 3bðz  2 4 3ð1 þ z Þb As an internal consistency check, we can operate to um also directly
2
d3z2 r2 d3z2 r2 2 1
ðz  2 Þ a 3z2 a 3
a2 2
4 with r2 , which in the end requires just d Rm ðrÞ=dr 2 , but gives
otherwise similar integration. Comparing the numerical results
from this and the two versions of Eqs. (147) and (148) give way
to estimate the accuracy of the numerical integration.
question. Like in Appendix B, we choose ul to be at the origin, and Note that the potential in Eq. (147) diverges as r ! R^z, and the
um to be at R ¼ R^z. potential in Eq. (147) diverges as r ! 0. For this reason we use
Explicitly, two-center polar grid, centered at the origin and at R^z, where the
Z two grids are divided by a plane parallel to xy-plane, and intersect-
3


Slm ðsÞ ¼ d r Rl ðr 1 ÞHs1 ðh1 ÞUs1 ðu1 Þ  Rm ðr 2 ÞHs2 ðh2 ÞUs2 ðu2 Þ ; ing with the ^z-axis at 12 R  ^z.

ð140Þ
References
where r 1 ¼ r and r 2 ¼ r  R^z. Switching to cylindrical coordinates
[1] A.H. Castro Neto, F. Guinea, N.M.R. Peres, K.S. Novoselov, A.K. Geim, Rev. Mod.
we get
Phys. 81 (2009) 109.
Z Z [2] D. Porezag, T. Frauenheim, T. Köhler, G. Seifert, R. Kaschner, Phys. Rev. B 51
Slm ðsÞ ¼ dzqdqRl ðr 1 ÞRm ðr 2 Þ  Hs1 ðh1 ÞHs2 ðh2 Þ (1995) 12947.
[3] M. Elstner, D. Porezag, G. Jungnickel, J. Elstner, M. Haugk, T. Frauenheim, S.
Z 2p Suhai, G. Seifert, Phys. Rev. B 58 (1998) 7260.
 Us1 ðu1 ÞUs2 ðu2 Þdu; ð141Þ [4] M. Elstner, T. Frauenheim, E. Kaxiras, G. Seifert, S. Suhai, Phys. Stat. Sol. (b) 217
u¼0 (2000) 357.
[5] T. Frauenheim, G. Seifert, M. Elstner, Z. Hajnal, G. Jungnickel, D. Porezag, S.
and we see that since ^z-axis remains the symmetry axis, the u-inte- Suhai, R. Scholz, Phys. Stat. Sol. (b) 217 (2000) 41.
gration can be done analytically. The last two factors in Eq. (141) [6] P. Koskinen, H. Häkkinen, B. Huber, B. von Issendorff, M. Moseler, Phys. Rev.
Lett. 98 (2007) 015701.
become an analytical expression /s ðh1 ; h2 Þ, and is given in Table 3. [7] K.A. Jackson, M. Horoi, I. Chaudhuri, T. Frauenheim, A.A. Shvartsburg, Phys. Rev.
Note here that r is a spherical coordinate, whereas q is the distance Lett. 93 (2004) 013401.
from a ^z-axis in cylindrical coordinates. We are left with [8] P. Koskinen, H. Häkkinen, G. Seifert, S. Sanna, T. Frauenheim, M. Moseler, New
Z Z J. Phys. 8 (2006) 9.
[9] P. Koskinen, S. Malola, H. Häkkinen, Phys. Rev. Lett. 101 (2008) 115502.
Slm ðsÞ ¼ dzqdqRl ðr 1 ÞRm ðr 2 Þ/s ðh1 ; h2 Þ; ð142Þ [10] E. Bitzek, P. Koskinen, F. Gähler, M. Moseler, P. Gumbsh, Phys. Rev. Lett. 97
(2006) 170201.
[11] S. Malola, H. Häkkinen, P. Koskinen, Phys. Rev. B 77 (2008) 155412.
a two-dimensional integral to be integrated numerically.
[12] C. Köhler, G. Seifert, T. Frauenheim, Chem. Phys. 309 (2005) 23.
[13] T. Frauenheim, G. Seifert, M. Elstner, T. Niehaus, C. Köhler, M. Amkreutz, M.S.Z.
C.2. Hamiltonian integrals Hajnal, A.D. Carlo, S. Suhai, J. Phys.: Condens. Matter 14 (2002) 3015–3047.
[14] G. Seifert, J. Chem. Phys. A 111 (2007) 5609.
[15] N. Otte, M. Scholten, W. Thiel, J. Phys. Chem. A 111 (2007) 5751.
The calculation of the Slater–Koster Hamiltonian matrix [16] M. Elstner, J. Phys. Chem. A 111 (2007) 5614.
elements [17] G. Seifert, H. Eschrig, W. Bieger, Z. Phys. Chem. 267 (1986) 529–539.
 P. Koskinen, V. Mäkinen / Computational Materials Science 47 (2009) 237–253 253

[18] W.M.C. Foulkes, R. Haydock, Phys. Rev. B 39 (1989) 12520. [42] J.J. Mortensen, L.B. Hansen, K.W. Jacobsen, Phys. Rev. B 71 (2005) 035109.
[19] O.F. Sankey, D.J. Niklewski, Phys. Rev. B 40 (1989) 3979. [43] GPAW wiki at <https://wiki.fysik.dtu.dk/gpaw>.
[20] T. Frauenheim, F. Weich, T. Köhler, S. Uhlmann, D. Porezag, G. Seifert, Phys. [44] J.P. Perdew, K. Burke, M. Ernzerhof, Phys. Rev. Lett. 77 (1996) 3865.
Rev. B 52 (1995) 11492. [45] H. Rydberg, M. Dion, N. Jacobson, E. Scroder, P. Hyldgaard, S.I. Simak, D.C.
[21] G. Seifert, D. Porezag, T. Frauenheim, Int. J. Quant. Chem. 58 (1996) 185–192. Langreth, B.I. Lundqvist, Phys. Rev. Lett. 91 (2003) 126402.
[22] J. Widany, T. Frauenheim, T. Köhler, M. Sternberg, D. Porezag, G. Jungnickel, G. [46] T.A. Halgren, J. Am. Chem. Soc. 114 (1992) 7827.
Seifert, Phys. Rev. B 53 (1996) 4443. [47] M. Elstner, P. Hobza, T. Frauenheim, S. Suhai, E. Kaxiras, J. Chem. Phys. 114
[23] F. Liu, Phys. Rev. B 52 (1995) 10677. (2001) 5149.
[24] T.N. Todorov, J. Phys.: Condens. Matter 13 (2001) 10125–10148. [48] P. Koskinen, L. Sapienza, M. Manninen, Phys. Scripta 68 (2003) 74–78.
[25] B. Torralva, T.A. Niehaus, M. Elstner, S. Suhai, T. Frauenheim, R.E. Allen, Phys. [49] N.W. Ashcroft, N.D. Mermin, Solid State Physics, Saunders College Publishing,
Rev. B 64 (2001) 153105. 1976.
[26] T.A. Niehaus, D. Heringer, B. Torralva, T. Frauenheim, Eur. Phys. J. D 35 (2005) [50] H.J. Monkhorst, J.D. Pack, Phys. Rev. B 13 (1976) 5188.
467–477. [51] D. Frenkel, B. Smit, Understanding molecular simulation, From Algorithms to
[27] J.R. Reimers, G.C. Solomon, A. Gagliardi, A. Bilic, N.S. Hush, T. Frauenheim, A.D. Applications, Academic Press, 2002.
Carlo, A. Pecchia, J. Phys. Chem. A 111 (2007) 5692. [52] S. Malola, H. Häkkinen, P. Koskinen, Phys. Rev. B 78 (2008) 153409.
[28] T.A. Niehaus, S. Suhai, F.D. Sala, P. Lugli, M. Elstner, G. Seifert, T. Frauenheim, [53] C.P. Liu, J.W. Ding, J. Phys.:Condens. Matter 18 (2006) 4077.
Phys. Rev. B 63 (2001) 085108. [54] V.N. Popov, New J. Phys. 6 (2004) 17.
[29] Hotbit wiki at <https://trac.cc.jyu.fi/projects/hotbit>. [55] P.W. Atkins, R.S. Friedman, Molecular Quantum Mechanics, Oxford University
[30] <http://www.gnu.org/licenses/gpl.html>. Press, 2000.
[31] ASE wiki at <https://wiki.fysik.dtu.dk/ase/>. [56] O. Kit, P. Koskinen, in preparation.
[32] R.G. Parr, W. Yang, Density-Functional Theory of Atoms and Molecules, Oxford [57] I. Mayer, Simple Theorems, Proofs, and Derivations in Quantum Chemistry,
University Press, 1994. Springer, 2003.
[33] N. Bernstein, M.J. Mehl, D.A. Papaconstantopoulos, Phys. Rev. B 66 (2002) [58] B. Yoon, P. Koskinen, B. Huber, O. Kostko, B. von Issendorff, H. Häkkinen, M.
075212. Moseler, U. Landman, Chem. Phys. Chem. 8 (2006) 157–161.
[34] R.S. Mulliken, J. Chem. Phys. 23 (1955) 1833. [59] F. Yin, J. Akola, P. Koskinen, M. Manninen, R.E. Palmer, Phys. Rev. Lett. 102
[35] J. Junquera, O. Paz, D. Sanchez-Portal, E. Artacho, Phys. Rev. B 64 (2001) (2009) 106102.
235111. [60] A.J. Bridgeman, G. Cavigliasso, L.R. Ireland, J. Rothery, J. Chem. Soc., Dalton
[36] J.C. Slater, G.F. Koster, Phys. Rev. 94 (1954) 1498. Trans. 2001 (2001) 2095–2108.
[37] J.A. Pople, J. Chem. Phys. 37 (1962) 53. [61] N. Börnsen, B. Meyer, O. Grotheer, M. Fähnle, J. Phys.: Condens. Matter 11
[38] T.B. Boykin, R.C. Bowen, G. Klimeck, Phys. Rev. B 63 (2001) 245314. (1999) L287.
[39] M. Graf, P. Vogl, Phys. Rev. B 51 (1995) 4940. [62] J.P. Perdew, Y. Wang, Phys. Rev. B 45 (1992) 13244.
[40] I.N. Bronshtein, K.A. Semendyayev, G. Musiol, H. Muehlig, Handbook of [63] D.D. Koelling, B.N. Harmon, J. Phys. C 10 (1977) 3107.
Mathematics, Springer, 2004. [64] R.M. Martin, Electronic Structure Basic Theory and Practical Methods,
[41] <http://www.dftb.org>. Cambridge University Press, 2004.