Tight Density Functional Theory PHD Thesis

Improvements to the Density-Functional
Tight-Binding method: new, efficient

parametrization schemes and prospects
of a more precise self-consistency
PhD thesis
Zoltan Bodrog
Universitat Bremen
2012
Dem Promotionsausschuss an der Universitat Bremen, im

Bremen Center for Computational Materials Science
als Dissertation zur Erlangung des akademischen Grades
Doktor der Naturwissenschaften (doctor rerum naturalium)
vorgelegt.
(Tippfehler und
ahnliche Fehler wurden korrigiert.)
Tag der Einreichung 10. Januar 2012

Tag der m
undlichen Pr
ufung 23. M
arz 2012
Promotionskommission
Erstgutachter
Zweitgutachter
Pr
ufer
Pr
ufer
Vertreter der wissenschaftl. Mitarbeiter
Vertreter der Studenten
Prof. Dr. rer. nat.

Prof. Dr. rer. nat.
Prof. Dr. rer. nat.
Prof. Dr. rer. nat.
Dr. B
alint Aradi
Michael Wehlau
Thomas Frauenheim
Gotthard Seifert
Thomas Heine
Peter Maa
Contents
0.1
0.2
0.3
0.4
Abstract . . . . . .
Kurzfassung . . . .
Summarium operis
Osszefoglal
as . . .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
3
4
5
6
1 An introduction to the DFTB method

1.1 What we shall know about DFT now . . . . . . .
1.1.1 What is density functional theory? . . .
1.1.2 The HohenbergKohn theorems . . . . . .
1.1.3 The KohnSham equations . . . . . . . .
1.2 DFTB, the tight-binding approximation of DFT
1.2.1 The non-self-consistent DFTB . . . . . .
1.2.2 SCC-DFTB . . . . . . . . . . . . . . . . .
1.2.3 Further second-order semi-empirical terms
1.2.4 The current practice of parametrization .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
7
7
7
8
10
13
13
16
22
25
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
27
27
27
29
30
34
34
35
35
36
37
41
41
44
55
55
58
63
63
66
67
72
2 Improvements to handling repulsive energy

2.1 Automated parametrization . . . . . . . . . . . . . . . . . .
2.1.1 Least-squares fitting of repulsive potentials . . . . .
2.1.2 Fitting to multiple fit system types and objectives .
2.1.3 Fitting to objectives other than energy . . . . . . . .
2.1.4 Weighting of fit targets . . . . . . . . . . . . . . . .
2.1.5 Basis function shapes . . . . . . . . . . . . . . . . .
2.1.6 Further automation in the parametrization workflow
2.1.7 Optimizing electronic parameters? . . . . . . . . . .
2.2 Applications of the parametrizer . . . . . . . . . . . . . . .
2.2.1 Carbohydrogens . . . . . . . . . . . . . . . . . . . .
2.2.2 Zinc-oxides . . . . . . . . . . . . . . . . . . . . . . .
2.2.3 Titanium-organics . . . . . . . . . . . . . . . . . . .
2.2.4 Halogen-organics . . . . . . . . . . . . . . . . . . . .
2.3 Some theoretical aspects of parametrization . . . . . . . . .
2.3.1 On the theoretical validity of parametrization . . . .
2.3.2 Handling dissociation energy level . . . . . . . . . .
2.4 Ab initio repulsives . . . . . . . . . . . . . . . . . . . . . . .
2.4.1 Theory . . . . . . . . . . . . . . . . . . . . . . . . .
2.4.2 Further technical details of repulsive calculation . .
2.4.3 Results . . . . . . . . . . . . . . . . . . . . . . . . .
2.4.4 Prospects of calculating repulsives . . . . . . . . . .
1
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
3 Enhancement proposals to SCC-DFTB

3.1 Multipole expansion of SCC energy . . . . . . . . . . . .
3.2 Semiempirical . . . . . . . . . . . . . . . . . . . . . . .
3.2.1 Semiempirical off-site part . . . . . . . . . . . . .
3.2.2 s in line with subshell hardness . . . . . . . . .
3.2.3 Fully resolved s . . . . . . . . . . . . . . . . . .
3.2.4 Further resolution of s as operators . . . . . . .
3.3 Supplementary material to SCC . . . . . . . . . . . . . .
3.3.1 Mulliken charge density operators . . . . . . . .
3.3.2 Properties of the generalized functions . . . . .
3.3.3 Transformations of tabulated matrices before use
3.4 Multipoles in TD-DFTB . . . . . . . . . . . . . . . . . .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
75
75
77
77
78
78
78
83
83
84
85
87
Summary
89
Acknowledgements
91
Bibliography
92
0.1
Abstract
The Density-Functional Tight-Binding (DFTB) method of quantum-chemical

calculations is a pairwise LCAO-based approximation to the KohnSham method
of Density Functional Theory (DFT). To put a reader who is not fully familiar
with DFT, DFT-like methods and DFTB into context, and to define the chief
parts of nomenclature used, a brief introduction to DFT and a more deep one
to DFTB is given in this thesis first, with emphasis on the grounds of enhancements in later sections. This initial section is also intended to be a very quick
introduction to DFTB, with the minimum necessary insight into KohnSham
DFT for a reader who wants to understand the basics of DFTB quickly.
The second part is about improving the treatment of the so-called repulsive
part of DFTB energy. It is the quasi-classical energy part specific to DFTB,
it is calculated as a sum of pairwise effective potentials, the so-called repulsive
potentials, and it is added to the quantum-mechanically calculated KohnSham
DFT electronic energy to give the total energy. Parametrization is optimizing
these pairwise potentials for every type of atomic pairs, and it has been a lengthy
work needing vast amounts of human effort, that has been posing a high barrier
of entry into DFTB calculations with atom types not used before. I present a
thorough description of a parametrizer automaton to take over this task. The
automaton was developed in the course of doctoral research basing this thesis.
Intertwined with the description, there are several additional points of possible enhancement revealed, e.g. implementation of new energetical objectives of
fitting, and a perspective of possible optimization of electronic parameters. A
few examples are also shown to illustrate the capabilities of the parametrizer
automaton (of which the halogen-organic set is yet to be published [23]). This
section is based on the publication of this part of work, see [5].
As a completely different answer to the parametrization problem, fully derived formulae of calculating the DFTB pairwise repulsive energy, instead of
fitting, are presented at the end of the second part. Directly calculated repulsive
potentials may be at least a good starting guess in cases when quick calculations
of bigger error tolerance are needed. These results will be published soon in a
separate article [6].
In the third part of this thesis, several pieces of enhancement are proposed
to the electronic part of DFTB method itself. They affect the so-called SelfConsistent-Charges part of second-order DFTB, which implements iterative selfconsistent KohnSham calculations in DFTB, or more accurately an approximation of it based on the point-like representation of atomic charge fluctuations.
Improvements of SCC are centered around making the effective interaction between point-like atomic excess charges more precise. First, a multipole expansion is suggested instead of the point-like charges. Second, the effective interaction profiles between atomic charge centres can be improved by semiempirical calculations, instead of the present heuristic interpolation formulas between
chemical hardness and long-range Coulombic interaction. Most of the proposals
described here were published in [4].
0.2
Kurzfassung
Die density-functional tight-binding (DFTB) Methode der quantenchemischen

Kalkulationen ist eine paarorientierte Approximation der KohnSham-Methode
der Dichtefunktionaltheorie (DFT). F
ur die, die DFT, bzw. DFTB nicht ganz
detailweise kennen, und auch um die weitere Entwicklungen vorzubereiten, beginnt eine kurze Einleitung zu DFT und eine mehr tiefgreifende Einleitung zu
DFTB diese Dissertation. Dieses erste Teil ist vorgehabt auch eine sehr kurze
Einleitung f
ur die zu sein, die DFTB mit dem kleinsten Eingriff zu DFT oder
mit DFT-Hintergrund schnell verstehen mochten.
In dem zweiten Teil werden Verbesserungen der Handlung von dem sogenannten repulsiven Anteil der DFTB-Energie beschrieben. Es ist das quasiklassische Energieanteil, das spezifisch zu DFTB ist, und es wird als eine Summe
von Paarpotentialen die sogenannten Repulsivpotentialen neben dem elektronischen Energieanteil gerechnet. Die Fittung von dieser Paarpotentialen
f
ur jedem Typ von Atompaaren, die Parametrisierung heisst, ist ohne Automatisierung ein sehr langer Prozess, der intensive Menschliche Arbeit nimmt,
und deshalb eine hoche Barrierie gegen neue DFTB-Kalkulationen anzufangen gibt. Hier wird, aber, ein Automat, das die Parametrisierung machen
kann, vorgeschlagen und detailreich beschrieben. Die Parametrisierungsautomat
wurde im Ramen der Doktorandforschung, die auch das Material dieser Dissertation gegeben hat, entwickelt. In der Beschreibung des Automaten werden auch
m
ogliche Richtungen der Weiterentwicklung geschrieben, z. B. Implementation
neuerer energetischen Zwecke der Parametrisierung, oder Optimierung elektronischer Parameter. Auch in diesem Teil werden einige Beispile der N
utzlichkeit
des Automaten dargestellt (von denen das Beispiel mit halogen-organischen Parametern wird in K
urze publiziert [23]). Dieses Teil ist auf [5] gegr
undet.
Am Ende des Kapitels wird ein komplett unterschiedliche Losung zum Parametrisierungsproblem vorgeschlagen. Vollig herabgeleitete Formeln der direkten
Kalkulation der Repulsivenergie anstatt der gefitteten Repulsivprofilen stehen
hier. Die direkt kalkulierte Repulsivenergie kann ein guter Tipp am Anfang
einer DFTB-Kalkulation sein, wann man schnelle Resultaten, aber keine hohere
Pr
azisit
at braucht. Diese Methode und Resultaten sind auch in K
urze zu publizieren [6].
In dem dritten Teil dieser Dissertation werden mehrere methodologische
Verbesserungen des elektronischen Anteils der DFTB Methode geschrieben. Sie
gehen um dem sogenannten Self-Consistent-Charges (selbstkonsistente Ladungen) Teil der DFTB von zweiter Ordnung, das approximative Selbstkonsistenz in
DFTB als eine KohnSham Methode implementiert (durch ein Bild von Punktladungen, die atomare Ladungsfluktuation representieren). SCC-Entwicklung
enth
alt zwei gr
oere Hauptgebiete. Das erste geht u
ber das Punktladungsbild
bei Ladungsfluktuationen, und schlagt eine Entwicklung der Ladung-Ladung
Wechselwirkung durch Multipolenreihen vor. Das zweite schlagt verbesserungen von dem Wechselwirkungsprofil zwischen Ladungen vor, von denen eine ist
die, die diese Profile semiempirisch anstatt der heutigen heuristischen Interpolation zwischen der atomaren H
arte und der coulombischen Grenzwert kalkuliert.
Fast alle Vorschl
age geschrieben in diesem Teil wurden in [4] publiziert.
0.3
Summarium operis
Methodus nomine Density-Functional Tight-Binding (DFTB) in chemia quantica est approximatio calculationum in theoria functionalis densitatis (DFT)
algorithmo KohnSham factarum, quae fundamenta in integralibus pro geminis
atomorum calculandis et basi minimali utendo habet. Contextum construendi
causa eis, qui non maxime periti in DFT, methodis similibus atque in DFTB
sint, nec non ad nomenclaturam sectionibus sequentibus scriptam introducendam habetur introductio in prima thesi, quae altitudines DFT minime, DFTB
autem maxime tractat. Introductio haec est scripta etiam eis, qui in DFT iam
periti aut DFT tangere non volentes methodum DFTB ipsam celerrime intellegere volunt.
In media thesi tractatur pars energiae in DFTB nomine energia repulsiva
appellata atque investigantur producenda parametra. Producenda parametra
est optimanda pars energiae quasi-classica, ad DFTB specifica, quae summa
potentialium effectivorum pro geminis atomorum, vulgo potentialium repulsivorum, calculatur et ad energiam methodo KohnSham calculatam additur. Haec
potentialia repulsiva omnibus geminis typorum atomorum adiustenda fuerunt
tamen labor immensis vires humanas maxime dissipans. Situatio haec obstitit
novis calculationibus methodo DFTB cum typis atomorum prior non praesentibus. Automaton problema hoc solvens propono et describo, quod in cursu
quaestionis doctoralis materiam huius dissertationis quoque producentis creatum est. Complexae in descriptione habentur etiam partes automati meliorem
faciendi, e.g. implementatio obiectivorum ad adiustanda parametra novorum
energeticalium et perspectiva optimandorum parametrorum electronicorum. In
hoc capite sunt etiam exempla nonnulla parametrisatorem automaticum multum bene functum illustrantia scripta (quorum parametra inter halogenia et
elementa organica optimata brevi tempore publicanda sunt [23]). Hoc caput
contenta maxime ab articulo has res describenti [5] ducit.
Altera solutio problematis parametrorum est calculatio potentialium repulsivorum arte ab initio facta, cuius formulae perfecte derivatae in ultimo capite
sunt scripta. Expectantur repulsiva sic calculata utilia esse in situationibus,
quibus calculationes celeriter currere oportet, errores autem magis tolerandi
sunt. Repulsiva hoc modo calculata sunt etiam publicanda [6].
In parte dissertationis tertia mutationes ad methodum DFTB ipsam et
partem electronicam meliorem fiendam propono. Has mutationes partem SelfConsistent-Charges (onera electrica sibi consistentia) methodi DFTB, implementantem calculationes sibi consistentes iterativas ad modum KohnSham
faciendas, sive praecisius approximationem earum oneribus electricis punctis
formatis utendam, afficiunt. SCC-DFTB possit melior fieri a melioribus interactionibus inter illas fluctuationes onerum electricorum. Ad primum continet
hoc thema expansionem multipolis constructam loco onerum punctis similium.
Ad secundum habentur hic interactionum functiones semiempirice calculatae,
loco interpolationis multum heuristicae. Plurimae harum theseon sunt etiam
in [4] publicatae.
0.4
Osszefoglal
as
A DFTB (Density-Functional Tight-Binding) modszer egy olyan kvantumkemiai

sz
amt
ogepes elj
ar
as, mely a s
ur
usegfunkcionalelmelet (Density-Functional Theory, DFT) KohnSham-fele megfogalmazasanak atomi parokon alapulo minimalis b
azis
u k
ozeltese. A DFT-ben, DFT-szer
u eljarasokban, illetve a DFTB-ben
melyen nem j
artas olvas
ok sz
am
ara, valamint a kesobbi fejlesztesek megalapoz
as
aul e disszert
aci
o elejen egy igen rovid bevezetes olvashato a DFT-be,
majd a DFTB-nek egy nagyon reszletes lerasa kovetkezik. Ezt a bevezeto
reszt sz
andekom szerint tank
onyveszer
uen is olvashatjak azok, akik a DFT
alapjainak minim
alis erintesevel, vagy kesz DFT-s hatterrel szeretnek gyorsan
megerteni a DFTB m
uk
odeset.
A m
asodik resz a DFTB-energi
anak az u
gynevezett repulzv reszet, e repulzv resz kezelesenek javt
as
at, illetve a DFTB parametrizalasat jarja kor
ul.
Maga a parametriz
al
as azoknak a parpotencialoknak (az u
gynevezett repulzv
potenci
aloknak) az optimaliz
al
as
at jelenti, amelyeknek az osszegekent a DFTBre jellemz
o kv
aziklasszikus energiareszt, a repulzv energiat szamoljuk az elektronok KohnSham-energi
aj
aval egy
utt ez adja a DFTB-ben szamolt teljes energi
at. Ezeknek a p
arpotenci
aloknak az optimalizalasa minden atompartpusra
egy rendkv
ul hosszadalmas es sok emberi erofesztest igenylo feladat, mely igen
nehezze teszi azt, hogy a DFTB-vel eddig nem szamolt atomtpusok kemiajat
sz
amoljuk. Javaslatot tesz
unk azonban egy automatikus parametrizaloeljarasra,
mely ezt a munk
at k
onnyebbe teszi, illetve nagyreszt elvegzi. A parametrizaloelj
ar
ast az e disszert
aci
o alapj
at kepezo doktori kutatasok kereteben fejlesztettem
ki. A parametriz
al
oautomata ler
as
aba szove annak lehetseges tovabbfejlesztesi
ir
anyait is felv
azolom, pelda ilyenre a jelenleg meglevokon fel
uli energetikai
celf
uggvenyek defini
al
asa, vagy az elektronikus parameterek optimalizacioja. E
fejezetben neh
any peld
aval illusztr
alom is az automata kepessegeit (a peldak
k
oz
ul a halogen-szerves parameterek onallo publikaciokent is meg fognak jelenni [23]). A fejezet eredmenyei publikaciokent meg is jelentek [5].
A parametriz
al
as problem
aj
anak eddigiektol teljesen eltero megoldasat adja
a repulzv energia k
ozvetlen, nem felempirikusan illesztett szamtasa. Ennek
kepleteit minden reszletre kiterjed
oen szinten kozlom. A szamtott repulzv
energia j
o kiindul
opont lehet olyan szamtasok eseten, ahol gyorsan kell eredmenyeket produk
alni, azonban az eredmenyek nagy pontossaga nem feltetlen
ul
sz
ukseges. E megold
as r
ovidesen publikaciokent is megjelenik [6].
A disszert
aci
o harmadik reszeben javaslatot teszek maganak a DFTB elj
ar
asnak, azaz az elektronikus resznek a fejlesztesere is. Ezek a fejlesztesi
javaslatok a DFTB-n bel
ul az SCC (Self-Consistent-Charges, azaz onkonzisztens t
oltesek) nev
u reszt erintik, mely a DFTB-n bel
ul masodrendben teszi
lehet
ove az
onkonzisztens KohnSham-fele szamtast, pontosabban annak egy
k
ozelteset, mely az atomi t
oltesfluktuaciok ponttoltesszer
u modelljen alapul.
Az SCC fejlesztese ezen t
oltesfluktu
aciok kolcsonhatasait teszi pontosabba. Ehhez az els
o m
odszer a pontt
oltesek helyett multipolussorok hasznalata. A masodik m
odszer az SCC pontost
as
ara a toltesek kozti effektv kolcsonhatasi profilok
pontosabb meghat
aroz
as
at t
uzi ki celul, ezen bel
ul is felempirikus modszerrel,
nem pedig az eddig jellemz
o heurisztikus interpolacioval, mely a kemiai kemenyseget a t
avoli Coulomb-k
olcs
onhat
assal kototte ossze. Az ebben a szakaszban
szerepl
o javaslatokat szinten publik
altuk [4].
Chapter 1
An introduction to the
DFTB method
1.1
What we shall know about DFT now
1.1.1
What is density functional theory?1
Although quantum mechanics brought a deeper-than-ever-before understanding

of microscopic world, and a revolution in physics connected therewith, practical
computations with the quantum mechanical machinery are impossibly complicated even with a perturbative treatment of fairly small systems, e.g. small
molecules. Solution of the
| = H|
(1.1)
stationary Schr
odinger equation for the ground state of the system is far from
trivial when is a multi-particle wavefunction with N particles. The spatial
uses
representation of in this case uses 3N spatial coordinates, while that of H
6N , and because of electron-electron interactions, all of these are inseparably
connected in the integro-differential Schrodinger equation.
However, quantum chemistry, as a prevalent, killer application of quantum
mechanics, quickly came out with the most successful simplification scheme of
the above problem, the density-based calculation of electronic configurations.
This means changing the Schrodinger equation or the equivalent
E() = |H|
= min., | = 1
(1.2)
variational problem to another variational problem, but the new problem searches
for the ground-state electronic density, not the ground-state wavefunction.
E() = min.,
(x)dx = N
(1.3)
x
In the presence of a Vext external electric potential, the above energy func1
For a thorough foundation and description of this topics, see [29].
CHAPTER 1. AN INTRODUCTION TO THE DFTB METHOD

tionals of wavefunction or density split to two parts,
1
E() = . . . (x1 , . . . , xN )
+
+
Vext (xi )
2
|x
xj |
i
i
i
i=j
(x1 , . . . , xN )dx1 . . . dxN = F () + Vext (x)(x)dx
= E() =F () + Vext (x)(x)dx.

(1.4)
Note that the last term in E() is trivially equal to Vext (x)(x)dx, and therefore the above equation is an implicit proof that the existence and validity of
density-based treatment equals to the existence of the F () = F () functional
which comprises kinetic energy as well as all terms of interaction between the N
electrons. This F () functional is called the universal part of energy functional,
since it is the same for every system containing N electrons.
The first attempt at such a scheme was the ThomasFermi method [15, 35]
with its simple kinetic energy functional and total neglection of exchange and
correlational energies (a term approximating exchange energy was later added
by Dirac [9] and the kinetic energy part was improved somewhat by Weizsacker
[37]), that is a very simplistic universal functional. This approximate densitybased method is far from fulfilling the (1.4) equivalence precisely, and serves
barely more than as a historical example today.
The very success of density-based quantum-chemical methods had begun
with the theoretically strict and exact formulation of density-functional theory
based on the HohenbergKohn theorems that establish the validity of handling
electronic density instead of wavefunction as the basic variable of calculation.
Validity is ensured by the existence of a F () universal functional that makes
the (1.4) equivalence exact. However, the actual form of this F () functional is
unknown, and thus it can only be approximated in practical calculations. Based
on the HohenbergKohn theorems, Kohn and Sham founded the most efficient
calculation method up to date that calculates density as that of non-interacting
electron-like quasiparticles.
1.1.2
The HohenbergKohn theorems
The first HohenbergKohn theorem establishes that different problems necessarily have different ground state electron densities, or equivalently, no two
different external potentials can lead to the same ground-state electronic density with the same number of electrons. The proof is simple: if both the V1 (x)
and V2 (x) external potentials give (x) as ground-state density and 1 , 2 are
the two ground-state wavefunctions for the two potentials, respectively, then
F (1 ) + (x)V1 (x)dx F (2 ) + (x)V1 (x)dx

(1.5)
since 1 is the ground state of the system with V1 . In the same way
F (2 ) + (x)V2 (x)dx F (1 ) + (x)V2 (x)dx

8
(1.6)
1.1. WHAT WE SHALL KNOW ABOUT DFT NOW

since 2 is the ground state of system 2 (here and in (1.4) we denoted the universal part of wavefunction-dependent energy functional with F , but although
it is not the same as the F density functional, this leads to no ambiguity). The
two inequalities together say that F (1 ) = F (2 ) and therefore the value of F
may be assigned to . With this, in a non-degenerate case
E1 = F () + (x)V1 (x) > E2 + (x) (V1 (x) V2 (x))

(1.7)
if we denote the ground-state energies of problem 1 and 2 with E1 and E2 , respectively. The inequality
holds because is the ground-state density of problem
2, and therefore F () + V2 = E2 . Similarly
E2 = F () + (x)V2 (x) > E1 (x) (V1 (x) V2 (x))

(1.8)
because is the ground-state energy of problem 1. If we add (1.7) and (1.8)
together, we get
E1 + E2 > E1 + E2
(1.9)
which is clearly a contradiction. We silently omitted treating degenerate cases,
when the sides of above inequalities become equal. The theorem also holds in
this more complicated case, but the proof is not that easy. As a plausibilization
instead of the proof, it must be noted that for the degenerate equations to hold,
it would need a very quirky situation: two ground states which are degenerate
ground states of two distinct Hamiltonians in parallel, but with the same density.
It is true therefore that a ground-state density can not belong to more than
one external potential, and thus ground-state electronic density determines the
external potential of the problem.
A corollary of the first HohenbergKohn theorem is that since the ground
state determines the potential, it determines the Hamiltonian as well, because
the universal part is a constant operator. Through the Hamiltonian, the groundstate density of a quantum-chemical problem determines the whole spectrum,
the ground-state wavefunctions and all of other eigenstates of it.
Another corollary of the above theorem can be the first stage of the proof:
the universal part of energy functional exists and it is well-defined. Its actual
form is not known, and is of unknown complexity, though.
The second HohenbergKohn theorem states that the ground state density
and energy of a quantum-chemical system can be calculated minimizing the
E() = F () + (x)Vext (x)dx

(1.10)
energy functional. The proof can be regarded as a third corollary of the first theorem, it is based on having a ground-state wavefunction determined by Vext and
N . Since this ground-state wavefunction is at the minimum of the wavefunctionbased energy functional, and E() = E( ), E() must take its minimum at
, the charge density determined by .
A further corollary of the HohenbergKohn theorems (the Levy constrainedsearch formulation) is that if the actual form of F () were known, the groundstate of a system would be found minimizing the density-dependent energy
functional first, then getting the 0 ground-state density, and finally minimizing
9

F () over the wavefunctions which have 0 as their density. This corollary is
again based on the energy-minimum property of ground state, and it introduces
a theoretical implementation of F ():
F () = min F ().
=
1.1.3
(1.11)
The KohnSham equations
The invention of Kohn and Sham about calculating the F () universal functional was calculating it with a method that is highly accurate compared to
earlier attempts, but circumvents the direct construction of an F () functional.
Instead, the density is calculated as a
(x) =
occ.
i(x)i(x)
(1.12)
density sum of independent, non-interacting electronic orbitals, the so-called

KohnSham orbitals, and being the KohnSham electrons non-correlated, the
kinetic energy part of universal functional is simply:
TKS () =
occ.
i(x) i(x)dx.
2
(1.13)
In the KohnSham machinery, the number of used orbitals equals to the number
of electrons, i.e. there are no unoccupied KohnSham orbitals in the original
scheme. Formulae can be upgraded, however, to the non-zero-temperature case
occ.
by changing
to
ni where ni are the appropriate occupation numbers.
i
To provide the collection of i KohnSham orbitals with a Schrodinger-like system of equations to calculate them, one places them into the noninteracting
reference system of the original problem. This reference system has an external
potential composed of the original external potential and the effective potential
of electrons, called the KohnSham potential:
VKS (x) = Vext (x) + VC (x) + Vxc (x)
where VC is the Coulomb potential of electrons
(y)
VC (x) =
dy,
|x y|
(1.14)
(1.15)
and Vxc is an effective potential arising from everything else, namely

the error of (1.13) KohnSham kinetic energy with respect to the true
one,
Coulomb exchange terms not present in the (1.15) classical Coulomb field,
correlational terms not present in the calculations with the multi-electron
KohnSham wavefunction, that is a determinant made of KohnSham
orbitals.
10
1.1. WHAT WE SHALL KNOW ABOUT DFT NOW

The KohnSham orbitals of noninteracting reference system obey a Schrodinger
equation of the form
KS |i = (TKS + VKS )|i = i |i
H
KS KohnSham Hamiltonian. This is equivalent to the
defining the H
occ
occ
KS |i +
0=
i|H
i i|i
i
(1.16)
(1.17)
variational problem if we use a zero-temperature occupation of states as the

simplest case. Here, the (1.12) condition of mimicking the electronic density of
the real system is included by Lagrangian multipliers.
To derive the KohnSham effective potential, we take total energy into parts
first.
E() = T + J() + K + Vext
= TKS () + (T TKS ()) + J() + K + Vext (1.18)

where quantities written as depending on are known to be a functional of
it, either in a known or a not known functional form, and K is the amount of
quantum mechanical Coulomb energy deviation from the classical one (Coulombic exchange and correlational terms). Making the implicit definition of Exc as
everything else more explicit,
Exc = E() J() TKS (),
(1.19)
whence it will be also evindent that Exc is also a well-defined functional of .

Knowing that all relevant parts of the energy are valid functionals of , we
can define the effective potential as
VKS (x) =
J() Exc ()
(E() TKS ()) =

+
.
(x)
(x)
(x)
(1.20)
The exchange-correlational part of KohnSham potential is the most challenging part of the whole construction. While all the other parts are quantummechanically motivated and well underpinned, any computationally reasonable
Vxc () exchange-correlational potential is far from exact and physically correct.
Nonetheless, approximations to it proved to be quite precise in quantum chemical calculations, although they are determined with surprisingly simple constructions. The simplest one is the local-density approximation (LDA), which
formulates the exchange-correlational energy as
Exc = ((x))(x)dx,
(1.21)
x
where : R R is an ordinary scalar to scalar function, thus making Vxc a

local quantity. The next sophisticated method of approximating Exc is the
generalized-gradient method (GGA), which uses almost the same structure,
11

but depends on not only itself, but also on some of its spatial derivatives. Since GGAs construction is able to take some nonlocality of and
therefore Vxc in a differential neighbourhood of x into account, it is sometimes
called a semilocal approximation of Exc . Further on this way, there are fully
nonlocal
approximations of exchange-correlational DFT energy available (i.e.
Exc =
(x)((x), (y))(y)dxdy, see e.g. [38]), but they are computationally very expensive, and thus they are not that widely applied as the former two
ones. Nonetheless, nonlocal xc energy kernels are promising tools to describe
e.g. the van der Waals interaction, that local and semilocal xc cannot do, and
this remains a strong motivation for further research on them.
In addition to the above scheme of non-locality, Exc can be also calculated
as a non-local quantity from the KohnSham orbitals (next to kinetic energy,
which is always calculated so in the leading order) as a sum of HartreeFock-type
exchange terms. In this case, all pairwise combinations of integrals

i(x)j(y) 1 i(y)j(x)dxdy
(1.22)
|x y|
y x
are summed up with the appropriate weights to give an exchange energy. As

DFT-based KohnSham exchangecorrelational energies and HartreeFock-based
ones tend to generate opposite errors in energetics, they are often combined with
some weighting to give a viable compromise. These blends of KohnSham and
HartreeFock exchange-correlational functionals are called hybrid xc functionals.
It can be clearly seen, that the KohnSham method is not a direct densitybased method as it would be imagined according to the HohenbergKohn theorems. This method employs the KohnSham orbitals as primary quantities to
express (x), and in addition to this, TKS () is also calculated from the orbitals,
in an algorithmically independent calculation (hybrid xc functionals calculate
even the xc energy partly with KS orbitals directly). With one can compute
the effective KohnSham potential more precisely, and the cycle begins with
calculating new KohnSham orbitals again. This algorithm renders the method
an iterative self-consistent-field approximation. Nevertheless, KohnSham DFT
proved to be by far the most efficient way of density-based quantum-chemical
calculations, when the trade-off between computational expense and precision
comes into question.
As a theoretical aspect of KohnSham method, one can notice that an intermediary energy term of it, the so-called KohnSham energy is
EKS = Tr(HKS P ) = Tr TKS P + VKS (x)(x)dx = Tr

P
(1.23)
P
x
where P is the one-electron density operator of the KohnSham system, and TKS
is the KohnSham kinetic energy operator. Notethat the kineticenergy operator is a constant with respect to P , this is why Tr P P Tr(TKS P ) = Tr(TKS P ).
This Legendre-transformation-like construction, as a thermodynamical analogy,
makes going over from the external potential to the ground-state density as
the basic determinative variable of the problem analogous to changing the N
particle number to the
E
=
(1.24)
N
12
1.2. DFTB, THE TIGHT-BINDING APPROXIMATION OF DFT

chemical potential as an independent variable in classical thermodynamics, that
E
term from the Helmholtz free energy to get the
implies distracting an N N
grand canonical potential. About this thermodynamical analogue and other
aspects of this analogy (e.g. even a classical DFT method!), see [3].
The above KohnSham energy is clearly different from the total energy of
the interacting system because the effective KohnSham potential contains energy terms arising from the electrons themselves, and therefore the VKS term
contains double-counting. After a KohnSham calculation of orbitals, one determines the total energy by distracting these double-counting terms, resulting
in
1
(x)V (x)dx (x)Vxc (x)dx + Exc ().
(1.25)
E = EKS
2
x
1.2
1.2.1
DFTB, the tight-binding approximation of

KohnSham DFT
The non-self-consistent DFTB
The original, non-self-consistent-charges (non-SCC) DFTB [33] is a tight-binding

DFT method to calculate electronic structures of chemical systems. It solves
the (1.16) KohnSham equations in an LCAO Hilbert space of wavefunctions
that is spanned by the electronic orbitals of the participating atoms.
KS |i = [T + VKS ] |i = i |i
H
(1.26)
where the VKS (x) KohnSham effective potential behind the VKS operator is the
(1.14) sum of the VC (x) Coulomb potential generated by the (x) charge density,
the Vxc (x) functional derivative of exchange-correlation energy with respect to
(x), and the Vext (x) external (non-electronic) potential of the atomic nuclei:
(y)
Exc [(x)]
VKS (x) =
dy +
+ Vext (x).
(1.27)
|x y|
(x)
y
Expressing all the |i KohnSham electronic orbitals in the atomic orbital basis,
|i =
ci, |
(1.28)
where | runs over all atomic orbitals around all atomic centres in the entire
basis. The stationary Schrodinger equation in (1.26) thus numerically becomes
a generalized eigenvalue problem:
H ci, = i
S ci,
(1.29)
KS | is the Hamiltonian matrix and S = | is the

where H = |H
overlap matrix (the matrix of the unit operator in this basis of atomic orbitals).
Strictly speaking, the basis does not consist of atomic, but pseudoatomic [14]
orbitals coming from pseudoatoms with a harmonic contractive potential in
addition to the Coulombic nuclear potential. This means
n
n
Z
r
r
= +
(1.30)
v (psat) (r) = v (at) (r) +
r0
r
r0
13

instead of the v (at) (r) atomic Coulomb potential (note that the pseudoatomic
potentials are used only at determining (pseudo)atomic quantities, like potentials and orbitals, Vext (x) is built of atomic ones). This construction makes our
atomic orbitals more compressed and thus much more apt for realistic LCAO
calculations.
In order to have more degrees of freedom with tuning electronic properties,
the above pseudoatomic potential can contain different compressive potentials
on top of the atomic Coulomb potential within one course of atomic ab initio
calculations that determine atomic densities, total atomic potentials and atomic
orbitals. In most cases, the first stages (determining overall electronic density
and potentials) of atomic ab initio calculations are carried out with the socalled density compression radius in the place of r0 , while calculation of actual
electronic orbitals used later in LCAO calculations is made with the so-called
wavefunction compression radius (denoted with r1 ) in the same place. Note
that the potentials calculated with the density compression radius are kept
and used in the second part, thus this is a bit of inconsistency brought in to
get more electronic parameters to be able to tune. Originally, the ability of
having multiple compression radii was invented to cancel the possible errors of
the so-called density superposition construction of molecular potential (see the
end of this section, after (1.36)) semiempirically, but has been being used very
often deliberately as another free parameter to improve electronic properties
since then. Further tunable parameters can be introduced by applying different
compression radii to orbitals of different angular momenta, but this is very rare,
though. The usual value of r0 is about five times the covalent radius of the
atom, while a typical r1 is twice the covalent radius. The whole construction
does not depend strongly on the value of n. This arbitrary n used to be 2 in
most cases.
With the above compressed pseudoatomic orbitals, it is even possible to get
good results without a self-consistent solution of equation (1.26). A Harrisfunctional-like [18] use of the KohnSham Hamiltonian leads to fairly good
energy values in one step, and thus, the non-SCC DFTB is a non-iterative
DFT-based method to calculate molecular electronics and chemical energies.
The total energy of the calculated chemical system is also calculated like in
KohnSham DFT:

KS 1 VC Vxc
+ Exc + Enuc
(1.31)
E = Tr P H
2
where P = i ni |i i| is the density operator of the electrons, ni being the

occupation number of state |i, Exc is the total exchange-correlational energy
and
nuclei
1 ZA ZB
Enuc =
(1.32)
2
rAB
A,B
is the energy of nuclear Coulomb repulsion.

The above total energy expression is further simplified in DFTB. First, the
atomic core electrons are frozen and they are not accounted for in the quantummechanically calculated total electronic energy, they only affect it via the contribution of their Pcores and cores (x) density to the VC and Vxc potentials (naturally, P = Pcores + Pvalence and = cores + valence ). Second, the energetics
14

KS ] is not calof atomic cores as well as the part of (1.31) after EKS = Tr[P H
culated with DFT methods, it is substituted with a set of pairwise effective
quasiclassical potentials (the repulsive potentials) acting between atoms:
Erep
= Tr Pvalence VC Vxc
2
+ Tr Pcores HKS VC Vxc

+ Exc + Enuc
2
nuclei
A=B
1
UZ Z (rAB )
2 A B
(1.33)
(we drop the valence index in the sequel and unindexed electronic densities will
denote those of the valence electrons from here on). The UZA ZB (rAB ) potentials
which depend on the type of atom pair AB and the rAB distance between
them are fitted preliminarily on fit systems against ab initio or experimental
energetics, in the process of parametrization.
Third, the KohnSham electronic energy part of total DFTB energy is also
simplified. When it is expressed with the basis wavefunctions
KS ] =
EKS = Tr[P H
KS |i =
ni i|H
ni ci, | T + VKS | ci, ,
(1.34)
i,,
the matrix elements of KohnSham potential in the sum break up to sums of

three-center terms
atoms
|VKS | =
| VB |
(1.35)
B
since VKS can also be regarded as a sum of VB atomic contributions2 . From

the above sum, proper three-center terms (A|B|C-like terms, where {A},
{C} i.e. and are centered at atoms A and C, respectively, and A = B,
B = C, C = A) are dropped because they are only very small corrections to the
energy, yet they are computationally expensive.
Two-center terms of the form A|B|A (the so-called crystal-field terms) are
also neglected due to their magnitude and two further reasons. First, A|A|Atype terms are treated in a special way: H = |T + VA | terms where ,
{A} are substituted with their respective free-atom values (matrix elements of
atomic potentials with atomic orbitals, where is the Kronecker delta
and is the atomic, not the pseudoatomic orbital energy) in order to recover
free-atom energies in the dissociation limit. Second, crystal-field terms (just like
three-center terms, if there were any) would mix valence orbitals with the frozen
core electronic states [33] and make a = cores +valence or P = Pcores + Pvalence
breakdown of density to core and valence contributions far from trivial (core
states not orthogonal to valence ones would imply a pseudopotential treatment).
These reasons leave no room for A|B|A-type contributions in A| . . . |A.
2 This is exactly true when the so-called potential superposition is used (V
molecular =
atoms
V (A )), but it is only almost true when density superposition is used (Vmolecular =
A
V ( atoms
A )) since exchange-correlational energy expressions are not linear. In the latter
A
case, however, the construction of pairwise potentials in (1.36) will handle the problem.
15

The remaining A|A|B and A|B|B-type terms are calculated as A|A + B|B
and they are calculated together with kinetic energy integrals
H (rAB ) = | T + VA + VB | ,
(1.36)
and they are tabulated as functions of atomic distances. These integral tables
are used by DFTB to calculate the KohnSham energy without actually calculating any integrals run-time. Building A + B-type Hamiltonian terms also
handles the problem of density superposition (see footnote 2). In this case,
V (A + B ) is used instead of VA + VB , and for the Hamiltonian matrix
atomselement
in question the difference between |V (A + B )| and |V ( A
A )|
can be approximated with the sum
|C (x)dx
(1.37)
| Vxc ()
d(x)
=A (x)+B (x)
x
C=A,C=B
which is a sum of three-center integrals (characterized by the A, B, C atomic

centres) and thus belongs to negligible order.
Shortly summarized, the following cases and types are distinguished calculating the elements of the KohnSham Hamilton matrix (using the [] notation
for the atomic centre of orbital ):
if [] = [],

= | T + VA + VB | if {A} and {B},
0
otherwise.
(1.38)
Thus the total non-SCC DFTB energy is a KohnSham energy calculated

from tabulated pairwise data perturbatively in one non-self-consistent step,
with the quasi-classical repulsive energy added. Being a one-step perturbative
method, it does not take the energetical effects of charge fluctuations between
atoms precisely into account. This behaviour makes non-SCC DFTB reliable
with chemical systems containing atoms of nearly the same electronegativity.
1.2.2
SCC-DFTB
To handle the energetics of charge fluctuations, the DFTB method was extended
to an iteratively solved self-consistent DFT-like method, called SCC-DFTB [12].
The iteration is done by constructing a new KohnSham Hamiltonian with the
new (x) charge density from the last step instead of the initial (0) (x) and
calculating KohnSham energies, LCAO KohnSham electronic states as well
as a more precise charge density in the next step, and so on. This iteration is
made until the KohnSham states converge. The actual realisation is detailed
below.
We write electronic charge density as a sum (x) = (0) (x) + (x) of
the starting charge density (the mere superposition of pseudoatomic electron
densities without interatomic interactions) and the fluctuations, and likewise
the density operator as P = P (0) + P .
As the charge density is in a very close linear relation with the density op KS [(x)] KohnSham Hamiltonian
erator, we will take the dependence of the H
16

KS as a Taylor series of the
on the charge density into account by expanding H
P density operator up to the first order in P . Note that since

KS = E ,
H
P
(1.39)
this is equivalent to examining the second-order term in the Taylor series of E

with respect to P
KS
2 E
H
=
(1.40)
P
P P
(where we silently did not define all the straightforward details about these
Frechet derivatives with respect to the valence density operator).
Since

1
E = Tr P T + VC (P ) + Vext
+ Exc (P ) + Ecores + Enuc
(1.41)
2
(where Ecores is the Tr[Pcores . . .] term in the (1.33) definition of repulsive energy), the second derivative of E with respect to P comes only from the Coulomb
and the exchange-correlation energy, the other parts either do not depend on P
or are only linear in it.
2 E
2 EC
2 Exc
=
+
P P
P P
P P
(1.42)
where EC = 21 Tr[P VC ].
We investigate the Coulomb part of the above derivative first. For this, (x),
the electronic charge density can be calculated as
(x) = x| P |x
(1.43)
where |x is the x position eigenstate. Furthermore

x| P |x = x|
ni |ii| |x =
ni |i(x)|2 ,
(1.44)
or expanding the i states according to (1.28) in the atomic orbital basis
x|
ni |ii| |x
ni ci, ci, x||x
i,,
ni ci, ci, (x)(x)
(1.45)
i,,
where i(x), (x) and (x) are position-space wavefunction values of the i Kohn
Sham orbital and the , basis atomic orbitals, respectively.
17

The Coulomb potential operator VC belonging to the (x) charge distribution is
VC = VC (x) |x x| dx
x
(y)
|x x| dydx
|x y|
x,y
y| P |y
|x x| dydx.
|x y|
(1.46)
x,y
The Coulomb part of the total energy is
y| P |y x| P |x
1
1
dxdy.
EC = Tr(P VC ) =
2
2
|x y|
(1.47)
x,y
Using (1.45) for x| P |x and y| P |y

EC =
1
2
ni nj cj, ci, ci, cj,
i,j,,,,
1
(y)(y)(x)dxdy.
(x)
|x y|
(1.48)
x,y
To step further, we rearrange the summation in (1.48):

EC =
1
2
A,B
ni nj cj, ci, ci, cj,
i,j,,
{A},{B}
1
(y)(y)(x)dxdy = 1
(x)
IAB
|x y|
2
(1.49)
A,B
x,y
defining the pairwise IAB subsums. With different but similar rearrangements,
we can also define other pairwise subsums: a constraint {A}, {B} in the
summation yields JAB , {A}, {B} yields KAB and {A}, {B}
yields LAB . With them
1
1
1
1
IAB =
JAB =
KAB =
LAB .
(1.50)
EC =
2
2
2
2
A,B
A,B
So
EC =
A,B
A,B
1
{IAB + JAB + KAB + LAB } .
8
(1.51)
A,B
The above summation can be taken apart to an on-site and an off-site part:
(on)
EC = EC
(off)
+ EC
1
{IAA + JAA + KAA + LAA }
8
A
1
+
{IAB + JAB + KAB + LAB } .
8
A=B
18
(1.52)

Off-site SCC energy
If, as an approximation, we suppose that all wavefunctions can be regarded
point-like relative to the distances of atoms, the off-site part of (1.52) goes into
(off)
EC
1
=
8
A=B
ni nj cj, cj, ci, ci,
i,j,,
{A},{B}
1
S S + . . .
|xA xB |
(1.53)
where xA is the position of the Ath atomic centre and S is the overlap matrix
(the point-like limit of IAB is presented in the equation, the limit of JAB +
KAB + LAB is symbolised by the dots). Defining
qA =
{A},,i
1
ci, ci, S + ci, ci, S )
ni (
2
(1.54)
as the Mulliken charge of atom A, the (1.53) form of off-site Coulomb energy
part in this point-like approximation easily becomes
(off)
EC
1 qA qB
2
rAB
(1.55)
A=B
where rAB = |xA xB | is the distance between atoms A and B.

The 1r Coulomb interaction profile of point-like charges must be modified to
take into account that these charges cannot be treated point-like approaching
each other at the typical molecular neighbour distances. The modified interaction profile is calculated as the interaction energy of two spherically symmetric
charge clouds with eA r and eB r radial distributions for atoms A and B.
This gives
1
(off)
EC =
AB qA qB
(1.56)
2
A=B
where
AB = (A , B , rAB )
(1.57)
1
and AB rAB
0 with rAB . The actual shape of AB at short distances
is determined by the A and B parameters [13] that depend only on the types
of atoms A and B, respectively.
On-site SCC energy

The on-site part is a bit different, as in this case the electronic orbitals cannot
be treated point-like in any aspect. To understand what the on-site part is, we
construct atomic Mulliken density operators first that give the foundation of
the (1.54) Mulliken charges (we use them only briefly here to give an insight on
the structure of the on-site energy; for more details, see 3.3.1):
1
PA =
2
ni (
ci, ci, || + ci, ci, ||) .
{A},,i
19
(1.58)

With them,
qA =
x| PA |x dx =
A (x)dx
(1.59)
Mulliken charge density of the Ath atom (of course,

A (x) being the atomic
A PA = P and
A A = ). With these partial atomic density operators we
write the on-site part
(on)
EC
1
=
2
A x,y
y| PA |y x| PA |x
dxdy
|x y|
=
1
2
A x,y
A (y)A (x)
dxdy.
|x y|
(1.60)
Thus we see that, just like the off-site part, the on-site part of Coulomb energy
is quadratic in PA or A .
The above on-site part of the Coulomb energy cannot be further simplified
and an approximative formula, unlike the off-site case, cannot be derived. To
circumvent this, one determines its second derivatives with respect to atomic
charges semiempirically. If the energy of a lone atom A with qA electrons is EA ,
then according to Janaks theorem [19], its first derivative is
EA
= A
qA
(1.61)
where A is the energy of the highest occupied atomic orbit of the free atom.
Furthermore,
2 EA
A
=
= UA
(1.62)
2
qA
qA
that is nothing else than the Hubbard parameter or chemical hardness of element
ZA . If we set the parameters in (1.57) so that
AA = UA ,
(1.63)
the second-order change of EC with respect to the electronic density can be

written
2 E =
1
1
1
2
AB qA qB +
UA qA
=
AB qA qB .
2
2
2
A=B
(1.64)
A,B
Note that the above semi-empirical on-site energy derivative is based on the
total energy, not only the Coulomb energy of electrons. This means that not
only the second-order change of Coulomb energy in (1.60) is covered by it, but
both relevant parts of (1.42) are included (Coulomb and xc). In addition to
this, the analytical interpolation [13] that is realized by the s between the
on-site chemical hardness and the long-range off-site Coulomb interaction can
be regarded also as a coarse description of the gradually vanishing second-order
xc energy between two neighbouring atoms.
20

The Hamiltonian
In addition to 2 E, it is also important for the DFTB machinery to construct
KS SCC correction of the Hamiltonian. We will construct the Kohn
the H
Sham Hamiltonian as the (1.39) derivative of total energy with respect to the
density operator. We start with a total energy expression in the SCC fashion
E=
1
AB qA qB .
2
(1.65)
A,B
Then
KS
H
E
P
AB
A,B
qA
P
qB =
AB
A,B
qA
qB
P
(1.66)
where the last step makes matrix elements of the operator derivative explicit.
Since
ni ci, ci, ,
(1.67)
P =
i
(1.54) can be rewritten as

qA =
1
2
P S + P S
{A},
1
2
(P S + P S )
(1.68)
{A},
(the last equality coming from the self-adjointness of P ). With all this
1
2 S if {A},
1S
if {A},
qA
= 2
S
P
if , {A},
0
otherwise,
(1.69)
and finally
KS
H
1
= S
([]B + []B )qB
2
(1.70)
where [] and [] denote the respective atomic centres of the orbitals. Because
the above Hamiltonian is linear in the Mulliken charges, its change with respect
to the non-interacting atoms is
1
KS
H = H
= S
[]B + []B qB .
(1.71)
2
Refinement of the on-site SCC energy

Having only one general chemical hardness value for every atom type is quite
satisfactory for electronic systems with s or s and p electrons on the valence shell
because the respective second energy derivatives for those systems are quite close
to each other, regardless of the s or p orbit which the electrons are fluctuating
21

to or from. With d orbits, however, a refinement of the above U s is necessary.
If we split up atomic Mulliken charges to subshell Mulliken charges (or more
exactly, if we use subshellwise Pl , l and ql Mulliken quantities in (1.58), (1.59)
and (1.60)), we arrive at calculating a Hubbard matrix
Ull =
l
ql
(1.72)
for every atom type (atom type indices were suppressed for brevity, l and l
denote angular momenta of subshells). With the Hubbard matrix, the secondorder on-site energy change for every atom can be calculated as
1
(on)
ql Ull ql
(1.73)
Eone atom =
2
l,l
taking into account that electron fluctuation is not restricted to only one atomic
subshell.
We are aware of only one publicly available DFTB implementation using
different Hubbard parameters for different subshell interactions [2]. But even
there, the off-diagonal Hubbard matrix elements are approximated with means
of diagonal elements instead
Ull
Ull + Ul l
,
2
(1.74)
following the idea of K

ohler et al. [24]. The simplistic, semidiagonal U matrix is
a good approximation for elements up to the third period, but it is significantly
wrong above (subshell hardness errors going well beyond 10%, see Table 1.1)
if d orbitals are present in the atomic basis. The bare magnitude of error
in Ull increases as one goes to the halogens within a period of the periodic
table, but the importance of d orbitals decreases in parallel. The resulting
overall importance of off-diagonal Hubbard matrix element errors as well as their
effective importance at molecular energetics is still a question. Nevertheless, it
affects the above semidiagonal Hubbard matrix, that is of great importance
according to K
ohler et al., e.g. with calculating magnetic properties of iron
clusters.
In our treatment of SCC, we will continue using atomwise chemical hardnesses in derivations for the sake of brevity and understandability unless we
especially treat refined matrices of chemical hardnesses. The results thus written can be easily extended to matrices of chemical hardnesses.
1.2.3
Further second-order semi-empirical terms
Spin-spin interaction
In the collinear spin model of spin-dependent DFTB [25], the density matrix of
electrons depends on the z-aligned spin values too. Thus traces in the spatial
representation are summations not only along spatial coordinates, but along
spin indices too. This makes a second-order correction to energy like this:
1
2 E
2 E =
(x, )(y, )dxdy
(1.75)
2
(x, )(y, )
, x,y
22
Table 1.1: Subshell hardnesses calculated as approximative averages of diagonal

Hubbard matrix elements in (1.74) or proper numerical derivatives according to
(1.73).
Upp +Udd
atom
Upp
Udd
Upd
relative error
2
Fe
0.237 0.447 0.303
0.280
0.084
Zn
0.266 0.523 0.345
0.317
0.091
Ga
0.277 0.577 0.365
0.326
0.119
Br
0.376 0.825 0.502
0.428
0.171
where (x, ) = x, |P |x, , and , =, are spin indices (cf. this expres
2 E
sion with 21 (x)(y)
(x)(y)dxdy of the spin-independent case). After
introducing (x) = (x, ) (x, ), we can write the above energy correction
as
1
E=
2
2
2 E
2 E
(x)(y) +
(x)(y)
(x)(y)
(x)(y)
x,y
2 E
2 E
(x)(y) +
(x)(y)dxdy
(x)(y)
(x)(y)
through (x,)
and (x,)
.
= 12 (x)
+ (x)
= 12 (x)
(x)
+
(1.76)
E
After realizing that (x)(y)
must be zero due to symmetry reasons3 , the
second-order energy correction stemming from magnetization (spin polarization)
of electrons can be simplified. The second-order energy change thus becomes
1
2 E
2 E
2
E=
(x)(y) +
(x)(y)dxdy (1.77)
2
(x)(y)
(x)(y)
x,y
in which the first term is the spin-independent second-order correction constructed with the energy derivatives in (1.42). Its reformulation within the
DFTB formalism was treated in 1.2.2. The second parts simplification and inclusion into the DFTB total energy expression goes in the same way, but due to
the more localized nature of spin-spin interaction, the sum of interacting spin
populations is strictly restricted to on-site terms. The DFTB spin-spin term
similar to (1.64) is
1
WA p2A
(1.78)
2
A
in a coarse atomwise approximation, and

1
WA,ll pA,l pA,l
2
(1.79)
A,l,l
in a finer, subshellwise breakdown of spin populations, where pA = q,A q,A

3 Because the theory is symmetric with respect to reverting the z direction, energy must
be an even function of magnetization. This symmetry is easily broken spontaneously by a
spin-polarized zeroth-order state of the system, however, but we do not go into details now.
23

and pA,l = q,A,l q,A,l are atomic and subshellwise Mulliken spin populations4 .
The W constants are the appropriate semi-empirical energy derivatives with
respect to the ps (or, as an alternative approach, W s are the DFTB approximations of the spin-dependent part of two-electron integrals in a second-order
expansion of total energy, just like U s follow a similar concept in the spinindependent part). The above energy expressions are parallel to the on-site
part of (1.64) in the atomwise and (1.73) in the subshellwise case, respectively.
WA,ll spin constants are calculated as a spin-dependent part of a perturbative response to certain electron population changes, quite similar to (1.72)
with chemical hardness
l,
1 l,
.
(1.80)
WA,ll =
2 nl
nl
With the above DFTB spin-spin interaction energy, e.g. magnetization properties of iron clusters can be computed quite accurately with DFTB [21]. The
need of a DFTB code calculating these magnetic properties accurately lead to
a subshellwise calculation of Mulliken populations and the respective almost
complete subshellwise implementation of atomic hardness mentioned in 1.2.2.
The simple collinear spin model can be generalized to non-collinear spins,
it involves extending wavefunctions to two-component non-relativistic spinors,
and substituting (1.79) with
1
WA,ll pA,l pA,l
2
(1.81)
A,l,l
where Mulliken charges first become 2 2 matrices with spinor indices (these
come from the orbital-times-orbital part of Mulliken density operator) and a
decomposition of them to a q + px x + py y + pz z linear combination of Pauli
matrices leads to the Mulliken charge itself along with spin populations in three
spatial dimensions. For an accurate derivation see [24].
Spin-orbit coupling
In the non-collinear treatment of spin effects in the DFTB+ code [2], the spin 2 2 (spinor)
orbit coupling Hamiltonian can be easily written as the L
is the angular momentum operator and is the vector of Pauli
matrix where L
matrices. For a bit deeper insight, see also [24].
Possibly, the non-collinear spin-orbit coupling scheme can also be simplified
to a collinear one by collapsing the interaction part into the z direction. This
leaves
m q p
(1.82)
that is a straightforward realization of collinear spin-orbit coupling in the DFTB

philosophy (m is the z-direction magnetic quantum number of ). Collinear
spin-orbit coupling does not seem to be implemented in any DFTB code, but
4 We use p instead of p as the spin population increment with respect to p(0) because
DFTB systems used to be built of spin-unpolarized atoms as the zeroth step of calculations
(this is in a strong correlation with the fact reflected in the previous footnote, that the treatment of spin-polarized zeroth-order states has a limited importance in DFTB traditionally).
24

it might be a good idea to try this regime as a simpler alternative to the noncollinear case.
Spin-orbit coupling constants are not calculated by atomic DFT codes related to DFTB+ in the way as Hubbard parameters or spin constants, but are
imported from other calculations or estimates.
1.2.4
The current practice of parametrization
Parametrization with pair potentials

According to equation (1.33), the repulsive energy is broken down to pairwise
potentials:
1
Erep ({ratoms }) =
UZA ZB (rAB )
(1.83)
2
A=B
where A and B both run over the atoms in the system and ZA ZB indicates the
type of atom pair AB. Naturally, UZA ZB (r) = UZB ZA (r).
The parametrization process optimizes these UZA ZB (r) pair potentials to
cover the difference between the reference energies of certain fit systems and the
corresponding electronic DFTB energy. The reference energies may be taken
from experimental data or ab initio calculations. We prefer the latter as it
allows a versatile reference data generation. The best parametrization can be
viewed as the one where the set of pair potentials minimizes the error:
2
UZA ZB (rAB ) (Eref EKS ) = min.

2
R=
(Eref EDFTB ) =
A=B
(1.84)
Due to the pair approximation of Erep and the otherwise approximative
nature of DFTB, parametrizations lack universal transferability, but as the cases
of successful parametrizations show, the validity of a good parameter set can
extend to a wide range of problems.
Hand-made repulsive potentials
In its usual course, parametrization for a bond type (e.g. the carbon-carbon
bond) begins with stretching one bond of that kind in an appropriate molecule,
as the simplest case, and creating a UCC (r) curve based on the energy difference
between the DFT reference and DFTB:
UCC (r) = EDFT (r) EKS (r) + const
(1.85)
with r being the length of the stretched bond. The constant term covers the
limit of EDFT (r) EKS (r) at r (this limit contains e. g. the repulsive
contributions from non-varying bonds, that may not be known in detail at all) in
order to ensure a zero limit for UCC (r), r . Of course, one chooses stretched
molecules so that stretching affects only one bond (or maybe several bonds, but
in a totally equivalent way), all the other pairs of atoms with changing distances
remain outside the ranges of their respective repulsives.
One can construct a reasonable curve for a given atom pair type by merging
curve sections created for different molecules which represent different chemical
25

bonds between the considered elements. For example, a carbon-carbon pair
potential can be constructed by taking the sections near 1.2
A, 1.34
A and
1.54
A from ethyne, ethene and ethane, respectively, in order to take single,
double and triple carbon bonds into account. The resulting compound curves
can then be heuristically improved by comparing DFTB results on some test
systems to DFT data, and fine-tuning them by hand. Unfortunately, the finetuning involves a tremendous amount of human work, making the fast extension
of a given set or creating a new set from scratch rather difficult.
26
Chapter 2
Improvements to handling
repulsive energy
2.1
Automatization of the parametrization process
In order to reduce the work involved in creating repulsive potentials, we propose

an automatic algorithm based on least-squares fitting of repulsive potentials
to reference energy values. During our early automatic fitting attempts we
had begun experimenting with genetic algorithms based on the early steps of
Knaup et al. [20], but the simpler least-squares fits, which were being used
originally only for checking purposes, turned out to be easier to handle and far
less resource-hungry while delivering results of the same or even better quality.
Knowing that the human work needed for parametrization consists of not only
the fitting itself, the process to be described below is not limited to the bare
fitting of repulsive potentials UZZ (r). It also helps selecting and producing fit
systems and fit data, tuning the priorities of different systems or properties etc.,
making the whole parametrization process largely automatic.
2.1.1
Least-squares fitting of repulsive potentials
In order to make a least-squares fitting for the pairwise repulsive potentials

possible, we express them in terms of some arbitrary basis functions as
UZZ (rAB ) =
ZZ , fZZ , (rAB ),
(2.1)
where ZZ is the type of atomic pair AB. Substituting this into the pair potential structure of Erep from equation (1.83), the total repulsive energy for a
given system becomes a linear combination
Erep =
ZZ , XZZ ,
Z,Z ,
27
(2.2)
CHAPTER 2. IMPROVEMENTS TO HANDLING REPULSIVE ENERGY

of the structure-dependent quantities
XZZ , =
1
2
fZZ , (rAB ).
(2.3)
ZA ZB =ZZ
A=B
The sum runs over all possible atom pairs where the pair AB belongs to pair
type ZZ .
Using the above, the best ZZ , coefficients may be easily approximated by
a least-squares fit to energy values of several different distortions of a chemical
system as a function of the changing X values. Due to the linearity of the energy
as a function of XZZ , s, this fitting is a multidimensional linear regression.
Running over a sequence of distortions denoted by s of the same system (we
will call this sequence a fit path and the distortions fit steps), the least-squares
fit minimizes the overall error
R=
all
steps
2
(s)
(s)
(s)
Erep
(Eref EKS )
all
steps
(s)
ZZ , XZZ ,
(s)
(Eref
(s)
EKS )
(2.4)
Z,Z ,
where EKS is the KohnSham energy of DFTB, that is the KohnSham energy
of valence electrons. With expression (2.2) of the total repulsive energy, the
stationary condition
R[ZZ , ] = 0
(2.5)
ZZ ,
of the above error leads to a matrix expression of the coefficients ZZ , :
1
A = XXT
XE.
(2.6)
The matrices E, X and A are constructed from the above energies, X structural
quantities and s in the following way:
(1)
(1)
Eref EKS
(2)
(2)
E EKS
E=
ref
..
.
(2)
XCH,1
. . . . . . . . . . .
(1)
(2.7a)
(2)
X
HH,1
(1)
XHH,2
..
.
X=
X (1)
CH,1
(1)
XCH,2
..
.
XHH,1
(2)
XHH,2
28
(2.7b)
2.1. AUTOMATED PARAMETRIZATION
HH,1
HH,2
..
.
A=
CH,1
CH,2
..
.
(2.7c)
where, as an example, we assumed that the enumeration of the investigated

atomic pairs begins with HH and contains CH.
As an example, a fit path could be built from a propane molecule with its
middle carbon atom being shifted by 40 small random displacements around
its equilibrium position. Each movement as well as the original configuration
is a different fit step. The energy and structure data of these 41 steps would
then give enough input to fit UCC (r) and UCH (r)1 provided the number of
independent fitting parameters ZZ , is well below 40, i.e. in this specific case
the number of basis functions used to describe one pairwise repulsive is well
below 20. This criterion is normally fulfilled, but if not, increasing the amount
of steps is always a straightforward remedy.
2.1.2
Fitting to multiple fit system types and objectives
An important expectation towards repulsive potentials is their transferability to

a broad range of different systems. Usually this requires compromises; transferability can be reached via a trade-off between individual systems. Our automatic
parametrization scheme enables the optimization of this trade-off by making the
fit possible on multiple test systems (multiple fit paths) at the same time. Staying with the example of the C-C and C-H repulsive fitting, by taking several
different carbohydrogen molecules and distorting them, one can generate several
molecular fit paths for the fit. Additionally, taking bulk diamond (with various
deformations) as an additional fit path, one can tune transferability towards the
description of crystalline systems as well.
The goal of fit becomes the minimization of overall error along all fit paths,
modifying (2.4) to
R=
all
paths
p
2
(ps)
(ps)
(ps)
Erep
(Eref EKS )
(2.8)
sp
with p enumerating the paths and Erep written as a function of s and Xes
within each path in the same way as the one-path case. It should be noted
here that the number of ZZ , parameters does not depend on the number of
fit paths (nor on the number of steps in the individual fit paths). Its value is
determined only by the choice of the basis functions to describe the repulsive
interactions in question.
The E column vector of target repulsive energies for the multiple-path fitting
1 Displacements of a carbon atom in a carbohydrogen molecule alter the CC and CH
bonds (but not the HH bonds), and thus enables fitting to these repulsive potentials only.
29

is created by putting the E(p) vectors for each path on the top of each other:
(11)
(11)
Eref EKS
(12)
(12)
(1) Eref EKS
E
..
E(2)
E=
= (21)
(21)
..
EKS
Eref
.
E (22) E (22)
ref
KS
..
.
(2.9a)
while the geometry matrix is created by putting single-path geometry matrices

together in a similar way:
X=
X(1)
X(2)
..
.
(2.9b)
This all gives back matrix equation (2.6) on A for the multiple-path fit.
2.1.3
Fitting to objectives other than energy
The scheme proposed here is not restricted to obtain repulsive potentials by

fitting to energy differences between ab initio DFT and DFTB calculations.
One can naturally extend it to interatomic forces or even Hessians as targets.
This way one gains the possibility of not just choosing the transferability range
by selecting various systems for the fitting procedure, but also of being able
to select the properties which are required to be transferable to the maximum
possible amount over those systems. Furthermore, by using energy differences
between successive steps as a target instead of the absolute energies, fitting on
molecular dynamics (MD) trajectories is also made efficient. Details for these
three target extensions (force, Hessian and energy difference), as well as an
approximation to the Hessian fit, and fitting to the stress tensor are given in
this section. Of the objectives introduced here, three, the force, energy difference
and force difference objectives are fully implemented and thoroughly tested in
our automaton.
Fitting to forces
The force objective from the repulsive interaction is repulsive force FA acting
on atom A projected onto a unit vector (a direction) u
FA,u = FA u =
ZA ZB , f ZA ZB , (rAB )
B=A,
30
rAB
u.
rAB
(2.10)

(Note that B = A in the summation goes over B values that do not equal to
A; A is not a running index.) This can be, similarly to the energy expression
(2.2), decomposed into a linear combination
(i,u)
FA,u =
ZZ , XZZ ,
(2.11)
ZZ ,
with coefficients ZZ , , built of geometry-dependent factors
rAB
(A,u)
XZZ , =
f ZZ , (rAB )
u
rAB
(2.12)
ZA ZB =ZZ
B=A
containing the first derivatives of basis functions fZZ

, (r). Because the coefficients in these force components are the same as the ones used for the energy
fitting, fitting to energies and forces can be unified when both are required. If
(Fref FKS )A,u takes the place of Eref EKS and the above new Xes are used
as independent variables, fitting to force components can be simply regarded as
a set of additional new fit paths (of course, Fref and FKS are the corresponding
derivatives of Eref and EKS ). The matrices E and X can then be extended in
the same way as in equation (2.9a) and (2.9b).
Fitting to MD trajectory energies

A problem that often compromises fitting to MD trajectories (or to large molecules
where only a tiny part is distorted) is the fact that equilibrium bond lengths are
heavily overweighted by their overwhelming presence in the sample fit paths.
This can make efficient fitting to ranges of bond lengths other than the covalent
equilibrium impossible with the original energy target described above.
A remedy of this problem can be found by fitting to energy differences between subsequent fit steps instead of energies of each fit step. As
(s)
(s+1)
(s)
(s+1)
(s)
Erep
= Erep
Erep
=
ZZ , XZZ , XZZ ,
(2.13)
ZZ ,
is a linear combination of structural quantities of type X (s+1) X (s) , it is a valid

target in our least-squares fit scheme. This modified energy target, however,
contains virtually nothing arising from those bonds which do not change over
the fit path, thus the overweighting of unchanged bonds is avoided. Even if
the MD trajectory does not contain many fixed bonds (e.g. with fluids), the
distribution of bond lengths is almost always concentrated around equilibrium
values and therefore the differences in (2.13) will contain sample data from
equilibrium bond lengths with a naturally very low weight. Of course, if fitting
to absolute energy values at molecular equilibrium bond lengths is required, it
can be brought back by an appropriate weighting between the original energy
objective and the current one, or by defining additional molecular fit paths.
As an alternative use, the fit target based on energy differences can also
be used in cases where retrieving force data from a DFT reference is for some
reasons problematic or meaningless (e.g. with symmetric distortions of symmetric systems atomwise total forces on some symmetrically positioned atoms are
constant zero). Using small distortion steps and the energy difference fit target,
one automatically obtains a fit mimicking the fit on certain force or stress tensor
components.
31

Fitting to Hessians
Similar to the forces, also the repulsive contribution to the Hessian matrix of a
chemical system can be projected onto unit vectors u and v (these unit vectors
can be regarded as virtual displacements of atoms). When both u and v are on
the same Ath atom (i.e. we examine the Ath on-site 3 3 hyperdiagonal block
of the 3N 3N collective molecular Hessian),
HA,uv = uHv =
um
B=A,mn
2 fZA ZB , (rAB )
vn
xA,m xA,n
B=A,mn,
1 f
1 f
1 2 f
(u r)(v r) +
(u v).
ZA ZB ,
2
2
3
r r
r r
r r
ZZ ZA ZB =ZZ
B=A,
2 UZA ZB (rAB )
vn
xA,m xA,n
ZA ZB , um
(2.14)
(with rAB = (xB,1 xA,1 , xB,2 xA,2 , xB,3 xA,3 ) at the beginning and f =
fZA ZB , (rAB ), r = rAB , r = rAB in the last step). So, with
(A,uv)
XZZ ,
ZA ZB =ZZ
B=A
1 f
1 2 f
3
2
2
r r
r r
(u r)(v r) +
1 f
(u v), (2.15)
r r
the usual linear combinations can be written again
(A,uv)
HA,uv =
ZZ , XZZ , .
(2.16)
ZZ ,
With the E vector composed of (Href HKS )A,uv s and the X matrix composed
of the above Xes, the fitting of the Hessian can be included as an additional
path into the fitting scheme.
When u and v are on the Ath and B th atoms, respectively,
HAB,uv = uHv = um
2 UZA ZB (rAB )
2 UZA ZB (rAB )
vn = um
vn .
xA,m xB,n
xA,m xA,n
(2.17)
Therefore a similar construction applies to off-site Hessian parts, but with the
opposite sign and without the summation over B.
As a linear combination of the above u and v atomic virtual displacements,
every collective distortion of a molecule can be constructed. This knowledge
can be used to fit to Hessians of DFT reference algorithms that give no detailed
Hessian matrix but only vibrational modes and frequencies in their output. If
e is a (3N -component) collective eigenmode of the molecular Hessian with
frequency,
2 Me = He = (HKS + Hrep )e
(2.18)
where M is the diagonal mass matrix. The vector of equations contained in
2 Me HKS e = Hrep e
32
(2.19)

can then be used as a new fit path with the left hand side as a vector of E
values and the right hand side as the usual linear combinations coming from the
repulsives and using s as coefficients2 . Note that the last equation contains
explicit Hessian data from DFTB only.
Fitting to force differences as an easy approximation for certain parts
of Hessian
Using differences of target quantities as new targets is also possible with repulsive forces. On the analogy of (2.13), differences of Frep,A,u = (Fref FKS )A,u
repulsive forces can be approximated by differences of right hand sides of (2.10).
Thus
(s+1)
(s)
(s+1;i,u)
(s;i,u)
Frep,A,u Frep,A,u =
ZZ , XZZ ,
XZZ ,
(2.20)
ZZ ,
defines a set of new targets for steps s and s + 1.

The above target can be used efficiently in two ways. First, it is an easyto-implement double (some sort of numerical derivative) of the Hessian target.
With dimers, it is able even to fit vibrational frequencies practically at an arbitrary precision. Second, in certain cases it makes a fit to forces more stable
when it is used along the vanilla force target.
Fitting to the stress tensor
The repulsive part of the stress tensor in periodical systems is calculated as
mn =
rep
1 E
1
=
V mn
V
=
FAB,m rAB,n
Aa cell
B
1
V
Aa cell
B
1
V
1 U (rAB )
rAB,m rAB,n
rAB rAB
ZA ZB , f (rAB )
Aa cell
B
rAB,m rAB,n
rAB
(2.21)
rep is the cellwise

where mn is the strain tensor, V is the unit cell volume, E
repulsive energy, rAB,m is a component of the relative position vector rAB from
the Ath atom to the B th and rAB is the length of it. A double projection of
mn onto unit vectors u and v can be written
(uv)
uv =
mn um vn =
ZZ , XZZ ,
(2.22)
ZZ ,
mn
2 An important issue is whether we must compare DFT-equilibrium Hessians to DFTBequilibrium Hessians or we must compare Hessians of the very same geometry (practically
the DFT equilibrium geometry, as DFT calculators tend to compute (real) eigenvalues and
eigenmodes instead of outputting raw Hessian matrices). From a theoretical point of view,
the latter comparison is more valid while from a semiempirical (practical) point of view, the
former comparison is much more justified.
33

if our structural quantities are
(uv)
XZZ , =
f (rAB )
ZA ZB =ZZ
Aa cell
B
(rAB u)(rAB v)
.
rAB
(2.23)
So with the above Xes, uv can be another valid target of our repulsive fitting
algorithm.
2.1.4
Weighting of fit targets
In the formalism described so far, every fit step contributes to the R overall
error with the same weight. As this may not always be the desired behaviour,
we allow each step in each path to have an individual weight for its contribution
to the total error. If the fit is done for multiple physical properties (e.g. energies
and forces) each property can also be weighted differently.
The weighting issues come to play mainly in two areas. First, one typically
would overweight near-equilibrium geometries to assure a higher precision at
near-equilibrium bond lengths at the cost of less precise description for strongly
distorted geometries. Furthermore, weighting becomes a key issue when multiple physical properties are invoked into the fit, since the numerical values of the
differences in the various properties (energy, force, etc.) must be converted to
the same scale. This requires some experimenting but it offers the possibility of
balancing the performance of repulsives for various physical quantities. For example, heavy weights for forces are usually necessary when the fitted repulsives
give poor results with geometry optimization otherwise.
2.1.5
Basis function shapes
We have experimented with several different f (r) basis functions for the repulsive fitting. The splines used in most of the current DFTB implementations
turned out to be inappropriate for a fitting procedure as they tend to give
very oscillatory behaviour. As a straightforward alternative, we decided for the
(2.24) cut-off polynomials first. They were used in the earliest DFTB implementation [31] and still retain popularity with making parametrization by hand,
since most of the parametrizations up to now have been able to be done with
them. The zeroth and first-degree terms are omitted from such a polynomial to
ensure a smooth decay at its cutoff distance r0 :
(r r0 )
if r < r0 and 2,
f (r) =
(2.24)
0
otherwise.
This representation showed to be successful at the relatively easy hydrocarbon
parametrization, among other examples.
To emulate spline-like behaviour in our scheme, we also tested bases containing the above cut-off polynomials, but having no universal cutoff value (these
bases can be regarded as sums of multiple single-cutoff bases). Bases with two
cutoff values are very efficient at improving polynomial repulsives, while more
cutoff values bring up the oscillatory nature of splines. Less successful but
still noteworthy examples of spline-like bases are wavelet bases, which we also
34

probed. It must be noted here that experiments have also begun with regularization schemes that refrain the oscillatory nature of above bases. Although
these experiments have not been going for a long time so far, regularization has
found another interesting usage: the level of overbinding needed may be reduced
by suppressing higher than first derivatives at the Sprungschanze region (see
2.3.2 and Figure 2.2), and thus making the resulting unphysical minimum in
the dissociation curve (see Figure 2.3) shallower or even away.
Another important basis was the family of exponential functions. ea r

( = 1, 2, 3, . . .) and their linear combinations, which seem to be a very natural
choice for a repulsive function basis. These exponential functions proved to be
a successful basis for our fittings with Ti and Zn. In these cases a fairly tiny
set of exponential basis functions (one to three of them) was quite enough to fit
remarkably good parameter sets.
2.1.6
Further automation in the parametrization workflow
Besides the automatized fitting process itself, there are three subprocesses of the
parametrization workflow in which our program substantially lowers the human
contribution.
The path building methods mentioned so far and some others are implemented to be executed automatically. They include bond stretchings,
displacing atoms, uniform volume changes, linear interpolations between
two configurations, and using predefined paths (e.g. MD trajectories or
reaction paths).
Instead of using fixed sets of metaparameters (input parameters determining the parametrization itself) for the fitting process, batch fits can
use intervals of them, e.g. instead of setting the longest chlorine-oxygen
polynomial cutoff to a fixed value of 2.5
angstroms, we can make a probethrough from 2.3 up to 2.8 by steps of 0.1
A. Scanning over all of these
values in all of these intervals in every combination spares a lot of try-andfail cycles for the user. At the end of the batch run, the set with lowest
total error on the targets (as defined in (2.8)) is picked as the fittest solution.
A module for defining test systems is built into our program too. It
tests the energetical and geometrical performance of fitted repulsives on
the specified test systems. This way it can give a first-glance feedback
about the performance of the fitted repulsive set on systems that were not
necessarily fit systems.
2.1.7
Optimizing electronic parameters?
Batch fitting of parameters outlined in the previous section is not only able
to run over large sets of metaparameters, but it is also capable of brute-force
choosing from a set of different electronic parameters. This means changing
electronic parameters (the pseudoatomic compression radii) defined in 1.2.1 step
by step, just like metaparameters of the fitting process, and then picking the
one that allowed the best fit. The best fit means the same as in the previous
section: the smallest accumulated error of repulsive targets defined in (2.8).
35

Of course, switching electronic parameters between batch steps needs more
than switching metaparameters. When electronic parameters change in the
batch, a full regeneration of SlaterKoster integral tables is necessary before
going on to the fitting step itself. This and the direct, batch-like probe-through
algorithm of comparing different electronic parameters instead of something
more sophisticated algorithm earns the brute-force title for the way of electronic parameter fitting proposed here.
Up to now, very little has been produced with the above electronic parameter fitting strategy in the traditional DFTB parametrization framework. As
results show so far, influence of electronic parameters on the quality of fitted
repulsives is less than expected, and this is in parallel with the contemporary
practice: electronic parameters are tuned almost always only to get better electronic properties, e.g. band structures of crystals resembling to reference ones,
improving energetics is handled by repulsive fitting later. Three changes in the
machinery can make electronic parameters more influencing, however. First,
if improvements of SCC-DFTB proposed in this thesis makes DFTB so precise that differences between repulsives on top of different electronics become
harsher, repulsives fitted with different electronics will become more distinguishable. Second, if the ab-initio-like calculation of (first-guess) repulsive potentials
prevails (see section 2.4), electronic parameters remain the only adjustable parameters of DFTB at the first, ab-initio-like level. This will make differences
in electronics more important than ever before. In the current state of this
project, optimization of electronic parameters seems to play a promising role,
see 2.4.3. Third, if a good method can be found to make electronic properties
(band structures, etc.) of DFTB-calculated systems an objective of the fitting
algorithm, the proposed brute-force electronic fitting will be able to optimize
these electronic properties with choosing optimal electronic parameters.
2.2
Applications of the parametrizer
Since the molecular reference ab initio calculations in the hand-made sets were
mostly done using the Gaussian [17] code, we also used it as a reference for
molecular systems. For the periodic systems, however, we have found the Siesta
[34] code far more stable (less prone to convergence failures) in our automatic
fitting environment, where distorted systems far from the equilibrium must be
calculated very often. Apart from stability issues, this choice is also a good
cross-calculator and cross-methodology (even between different xc functionals
in DFT references) consistency check of our algorithm and in general for the
DFTB parametrization philosophy. As it will be seen from the results, this
mixing of DFT references did not pose any problem.
The DFT calculations of the last application example, the half halogenorganic set were not made with the above ad hoc machinery, but with the
NWChem [36] code, the B3LYP exchange-correlation functional and a cc-pVTZ
basis. This setting was in line with the other part of the halogen-organic project,
carried by other authors.
The DFTB calculations were carried out using the DFTB+ package [2].
36
2.2. APPLICATIONS OF THE PARAMETRIZER
2.2.1
Carbohydrogens
The carbohydrogen case is a relatively easy case of parametrization in the sense

that quite useful parameter sets can be fitted to it even with a small effort.
Fitting to DFT references with the PBE [30] exchange-correlation functional,
the resulting parameter sets produce, according to our experiences, geometrical
errors typically within a few 102
A and atomization energy errors in the range
of a few 102 a. u. This quality, which is almost comparable with the hand
made mio set [12], is pretty easy to reach at an automatic fit with a non-trivial
handful of fit systems and a couple of hours working with them.
Adding more configurations, the results can be further improved. In order
to demonstrate the automatism in our procedure, we give the instructions used
to generate those configurations:
a methane molecule with its central carbon atom randomly displaced on
five shells within a sphere3 of diameter 0.75
A;
an ethane molecule with one carbon atom displaced on ten equidistant
shells within a 0.75
A sphere;
a butane molecule with its 1-2 carbon-carbon bond stretched in fifteen
0.1
A steps, from a shortening of 0.6
A to a lengthening of 0.9
A;
a benzene ring with one of its carbon atoms displaced on five equidistant
shells within a 0.75
A diameter sphere;
an ethene molecule with one carbon atom displaced on five shells in a
0.75
A diameter sphere;
a series of random displacements similar to the above with an ethyne
molecule;
a hydrogen molecule with its only bond shortened in eight and lengthened
in twelve 0.025
A steps;
an isobutane molecule with its central carbon atom displaced in a 1
A
diameter sphere.
As the mio set, the basis of comparison, was fitted to calculations with the
B3LYP xc functional and the 6-31G* basis, we also used this as a reference. The
force objective had a weight of three while energy had a weight of one, and each
path had its near-equilibrium steps (at most three steps away from equilibrium)
overweighted by five. For the diamond test system we used the CRYSTAL2003
code [32] (because of the problems with Gaussian mentioned above) with a 621G* [11] basis set and a k-space mesh of an 8 8 8 MonkhorstPack scheme.
During the fitting process the automaton was allowed to sweep over the
following metaparameters to search for the best fit:
3 Displaced on n shells in a sphere means that the atom is dislocated with a random
vector on n spherical shells around its original position; the n equidistant shells are defined
within the largest sphere, from radius 0 up to the largest radius. The random vectors are
generated isotropically, one with length zero, and at least four on each non-trivial shell. This
way, a path with an atom jumping around n times contains at least 4n steps, plus one for the
original configuration.
37

The cutoff of CC: 2.0
A 2.3
A,
The cutoff of CH: 1.3
A 2.1
A,
The cutoff of HH: 1.3
A 2.1
A,
The highest degree of polynomials: 10 12.
The best fit was achieved with values of 2.3
A, 2.1
A, 1.3
A and 11 for the
above metaparameters, respectively. For the sake of smoothness, the polynomials contained a minimal power of 4. Table 2.1 shows the performance of the
resulting repulsives in comparison to the mio set (columns mio and hom) on
the respective equilibrium structures. As our method aims at not only describing equilibrium properties as close to ab initio results as possible, but also to
provide a reasonable accuracy when dealing with structures out of equilibrium,
we calculated also the energy errors over all non-equilibrium configurations in
the fit paths. They remained generally within the error of several 102 hartrees
compared to DFT reference except some of the extremely distorted geometries.
Table 2.1: Molecular and crystalline data calculated with the three
parameter sets (the mio set [12] and the two automatically fitted
ones) compared to reference values. E means atomization energy
error relative to the reference in kcal/mol and ZZ atom pairs denote distances of the appropriate neighbouring atoms in
A. The
column hom contains a fit without dissociation energy correction,
inhom contains a fit with it. Values in parentheses indicate errors for the set with dissociation energy correction when used in a
DFTB implementation without this correction scheme. Italicized
names denote systems that were fit systems too, the other molecules
are the rest of the carbohydrogen part of the G2 [8] test set.
property
methane
E
CH
ethane
E
CC
CH
ethene
E
C=C
CH
ethyne
E
CC
CH
benzene
E
CC
CH
reference
mio
hom
inhom
0
1.093
7.3
1.089
-2.5
1.094
-0.1 (52.9)
1.080
0
1.531
1.096
17.7
1.501
1.098
-1.0
1.535
1.102
0.1 (94.8)
1.516
1.088
0
1.331
1.087
14.6
1.327
1.094
-2.4
1.327
1.099
-3.6 (68.6)
1.326
1.084
0
1.205
1.067
21.7
1.203
1.075
10.2
1.200
1.080
-3.5 (53.1)
1.204
1.066
0
1.397
1.087
52.9
1.396
1.098
-1.7
1.405
1.104
0.8 (170.8)
1.397
1.090
38
butane
E
(1, 2) CC
(2, 3) CC
(1) CH
isobutane
E
CC
CH
diamond
CC
cyclobutane
E
CC
CH
isobutene
E
CC
C=C
CH (in CH3 )
CH (in CH2 )
bicyclobutane
E
CC (edge)
CC (middle)
CH (in CH2 )
CH (in CH)
cyclobutene
E
CC
C=C
CH (in CH2 )
CH (in CH)
cyclopropane
E
CC
CH
propane
E
CC
CH (end)
CH (middle)
cyclopropene
E
CC
C=C
CH (opposite to C=C)
CH (neighbour to C=C)
0
1.547
1.536
1.097
38.7
1.519
1.518
1.097
2.8
1.555
1.552
1.102
1.6 (179.8)
1.537
1.534
1.088
0
1.535
1.097
38.0
1.518
1.098
1.9
1.552
1.102
1.5 (178.7)
1.534
1.088
1.555
1.540
1.575
1.558
0
1.557
1.095
40.0
1.539
1.102
7.3
1.569
1.107
13.4 (169.2)
1.534
1.094
0
1.509
1.337
1.099
1.087
36.2
1.493
1.341
1.100
1.093
0.5
1.524
1.34
1.104
1.099
-2.7 (153.0)
1.505
1.339
1.090
1.084
0
1.510
1.900
1.112
1.095
26.9
1.464
2.003
1.195
1.066
-2.1
1.549
2.112
1.161
1.021
10.3 (143.4)
1.486
1.980
1.158
1.065
0
1.573
1.519
1.097
1.087
29.5
1.569
1.524
1.104
1.097
-1.8
1.597
1.548
1.109
1.103
7.4 (140.6)
1.538
1.493
1.097
1.089
0
1.509
1.087
18.9
1.489
1.096
-9.6
1.523
1.100
-0.2 (113.8)
1.502
1.087
0
1.532
1.097
1.099
27.7
1.509
1.098
1.107
0.2
1.544
1.102
1.110
0.0 (136.6)
1.525
1.088
1.097
0
1.508
1.295
1.095
1.080
12.2
1.495
1.319
1.107
1.090
-13.8
1.528
1.319
1.109
1.095
-7.7 (83.7)
1.508
1.318
1.096
1.081
39
spiropentane
E
CC (radial)
CC (outer)
CH
methylene-cyclopropane
E
C=C
CC (radial)
CC (outer)
CH (in CH2 )
CH (on ring)
propadiene
E
C=C
CH
1,3-butadiene
E
CC
C=C
CH (middle)
CH (end)
2-butyne
E
CC
CC
CH
propyne
E
CC
CC
CH (in CH3 )
CH (in CH)
propene
E
CC
C=C
CH (in CH3 )
CH (in CH)
CH (in CH2 )
0
1.485
1.530
1.088
29.4
1.479
1.508
1.097
-18.9
1.508
1.547
1.102
-2.0 (173)
1.488
1.524
1.088
0
1.322
1.470
1.540
1.088
1.089
25.9
1.328
1.465
1.512
1.095
1.098
-10.6
1.327
1.491
1.551
1.101
1.102
-4.5 (128.7)
1.327
1.472
1.529
1.086
1.089
0
1.307
1.088
21.2
1.312
1.096
-2.6
1.312
1.102
-7.8 (83.6)
1.312
1.087
0
1.439
1.392
1.089
1.086
49.9
1.436
1.372
1.098
1.104
16.2
1.457
1.373
1.103
1.085
10.0 (143.2)
1.441
1.370
1.089
1.095
0
1.462
1.209
1.097
38.8
1.455
1.209
1.100
6.5
1.477
1.205
1.105
-1.2 (132.0)
1.461
1.209
1.091
0
1.460
1.207
1.097
1.066
30.3
1.453
1.206
1.100
1.074
8.3
1.475
1.203
1.104
1.079
1.2 (92.6)
1.459
1.207
1.090
1.066
0
1.502
1.333
1.098
1.091
1.087
24.9
1.485
1.334
1.100
1.102
1.093
-1.5
1.517
1.334
1.105
1.106
1.098
-3.7 (110.3)
1.497
1.333
1.091
1.092
1.084
Using one-body repulsive terms

With this carbohydrogen fit, we also experimented with using one-body terms
in the repulsive energy
1
Erep ({Rnuclei }) =
UZA ZB (rAB ) +
UZA
(2.25)
2
A=B
40

(see (2.45)). One-body terms are a special case of inhomogeneous or dissociation
energy terms: they represent a fixed, geometry-independent energy part as a
sum of atomwise parts, that does not come from the linear combinations of
pairwise basis functions, and that remains the asymptotic value of Erep at the
dissociation limit. One-body energies are the only mathematically correct means
of putting any correction to dissociation energy because only a sum of oneatomic dissociation energy terms behaves like an extensive quantity, i. e. is a
additive function of stoichiometry. This fact strongly encourages investigating
their use.
As the results in Table 2.1 illustrate (column inhom), one-body terms can
slightly improve geometry results via eliminating the need of trying to set absolute atomization energy levels by the pair potential profiles. The resulting
one-body terms were UC = 0.030633H and UH = 0.017967H for C and H atoms,
respectively. The optimal cutoff distances were (determined by a similar batch)
equal to those of the homogeneous case.
To maintain compatibility with the current DFTB implementations lacking
one-body repulsive parts, we also took another way of improving results by onebody terms into account. Using them only at the fitting process but dropping
them after it retains improved geometries and reaction energies calculated with
the produced set, yet leaves the pair potential structure of the repulsive energy
built in DFTB intact (deteriorated bare atomization energy values are shown
in parentheses in the appropriate column of Table 2.1).
2.2.2
Zinc-oxides
As a further demonstration for our fitting procedure, we attempted to create a

parametrization for the Zn-O interaction. A high-quality and well-tested parameter set had been recently created manually for the zinc-organic chemistry by
Moreira et al. [27] which should serve as an etalon for our Zn-O repulsive. For
the DFT references, the same settings had been used as for the hand-made
parametrization (PBE functional, double- polarized basis, norm-conserving
TroullierMartins pseudopotentials, 8 8 8 MonkhorstPack scheme for k
sampling). The fit paths were made with distortions applied to the test systems
(see table 2.2) in addition to Zn-Zn and Zn-O dimers with very low weights.
The distortions applied to crystalline paths were uniform volume scaling and
moving a Zn atom around. We show a comparison between the performance of
the two Zn-O sets in 2.2. As fit targets, we used the two energy targets (energy
and energy differences between steps weighted by 1:10), step weighting was by
10 and 2 in the immediate and in a wider neighbourhood of equilibria. Here, the
2
3
basis of repulsives consisted of exponential functions of type ea2 r and ea3 r ,
as these shapes offered good results quickly in situations where absolute energy
targets were not heavily weighted. As it can be seen in Table 2.2, the resulting
set is superior in reproducing crystalline geometries to the hand-made one.
2.2.3
Titanium-organics
Our next test of the fitting automaton was producing a titanium-oxygen set and
extending it to a titanium-organic set. For this parametrization a good-quality
hand-made set (tiorg) has been recently created by Dolgonos et al. [10]. We
41
Table 2.2: Reference data and its comparison with previous hand-made
parametrization [27] (znorg) and the automatically created one (auto) for
Zn and ZnO crystals. The values given here refer to the equilibrium structure
of each system. E denotes the atomization energy difference with respect to
the reference in kcal/mol. Atom pairs denote distances in
A.
property
reference znorg auto
Zn hcp
E (per Zn2 )
0
115.5
94.1
ZnZn (#1)
2.523
2.796 2.433
ZnZn (#2)
2.886
2.864 2.931
ZnZn (#3)
3.831
4.051 3.788
ZnZn (#4)
4.591
4.872 4.524
ZnO zincblende
E (per ZnO)
0
22.3
-1.1
ZnO
2.005
2.015 2.011
ZnZn
3.274
3.290 3.281
ZnO wurtzite
E (per Zn2 O2 )
0
46.5
-0.7
ZnO
2.017
2.015 2.018
ZnO
2.037
2.014 2.004
used the same reference structures and ab initio reference data (various molecular systems calculated with the B3LYP functional and with mixed SDD+ basis
set) augmented with crystalline reference systems. For the reference calculations
of the periodic systems, the PBE functional, double- plus polarized basis functions and norm-conserving TroullierMartins pseudopotentials had been used.
K-point sampling was set to an 8 8 8 MonkhorstPack scheme with both
Siesta and DFTB in this fit session.
In order to fit repulsive potentials for the Ti-Ti and Ti-O interactions, we
used a fit set including a titanium dimer (with a very low weight), a TiO2
molecule, a planar Ti2 O2 molecule, a tribridged Ti2 O4 , the bulk hcp titanium,
and the bulk anatase and rutile forms of TiO2 . The molecular fit paths were
created by stretching bonds and displacing titanium atoms while the crystalline
paths were created by uniformly changing the volume of the crystal lattices and
by using crystals with displaced titanium atoms. We used both energy and force
targets (generally weighted 1:2) in the fit.
In a fit session of a few days, we were able to produce a set of Ti-Ti and
Ti-O repulsive potentials which reproduce energy and geometrical data in the
same quality as the reference hand-made set. A detailed comparison is given in
2
Table 2.3. These results were obtained using ea1 r and ea2 r -type exponential
functions as basis functions for the fit because this analytical basis gave very
good results quickly with the Ti and Ti-O chemistry.
After creating the repulsives for the TiTi and TiO interactions, we extended the set to a complete Ti-organic set, still using exponential basis functions for expressing the repulsive potentials. The extension turned out to be
more difficult than expected, mainly due to the sudden cutoff in the hand-made
Ti-Ti repulsive giving a very stable 3
A TiTi distance in hcp titanium, titanium nitride and titanium carbide. This feature (shown in Figure 2.1) can be
42
Table 2.3: Titanium-oxygen compound reference data and its comparison with
previous hand-made parametrization [10] (tiorg) and the automatically created one (auto). The values given here refer to the equilibrium structure of
the various systems. E means atomization energy difference with respect to
the reference in kcal/mol. Atom pairs denote neighbour distances in
A, triples
denote angles in degrees (double distance values show artificially broken symmetries of DFTB-optimized lattices). Italicized system names denote fit systems.
property
reference tiorg auto
TiO
E
0
55.0
40.2
TiO
1.586
1.592 1.586
Ti2 O2 planar
E
0
87.7
69.4
TiTi
2.198
2.355 2.092
TiO
1.857
1.891 1.866
Ti2 O2 non-planar
E
0
68.0
57.6
TiTi
2.127
2.249 2.133
TiO
1.838
1.888 1.826
Ti2 O4 #1 (dibridged with end O atoms in cis position)
E
0
49.4
81.3
TiTi
2.716
2.800 2.635
bridging TiO
1.848
1.887 1.812
end TiO
1.622
1.606 1.589
OTiTi (ending O)
126.1
124.7 123.7
Ti2 O4 #2 (dibridged with end O atoms in trans position)
E
0
145.7 82.8
TiTi
2.709
2.726 2.709
bridging TiO
1.840
1.831 1.806
end TiO
1.625
1.608 1.590
OTiTi (ending O)
123.7
122.3 122.2
Ti2 O4 #3 (tribridged with an O atom at one end)
E
0
123.3 66.4
TiTi
2.394
2.540 2.399
bridging TiO (opposite to end O)
1.763
1.801 1.742
end TiO
1.628
1.606 1.586
Ti hcp
E (per Ti)
0
-40.6
5.4
TiTi
2.900
2.993 2.915
TiO2 anatase
E (per TiO2 )
0
22.5
-4.5
shortest TiTi
3.028
2.996 3.082
1.921
1.957
shortest TiO
1.933
1.995
1.958
TiO2 rutile
E (per TiO2 )
0
19.4
1.1
shortest TiTi
3.559
3.613 3.605
1.914
1.976
shortest TiO
1.974
1.992
1.995
43
Figure 2.1: Comparison of the tiorg (dashed) and the automatically generated
(solid line) Ti-Ti repulsives in the area of the sharp cutoff of the former.
repulsive energy [a. u.]
0.2
0.1
2.5
3.5
distance []
hardly reconstructed with analytical sets. Although this peculiar shape is the
numerically most convenient way to confine the range of the Ti-Ti repulsive well
below the second-neighbour distance as well as it gives good results for various
systems, it may be an interesting question for future investigations whether it
is a precise representation of the underlying physics.
Similar to the case of the titanium-oxygen fit, we used the same molecules
(TiH4 , Ti(CH3 )4 , Ti(NH2 )4 ) as during the hand-made parametrization [10] extended by crystalline fit systems (TiN and TiC). The fit paths were similar
constructions to the titanium-oxygen case: every molecule had two fit paths
with the relevant bonds stretched and the titanium atom displaced, while the
crystalline systems had a volume change path and a path with a titanium atom
displaced. We handled the relative difficulty of Ti-C and Ti-N fitting compared
to Ti-Ti and Ti-O by lowering the relative weights of the energy targets in the
Ti-C and Ti-N case. As it can be seen from the results (Table 2.4), this resulted
in fairly good geometries at the expense of less accuracy in energies (one-body
terms as a tool for resolving the conflict between energy and geometry accuracy
were not used here).
2.2.4
Halogen-organics
The halogen-organic project part summarized here4 includes producing ZH, Z

N and ZO repulsive potentials where Z means Cl, Br and I. As it was given in
the introduction, the DFT reference side of these fittings use the NWChem [36]
4 These fits, or their possible improved versions will constitute a part of a halogen-organic
series of fittings that extend the mio organic set and the fluor-organic parameters included
in the pbc set [22] to a complete halogen-organic set up to iodine. Production of electronic
tables and the carbon-halogen and heteroatomic halogen-halogen fits are being made by Tom
as
Kubar and Michael Gaus. The whole product will be published [23] and released for general
use.
44
Table 2.4: Reference data and its comparison with previous hand-made
parametrization [10] (tiorg) and the automatically created one (auto) for
various titanium compounds. The values given here refer to the equilibrium
structure of each system. E denotes the atomization energy difference with
respect to the reference in kcal/mol. Eb indicates the binding energy between the central Ti atom and the ligands compared to the reference value in
kcal/mol. Atom pairs denote neighbour distances in
A, triples denote angles in
degree. Italicized system names refer to fit systems.
property reference tiorg auto
Ti(CH3 )4
Eb
0
64.6 180.2
TiC
2.072
2.096 2.025
Ti(CH3 )2
Eb
0
-38.8
93.4
TiC
2.038
2.096 2.025
CTiC
113.7
110.2 109.9
crystalline TiC
E
0
111
91.7
TiC
2.141
2.159 2.170
TiTi
3.024
3.047 3.067
Ti(NH2 )4
Eb
0
30.6 287.4
TiN
1.899
1.902 1.853
H3 Ti(NH2 )
Eb
0
12.0
76.2
TiN
1.846
1.898 1.837
HN=Ti=NH
Eb
0
15.5 156.1
TiN
1.707
1.703 1.671
NTiN
114.8
114.7 113.7
crystalline TiN
E
0
196.6 192.2
TiN
2.094
2.159 2.115
TiTi
2.958
3.043 2.982
Ti2 H2 (dibridged planar)
E
0
123.5
131
TiTi
1.985
1.967 2.011
TiH
1.868
1.827 1.899
TiHTi
64.2
65.2
63.9
45

code with the B3LYP xc functional5 and a cc-pVTZ basis.
In order to give a novel systematic and reconstructible treatment of the problem of dissociation levels and overbindings written in 2.3.2, we use the following
scheme here. A preliminarily fixed (natural) one-body repulsive energy is calculated for every Z atom type by halving the binding energy of the heavily
stretched Z2 dimer. Heavy stretching means here a quasidissociated state with
a distance of ten atomic covalent radii between the two atoms (and therefore
its energy with respect to the two free atoms may be considered a numerical
corrective term rather than a physical binding energy). As the above one-body
energies would give a huge overbinding, they are reduced then to 20% of their
original amount. (For the sake of preciseness, we note that the latter reduction was achieved by generating the final parameter set as a 20%80% weighted
sum of dissociation-corrected and dissociation-uncorrected sets, but a new fit
with the reduced one-body terms would have resulted in roughly the same parameters.) These reduced one-body energies (see Table 2.5) provided enough
overbinding to eliminate the unphysical second minima of dissociation curves
(see Figure 2.3).
Table 2.5: One-body repulsives used for halogenorganics. Note
that pairwise sums of them give equivalent traditional pairwise
overbinding values.
atom
H
N
O
Cl
Br
I
one-body
natural
0.0525
1.8266
2.1182
0.8363
0.6911
0.5817
rep [eV]
reduced
0.0105
0.3653
0.4236
0.1673
0.1382
0.1163
The cutoff values of polynomials building the repulsives were partly allowed
to vary during fitting. The optimized values of these metaparameters are shown
in Table 2.6. The quality of parameters is illustrated by tables 2.72.16.
5 It would be worth investigating that fitting to a B3LYP reference is much more difficult
than fitting to PBE or PBE0 (with which these fits were made more successfully, but they were
rejected for a compatibility reason with mios B3LYP reference). The most straightforward
explanation is that DFTB itself uses integral tables produced with PBE, but there may or may
not be other reasons in the background too. Cf. this with the similar note at the beginning
of 2.2.1.
46
Table 2.6: Cutoff values in the repulsives (note that close multiple values come
from the parts of blended parameter sets, while more distant values come from
multiple-cutoff fits).
pair
ClH
ClO
ClN
ClCl
BrH
BrO
BrN
BrBr
IH
IO
IN
II
cutoff(s) [
A]
1.9, 2.0
1.5, 1.95, 2.8
2.8
2.8
1.9, 2.1
1.8, 2.1, 3.0
3.0
3.0
2.4
1.9, 2.0, 2.3, 3.2
3.1
3.4
47
Table 2.7: Results with chlorine compounds. E: atomization energy in eV

(E : calculated with one-body terms on the DFTB side), doubles are bond
lengths (distinguished with bond quality if needed) in
A, triples are angles in
degrees.
property
ref.
DFTB
error
rel. error.
Cl2
E
-2.3482 -2.6858 -0.3377
0.14379
E
-2.3482 -2.3513 -0.0031
0.00134
ClCl
2.02249 2.02366 0.00117
0.00058
Cl2 O
E
-4.1138 -6.9317 -2.8179
0.68500
E
-4.1138 -6.1736 -2.0598
0.50071
OCl
1.71496 1.64923 -0.0657
-0.0383
ClOCl
112.572 113.370 0.79738
0.00708
HCl
E
-4.5256 -4.5718 -0.0462
0.01022
E
-4.5256 -4.3941 0.13154
-0.0291
ClH
1.28337 1.28587 0.00250
0.00195
HClO
E
-7.0008 -8.7040 -1.7033
0.24330
E
-7.0008 -8.1026 -1.1019
0.15739
OCl
1.70813 1.63403 -0.0741
-0.0434
HOCl
102.955 106.920 3.96446
0.03851
HClO2
E
-8.3690 -10.127 -1.7580
0.21006
E
-8.3690 -9.1019 -0.7329
0.08758
ClO
1.72867 1.66153 -0.0671
-0.0388
Cl=O
1.52336 1.5859 0.06254
0.04105
OCl=O 112.334 98.5918 -13.742
-0.1223
HClO3
E
-10.810 -13.217 -2.4067
0.22263
E
-10.810 -11.768 -0.9580
0.08862
ClO
1.73379 1.75060 0.01682
0.0097
Cl=O
1.46299 1.47015 0.00716
0.00489
OCl=O 104.575 100.075 -4.5007
-0.0430
O=Cl=O 112.646 108.001 -4.6452
-0.0412
HClO4
E
-12.761 -17.304 -4.5427
0.35599
E
-12.761 -15.431 -2.6704
0.20927
ClO
1.68073 1.66307 -0.0177
-0.0105
Cl=O
1.43962 1.41986 -0.0198
-0.0137
OCl=O 100.908 104.382 3.47457
0.03443
O=Cl=O 115.073 115.582 0.50981
0.00443
48
Table 2.8: Results with bromine compounds. E: atomization energy in eV

degrees.
property
ref.
DFTB
error
rel. error.
Br2
E
-2.1199 -2.7969 -0.6770
0.31937
E
-2.1199 -2.5204 -0.4006
0.18896
BrBr
2.31469 2.32797 0.01328
0.00574
Br2 O
E
-3.8942 -6.2094 -2.3152
0.59453
E
-3.8942 -5.5093 -1.6151
0.41475
OBr
1.85118 1.79333 -0.0578
-0.0312
BrOBr 114.842 115.369 0.52648
0.00458
HBr
E
-3.9826 -4.0572 -0.0746
0.01872
E
-3.9826 -3.9085 0.07418
-0.0186
BrH
1.42494 1.42620 0.00126
0.00089
HBrO
E
-6.9001 -8.3665 -1.4664
0.21252
E
-6.9001 -7.7941 -0.8940
0.12957
OBr
1.84191 1.77491
-0.067
-0.0364
HOCl
103.350 106.842 3.49245
0.03379
HBrO2
E
-8.4140 -9.8850 -1.4711
0.17483
E
-8.4140
-8.889
-0.4750
0.05646
BrO
1.85360 1.80672 -0.0469
-0.0253
Br=O
1.66639 1.72212 0.05572
0.03344
OBr=O 109.957 96.5473 -13.410
-0.1220
HBrO3
E
-10.938 -12.864 -1.9254
0.17603
E
-10.938 -11.444 -0.5058
0.04624
BrO
1.83944 1.86485 0.02541
0.01381
Br=O
1.61552 1.63254 0.01702
0.01053
OBr=O 103.646 98.1828 -5.4635
-0.0527
O=Br=O 110.498 106.704 -3.7942
-0.0343
HBrO4
E
-12.433 -16.220 -3.7873
0.30462
E
-12.433 -14.377 -1.9440
0.15636
BrO
1.79985 1.7691
-0.0308
-0.0171
Br=O
1.60099 1.57973 -0.0213
-0.0133
OBr=O 100.362 103.664 3.30144
0.03290
O=Br=O 115.452 116.522 1.06988
0.00927
49
Table 2.9: Results with iodine compounds. E: atomization energy in eV

degrees.
property
ref.
DFTB
error
rel. error.
I2
E
-1.8621 -2.4641 -0.6020
0.32331
E
-1.8621 -2.2314 -0.3693
0.19834
II
2.70438 2.72377 0.01939
0.00717
I2 O
E
-3.8071
-5.803
-1.9959
0.52427
E
-3.8071 -5.1467 -1.3396
0.35187
OI
2.01955 1.96827 -0.0513
-0.0254
IOI
118.785 117.962 -0.8224
-0.0069
HI
E
-3.3792 -3.2881 0.09111
-0.0270
E
-3.3792 -3.1612 0.21797
-0.0645
IH
1.62167 1.61546 -0.0062
-0.0038
HIO
E
-6.8651 -8.2057 -1.3406
0.19528
E
-6.8651 -7.6552 -0.7901
0.11509
OI
2.01066 1.94386 -0.0668
-0.0332
HOI
104.630 107.344 2.71457
0.02594
HIO2
E
-8.7164 -9.8946 -1.1782
0.13517
E
-8.7164 -8.9205 -0.2041
0.02341
IO
2.00402 1.97942 -0.0246
-0.0123
I=O
1.83408 1.89457 0.06049
0.03298
OI=O
106.608 93.5438 -13.065
-0.1225
HIO3
E
-11.650 -13.029 -1.3792
0.11839
E
-11.650 -11.631 0.01857
-0.0016
IO
1.97381 1.98529 0.01148
0.00581
I=O
1.78899 1.81169 0.02270
0.01269
OI=O
102.414 95.4984 -6.9158
-0.0675
O=I=O 108.221 106.861 -1.3595
-0.0126
HIO4
E
-13.180 -16.474 -3.2941
0.24993
E
-13.180 -14.652 -1.4726
0.11173
IO
1.94152 1.90838 -0.0331
-0.0171
I=O
1.78190 1.75823 -0.0237
-0.0133
OI=O
99.5880 102.241 2.65285
0.02664
O=I=O 116.022 117.937 1.91501
0.01651
50
Table 2.10: Results with molecules containing chlorine and nitrogen. E: atomization energy in eV (E : calculated with one-body terms on the DFTB
side), doubles are bond lengths (distinguished with bond quality if needed) in
A, triples are angles in degrees.

property
ref.
DFTB
error
rel. error.
NCl3
E
-8.8834 -9.9016 -1.0182
0.11462
E
-8.8834 -9.0345 -0.1511
0.01701
NCl
1.78150 1.7784
-0.0031
-0.0017
ClNCl 107.82 108.165 0.34498
0.0032
NHCl2
E
-11.234 -11.890 -0.6555
0.05834
E
-11.234 -11.179 0.05489
-0.0049
NCl
1.77056 1.76391 -0.0066
-0.0038
ClNCl 110.671 110.655 -0.0161
-0.0001
NH2 Cl
E
-13.486 -13.988 -0.5021
0.03723
E
-13.486 -13.434 0.05151
-0.0038
NCl
1.76834 1.75424 -0.0141
-0.0080
NOCl
E
-11.168 -13.122 -1.9543
0.17499
E
-11.168 -12.166 -0.9981
0.08937
NCl
2.00064 2.06546 0.06482
0.0324
Table 2.11: Results with molecules containing bromine and nitrogen. E:

atomization energy in eV (E : calculated with one-body terms on the DFTB

property
ref.
DFTB
error
rel. error.
NBr3
E
-8.0083 -8.8055 -0.7972
0.09955
E
-8.0083 -8.0255 -0.0172
0.00215
NBr
1.94118 1.93520 -0.0060
-0.0031
BrNBr 108.829 109.159 0.33077
0.00304
NHBr2
E
-10.663 -11.218 -0.5556
0.05210
-10.663 -10.566 0.09672
-0.0091
E
NBr
1.92381 1.91180 -0.0120
-0.0062
BrNBr 112.141 112.255 0.11468
0.00102
NH2 Br
E
-13.218 -13.680 -0.4625
0.03499
E
-13.218 -13.155 0.06208
-0.0047
NBr
1.91561 1.89528 -0.0203
-0.0106
NOBr
E
-10.871 -12.873 -2.0020
0.18416
E
-10.871 -11.946 -1.0748
0.09887
NBr
2.15763 2.16755 0.00991
0.00459
51
Table 2.12: Results with molecules containing iodine and nitrogen. E: atomization energy in eV (E : calculated with one-body terms on the DFTB

property
ref.
DFTB
error
rel. error.
NI3
E
-7.1962 -7.9194 -0.7233
0.10050
E
-7.1962 -7.2051 -0.0089
0.00123
NI
2.13754 2.12806 -0.0095
-0.0044
INI
110.727 110.914 0.18694
0.00169
NHI2
E
-10.124 -10.714 -0.5891
0.05819
E
-10.124 -10.105 0.01945
-0.0019
NI
2.11405 2.09227 -0.0218
-0.0103
INI
114.775 115.006 0.23122
0.00201
NH2 I
E
-12.958 -13.463 -0.5052
0.03899
E
-12.958 -12.960 -0.0025
0.00019
NI
2.09975 2.06760 -0.0321
-0.0153
NOI
E
-10.487 -12.566 -2.0789
0.19823
E
-10.487 -11.661 -1.1736
0.11191
NI
2.36624 2.36639 0.00015
6.2E-05
Table 2.13: Results with chlorine compound ions. E: atomization energy in

eV (E : calculated with one-body terms on the DFTB side), doubles are bond
degrees.
property
ref.
DFTB
error
rel. error.
ClO
3
E
-10.083 -11.589 -1.5062
0.14938
E
-10.083 -10.151 -0.0680
0.00674
ClO
1.51392 1.57624 0.06232
0.04117
OClO 108.094 108.355 0.26105
0.00242
ClO
4
E
-12.907 -17.487 -4.5805
0.35489
E
-12.907 -15.625 -2.7187
0.21064
ClO
1.47527 1.45669 -0.0186
-0.0126
ClO
E
-4.5841 -5.4154 -0.8313
0.18133
E
-4.5841 -4.8245 -0.2404
0.05244
ClO
1.72102 1.73016 0.00914
0.00531
ClO
2
E
-6.7945 -7.9260 -1.1316
0.16654
-6.7945 -6.9115 -0.1171
0.01723
E
ClO
1.59790 1.66100 0.06310
0.03949
OClO 115.070 114.049 -1.0205
-0.0089
52
Table 2.14: Results with bromine compound ions. E: atomization energy in

eV (E : calculated with one-body terms on the DFTB side), doubles are bond
degrees.
property
ref.
DFTB
error
rel. error.
BrO
E
-4.5172 -5.0072 -0.4900
0.10848
E
-4.5172 -4.4454 0.07185
-0.0159
BrO
1.83584 1.89874 0.0629
0.03426
BrO
2
E
-6.7873 -7.4236 -0.6363
0.09374
E
-6.7873 -6.4381 0.34923
-0.0515
BrO
1.73371 1.80432 0.07061
0.04072
OBrO 113.866 114.923 1.05687
0.00928
BrO
3
E
-10.199 -11.000 -0.8003
0.07846
E
-10.199 -9.5905 0.60886
-0.0597
BrO
1.65927 1.72552 0.06625
0.03993
OBrO 107.415 107.868 0.45317
0.00422
BrO
4
E
-12.740 -16.136 -3.3960
0.26657
E
-12.740 -14.303 -1.5633
0.12271
BrO
1.62869 1.61748 -0.0112
-0.0069
Table 2.15: Results with iodine compound ions. E: atomization energy in eV
degrees.
property
ref.
DFTB
error
rel. error.
IO
E
-4.4664 -4.6689 -0.2025
0.04533
E
-4.4664 -4.1289 0.33751
-0.0756
IO
1.97194 2.07088 0.09893
0.05017
IO
2
E
-6.9908 -7.0284 -0.0376
0.00538
E
-6.9908 -6.0648 0.92603
-0.1325
IO
1.89059 1.98852 0.09793
0.05180
OIO
111.287 115.971 4.68429
0.04209
IO
3
E
-10.785 -10.701 0.08433
-0.0078
E
-10.785 -9.3134 1.47159
-0.1364
IO
1.82580 1.90422 0.07842
0.04295
OIO
106.672 107.914 1.24267
0.01165
IO
4
E
-13.450 -15.839 -2.3891
0.17763
E
-13.450 -14.028 -0.5782
0.04299
IO
1.80423 1.79588 -0.0084
-0.0046
OIO
109.471 109.472 0.00114
0.00001
53
NCl3
HClO4
molecule
Cl2
HCl
HClO
ref.
539.8
2945.5
3774.4
1267.8
740.8
3710.9
1297.2
1211.6
1184.5
1003.9
692.8
550.9
544.4
520.6
395.5
387.6
169.1
614.0
614.0
547.5
349.6
252.6
252.6
DFTB
547.6
2929.6
3665.3
1364.6
1038.4
3649.5
1350.0
1334.3
1175.2
1107.4
771.2
526.0
524.6
496.9
394.7
363.2
96.2
757.1
744.3
633.7
383.6
289.2
280.6
NBr3
HBrO4
molecule
Br2
HBr
HBrO
ref.
469.4 (325.3)
2330.0 (2649.0)
3812.7
1614.4
743.5
3724.8
1241.8
1030.4
1022.2
895.0
678.1
488.8
471.8
404.2
305.0
302.9
212.4
516.0
516.0
418.0
139.1
DFTB
298.3
2615.6
3690.7
1282.1
898.8
3683.3
1092.9
1045.3
1033.4
964.4
708.4
366.8
361.5
343.9
292.3
272.5
101.7
667.8
652.8
493.3
240.7
170.9
166.4
NI3
HIO4
molecule
I2
HI
HIO
ref.
218.1
2314.5
3803.2
1135.7
591.6
3732.9
1039.8
886.1
881.0
825.6
607.2
288.2
287.9
270.1
240.3
229.7
511.3
511.3
360.2
156.2
104.7
104.7
DFTB
202.4
2287.6
3715.5
1203.7
774.7
3737.5
1041.0
930.7
917.2
873.4
633.3
285.5
275.3
253.0
217.5
198.6
114.9
620.0
613.7
395.0
179.3
122.3
118.6
Table 2.16: Vibrational modes of molecules calculated with the reference method and the new parameters (frequencies in cm1 ). Data in
parentheses are experimental values from [1].
54
2.3. SOME THEORETICAL ASPECTS OF PARAMETRIZATION
2.3
2.3.1
Some theoretical aspects of parametrization

On the theoretical validity of parametrization
The representation of Erep in (1.33) as a sum of quasiclassical multibody potentials is not only a matter of convenience. Its validity is strictly underpinned
by the fact that Erep does not depend on molecular charge fluctuations in the
leading order. This means that Erep is a functional of the superposed (0) electronic density of atoms (or the corresponding P (0) density operator) up to a
satisfactory precision, therefore the repulsive energy of a certain chemical system is influenced by only its geometry (i.e. how (0) is composed from atomic
densities), not its particular electronic state. This enables using transferable
quasiclassical repulsive potentials.
Erep (P ) in general. Generally, the E total energy of a chemical system can
be approximated with its Taylor expansion by P , which we investigate here up
to the second order:

2 E
1
E
(0)
P .
(2.26)
+ Tr P
E = E + Tr P
P 0
2
P2 0
The KohnSham energy of the same system6 in DFT is then given by
E
EKS = Tr P
.
P
(2.27)
In the first-order machinery of non-SCC DFTB, the above KohnSham energy is approximated by
(1)
EKS
E
(0)
(0)
(1)
= Tr[P HKS ] = Tr (P + P)
P 0

E
(0) E
= Tr P
+ Tr P
,
P 0
P 0
(2.28)
(1)
and this is supplemented by the repulsive energy to give E = EKS +Erep . Finally
(1)
Erep
=E
(0)
1
2 E
(0) E
+
Tr
P
Tr P
P
.
P 0
2
P2 0
(2.29)
Clearly, regardless of any further detail on the structure of E, we can state

that in first-order DFTB the repulsive energy does not depend on P in the first
order. (This follows from the fact that Erep is the Legendre transform of E with
respect to P , this is a direct consequence of (1.23))
6 In this section, a thorough treatment of energetics would handle the fact that DFTB
KohnSham energy belongs to the valence electrons only, but since core electrons are frozen,
(0)
their energetics is always depending on cores only. This makes a precise treatment of core
electrons unnecessary here.
55

In the second-order regime of SCC-DFTB, a theoretical, completely consistent KohnSham treatment would imply a refinement of (2.28) to
(2)
EKS
2 E
E
(0)
P
= Tr (P + P)
+
P 0
P2
0

E
(0) E
+ Tr P
= Tr P
P 0
P 0
2
2 E
(0) E
+ Tr P
P
+
Tr
P
P
(2.30)
P2 0
P2 0
by the (1.34) definition of the KohnSham energy. The repulsive energy in this
case would be
(2)
(2)
Erep
= E EKS

1
2 E
(0)
(0) E
(0) E
= E Tr P
Tr
P
P
P
,
Tr
P
P 0
P2 0
2
P2 0
(2.31)
which is clearly not independent of P in the relevant (second) order. However,

in SCC-DFTB this second-order KohnSham energy expression is not used. One
uses the
E
2 E
HKS =
+
P
(2.32)
P 0
P2 0
KohnSham Hamiltonian to calculate electronic states, but total energy is calculated as its first-order DFTB approximation plus the second-order term in its
Taylor series:
1
2 E
(1)
(2)
E = EKS + Erep + Tr P
P .
(2.33)
2
P2 0
This construction propagates the usage of charge-fluctuation-independent firstorder repulsive potentials into SCC-DFTB. Furthermore, if we take a closer look
at what Erep in SCC-DFTB is
1
2 E
(1)
(2)
(0) E
Erep = E EKS Tr P
P = Tr P
,
(2.34)
2
P2 0
P 0
we notice that the repulsive potential of SCC-DFTB is independent of P up
to the second order.
In the following, we calculate the repulsive energy part in detail when the
total energy is of the well-known structure: E = EC + Exc , the sum of Coulomb
and xc energy (we omit kinetic energy
and energy from external potential as

T
= T ). Instead of the density-operator

they are trivial cases here, e.g. Tr P
P
formalism, we use (x) = x|P |x as the basic variable.7

7 The transition from density operator formalism to spatial charge density formalism is
for a general A
opersimply writing the traces in the spatial representation. Thus Tr[P A]
it is x|P |xx|A|xdx
ator will be x|P |yy|A|xdx,

and in the case of a local A,
=
x(x)x|A|xdx.
56

The Coulomb part of total energy is
1
EC =
2
((x)VC (x))dx
1
1
(0)
(0) (x)VC (x)dx +
(0) (x)VC (x)dx
=
2
2
1
1
(0)
+
(x)VC (x)dx +
(x)VC (x)dx,
2
2
(2.35)
if we expand charge density in the first order as = (0) + . Also using that
the Coulomb potential itself is a linear function of charge density (i.e. the total
Coulomb energy is quadratic in ), we get that
(0)
(0)
(x)VC (x)dx = (x)VC (x)dx
(2.36)
and therefore
1
1
(0)
(0)
(0)
(x)VC (x)dx + (x)VC (x)dx +
(x)VC (x)dx.
EC =
2
2
(2.37)
(1)
Then the Coulomb part of EKS is
EC
(1)
(0)
EKS,C = (x)
dx
=
((0) (x) + (x))VC (x)dx.
(2.38)
(x) 0
It implies a repulsive energy of the form
(1)
(1)
Erep,C = EC EKS,C

1
1
(0)
(0) (x)VC (x)dx +
(x)
(y)dxdy,
=
2
|x y|
(2.39)
which is a functional of only in the second order in the non-SCC case. In

the SCC case, the second-order part is also dropped from the repulsive energy
part.
The exchange-correlation energy part
approximation in the same system is
of total energy in the LDA DFT
Exc =
(x)((x))dx
x
=
x
(0) (0) + (0) + (0)
d
d
1 d2
(0)
dx,
d 0
d 0
2 d2 0
(2.40)
d
=
. In a semilocal
where (0) = ((0) ). In the above mentioned LDA case, d
d
xc
GGA case of xc functional, the part denoted by d (a part of the E
operator
P
derivative) is a differential operator in this integral representation.
57

The xc part of first-order DFTB KohnSham energy is
(1)
EKS,xc
Exc
(0)
= (x)
dx = ((0) (x) + (x))Vxc
(x)dx
(x) 0
Exc
Exc
(0)
= (x)
dx + (x)
dx
(x) 0
(x) 0

(0)
(0)
(0) d
(0)
(0) d
=
+
+ +
(x) dx.
d 0
d 0
(2.41)
The corresponding part of repulsive energy is
(0)
1
(1)
(1)
(0) d
dx,
Erep,xc = Exc EKS,xc = Wxc
2
d 0
(2.42)
where Wxc is the second derivative of Exc with respect to . It may be written
2
d
(0) d
Wxc = 2
+
,
(2.43)
d 0
d2 0
where differentials of are partial derivatives in LDA, and differential operators
in GGA, just like above. Again here, the investigated non-SCC DFTB repulsive
energy part depends on in the second order only, and even the second-order
term is dropped from repulsive energy in SCC-DFTB.
2.3.2
Handling dissociation energy level
Theoretically, if parametrization were a search for a rich-enough set of multibody

potentials, there would be no need for several molecular properties as fit targets
because a proper fit to forces would give good energies and vice versa, as well
as they both would result in good vibrational modes of DFTB molecules. In
this ideal world, a fit to any derivative of energy would effectively lead to a fit
to energy profiles themselves and the practically infinite flexibility of multibody
repulsives would handle all the possible different chemical situations8 .
As we know, however, the above ideal picture is far from practice. First,
reducing the form of repulsive energy to a sum of pairwise potentials makes
various chemical environments9 of a bond non-trivially different and therefore
leads to our tedious work of balancing between different fitsystem chemistries.
Second, a serious contradiction between energy, force and vibrational targets
is produced by restricting repulsive energy as a sum of short-range pairwise
potentials. The latter restriction prevents fitted repulsives from following any
difference between DFTB electronic energy and reference energy properly in
the middle and far distance region. This produces a conflict between absolute
repulsive energy levels and force profiles: bridging the gap between short-range
8 As there may be only a finite number of other atoms acting as influencing neighbours on
an atom, a theoretical fit with multibody potentials up to some high order of centres would be
able to describe chemistry in an extremely accurate way, just like calculated with quantummechanical methods, within a reasonable range of chemical situations (e.g. high ionization
levels and other extremistic states excluded).
9 With three-body and higher-order multibody parts neglected, the effects of all these
appear as different chemical environments for a certain type of atom pair.
58
Figure 2.2: The Cl-Cl repulsive potential fitted without any overbinding strategy
and the one with one-body energies (the upper and lower strong dashed lines,
respectively). Long-range Eref EKS differences (dashed thin line originally,
solid line shifted to zero dissociation) and the pairwise overbinding potential
(dotted) that is equivalent to the applied one-body terms for this dimer are also
shown. (Note, however, that one-body energies and pairwise overbinding potentials cannot be equivalent for several distinct stoichiometries; an equivalence
is valid only for one stoichiometry, e.g. Cl2 in this case. Note also the difference between the lower, dissociation-adjusted fitted repulsive and the reference,
resulting from the absence of middle-range repulsives.)
E [a.u.]
0.2
0.1
-0.1
r []
and long-range energies (i.e. a proper energy fit) but having a slope that bridges
this gap over a much longer range (a proper force fit) is an unsolvable conflict.
As it is also evident from Figure 2.2, covering repulsive energy with short-range
potentials is simply not enough. This is the problem that will be studied briefly
in this section.
Taking a closer look at Figure 2.2, one can see that discrepancy on this dimer
system with respect to the short-range pair potential picture is composed of two
main parts. First, the repulsive energy dies out in a much longer range than we
implicitly expect with our short-range potentials. Second, there is a remanent
discrepancy of repulsive energy in the dissociation limit (r ). While the
former difference in pair potential ranges is rather easy to understand (but not
to solve!), one may be shocked about the fact that dissociated molecules are
not converging energetically to sets composed of lone model atoms in DFTB
and the reference method, which means that for a molecule of stoichiometry
Zl Zm
. . . Zn ,
E (dissoc) (Zl Zm
. . . Zn ) = lE(Z) + mE(Z ) + nE(Z )
(2.44)
in the dissociation limit. The most plausible explanation of this phaenomenon

is that differences of electronic structures and energies of dissociated chemical
59

systems with respect to the sum of lone atoms are not the same10 . Still, if the
shifts were consistent with each other in DFTB and in reference methods, the
problem would be away (a consistent change in electronics and a consistent shift
in atomic energies would cancel even if dissociated systems are not converging
to sets of lone atoms in DFTB and the references), but this is ostensibly not
the case.
To address the dissociation discrepancies outlined above, traditional parametrization uses the so-called overbinding strategy. This means adding pairwise more or less regularized profiles to each fitted short-range repulsive pair
potential in order to kill their sudden slopes at the outer side of the physical
first-neighbour region (the so-called Sprungschanze, meaning ski jumping hill
in German, see the repulsive without overbinding in Figure 2.2). There are two
serious problems with this approach. First, the size of overbinding and its step
profile is chosen in a rather ad hoc way. It is set up without any guarantee that
it will never affect physical bonding or other interesting situations. The only requirement against it is to eliminate the Sprungschanze on fit systems and thus to
produce reasonable dissociation curves. Second, this solution means a pairwise
energy correction in the dissociation limit, that is mathematically not feasible
(see below). In addition to the two major problems, the shape and sometimes
even the magnitude of overbinding used to go rather poorly documented.
To eliminate the shortcomings of traditional overbinding, yet keeping in
line with the strict pairwise nature of DFTB and knowing that a fit to higherorder repulsive potentials would be practically impossible, we remain within the
framework of using no higher multibody potentials than pairwise ones. In this
framework, however, one can use one-body and two-body repulsive terms
Erep =
atoms
UZA
atoms
1
UZA ZB (rAB ).
+
2
(2.45)
A=B
Of the two problems arising from the short-range nature of fitted repulsives,
one may try to handle the middle-range problems with middle-range repulsives,
and the dissociation problem with one-body and pairwise long-range terms. At
least for the first look.
However, infinite-range (dissociation) limit of energy correction cannot be
treated with pairwise terms. Energy being an extensive quantity, it must be an
additive function of stoichiometry when separated systems are added together
10 The argumentation of this paragraph is rather sloppy, preparing the unprecise conclusion
that one-body repulsive energy terms are the virtually only correct solution to the dissociation
energy problem. Carefully stated, the remanent dissociation energy discrepancy is present in
the practical, not the actual dissociation limit, the state where atom pairs have reached distances from each other outside the ranges of respective repulsive potentials (this is a practical
r state from the point of view of repulsive fitting although all r values are finite, of
course.) Therefore in fact, using one-body energies is the only mathematically acceptable way
of handling this problem if one accepts that any additional information coming from atomic
pairs being outside of their respective designated repulsive pair potentials must not influence
the physics of DFTB systems. On the contrary, the basic problem of traditional overbinding
strategy lies in carrying a hidden knowledge about which atom pairs can be, will be or have
ever been in certain types of bonds in their history even beyond the distance of their respective
repulsive ranges. This footnote has been added to the text as a means of clarifying a possible
logical flaw after submitting the thesis for supervision.
60
Figure 2.3: The Cl-Cl dissociation curve calculated with the non-overbinding
and the one-body-corrected repulsive of Figure 2.2 (the dotted and the dashed
strong lines, respectively). The thin solid line is the DFT reference, and the
dashed one is what it seems to be without dissociation adjustment. (Note the
middle-range difference between the reference and the dashed strong line also
here.)
E [a.u.]
0.1
-0.1
-0.2
r []
(but kept separated)
E (dissoc) (Zl Zm
. . . Zn + Zp Zr . . . Zs )
= E (dissoc) (Zl Zm
. . . Zn ) + E (dissoc) (Zp Zr . . . Zs )
= (l + p)E(Z) + (m + r)E(Z ) + . . . (n + s)E(Z ).
(2.46)
It can be easily derived from the above formula, that only one-body repulsive energy terms fulfill the additivity criterion and can make the two sides of (2.44) if
not equal, consistent between reference and DFTB. Therefore, a correct, consistent replacement of traditional overbinding will use one-body repulsive11 energy
parts.
On top of one-body regime, the non-dissociation (short and middle-range)
repulsive part in question must be fitted as a sum of pairwise profiles. The
fact that pairwise repulsives are not only a short-range phaenomenon can also
be seen in Figure 2.2 and Figure 2.3. Unfortunately, the attempts made at
fitting middle-range repulsive pair potentials are highly limited by a rather
practical circumstance: DFTB and most DFT methods tend not to converge
when atoms leave each others region of influence. The chlorine-chlorine dimer
used to draw Figure 2.2 and Figure 2.3 was a lucky example to present a longrange behaviour of pair energy, but the usual course of these calculations leads
to heavy convergence problems even with simple molecules of chlorine and other
11 Although it is meaningless to call a one-body energy repulsive, we keep this terminology
as a name for elements of the parametrized, quasi-classical part of DFTB total energy.
61

elements. One can try more complex molecules or reaction paths that keep all
the partner atoms of a pair bound throughout the whole fit as fit systems, but
this leads to an overwhelmingly complicated system of middle-range interactions
between participating atoms. This is why repulsive potentials reaching beyond
the typical first-neighbour distance are generally rejected as an option when
possible ways of improving DFTB parametrizations comes into question.
Anyway, a practical full implementation of the one-body + middle-range
strategy outlined and proposed here will need further considerations on the
practical side. The one-body part of it can be regarded working and efficient,
however.
Fitting one-body repulsives
Having one-body repulsives modifies the (2.2) composition of Erep only slightly;
Erep =
ZZ , XZZ , +
UZ NZ
(2.47)
ZZ ,
where UZ is the one-body repulsive for atom type Z and NZ is the number
of such atoms in the fit system in question. As it can be seen from the above
composition, one-body terms behave just like two-body terms in the fit if UZ s
go along s and NZ s along Xes. Thus the (2.7) X and A matrices of fitting
are modified so:
(1)
(2)
XHH,1 XHH,1
(1)
(2)
XHH,2 XHH,2

.
..
(2)
(1)
XCH,1 XCH,1
(1)
X=
(2.48a)
XCH,2 . . . . . . . . . . .
.
.
.
(1)
(2)
NH

NH
(1)
(2)
NC
NC

..
.
HH,1
HH,2
..
.
CH,1
A = CH,2
(2.48b)
..
.
UH
UC
..
.
where we remained at the example composition of systems containing H and
C. With the above modifications, the fitting algorithm can fit for one-body
terms in the way presented for two-body terms. Of course, fit targets that are
homogeneous, in the sense that they do not depend on the absolute value of
62
2.4. AB INITIO REPULSIVES

Erep , cannot interfere with one-body terms. Thus with energy derivative (e.g.
force) and energy difference targets, the corresponding X sections contain zeros
only. If a fit contains no inhomogeneous energy target, a fit for one-body terms
is naturally impossible.
An important detail of one-body repulsive fitting is that for the matrices to
have the proper rank, the rows of X containing N s must be linearly independent. If this is not the case (e.g. having NCl3 , NHCl2 and NH2 Cl as fit systems,
the three rows of N, Cl and H are linearly interdependent), aggregates of onebody energies for each fit system may be fitted as a circumvention. In this case,
dissociation energies of single atomic species cannot be calculated. One builds
the fit systems stoichiometries from the aggregates, and the X matrix has the
rows corresponding to this. In the former example, these aggregates would be
NCl3 and NH2 Cl (then 2NHCl2 = NCl3 + NH2 Cl) or NCl3 and HCl1 (in the
latter case, NHCl2 = NCl3 + HCl1 and NH2 Cl = NCl3 + 2HCl1 ).
2.4
Ab initio calculation of repulsive energy
In this section, we try to investigate the possibilities of replacing the fitted

repulsive potentials of DFTB with directly calculated ones. While this leads
to losing flexibility, reasonable first-guess repulsive sets without fitting would
mean a tremendous improvement at getting DFTB work with new chemical
elements. It must be also noted that being a fit for middle-range repulsive pair
potentials problematic (see 2.3.2), calculated repulsives seem to give the only
opportunity to produce physical repulsives in this region. In addition to this,
the dissociation limit of repulsive energy is also handled correctly with them by
nature. [26] can be regarded a very simple first step in this direction, trying to
cover pair repulsives with Coulomb interactions of compressed atomic densities
with freely optimized compression radii.
2.4.1
Theory
As it was seen in (1.33), the repulsive energy is the sum of the so-called counterterms of double-counting terms in KohnSham energy (for valence electrons),
as well as the total energy of core electrons. These together formed
1
Erep = Tr Pvalence VC Vxc
2
KS 1 VC Vxc
+ Exc + Enuc
+ Tr Pcores H
2
63
(2.49)

that may be written using integrals instead of operator traces and dissolving
KS to its parts
H
Erep =Tr[Pcores T]

1 (x)((y) + (y))
+
dy
2
|x y|
x
(x)Vxc ((x) + (x))
(x)((y) + (y))
+
dy
|x y|
y
+(x)Vext (x)
(2.50)
+(x)Vxc ((x) + (x))
1 (x)((y) + (y))
dy
2
|x y|
y
(x)Vxc ((x) + (x))
+((x) + (x))((x) + (x)) dx

+Enuc
where we denote the core and valence electronic charge densities by and .
Rearranging the above equation,
Erep = Tr[Pcores T] +
1 (x)(y)
dy +
2 |x y|
1 (x)(y)
dy
2 |x y|
+ ((x) + (x))((x) + (x)) (x)Vxc ((x) + (x)) + (x)Vext (x) dx

atoms
1 ZA ZB
+
.
2
rAB
(2.51)
A=B
The terms of the above sum are nothing else than

the kinetic energy of frozen core electrons,
minus the Coulomb energy of valence electrons,
plus the Coulomb energy of core electrons,
plus Exc ,
minus the KohnSham xc term for valence electrons,
plus the potential energy of core electrons in the nuclear potentials,
plus internuclear Coulomb repulsion.
64

Within the validity range of DFTB approximation, the above integrals can be
calculated in pairwise calculations along the Hamiltonian, overlap and position
matrix elements. It is an important detail of these calculations that intraatomic
Coulomb interactions can be omitted from the sum since they are also constant
just like the core kinetic term. Therefore, the terms in (2.51) can be calculated
with an AB atom pair as AB and BA interaction energies, AA and BBtype terms may be omitted. Thus for an AB pair of atoms,

Erep (AB) =
A B
+
r
A B
r
(A + B + A + B )(A + B + A + B )
(A + B )Vxc (A + B + A + B )
ZA ZB
+ A Vext,B + B Vext,A +
rAB
(2.52)
(with spatial arguments of the integrands omitted, and |x y| of (2.51) substituted with r). In order to maintain the zero dissociation limit of repulsive
pair profiles (see 2.3.2), we must subtract its dissociation limit from the pair
repulsive energy. Thus
UAB = Erep (AB)
(A + A )(A + A ) + A Vxc (A + A )
(B + B )(B + B ) + B Vxc (B + B ).
(2.53)
With this definition of UAB it remains true within the precision of negligible
three-center terms (see the argumentation about density superposition at the
end of 1.2.1) that for any chemical system
Erep =
1
UAB (rAB )
2
(2.54)
A=B
up to a stoichiometry-dependent constant, which is the sum of one-body repulsive energies:

(dissoc)
Erep
=
UA =
(A + A )(A + A ) A Vxc (A + A ) .
A
(2.55)
If one pursues absolute energy level compatibility (and comparability) with
all-electron DFT calculations, the constant terms should contain all the oneatom contributions in the above sum. The terms that constitute the atomic
part of the full repulsive energy and must be used instead of the UA s in (2.55)
are
The kinetic energy of core electrons:
T(A) = Tr(Pcore,A T) =
core of A
65
|T |
(2.56)

The Coulomb energy of core electrons:

1
A (x)A (y)
(A)
EC, core =
dxdy
2
|x y|
(2.57)
x y
The atomic Coulombic counterterms of valence electrons:

A (x)A (y)
1
(A)
EC, valence =
dxdy
2
|x y|
(2.58)
x y
The one-atomic xc counterterm in (2.55) for A:
(A)
(A)
Exc
EKS,xc = (A + A )(A + A ) A Vxc (A + A )
(2.59)
The potential energy of core electrons in their nuclear Coulomb potential:
Eext,core = A Vext,A
2.4.2
Further technical details of repulsive calculation
The specific structure of DFTB machinery gives two more particular aspects
which must be considered at repulsive calculations. First, the calculated repulsive energy uses atomic quantities that need to come from distinct pseudoatomic
calculations. Exchange-correlational potentials and energy densities come from
pseudoatoms compressed with the density compression radius. Densities used
as the source of Coulombic potentials also come from these pseudoatoms. Densities, however, which counterbalance P terms in counterterms, must be compressed with the wavefunction compression radius. Core densities that are not
used to build potentials do not play a role in the Hamiltonian and the Kohn
Sham energy, and they can be attributed a quite arbitrary compression. To
avoid further complication of the structure of repulsive calculations, we make
them behave always parallel to valence densities. Thus the pair repulsive becomes
UAB =
1

A B
1
A B
1
A B
1
A B
+
+
r
2
r
2
r
2
r
+ (A
+ B
+ A + B )(A
+ B
+ A + B )
(A + B )Vxc (A
+ B
+ A + B )
ZA ZB
+ A
Vext,B + B
Vext,A +
rAB
(A
+ A )(A
+ A ) + A Vxc (A
+ A )
(B
+ B )(B
+ B ) + B Vxc (B
+ B ) (2.60)
66

and the one-body term for atom A is

A
A
1
A A
UA =
r
2
r
core of A
+ (A + A )(A + A ) A Vxc (A + A ) + A
Vext,A ,
1
|T | +
2
(2.61)
where quantities marked with come from an atomic calculation with the socalled density compression radius (which would be better called potential compression radius), and those marked with come from atomic calculations compressed with the so-called wavefunction compression radius. Later, it may be
needed to bring in more tunable parameters in calculated repulsives. To fulfill this need, we will let the potential and wavefunction compression radii of
repulsive calculation (r2 and r3 ) be other than the same radii of Hamiltonian
construction (r0 and r1 as it was introduced in 1.2.1).
The second technical detail that needs consideration here, is how potential
and density superposition affects calculating repulsive potentials. Clearly, this
question affects the
(A +B +A +B )(A +B +A +B ) (A +B )Vxc (A +B +A +B )
(2.62)
part of repulsives. In density superposition, the formulae to calculate the xc part
of repulsives are literally as above, densities are superposed in the arguments
of xc-density-type quantities ( and Vxc ). In potential superposition, the dimer
xc potential is the Vxc (A + A ) + Vxc (B + B ) superposition of atomic ones,
while the dimer can be two different ways. Either we leave the above (A +
B + A + B ) xc energy density because we do not want to tinker with the
inner structure of , or we can use the
(A + B + A + B ) =
(A + A )(A + A ) + (B + B )(B + B )
(2.63)
A + A + B + B
weighted sum of atomic xc energy densities (which gives back the potentialsuperposed Vxc as its derivative with respect to , as well as it gives back the
dimer xc energy as the sum of atomic xc energies). In our attempts made so far,
we have not tried and experimented with the potential superposition schemes
proposed here.
2.4.3
Results
As a first example of CC repulsive calculation, we took the electronic parameters of the mio set (r0 = 7 a.u., r1 = 2.7 a.u.), and optimized the wavefunction compression radius used for the repulsive, independent of the wavefunction
compression radius used for the Hamiltonian (the potential compression used
for repulsive remained equal with the potential compression used for electronic
tables, r2 = r0 ). This choice was an ad hoc first trial to calculate repulsives
after realizing that letting both compression radii for making repulsives equal to
the electronic tabulation counterpart gives unusable results (this problem may
be an effect of choosing compression radii rather freely in the world of fitted
67

repulsives, and may be overcome with an extensive reoptimization of electronic
parameters).
The optimization, which was carried out by the brute-force electronic parameter search12 described in 2.1.7, produced a very unusual result: the optimal wavefunction compression radius for repulsive calculation was as small as
possible, r3 = 0.01 a.u. in this optimization (recall that this radius was made
independent of the wavefunction compression radius used for integral tables,
and therefore this optimization had no effect on the electronic properties of
molecules). Despite of the surprisingly small value, the quality of the resulting
dissociation curve of CC dimer (Figure 2.4) encourages further investigation.
To assess the method in real-world circumstances, we optimized a set for carbohydrogens. The hydrogen parameters optimized on the CH bond behaved
similar to the CC case: an optimal independent wavefunction compression
radius of H for repulsives turned out to be the smallest value allowed by the
optimization range, 0.001 a.u. in this case. It must be noted that optimizing
the CC repulsive on carbohydrogen systems gave the same optimal parameter
(note that only r3 was optimized here) as the optimization on the CC dimer.
The geometrical and energetical test results of the resulting CH parameter set
can be seen in Table 2.17 (the ab initio reference for comparison with calculated
repulsives was similarly B3LYP as in 2.2.1). It must be also noted that vibrational frequencies are systematically underestimated with these repulsives by
1015%, which is an effect of the extremely compressed nature of the wavefunctions used to calculate repulsives. We expect huge improvement from a more
thorough reoptimization of all electronic parameters.
Another test of the calculated repulsives was a HCO set (an extension of
the above set with oxygen). Producing it was a tougher process, as traditional
electronic parameters of oxygen turned out to be unusable with our calculations
(it was impossible to optimize a reasonable pair of compression radii for repulsives on top of the mio or the matsci [16] electronic tables). After a heuristic
search for usable ranges of new tabulation electronic parameters for oxygen, a
more systematic optimization of a scheme seen at C and H has been made: a
single value for both tabulation and repulsive potential compression radii was
optimized with an optimization of the tabulation wavefunction compression in
parallel, with the repulsive wavefunction compression fixed to r3 = 0.01 a.u.
(this was a guess based on the CH search). The fit systems for this brute-force
fit were CO2 , H2 O, O2 and C2 H5 OH. This resulted in the optimized potential
compression radius for electronic tables and repulsives being r0 = r2 = 3.4 a.u.,
and the tabulation wavefunction compression being r1 = 1.5 a.u.. The performance of the resulting repulsives is illustrated in Table 2.18. It must be noted
here, that this reoptimization of electronic parameters is a foreshadow of a possible reoptimization of all electronic parameters, based on energetical objectives,
if calculation of repulsives becomes a successful, method and common practice
(note the good absolute energy values of CO and HO bonds with the more
deeply reoptimized oxygen part in Table 2.18).
12 Of course, the brute-force optimization process was a little bit altered to use calculated
repulsives instead of fitted ones: the fitted part was a dummy set while the electronic part
contained all the information about calculated repulsives.
68
Figure 2.4: The CC dimer dissociation curve. The solid line is the DFT reference, and the dashed line is the one produced with DFTB and a calculated repulsive where electronic parameters of the mio set were used for electronic tables,
and the potential compression radius for repulsive was equal to the tabulation
one (r2 = r0 ). The wavefunction compression radius for repulsive calculation
was 0.2 a.u. in the top graph (for a better comparison of curve shapes, both
curves were shifted to an arbitrary absolute level determined by the one-body
repulsives used here). The same parameter of the bottom graph was optimized
further to 0.01 a.u. to give better geometries (this latter parametrization was
used for test calculations).
E [a.u.]
-150
-150.5
-151
1.5
2.5
r []
E [a.u.]
1.5
0.5
1.5
r []
69

Table 2.17: Test data of first-glance ab initio repulsives for carbohydrogen systems (the G2 [8] test set and diamond). Es are
atomization energies in eV, doubles are bond lengths in
A, triples
are angles in degrees.
property
methane
E
CH
ethane
E
CC
CH
ethene
E
C=C
CH
ethyne
E
CC
CH
benzene
E
CC
CH
butane
E
CC (end)
CC (middle)
CH (end)
isobutane
E
CC
CH
diamond
CC
cyclobutane
E
CC
CH
isobutene
E
CC
C=C
CH in (CH3 )
CH in (CH2 )
ref.
DFTB
error
rel. error
-18.266
1.09302
-15.344
1.10835
2.92274
0.01533
-0.1600
0.01403
-30.881
1.53061
1.09615
-26.469
1.58222
1.11910
4.41151
0.05162
0.02295
-0.1429
0.03372
0.02093
-24.437
1.33103
1.08743
-21.201
1.31039
1.11682
3.23639
-0.0206
0.02940
-0.1324
-0.0155
0.02703
-17.485
1.205
1.06653
-16.442
1.15838
1.08725
1.04324
-0.0466
0.02073
-0.0597
-0.0387
0.01944
-59.083
1.39663
1.08699
-53.143
1.41683
1.12142
5.94045
0.02021
0.03444
-0.1005
0.01447
0.03168
-55.967
1.54680
1.53643
1.09667
-48.644
1.60341
1.59862
1.11864
7.32247
0.05661
0.06219
0.02196
-0.1308
0.03660
0.04048
0.02003
-56.235
1.53530
1.09665
-48.835
1.59850
1.11919
7.40019
0.06320
0.02254
-0.1316
0.04116
0.02055
1.555
1.62533
0.08533
0.05541
-49.524
1.55728
1.09484
-43.599
1.62048
1.12521
5.92536
0.06319
0.03037
-0.1196
0.04058
0.02774
-50.066
1.50876
1.33658
1.09898
1.08736
-43.821
1.57347
1.3234
1.12092
1.11584
6.24536
0.06471
-0.0132
0.02194
0.02848
-0.1247
0.04289
-0.0099
0.01996
0.02619
70
bicyclobutane
E
CC (edge)
CC (middle)
CH (in CH2 )
CH (in CH)
cyclobutene
E
CC
C=C
CH (in CH2 )
CH (in CH)
cyclopropane
E
CC
CH
propane
E
CC
CH (end)
CH (middle)
cyclopropene
E
CC
C=C
CH (opposite to C=C)
CH (neighbour to C=C)
spiropentane
E
CC (radial)
CC (outer)
CH
methylene-cyclopropane
E
C=C
CC (radial)
CC (outer)
CH (in CH2 )
CH (on ring)
propadiene
E
C=C
CH
-39.222
1.51018
1.90017
1.11157
1.09512
-34.175
1.57435
2.08801
1.15735
1.12634
5.04709
0.06418
0.18785
0.04578
0.03122
-0.1287
0.04250
0.09886
0.04119
0.02851
-43.155
1.57318
1.51869
1.09709
1.08716
-38.052
1.64547
1.59169
1.12671
1.12103
5.10254
0.07229
0.07299
0.02962
0.03387
-0.1182
0.04595
0.04806
0.02700
0.03115
-36.872
1.50922
1.08659
-31.570
1.57198
1.11653
5.302
0.06276
0.02994
-0.1438
0.04158
0.02755
-43.548
1.53199
1.09717
1.09854
-37.630
1.59033
1.11925
1.12985
5.91752
0.05834
0.02208
0.03131
-0.1359
0.03808
0.02013
0.02850
-29.452
1.50826
1.29488
1.09518
1.07983
-25.360
1.57681
1.30087
1.12367
1.11128
4.09168
0.06855
0.00599
0.02849
0.03145
-0.1389
0.04545
0.00463
0.02601
0.02913
-55.408
1.48491
1.52990
1.08815
-47.655
1.55222
1.60192
1.11817
7.75290
0.06731
0.07202
0.03002
-0.1399
0.04533
0.04707
0.02758
-42.934
1.32232
1.46992
1.53954
1.08776
1.08928
-37.323
1.30854
1.53619
1.60484
1.11987
1.11842
5.61093
-0.0138
0.06627
0.06530
0.03210
0.02914
-0.1307
-0.0104
0.04508
0.04242
0.02951
0.02675
-30.567
1.30686
1.08824
-27.153
1.29116
1.12033
3.41345
-0.0157
0.03209
-0.1117
-0.0120
0.02949
71
1,3-butadiene
E
CC
C=C
CH (middle)
CH (end)
2-butyne
CC
CC
CH
propyne
E
CC
CC
CH (in CH3 )
CH (in CH)
propene
E
CC
C=C
CH (in CH3 )
CH (in CH2 )
CH (in CH)
2.4.4
-39.068
1.43937
1.39158
1.08870
1.08596
-34.610
1.49511
1.36317
1.12067
1.12845
4.45754
0.05574
-0.0284
0.03197
0.04249
-0.1141
0.03872
-0.0204
0.02937
0.03913
1.46169
1.20920
1.09707
1.51404
1.16399
1.12084
0.05235
-0.0452
0.02377
0.03582
-0.0374
0.02167
-30.483
1.46027
1.20715
1.09689
1.06610
-27.785
1.51272
1.16164
1.12004
1.08671
2.69794
0.05245
-0.0455
0.02314
0.02061
-0.0885
0.03592
-0.0377
0.02110
0.01933
-37.260
1.50206
1.33336
1.09846
1.08675
1.09114
-32.486
1.56603
1.31711
1.12126
1.11577
1.12569
4.77352
0.06397
-0.0163
0.02280
0.02902
0.03455
-0.1281
0.04259
-0.0122
0.02076
0.02670
0.03166
Prospects of calculating repulsives
The results shown here illustrate that the proposed machinery of calculating
repulsive potentials instead of fitting them is able to produce parametrizations
that give less precise results than the best fitted repulsive sets, but they still
ensure reasonable, qualitatively good results right at the first glance.
The calculated repulsives still contain parameters that allow semiempirical
tuning of them to ensure the least difference with respect to reference methods,
but optimization of them means a huge simplification of the parametrization process compared to fitting pair potentials. First, they are two new real numbers
per atom, instead of the pair potentials, that have infinite degrees of freedom;
and second, they are atomic, not pairwise parameters. Once the transferability
of these new parameters is proved (this is foreshadowed by the fact that CC
repulsives got from dimer optimization were the same as the CC parameters
from the CH optimization), their atomwise nature ensures that extending a
parameter set to a new atomic species means tuning only one set of atomic parameters (two compression radii), not a complete set of new pairwise repulsives
between the new species and all the previous members of the parameter set.
The electronic parameters for calculating HCO systems optimized in our
first usages of the repulsive calculation scheme seemingly leave more space for
further optimization. Besides the performance of the calculated parameters,
that is surprisingly good compared with the simplicity of theory, but still lag
72
Table 2.18: Results with HCO compounds. E: atomization energy in eV,

doubles are bond lengths in
property
ref.
DFTB
error
rel. error
CO2
E
-16.851 -16.538 0.31332
-0.0186
O=C
1.16005 1.19163 0.03158 0.02722
H2 C=O
CH
1.10671 1.17340 0.06669 0.06026
C=O
1.19896 1.21711 0.01815 0.01513
HC=O
122.115 122.492 0.37638 0.00308
H2 O
E
-9.9215 -10.691 -0.7700
0.07760
OH
0.96107 0.97234 0.01126 0.01172
HCOOH
E
-21.446 -21.038 0.40824
-0.0190
C=O
1.19023 1.21953 0.02931 0.02462
CH
1.10330 1.17803 0.07473 0.06773
CO
1.35071 1.42625 0.07554 0.05593
OH
0.96462 0.98630 0.02168 0.02247
O=CH
124.024 126.909 2.88527 0.02326
COH
109.933 109.365 -0.5684
-0.0052
C2 H5 OH
E
-34.971 -32.068 2.90355
-0.0830
CC
1.51430 1.58271 0.06841 0.04517
CO
1.42674 1.47539 0.04865 0.03410
CCO
107.982 107.458 -0.5239
-0.0049
CO
E
-11.087 -11.486 -0.3994
0.03602
CO
1.12555 1.12138 -0.0042
-0.0037
CH3 OCH3
E
-33.057 -29.072 3.98538
-0.1206
CO
1.37698 1.42937 0.05239 0.03805
CH
1.09697 1.14405 0.04708 0.04292
CH3 OH
E
-22.149 -20.769 1.38017
-0.0623
CH
1.08850 1.12324 0.03474 0.03192
CO
1.42018 1.47134 0.05116 0.03603
OH
0.96029 0.98862 0.02833 0.02950
H2 O2
E
-11.228 -14.415 -3.1875
0.28389
HO
0.96560 0.99593 0.03033 0.03141
OO
1.46043 1.36649 -0.0939
-0.0643
HOO
104.931 104.601 -0.3295
-0.0031
O2
OO
1.21
1.17507 -0.0349
-0.0289
73

behind that of fitted ones, the actual values of optimized parameters give a
hint about their imperfection: the extremely low r3 values suggest that they
try to compensate for a too large r2 . We have tried to re-optimize the CH set
with r2 = r3 , in order to keep the number of independent parameters as few as
possible, but the resulting parameters were not as good as those presented above.
A thorough re-optimization of electronic compression radii may be the solution,
however. This new optimization may also be an opportunity to optimize the
electronic properties of atoms (through r0 and r1 ) along energetical objectives.
74
Chapter 3
Enhancement proposals to
SCC-DFTB
3.1
Multipole expansion of SCC energy
The off-site energy of interacting atomic charges presented in (1.56) is a very

simple approximation; the SCC mechanism based on it considers every atomic
excess charge point-like. Remaining within the computational complexity level
of present SCC, we can make SCC potentials more precise, however. This is
done via a multipole extension of the current SCC scheme.
To derive an extension of SCC with multipole terms, we start with the
1
(y)(y)(x)dxdy
(x)
(3.1)
|x y|
x,y
integrals in (1.48) and expand them as multipole series. A multipole expansion

of an expression like
(x)(y)
dxdy
(3.2)
|x y|
x,y
around an x0 , y0 pair of base points is

1
m!n!
m,n
m n
1
xm yn |x y| x0 ,y0
m
(x x0 ) (x)dx (y y0 )n (y)dy
x
(3.3)
1
kernel. After rearranging
by simply using its Taylor series instead of the |xy|
the off-site part of the (1.48) sum (quite similarly to the rearrangement that led
to the (1.55) original off-site SCC formula) and using multipole-series expansions
instead of the (3.1) integrals, we get a multipole extension of (1.55). It is the
same as if we took the atomic Mulliken densities defined in (1.59) as and ,
and used the multipole expansion of their interaction integrals. Of them, the
75
CHAPTER 3. ENHANCEMENT PROPOSALS TO SCC-DFTB

n = 0, m = 0 term is a monopole-monopole interaction, the n = 1, m = 0 term
is a dipole-monopole interaction, etc.
In the course of the above rearrangement, retaining all multipole levels in the
series of (1.48) integrals, the rearranged sums of (1.50) become not as simple
as in (1.53) but containing a series in the place of
1
S S ,
|xA xB |
{A}, {B}.
(3.4)
Plausibly we choose atomic centres A and B as base points of the Taylor series
in the multipole expansion, and so the first terms of the series taking the place
of (3.4) are
1
|xA xB | S S
(A)
(B)
xA xB
|xA xB |3 S X
(A) (B)
A xB )(xA xB )
|xA xB | 2(x
X X + . . .
|xA xB |5
xB
|xxAAx
3 X S +
B|
(3.5)
(A)
where the X position operator matrix is quite similar to the S overlap
matrix except that it calculates the first moment of the (x)(x)

product with
respect to the Ath atomic centre:
(A)
X = |X XA |x = (x)(x
xA )(x)dx.
(3.6)
x
We keep the SCC philosophy of including smearing and xc energy via

profiles instead of the classical Coulomb profile, therefore our new SCC energy
will employ the derivatives of s1 instead of those of 1r (which can be regarded
the long-distance limit) in (3.5). With this all, the total SCC energy (formerly
(1.64) in the monopole-monopole approximation) becomes
2 E =
1 (00)
(10)
AB qA qB + AB dA qB
2
A,B
(01)
(11)
+ AB qA dB + AB dA dB +
where
(mn)
AB
1 m n
n AB ,
m!n! xm
A xB
and we introduced the dA Mulliken dipole moment of the Ath atom:
1
(A)
(A)
dA =
ni ci, ci, X + ci, ci, X .
2
(3.7)
(3.8)
(3.9)
{A},,i
d is the change of the Mulliken dipole with respect to the non-interacting

atoms.
1 Instead of using derivatives of the old s as multipole-multipole interaction profiles,
theoretically we could use semiempirical interaction profiles for each level of interaction independent of each other. This is a key feature of improving results with constructions similar to
multipole series, for example in neglect of differential overlap methods and their derivatives
(see e.g. [7]), but in the case of DFTB, increasing the number of parameters is highly undesired. We also want to make the first steps of enhancing SCC energy relatively simple and
theoretically not complicated.
76
3.2. SEMIEMPIRICAL
In the following, we will use the very first terms above the original SCC level,
the monopole-dipole ones. The construction can be extended to higher terms
as well. Up to the desired monopole-dipole level, the SCC Hamiltonian matrix
becomes (cf. (1.71) in the monopole-monopole approximation)
(00)
1
(00)
S
[]B + []B qB
2
B
(01)
1
(01)
+ S
[]B + []B dB
2
B
1 [] (10)
[] (10)
+
X []B + X []B qB .
2
H =
(3.10)
The transformation rules needed to relate the (3.6) position matrix elements
in the molecular coordinate system to the tabulated ones are presented in 3.3.3.
3.2
3.2.1
Semiempirical full chemical hardness matrices and profiles

Semiempirical off-site part
As described above in section 1.2.2, the current SCC-DFTB method uses a

interaction profile in the off-site part of SCC that is an interpolated shape
between the Coulombic 1r at r and the atomic chemical hardness at r = 0.
The heuristic derivation of this interpolated shape is based on the classical
Coulombic interaction of two exponentially smeared diffuse charges, and then it
is renormed to give the chemical hardness at zero distance. The latter correction
is considered to take the xc energy into account up to a relevant precision.
Extending the SCC-DFTB philosophy as a pairwise approximation to LCAO
DFT with semi-empirical elements and a simplified self-consistency, we can determine the profiles (that are semiempirical approximations of two-electron
integrals in a proper expansion of DFTB energy, see 3.2.4) with pairwise energy
calculations. These calculations use atom pairs, but otherwise they are just like
those with lone atoms when chemical hardnesses are determined. As in (1.62),
the strength of second-order response is the derivative of the energy of a highest
occupied state with respect to a population, but the occupied state and the
population belong to two different atoms:
A
(3.11)
(ZA , ZB , r) =
qB r
where atoms A and B, being of types ZA and ZB , are in a distance of r, A
is the highest occupied electronic energy level of the atom A, and qB is the
total electronic population of atom B. Similarly to the calculation of chemical
hardness in traditional SCC scheme, an excess charge on atom B is created by
varying the population on its highest occupied orbital.
Because in a self-consistent diatomic calculation the energy levels cannot
be attributed to specific atoms, the derivative in (3.11) can be calculated in
a perturbative way. If A is the atomic state behind A and VB is the total
77

potential that electrons feel around atom B,
A |VB |A
A
=
.
qB r
qB
3.2.2
(3.12)
s in line with subshell hardness
A further improvement of SCC, but independent of the above proposals, could

be the usage of the full Hubbard matrix as described in section 1.2.2 as a technical enhancement of the on-site limit of with respect to the current DFTB
implementations. In line with this, the above semi-empirical profiles could also
be enriched with a similar matrix structure. Thus (3.11) and (3.12) is modified
to
1 A,l,m |VB |A,l,m

A,l
=
(3.13)
ll (ZA , ZB , r) =
qB,l
2l + 1
qB,l
r
where l and l are the angular momentum indices of the interacting orbitals of
atoms A and B, respectively, and m is the magnetic quantum number of the
participating A orbitals, over which we make a mean because we want to have
a general l-l interaction profile (angular dependence of SCC interaction which
arises due to uneven filling of different-m orbitals is covered by multipole SCC).
We do not have to average over magnetic quantum numbers of atom B since
VB (x) is calculated with a spherically symmetric atom.
3.2.3
Fully resolved s
At the end of the way described up to now stands a complete refinement of s.

It goes beyond the Hubbard matrix of (1.72) and the of (3.13) and it resolves
all the m quantum number cases. Having all the electronic quantum number
cases resolved, we can index s with orbitals. In this way
=
|V[] |
=
q
q
(3.14)
where, of course, is the atomic energy eigenvalue of the orbital and q is

the population sitting on , and [] is the atom owning (naturally, there are
also U Hubbard matrices as the homonuclear, r = 0 cases of the above s).
We do not designate atom types, quantum numbers and the distance of [] and
[] explicitly here, as these are all obvious from the and elements of the
LCAO basis.
Since the orbitals used for tabulation will not be the same as the ones in the
molecular system, the above orbital-indexed s have to be built of stored s
in a similar way to what applies to the other tabulated matrices. The rules are
described in 3.3.3.
3.2.4
Further resolution of s as operators2
Although there is a strong competition between computational costs and calculational precision, that must result in a fair compromise, it is not worthless
2 This section takes motivation from discussions with Adriel Dominguez and Thomas
Niehaus.
78
3.2. SEMIEMPIRICAL
investigating where the further resolution of s can theoretically lead to. Some
aspects of such a further resolution are presented in the following.
Semioperator s
One step beyond the Mulliken-analysis-based improvement of more and more
resolved s, one can begin resolving them at least on one side to get fullyfledged orbital-indexed quantities, that is, operators written in the LCAO space,
not only the diagonal elements described in 3.2.3. This resolution beyond every
limit refines (3.14) on one side to a well-behaved orbital-indexed matrix:
; =
|V[] |
H
=
.
q
q
(3.15)
These s can be used in the calculation of total energy and self-consistent

Hamiltonian matrix much easier than former ones. In total energy
2 E =
1
P ; q
2
(3.16)
i,,,
would represent the SCC term, where P is the change of density matrix
(0)
element i ni ci, ci, with respect to P , and
; q
(3.17)
is an effective kernel that represents the effective semi-classical interaction of

orbital charges with point-like charges accounted for orbitalwise (represented
with the same q quantities as in 3.2.3). In the Hamiltonian matrix, the corresponding part is nothing else than the above kernel
KS
H
=
; q .
(3.18)
Of course, the advantage of possible refinement of calculation with these new

s may be less than the disadvantage of tabulating and using |{A}| |{A}|
|{B}| of them3 for every distance of every atom type pair, where |{A}| denotes
the number of orbitals centered around A. To remedy the problem of extensive
tables, and taking into account that the same refinement may not be needed
at either sides of , one can begin making the source side more accumulated.
This way, even an orbital-center-type set of s can be defined:
;B =
H
|VB |
=
.
qB
qB
(3.19)
The number of the above s is only |{A}| |{A}| per distance step per atom
pair type (this means practically |{A}| |{A}| + |{B}| |{B}| entries in an
A-B-type table row).
Besides the (3.15) arrangement where the three indices of came from centres AA; B, we must shortly define all the other possible combinations of centres
3 This number of
; values is valid when cases with [] = [] (i.e. three-center-type
ones) are neglected.
79

in one or two-center-like s. As long as we remain at point-like charges at the
source (nonoperator) side, calculation of s goes as in (3.15) or (3.19). The
complexity of dealing with multipole contributions is detailed in the following.
Accumulation of the Mulliken side needs further consideration in the non-AA; B
case, however. Keeping in line with the importance of Mulliken-side resolution
in different cases, an adequate minimization of necessary table sizes may be
realized by using a nonuniform resolution of the nonoperator (Mulliken) side;
e.g. for off-site s an atomwise resolution suffices while on-site s (chemicalhardness-type interaction intensities with AA; A-type and maybe AB; A-type
indices) need a subshell or orbitalwise resolution. Note that, depending on the
resolution of the quasiclassical side and the number of relevant nonequivalent
arrangements of the two atomic centres involved (e.g. AA; B and AB; A in the
heteroatomic case of semioperator s), the number of integral table entries
mentioned above must be multiplied with a small integer, but we will silently
omit treating this in detail, as we are only interested in magnitudes of table
sizes now.
Another aspect that must be clarified about the semioperator s in the
future is if some HartreeFock-type (HF) exchange energy can be expressed with
them, and if yes, how this can go. Certainly, fully operator-like s facilitate
calculating HF exchange, but it raises two major questions. First, calculating
exchange energy implies using a two-electron density operator (or a construction
alike), not only the one-electron one used in DFTB and in this thesis. Second, xc
parts are also expressed in profiles semiempirically, and it is a serious question
if we are able to balance between a reduced semiempirical xc energy and a new
HF exchange energy correctly.
Handling multipole contributions with semioperator s
An above-defined one-sided resolution of operator-level s affects the multipole
expansion too. Since the operator-like side does not need approximations like
multipole expansion any more, only the other side may conserve it. This goes
by defining the same Mulliken charges and dipole moments as with the atomresolved (i.e. non-resolved) or any other finer Mulliken case, but using s and
energy expressions like in (3.18). In this way, a matrix element of the SCC
Hamiltonian looks
(0)
(1)
KS
H
=
;B qB + ;B dB .
(3.20)
Note that higher-multipole s are produced also here as derivatives, but

this construction affects only one side of that is marked by only one derivative index. Being the interacting quantities (atomic potentials and products of
electronic orbitals) at calculations not spherically symmetrical at both participant centre any more, the derivatives of these s cannot be constructed
from the bare derivative with respect to the distance between centres. Different
derivatives must be calculated in all the three directions of space. This means
additionally calculated and stored quantities in a table when one wants multipole moments, that negatively affects the usability of semi-operator s. In this
(1)
(1)
(1)
case, one would store (0) , x , y and z .
The on-site (one-center) case is a bit more difficult because one has to put
(1)
dipole moments in atomic DFT calculations by hand to calculate ;B where
80
3.2. SEMIEMPIRICAL
, {B} because there is no distance parameter in chemical hardness with
respect to which this profile could be derived, as well as there are no nuclei able
to be moved with respect to each other. ;B when {B} but
/ {B} is a
similar complicated case: the problem of handling of higher multipole order s
also applies here because {B}.
There is, however, a possible way to circumvent the calculation of multipole
contributions with semi-operator s. Because of the symmetry between the two
sides of (recall that they are nothing else than a complete set of two-electron
integrals, or in other words, elements of a kernel of a second derivative of energy),
a total sum of all possible multipole contributions above monopole-monopole
level in SCC-DFTB can be approximated (exactly at the monopole-dipole level,
and with gradually increasing relative errors at higher levels) by
2 Eall multipoles = 2(2 Esemi-operator 2 Epoint charges ),
(3.21)
where the first term is calculated without any attempt to handle multipole
corrections, and the last term is a point-like-charges SCC energy calculated
with adequately refined s (in (3.11), (3.13) or (3.14)). This makes handling
multipoles even with semi-operator s unnecessary, because
2 E = 2 Epoint charges + 2 Eall multipoles
= 22 Esemi-operator 2 Epoint charges ,
(3.22)
containing the old SCC energy in addition the the semi-operator one, which can
be calculated far more easily than a multipole extension of semi-operator s.
Two-sided exact operators and their approximations implemented
by other schemes
Naturally, one can proceed to refine both sides of . This is nothing else than
giving back the (1.48) two-electron integrals faithfully in the SCC part. While
this means an extremely populous set of ; parameter tables, one can consider neglecting three-center terms also here, arriving at the still very large
|{A}| |{A}| |{B}| |{B}| number of table entries per atom pair per distance
value (this is only an estimation of magnitude, also here). Of course, this regime
makes multipole calculations away, and it gives a possibility of calculating HF
exchange in a more perfect way.
A useful byproduct of the investigation of operator-like s made so far is
that we are now able to compare them with Mulliken-level s and clarify what
and how is simplified in point-like-charges SCC, and what and how is improved
in the multipole extension of it. The two-sided operator-like s represent the
(1.48) four-orbital integrals without any approximation beyond the basic DFTB
rules, but with the second derivative of xc energy added to the Coulomb kernel,
of course:

2 Exc
1
+
(y)(y)dydx.
(3.23)
; =
(x)(x)
|x y| (x)(y)
x y
The fully-resolved, but not operator-like s in 3.2.3 simplify the above plenty
of integrals to

1
2 Exc
=
(x)(x)
+
(y)(y)dydx,
(3.24)
|x y| (x)(y)
x y
81

but considering the structure of the SCC Hamiltonian, this gives also an approximation to the four-indexed :
;
1
S ( + + + )S .
4
(3.25)
In the subshellwise SCC Hamiltonian, the above structure is maintained while

the s inside are replaced with an average of e.g. s resulting in ll . In
the original, atomwise SCC, ll s for every atom pair are replaced with the
one describing interaction between the respective highest occupied subshells of
atoms in the pair.
The multipole expansion makes the above approximation more precise by the
multipole regime described in (3.3). This is built of the derivatives of the kernel
and the multipole moments of the two sides. Up to the first (monopole-dipole)
level, the multipole approximation of the four-indexed reads
1
(00)
(00)
(00)
(00)
S ( + + +
)S
4
1 [] (10)
[] (10)
[] (10)
[] (10)
)S
+ (X + X + X + X
4
1
(01) []
(01) []
(01) []
(01) []
X ),
+ S ( X + X + X +
4
(3.26)
()
where s can be changed to whichever lower-resolution scheme is used. The

above formula gives a clear impression about the structural change that multipole extension means to SCC.
It may also be worth taking a look at what approximations the semioperator
s to operator ones mean. The (3.16) energy expression and the respective
(3.18) Hamiltonian implies an approximation
;
1
(; + ; ) S .
2
(3.27)
The (3.22) workaround of computing multipole contributions with semioperator

s goes beyond the above approximation by doubling it with the operator and
quasiclassical sides symmetrized, and then subtracting the fully quasi-classical
energy:
1
(; S + ; S + S ; + S ; )
2
1
S ( + + + )S
4
where the definition of ; s is self-evident. With a

1
2 Exc
; =
(x)(x)
+
(y)(y)dydx
|x y| (x)(y)
(3.28)
(3.29)
x y
semi-simplified two-electron integral representation of the semioperator s, the

(3.28) approximation of the exact s can be seen quite accurate and able to
produce even HF exchange, just like exact ones, yet at a much lower cost (much
smaller, but still large integral tables). The quasiclassical side can be even more
simplified, an l-indexed resolution seems to be more than satisfactory.
82
3.3. SUPPLEMENTARY MATERIAL TO SCC
3.3
3.3.1
Supplementary material to SCC improvements

Mulliken charge density operators
At the investigation of on-site Coulombic energy, it was worth introducing the

1
PA =
ni (
ci, ci, || + ci, ci, ||)
(3.30)
2
{A},,i
atomwise Mulliken density operator. It further implied the A atomic Mulliken

charge density and the qA atomic Mulliken charge such that
(3.31)
qA = TrPA = x| PA |x dx = A (x)dx
x
With the atomic Mulliken density operators and using

1
P =
P + P =
PA ,
2
(3.32)
derivation of off-site SCC Coulomb energy with atomwise Mulliken density operators would also have been a bit shorter in 1.2.2, but it would have hidden
some intuitive details.
To be seen in a broader view, the Mulliken density operators give a decomposition of the P density operator using a complete set of projections in our
LCAO space.
1
)
PA =
( P + P
(3.33)
2
{A}
where
=
|| S 1
(3.34)
and
= 1.
(3.35)
Similarly to the general refinement of SCC equations to orbital hardnesses

and up to orbitally resolved s, the Mulliken density operators can also be
refined to reflect this structure. An orbitalwise density operator is
1
1

P = (
ni (
ci, ci, || + ci, ci, ||) .
(3.36)
P + P ) =
2
2 ,i
With this,
PA =
P .
(3.37)
{A}
Orbitalwise Mulliken density operators imply orbitalwise Mulliken charge densities and charges, and they may be useful at the investigation of the totally
resolved s in section 3.2.3.
As an intermediate step between the atomwise and the totally resolved
regime, there can be Mulliken densities defined which are resolved only up to
the angular momenta of orbitals. They are instruments of studying the almosttotally refined SCC in section 3.2.2.
83
3.3.2
Properties of the generalized functions
Shape
For the original shape, see [13]. As the new shapes that we propose in this
paper are semi-empirical, there must be analytical shapes (similar to the original
ones in [13] or not) fitted to the semi-empirical data in order to be able to carry
out derivations on them, if one decides to use them.
Either in the interpolated or in the semiempirical case, the multipole expansion of SCC energy requires derivatives of basic (00) as higher-order s
(00)
(see (3.8)). Since AB (ZA , ZB , xA , xB ) depends on its xA and xB arguments
as a spherically symmetric function of r = |r| = |xA xB | (for brevity we use
r instead of xA xB and we drop the atom-label arguments in the following
equations), the higher-order s are quite simple, e.g.
(10) (r) =
r (00) (r)
r
r
(3.38)
and the second derivative by components

(20)
(r) =
1 (00) (r) 1 r r 2 (00) (r)

+
.
2 r
r
2 r2
r2
(3.39)
The above formulae can be seen in (3.5) applied to the 1r long-range limit potential in the place of a proper (00) .
As we can see from the above too, higher-order s are matrices according to
their multipole indices. (mn) is a matrix of order m + n, e.g. (10) is a vector
(not mentioning the atomic centre indices of s when they are listed as AB ).
Identities
As every is a AB = (ZA , ZB , rAB ) = (ZA , ZB , |xA xB |), it must hold
that
AB
AB
=
.
(3.40)
xA
xB
From physical reasons (from the fact that nothing must depend on the listing
order of atoms or from the second-derivative nature of s) it is also obvious
that
(mn)
(nm)
AB = BA .
(3.41)
In the multipole expansion of SCC energy, the s were coefficients of a Taylor
series, therefore
(mn)
(m+1,n)
(m + 1)AB
AB
xA
(mn)
AB
xB
(m,n+1)
= (n + 1)AB
(3.42)
and hence
(mn)
AB
= (1)
mn (nm)
AB .
(3.43)
mn (mn)
BA .
(3.44)
From the last equation and (3.41)

(mn)
AB
= (1)
Based on the above identities, it is enough to calculate a for each value of

m + n, e.g. (m+n,0) , all the other s of the same level are of the same shape.
84
3.3. SUPPLEMENTARY MATERIAL TO SCC

Calculating molecular SCC energy derivatives
Higher-order s are not only needed for higher-order multipole terms, but also
for calculating the derivatives of SCC energy. As an example, we show how to
calculate a molecular SCC force up to the monopole-dipole order.
The molecular SCC force acting on atom B is the negative gradient of (3.7)
(now up to the monopole-dipole order) with respect to xB , and its part that
contains derivatives is
()
FB, = qB
(00)
AB
qA
xB,
A=B
dB
(10)
(01)
AB
AB
qA qB
dA .
xB,
xB,
(3.45)
A=B
A=B
Writing all the vector quantities by components and using the identities of s,
()
FB, = qB
(01)
AB, qA
A=B
dB,
(02)
AB, qA
A=B
qB
(02)
AB, dA, .
(3.46)
A=B
3.3.3
Transformations of tabulated matrices before use
In current SCC-DFTB implementations, there are two matrices of quantummechanical integrals tabulated for run-time use with DFTB; the Hamiltonian
matrix and the overlap matrix. Using them needs a transformation because the
coordinate system in which the tabulated integrals have been carried out is not
the same as the molecular coordinate system in which the DFTB calculation
takes place (usually, the system of integration is fixed to the line between the
participating two atoms, while the molecular one is quite arbitrary).
The overlap and the position matrix
In this section we recall the usage rules of overlap matrix tables in a formalized
way, and parallel with it, we construct the similar rules for the position operator
matrix.
,
Let us call the molecular basis , , etc. and the integration basis ,
etc.
From the relative orientations of the molecular and the integration coordinate
systems, a linear transformation can be set up between the two:

=
T
(3.47)
in which orbitals of different angular momenta do not mix, i.e. s, p, etc. orbitals
of the molecular basis are linear combinations of integration basis elements of
85

the same respective types. With this known, the elements of the molecular
overlap matrix can be constructed as
T ST
(3.48)
S =
,
similar to the transformation of canned Hamiltonian elements.

In the enhanced SCC scheme, elements of the position operator matrix in
(3.6) must be also stored besides the overlap matrix and the Hamiltonian. Taking them out of the box is a bit more difficult since they are not only depending
on the orbitals with whom they were integrated, but, as a vector-valued operator, they are also directly depending on the coordinate system chosen. If x s
are the molecular coordinates, x s are the integration coordinates and
x =
t x
(3.49)
is the transformation between the two, a position matrix element with origin a
in the molecular system is
(a)
X, =
(x)(x
a )(x)dx
T(x)
t (x b ) + b a (x)T
dx
,
,
x
(b)
+ S (b a )
Tt T
X,
(3.50)
,
,

if the coordinates of origin in the integration reference are b = , t b .

This is the transformation of the stored position matrix elements to the molecular system.
In addition to the transformation rules, we note that in the case of the
position matrix, there are vanishing elements, just like with the overlap matrix,
but the selection rules differ. In the case of the overlap matrix, S was zero
if the magnetic quantum numbers of and were different, non-vanishing
elements had m m = 0 in a coordinate system whose z axis was fixed to the
()
line between the two participating atoms. For X, a similar rule applies, but
the criterion of a nonzero element includes the coordinate; m + m m = 0
is the proper selection criterion with m being 0 for z, 1 for x and -1 for y (this
rule applies to real-valued orbitals, complex ones follow a bit more complicated,
yet similar set of rules).
The fully resolved s
The definition of orbitalwise in (3.14) contains | . . . | as the -dependent
part (the V[] potential and q are independent from the particular choice of ),
therefore it transforms with
=
T
(3.51)
T
86
3.4. MULTIPOLES IN TD-DFTB

of the integration system to the s of the
when we make the transition from s
molecular calculation. Transformation of the index is a bit more complicated
to construct, but since is nothing else than a second derivative of the total
energy with respect to q and q (or at least an approximation of it), =
and therefore the transformations in both indices must obey exactly the same
rules
T
=
TT
(3.52)
T
A more precise derivation of the above transformation rules would be built on

that the -dependent part of the (3.14) definition of is
V[]
=
q
V[] [] (x)
dx
[] (x) q
(3.53)
in which only q is explicitly depending on the particular choice of . Being

q =
1
(P S + P S ) ,
2
(3.54)
it shall follow the transformation rules of a product containing | and |,

similar to (3.51).
Note that the (3.52) transformation rules are not of a well-behaved tensor
of the LCAO space. A tensor transforming with the product of four transformation matrices should have four independent indices, and all of these four
indices should be meaningful. It is supposed, however, that the real four-index
s (cf. these four indices with the four orbitals in integrals of (1.48)) are
diagonal-dominated in their first and second two indices respectively and the
perturbative non-diagonality is approximated with the multipole structure of
SCC. The problem of this approximate transformation of is naturally away
with the operator-like s introduced in 3.2.4
3.4
Extensions of TD-DFTB in line with the

multipole expansion of SCC
In time-dependent DFTB (TD-DFTB, which is a tight-binding approximation

of time-dependent DFT) [28], an improvement of SCC nomenclature comes into
play in two aspects. Here, we show these differences very briefly, investigating
the details shall take place in further research.
The first aspect shown here is the coupling matrix of TD-DFTB4 . For singlet
transitions in the spin domain, the coupling matrix is nothing else than the
second derivative of total energy with respect to the electronic density operator
(which gives also the base of a theoretical non-approximated version of the SCC
Hamiltonian), except that it contains matrix elements between occupied (i, j)
and unoccupied (a, b) KohnSham orbitals:

1
2 Exc
S
j(y)b(y)dydx
Kia,jb =
i(x)a(x)
+
(3.55)
|x y| (x)(y)
x y
4
On the usage and construction of the coupling matrix, see section 3 in [28].
87

As it was shown in 3.2.4, the multipole expansion of two-electron integrals like
above is a valid tool of making their approximation better than the point-like
Mulliken one. The improvement that the multipole expansion of SCC energy
can give to the approximation of the above integrals is also detailed there. As
S
Kia,jb
=
ci, cj, ca, cb, ; ,
(3.56)
,,,
where |i =
ci, | etc. with the ; of (3.26), the structure of the
S
improved Kia,jb
shall be also clear from (3.26).
It must be noted here, however, that the coupling matrix of spin triplet
excitations cannot be improved by the multipole expansion introduced in 3.1.
As it was written in 1.2.3, spin-spin interaction terms in DFTB are confined
2 Exc
. In addition to this,
to on-site terms, because of the short range of (x)(y)
on-site multipole-multipole interactions are zero in our approximation above
the monopole-monopole level because the (1.62) semiempirical calculation of
on-site second-order energy change (the chemical hardness, and its refinements
introduced later) includes every level of multipole-multipole interactions. This
is also why our new higher-order s tend to zero when r = 0. Nonetheless, if one
wants to make improvements in this part of DFTB, there remains a possibility
of using fitted monopole-dipole, etc. interaction strengths, if the underlying
physics needs it. The more promising way of improving the description of triplet
excitations is including spin-orbit interaction in TD-DFTB. Using some kind of
operator-like s (see 3.2.4) would be a major improvement in this field, but at
rather high cost.
The second aspect of making TD-DFTB better by the multipole SCC which
we treat here is the improvement of its oscillator strength expression. The
oscillator strength of a transition from state i to a is proportional to the square
of a transitional dipole moment
2
2

fia
(3.57)
D(ia) =
i|X |a
is the component
where = x, y, z is the index of spatial coordinates, and X
of position operator. The above matrix element can be easily calculated with the
aid of the previously
stored position
operator integral table elements discussed

in 3.3.3. If i = ci, and a = ca, , the dipole matrix elements of (3.57)
become (the sign is only for conventional reasons)
D(ia) =
atoms
(A)
where R
atoms

{A}
atoms
,{A}
| =
ci, ca, |X
{B}
(A)
(A)
R
ci, ca, |(X
) + R
|+
atoms
atoms

{A} B=A {B}
R(A) ) + R(A) | (3.58)

ci, ca, |(X
is the component of the position of atom A. By the definition of

88
3.4. MULTIPOLES IN TD-DFTB

the tabulated overlap and position matrices (see 3.1), we finally get
(ia)
D(ia) =
atoms
1
2
A
(ia)
+ D
=
2
(A)
(A)
(A)
(A)
ci, ca, (X, + S R
) + ci, ca, (X, + S R
)=
{A},
atoms
(ia)
dA
(ia)
qA R(A)
(3.59)
where R(A) is the vector of the position of atom A, and the first, trivial step
is didactically included in order to symbolise that half of the sum over A, B
and , is made with a A B, relabelling of the summand. The
quantities defined in the last step are called transitional Mulliken charges and
dipole moments between states i and a.
89
Summary
In this dissertation, I treat two areas of improving the DFTB method. The
first area of interest is the repulsive part of DFTB energy. There has been
a parametrizer automaton developed that performs a least-squares fitting of
repulsive potentials as linear combinations of basis profiles through a linear
regression of coefficients. It makes fitting of repulsive profiles to multiple fitsystems incomparably easier than before; that is illustrated by several examples
as successful usages of the fitter. The automaton is also able to tune electronic
parameters. To make the fitting process truly automatic as seen from the user
side, the parametrizer does many auxiliary tasks automatically as well, e.g. fitsystem generation. Among the objectives of fits, not only energies can be used,
but several other energetical properties (mainly but not exclusively, derivatives
of energy). User-input weighting of fit targets ensures the most precise fitting
to the most important properties of most important fitsystems selected by the
user.
As a further question regarding repulsive energy, its structure is also treated
here briefly. It is quite clear now that in addition to the pairwise repulsive
profiles, one-body repulsive energy terms can be used too to tune the dissociation
behaviour of DFTB. In addition to short-range pair repulsives and one-body
terms, there must be middle-range pairwise repulsives also used, but fitting
them is not (yet?) feasible.
As a perfect way of building repulsives, calculation of them instead of fitting
is also proposed. After deriving the formulas needed for repulsive calculation and
implementing them, the first results suggest that this scheme is able to be used
to get at least first-glance repulsives for any set of elements that has not been
parametrized with DFTB before, if one optimizes atomic electronic parameters
successfully. Being able to make transferable parametrizations depending on
atomic, not pairwise parameters gives the prospect of producing and extending
parametrizations including large sets of elements with reasonable effort.
The second area of DFTB improvement that is investigated here is making
its SCC part better. A proposition of a multipole expansion of SCC energy
instead of the currently used point-like charges is made as well as a proposition to calculate the effective interaction profiles between atomic excess charges
semiempirically, instead of the current heuristical interpolation. Fully derived
formulae of the above improvements are ready to implement in SCC. Furthermore, improvements to time-dependent DFTB are also introduced briefly.
90
The above enhancements, enhancement proposals and descriptions of them

are introduced by a thorough derivation of current DFTB method in the very beginning of the thesis, which prepares the methodological improvements, among
others, by a formulation of the DFTB equations based on the density matrix.
Through treating improvements to DFTB, a deeper understanding of the nature
of it, including its approximative self-consistency (SCC) could be achieved too.
91
Acknowledgements
It is hardly possible to write down how deeply I am grateful to the colleagues
who made my doctoral research possible and who helped it. First, I thank to
Professor Thomas Frauenheim to give me the opportunity of working at the
Bremen Center of Materials Science as a PhD student. As large parts of the
research has been carried out in working together with Balint Aradi, I must
acknowledge that my results would have been incomparably less without the
help of him.
In addition to the above, cooperations and other discussions canalized a lot of
helpful stimuli and other useful impressions from Gotthard Seifert, Martin Persson, Karl Jalkanen, Grygoriy Dolgonos, Thomas Niehaus, Adriel Dominguez,
Michael Gaus and Tom
as Kubar into my research.
Many thanks shall go to Christof Kohler and Zoltan Hajnal for facilities I
was able to use in Bremen and temporarily in Budapest, respectively.
Last but not least, my thanks go to my family, as well as to my friends in
Bremen.
92
Bibliography
[1] NIST computational chemistry comparison and benchmark database. NIST
Standard Reference Database Number 101, http://cccbdb.nist.gov. Release 15b, August 2011, Editor: Russell D. Johnson III.
[2] B. Aradi, B. Hourahine, and Th. Frauenheim. DFTB+, a sparse-matrixbased implementation of the DFTB method. J. Phys. Chem. A, 111:5678
5684, 2007.
[3] Nathan Argaman and Guy Makov. Density functional theory an introduction. 1999. arXiv:physics/9806013v2.
[4] Zolt
an Bodrog and Balint Aradi. Possible improvements to the selfconsistent-charges density-functional tight-binding method within the second order. Phys. Status Solidi (b), 249:259269, 2012.
[5] Zolt
an Bodrog, B
alint Aradi, and Thomas Frauenheim. Automated repulsive parametrization for the DFTB method. J. Chem. Theory Comput.,
7:26542664, 2011.
[6] Zolt
an Bodrog, B
alint Aradi, and Thomas Frauenheim. Atomic-parameterbased repulsive potentials for the density-functional tight-binding method.
2012. To be published soon.
zek. A new method of approximative calculation of polycentric inte[7] Jir C
grals used in the quantum mechanical study of molecular structure. Mol.
Phys., 6:1931, 1963.
[8] Larry A. Curtiss, Krishnan Raghavachari, Paul C. Redfern, and John A.
Pople. Assessment of Gaussian-2 and density functional theories for the
computation of enthalpies of formation. J. Chem. Phys., 106:10631079,
1997.
[9] P. A. M. Dirac. Note on exchange phenomena in the Thomas atom. Proc.
Cambridge Phil. Soc, 26:376385, 1930.
[10] Grygoriy Dolgonos, Balint Aradi, Ney H. Moreira, and Thomas Frauenheim. An improved self-consistent-charges density-functional tight-binding
(SCC-DFTB) set of parameters for simulation of bulk and molecular systems involving titanium. J. Chem. Theory Comput., 6:266278, 2010.
[11] R. Dovesi, M. Causa, R. Orlando, C. Roetti, and V.R. Saunders. Ab-initio
approach to molecular-crystals a periodic hartree-fock study of crystalline
urea. J. Chem. Phys., 92:74027411, 1990.
93
[12] M. Elstner. SCC-DFTB: What is the proper degree of self-consistency. J.

Phys. Chem. A, 111:56145621, 2007.
[13] M. Elstner, Dirk Porezag, G. Jungnickel, Joachim Elstner, M. Haugk,
Thomas Frauenheim, S
andor Suhai, and Gotthard Seifert. Self-consistentcharge density-functional tight-binding method for simulations of complex
materials prooperties. Phys. Rev. B, 58:72607268, 1998.
[14] H. Eschrig and I. Bergert. An optimized LCAO version for band structure
calculations. Phys. Status Solidi (b), 90:621628, 1978.
[15] Enrico Fermi. Un metodo statistico per la determinazione di alcune propriet`
a dellatomo. Rend. Accad. Naz. Lincei, 6:602607, 1927.
[16] J. Frenzel, A. F. Oliveira, N. Jardillier, T. Heine, and G. Seifert.
Semi-relativistic, self-consistent-charges Slater-Koster tables for densityfunctional-based tight-binding (DFTB) for materials science simulations.
Technical report, TU Dresden, 20042009.
[17] M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Rob,
J. R. Cheeseman, J. A. Montgomery Jr., T. Vreven, K. N. Kudin, J. C.
Burant, J. M. Millam, S. S. Iyengar, J. Tomasi, V. Barone, B. Mennucci,
M. Cossi, G. Scalmani, N. Rega, G. A. Petersson, H. Nakatsuji, M. Hada,
M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima,
Y. Honda, O. Kitao, H. Nakai, M. Klene, X. Li, J. E. Knox, H. P. Hratchian,
J. B. Cross, V. Bakken, C. Adamo, J. Jaramillo, R. Gomperts, R. E. Stratmann, O. Yazyev, A. J. Austin, R. Cammi, C. Pomelli, J. W. Ochterski,
P. Y. Ayala, K. Morokuma, G. A. Voth, P. Salvador, J. J. Dannenberg,
V. G. Zakrzewski, S. Dapprich, A. D. Daniels, M. C. Strain, O. Farkas,
D. K. Malick, A. D. Rabuck, K. Raghavachari, J. B. Foresman, J. V. Ortiz,
Q. Cui, A. G. Baboul, S. Clifford, J. Cioslowski, B. B. Stefanov, G. Liu,
A. Liashenko, P. Piskorz, I. Komaromi, R. L. Martin, D. J. Fox, T. Keith,
M. A. Al-Laham, C. Y. Peng, A. Nanayakkara, M. Challacombe, P. M. W.
Gill, B. Johnson, W. Chen, M. W. Wong, C. Gonzalez, and J. A. Pople.
Gaussian 03. Gaussian, Inc., Wallingford, CT, 2003.
[18] J. Harris. Simplified method for calculating the energy of weakly interacting
fragments. Phys. Rev. B, 31:1770, 1985.
[19] J. F. Janak. Proof that
B, 18:71657168, 1978.
E
ni
= i in density-functional theory. Phys. Rev.
[20] J. M. Knaup, B. Hourahine, and Th. Frauenheim. Initial steps toward

automating the fitting of DFTB Erep (r). J. Phys. Chem. A, 111:5637
5641, 2007.
[21] Christof K
ohler. Ber
ucksichtigung von Spinpolarisationseffekten in einem
dichtefunktionalbasierten Ansatz. PhD thesis, Universitat Paderborn, 2003.
[22] Christof K
ohler and Thomas Frauenheim. Molecular dynamics simulations
of Cfx (x = 2, 3) molecules at Si3 N4 and SiO2 surfaces. Surf. Sci., 600:453
460, 2006.
94
[23] Tom
as Kubar, Zolt
an Bodrog, and Michael Gaus. Density-functional tightbinding parameters for halogen-organic chemistry. 2012. To be published
soon.
[24] Christof K
ohler, Thomas Frauenheim, Ben Hourahine, Gotthard Seifert,
and Michael Sternberg. Treatment of collinear and noncollinear electron
spin within an approximate density functional based method. J. Phys.
Chem. A, 111:56225629, 2007.
[25] Christof K
ohler, Gotthard Seifert, and Thomas Frauenheim. Densityfunctional-based calculations for Fe(n), (n 32). Chem. Phys., 309:2331,
2005.
[26] Andre Mirtschink. Berechnungen von Bindungsenergien zweiatomiger
Molek
ule der Elemente der ersten und zweiten Periode im Rahmen der
DFTB-Methode. Masters thesis, Technische Universitat Dresden, 2010.
[27] Ney H. Moreira, Grygoriy Dolgonos, Balint Aradi, Andreia L. da Rosa,
and Thomas Frauenheim. Toward an accurate density-functional tightbinding description of zinc-containing compounds. J. Chem. Theory Comput., 5:605614, 2009.
[28] Thomas Niehaus. Approximate time-dependent density functional theory.
J. Mol. Struct. THEOCHEM, 914:3849, 2009.
[29] Robert G. Parr and Weitao Yang. Density-Functional Theory of Atoms
and Molecules. Oxford University Press, 1994.
[30] J. P. Perdew, K. Burke, and M. Ernzerhof. Generalized gradient approximation made simple. Phys. Rev. Lett., 77:38653868, 1996.
[31] Dirk Porezag, Thomas Frauenheim, Thomas Kohler, Gotthard Seifert, and
R. Kaschner. Construction of tight-binding-like potentials on the basis of
density-functional theory: Application to carbon. Phys. Rev. B, 51:12947,
1995.
[32] V.R. Saunders, R. Dovesi, C. Roetti, R. Orlando, C. M. Zicovich-Wilson,
N.M. Harrson, K. Doll, B. Civalleri, I. Bush, Ph. DArco, and M. Llunell.
CRYSTAL2003 Users Manual. University of Torino, Torino, 2003.
[33] G. Seifert. Tight-binding density functional theory: An approximate Kohn
Sham DFT scheme. J. Phys. Chem. A, 111:56095613, 2007.
[34] Jose M. Soler, Emilio Artacho, Julian D. Gale, Alberto Garca, Javier Junquera, Pablo Ordejon, and Daniel Sanchez-Portal. The SIESTA method
for ab initio order-N materials simulation. J. Phys., 14:2745, 2002.
[35] L. H. Thomas. The calculation of atomic fields. Proc. Cambridge Phil.
Soc., 23:542548, 1927.
[36] M. Valiev, E.J. Bylaska, N. Govind, K. Kowalski, T.P. Straatsma, H.J.J.
van Dam, D. Wang, J. Nieplocha, E. Apra, T.L. Windus, and W.A. de Jong.
NWChem: a comprehensive and scalable open-source solution for large
scale molecular simulations. Comput. Phys. Commun., 181:14771489,
2010.
95
[37] Carl Friedrich von Weizs

acker. Zur Theorie der Kernmassen. Zeitschrift
f
ur Physik, 96:431458, 1935.
[38] Oleg A. Vydrov and Troy Van Voorhis. Nonlocal van der Waals density
functional: The simpler the better. J. Chem. Phys., 133:244103, 2010.
96

Tight Density Functional Theory PHD Thesis

Uploaded by

Copyright:

Available Formats

Tight Density Functional Theory PHD Thesis

Uploaded by

Document Information

Original Description:

Original Title

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

Tight Density Functional Theory PHD Thesis

Uploaded by

Copyright:

Available Formats

Improvements to the Density-Functional

Tight-Binding method: new, efficient

Dem Promotionsausschuss an der Universitat Bremen, im

Tag der Einreichung 10. Januar 2012

Prof. Dr. rer. nat.

1 An introduction to the DFTB method

2 Improvements to handling repulsive energy

3 Enhancement proposals to SCC-DFTB

The Density-Functional Tight-Binding (DFTB) method of quantum-chemical

Die density-functional tight-binding (DFTB) Methode der quantenchemischen

A DFTB (Density-Functional Tight-Binding) modszer egy olyan kvantumkemiai

What we shall know about DFT now

What is density functional theory?1

Although quantum mechanics brought a deeper-than-ever-before understanding

For a thorough foundation and description of this topics, see [29].

CHAPTER 1. AN INTRODUCTION TO THE DFTB METHOD

(x1 , . . . , xN )dx1 . . . dxN = F () + Vext (x)(x)dx

= E() =F () + Vext (x)(x)dx.

The HohenbergKohn theorems

F (1 ) + (x)V1 (x)dx F (2 ) + (x)V1 (x)dx

F (2 ) + (x)V2 (x)dx F (1 ) + (x)V2 (x)dx

1.1. WHAT WE SHALL KNOW ABOUT DFT NOW

E1 = F () + (x)V1 (x) > E2 + (x) (V1 (x) V2 (x))

E2 = F () + (x)V2 (x) > E1 (x) (V1 (x) V2 (x))

E() = F () + (x)Vext (x)dx

CHAPTER 1. AN INTRODUCTION TO THE DFTB METHOD

The KohnSham equations

density sum of independent, non-interacting electronic orbitals, the so-called

and Vxc is an effective potential arising from everything else, namely

1.1. WHAT WE SHALL KNOW ABOUT DFT NOW

variational problem if we use a zero-temperature occupation of states as the

E() = T + J() + K + Vext

= TKS () + (T TKS ()) + J() + K + Vext (1.18)

whence it will be also evindent that Exc is also a well-defined functional of .

(E() TKS ()) =

where : R R is an ordinary scalar to scalar function, thus making Vxc a

CHAPTER 1. AN INTRODUCTION TO THE DFTB METHOD

are summed up with the appropriate weights to give an exchange energy. As

EKS = Tr(HKS P ) = Tr TKS P + VKS (x)(x)dx = Tr

1.2. DFTB, THE TIGHT-BINDING APPROXIMATION OF DFT

DFTB, the tight-binding approximation of

The original, non-self-consistent-charges (non-SCC) DFTB [33] is a tight-binding

KS | is the Hamiltonian matrix and S = | is the

CHAPTER 1. AN INTRODUCTION TO THE DFTB METHOD

where P = i ni |i i| is the density operator of the electrons, ni being the

is the energy of nuclear Coulomb repulsion.

1.2. DFTB, THE TIGHT-BINDING APPROXIMATION OF DFT

+ Tr Pcores HKS VC Vxc

ni ci, | T + VKS | ci, ,

the matrix elements of KohnSham potential in the sum break up to sums of

since VKS can also be regarded as a sum of VB atomic contributions2 . From

CHAPTER 1. AN INTRODUCTION TO THE DFTB METHOD

which is a sum of three-center integrals (characterized by the A, B, C atomic

= | T + VA + VB | if {A} and {B},

Thus the total non-SCC DFTB energy is a KohnSham energy calculated

1.2. DFTB, THE TIGHT-BINDING APPROXIMATION OF DFT

P density operator up to the first order in P . Note that since

this is equivalent to examining the second-order term in the Taylor series of E