Praca dyplomowa doktorska
In silico studies of electronic structure and bonding
properties of selected thermoelectric materials from
ternary chalcogenides and silicon clathrates groups
Studia in silico struktury elektronowej oraz właściwości wiązań dla
wybranych materiałów termoelektrycznych z grup chalkogenków
oraz klatratów krzemowych
Autor:
Wojciech Szczypka
Dyscyplina naukowa:
Chemia
Opiekun pracy:
dr hab. inż. Andrzej Koleżyński, prof. AGH
Uprzedzony o odpowiedzialności karnej na podstawie art. 115 ust. 1 i 2 ustawy z dnia 4 lutego 1994 r. o prawie autorskim i prawach pokrewnych (t.j. Dz.U. z 2006 r. Nr 90, poz. 631 z późn. zm.): „Kto przywłaszcza sobie autorstwo albo wprowadza w błąd co do autorstwa całości lub części cudzego utworu albo artystycznego wykonania, podlega grzywnie, karze ograniczenia wolności albo pozbawienia wolności do lat 3. Tej samej karze podlega, kto rozpowszechnia bez podania nazwiska lub pseudonimu twórcy cudzy utwór w wersji oryginalnej albo w postaci opracowania, artystyczne wykonanie albo publicznie zniekształca taki utwór, artystyczne wykonanie, fonogram, wideogram lub nadanie.”, a także uprzedzony o odpowiedzialności dyscyplinarnej na podstawie art. 211 ust. 1 ustawy z dnia 27 lipca 2005 r. Prawo o szkolnictwie wyższym (t.j. Dz. U. z 2012 r. poz. 572, z późn. zm.) Za naruszenie przepisów obowiązujących w uczelni oraz za czyny uchybiające godności studenta student ponosi odpowiedzialność dyscyplinarną przed komisją dyscyplinarną albo przed sądem koleżeńskim samorządu studenckiego, zwanym dalej „sądem koleżeńskim”, oświadczam, że niniejszą pracę dyplomową wykonałem osobiście i samodzielnie i że nie korzystałem ze źródeł innych niż wymienione w pracy.
. . . .
Dziękuję mojemu wieloletniemu promotorowi, mentorowi, i opiekunowi naukowemu, profe-sorowi Andrzejowi Koleżyńskiemu, za jego ogromne zaangażowanie, otwieranie nowych horyzon-tów wiedzy i cierpliwość. Współpraca zawiązana jeszcze na pierwszym roku studiów inżyniers-kich umożliwiła zajście długotrwałego procesu rozwoju i kształtowania się we mnie dociekliwości naukowej. Jego wparcie i życiowe doświadczenie pozostanie dla mnie nieocenione.
Serdeczne podziękowania składam pani Monice Gałkiewicz, która wspierała mnie w odkry-waniu piękna chemicznego świata w czasie nauki w III Liceum Ogólnokształcącym im. A. Mick-iewicza w Katowicach. Będąc absolutnym wzorem nauczyciela pokazała mi i moim przyjaciołom chemikom, że “chcieć to móc”.
Pragnę podziękować moim Rodzicom za cały trud, który włożyli w moje wychowanie, za przekazane fundamenty i przepis na szczęśliwe życie. To dzięki Nim całe dzieciństwo czułem się ogromnie zmotywowany do ciągłego poszerzania wiedzy i umiejętności, zawsze poprzez docenianie najmniejszych sukcesów, tak bardzo ważne w życiu dziecka, lecz nigdy nie wymagając niczego ponad naturalne minimum wychowawcze. Swoboda w rozwoju, którą dostałem, choć czasami trudna w konsekwencjach, pozwoliła mi dokonać dobrych wyborów i znaleźć się tu, gdzie jestem. Chyba najlepszym wyrazem podziękowania i pochwały dla Nich za to wszystko, jest stwierdzenie, że sam teraz pragnę w takich samych ramach wychować moje własne dzieci.
Chcę tutaj podziękować moim małym słoneczkom, Arturowi i Izabeli (i D.), za zrozumienie, że tata czasem musiał “zamknąć się w pokoju”. Przede wszystkim jednak jestem im wdzięczny za ich nieograniczoną i bezwarunkową miłość, która daje mi szczęście i jest siłą napędową do życia. Dziękuję moim przyjaciołom i innym ludziom, z którymi wspólne rozmowy, czasem wielo-godzinne, uczyły mnie i uczą nadal patrzenia na świat z innej perspektywy.
Na końcu dziękuję najważniejszej osobie w moim życiu, mojej żonie Monice, którą poznałem dzięki chemii, i to dosłownie. Jest dla mnie wsparciem w każdym moim działaniu i sprawia, że z każdym dniem robię kolejny krok w kierunku samospełnienia. Zarówno przed, jak i w trakcie mojego doktoratu wiele poświęciła, bym mógł osiągnąć to, czego się podjąłem. Trudno wyrazić słowami, jak cenna jest jej obecność we wszystkim, czym żyję.
Praca naukowa finansowana ze środków budżetowych na naukę w latach 2015-2019 jako projekt badawczy w ramach programu pod nazwą “Diamentowy Grant”. Praca została wykonana z wykorzystaniem Infrastruktury PL-Grid.
This research was financed by Polish Ministry of Science and Higher Education from the budget for science in the years 2015-2019, as a research project under the program "Diamond Grant" (grant no. DI2014 019144) and supported in part by PL-Grid Infrastructure.
Contents
Abbreviations ... 7
List of symbols ... 9
1. Introduction ... 13
2. Solid state computational chemistry and principles of thermoelectrics ... 15
2.1. Calculations from first principles ... 15
2.1.1. Born–Oppenheimer approximation ... 16
2.1.2. Self-consistent field theory ... 17
2.1.3. Hohenberg–Kohn theorems ... 18
2.1.4. Kohn–Sham equations ... 19
2.1.5. Exchange and correlation in density functional theory ... 20
2.1.6. Exchange-correlation functionals ... 21
2.2. Electronic structure calculations in periodic systems ... 24
2.2.1. Bloch’s theorem ... 24
2.2.2. Reciprocal space ... 25
2.2.3. Band structure of solids ... 25
2.2.4. Modelling of point defects in solids ... 29
2.3. Quantum theory of atoms in molecules ... 30
2.3.1. Principles of Bader’s approach ... 31
2.3.2. Topological critical points ... 32
2.3.3. Topological atoms ... 34
2.3.4. Chemical bond and its properties ... 35
2.4. Bond valence theory ... 37
2.4.1. The principles of the bond valence model ... 37
2.4.2. Distorted systems ... 38
2.5. Thermoelectricity ... 39
2.5.1. Thermoelectric phenomena ... 39
2.5.2. Microscopic view on thermoelectric effect ... 41
2.5.4. Maximisation of thermoelectric figure-of-merit ... 44
3. Object of study... 50
3.1. Silver chalcogenides ... 50
3.1.1. Binary and ternary chalcogenides ... 50
3.1.2. Structure of AgSbTe2 ... 52
3.1.3. Band gap of AgSbTe2 ... 54
3.1.4. Point defects in AgSbTe2 ... 54
3.2. Silicon clathrates ... 55
3.2.1. Structure of clathrates ... 57
3.2.2. Searching for new silicon clathrates ... 58
3.2.3. Type I silicon clathrates with zinc ... 59
4. Aim of the work ... 60
4.1. Silver chalcogenides ... 60
4.2. Silicon clathrates ... 61
5. Methodology ... 62
5.1. DFT calculations in a solid state ... 62
5.1.1. Full potential method ... 63
5.1.2. Parameters for calculations ... 63
5.1.3. Geometry optimisation ... 65
5.1.4. Calculation of electronic properties ... 67
5.2. Topology of total electron density ... 69
5.3. Bond valence model ... 69
5.4. Transport properties ... 70
6. Silver chalcogenides – results ... 71
6.1. Initial selection ... 71
6.1.1. Structural properties ... 71
6.1.2. Band structure ... 71
6.1.3. Electron density topology ... 72
6.1.4. Transport properties ... 75
6.2. Cation sublattice ordering ... 75
6.2.1. Structure stability ... 78
6.2.2. Electron density topology ... 81
6.2.3. Bond valence model ... 86
6.2.4. Carriers effective mass ... 86
6.2.5. Transport properties ... 88
6.3.1. Analysed structures ... 89
6.3.2. Defects ordering and formation energies ... 93
6.3.3. Bond valence analysis ... 95
6.3.4. Electron density topology ... 98
6.3.5. Transport properties ... 102
6.4. d-block metal doping ... 103
6.4.1. Defect formation ... 104
6.4.2. Band structure ... 105
6.4.3. Electron density topology ... 105
6.4.4. Transport properties ... 109
7. Silicon clathrates – results ... 115
7.1. Preliminary analysis ... 115
7.1.1. Electron density topology ... 115
7.1.2. Band structure ... 118
7.1.3. Transport properties ... 118
7.2. Partial exchange of host by Zn ... 120
7.2.1. Zn location in silicon framework ... 122
7.2.2. Electron density topology ... 123
7.2.3. Band structure ... 125
7.2.4. Transport properties ... 127
8. Summary and conclusions ... 132
8.1. Silver chalcogenides ... 132 8.2. Silicon clathrates ... 134 8.3. General conclusions ... 135 References ... 136 List of Figures ... 148 List of Tables ... 150
Abbreviations
a.u. atomic unit. 36, 64, 65
B3LYP Becke–three-parameter–Lee–Yang–Parr XC functional. 23
BCP bond critical point. 33–37, 72, 75, 81–84, 98, 99, 105, 108, 109, 115, 116, 118, 122, 123, 126, 134, 148, 150
BS band structure. 25, 27, 59, 67, 71, 72, 88, 105
BVM bond valence model. 37, 38, 65, 69, 86, 87, 96, 133, 149, 150
CB conduction band. 27, 29, 46, 47, 54, 72, 88, 105, 118, 134 CCP cage critical point. 33
CI configuration interaction. 20 CP critical point. 31–33, 69, 148
CPA coherent potential approximation. 30
DFT density functional theory. 15, 18, 20, 21, 23, 30, 47 DOS density of states. 67, 72
EDA energy decomposition analysis. 30 ELF electron localisation function. 30, 31, 36 ELI electron localizability index. 36
f.u. formula unit. 52, 53, 78, 80
FBZ first Brillouin zone. 25, 26, 54, 63–65, 67, 72, 88, 105, 118
GGA generalised gradient approximation. 22–24, 54, 67, 72
HF Hartree–Fock. 15, 17, 18, 20
HSE Heyd–Scuseria–Ernzerhof XC functional. 54, 72
IWS irreducible Wigner–Seitz. 69
LAPW linearised augmented plane wave. 62, 63, 148 LAST Lead-Antimony-Silver-Tellurium. 50–53, 74, 103 LCAO linear combination of atomic orbitals. 18 LDA local density approximation. 22–24, 54
mBJ modified Becke–Johnson XC potential. 23, 72, 74, 148 NBO natural bond orbital. 30
NCP nuclear critical point. 31, 33, 34
PBE Perdew–Burke–Ernzerhof XC functional. 22, 54, 118
PBEsol revised Perdew–Burke–Ernzerhof XC functional for solids. 22, 65, 72, 74, 127, 133 PES potential energy surface. 17
PGEC phonon glass–electron crystal. 48, 55, 57, 122 QCT quantum chemical topology. 31
QTAiM quantum theory of atoms in molecules. 30, 31, 34–36, 67, 69, 133, 134 RCP ring critical point. 33, 34
SCF self-consistent field. 18, 65, 67
TAGS Tellurium-Antimony-Germanium-Silver. 51, 53 TE thermoelectric. 13, 14, 48, 50, 102, 132, 134, 135 TEG thermoelectric generator. 44, 45
TPSS Tao–Perdew–Staroverov–Scuseria XC functional. 23 v.u. valence unit. 39, 86
VB valence band. 27, 29, 46, 47, 54, 55, 71, 72, 105, 111, 118, 127 VCA virtual crystal approximation. 30
VSEPR valence shell electron pair repulsion. 30, 31 WS Wigner–Seitz. 69
XC exchange-correlation. 20–24, 54, 65, 71, 72, 127 ZK rule Zintl–Klemm rule. 58, 59
List of symbols
Symbol Description Pages
BO bond order [35]
B bulk modulus [65, 80]
D global instability index calculated within BVM
the-ory [39, 86, 87, 96, 98, 149]
E[ρ] total energy as functional of ρ [19] EH[ρ] Hartree energy as functional of ρ [19, 20]
Exc[ρ] exchange-correlation energy as functional of ρ [19–21]
Eee[ρ] energy of electron-electron interaction as functional
of ρ [20]
E energy [15–17, 19, 23, 64–66,
78, 80, 86, 93–95, 122]
I electric current [39, 40]
L Lorentz number [45]
N (Ω) electron population of a topological atom [35] OS(Ω) formal oxidation state of a topological atom [35]
Q heat [39, 40]
R0 length of a bond of particular type with valence
equal to 1 – empirical BVM parameter [37, 69, 96, 98] RM T size of an atomic sphere (muffin tin) in LAPW
method [63–66]
Rij length of a bond between ith and j th ions [37]
S(Ω) surface of a basin [34]
T [ρ] total kinetic energy as functional of ρ [20] Ts[ρ] Kohn-Sham kinetic energy as functional of ρ [19, 20]
T temperature [39–41, 43–46, 48]
Vef. Kohn-Sham effective potential [19, 20]
V potential [18, 19, 22, 24]
Ω basin (subregion of space) of topological atom [34, 35, 76, 82, 84, 85, 100, 110, 117, 148–150] Πab differential Peltier coefficient between materials a
and b [39, 40]
Π Peltier coefficient [43, 44]
αab differential Seebeck coefficient between materials a
Symbol Description Pages
α Seebeck coefficient (thermopower) [39, 40, 43–46, 70, 75– 77, 88, 90, 92, 102, 109, 111, 112, 114, 118, 120, 121, 127, 129, 131, 148, 149]
β Thomson coefficient [40]
χ elementary basis set function [18, 63]
δk strain of kth bond [38]
δ global strain factor calculated within BVM theory [38, 86, 87, 96, 98, 149]
F Fermi level [27, 53, 105]
g size of a band gap [23, 29, 118, 127]
energy in a meaning of an energy level [29, 43, 54, 63, 71, 105] ~ reduced Planck constant, ~ = 2πh =
1.054571800(13) · 10−34[J s] = 6.582119514(40) · 10−16[eV s]
[16, 19, 29, 36, 37, 42]
κe electronic part of thermal conductivity [43–45, 70]
κl lattice part of thermal conductivity [44, 45]
κ thermal conductivity [40, 44, 45]
G(r) gradient kinetic energy density [36]
He(r) total electronic energy density [36, 37, 75, 82, 84, 115, 116, 123, 126, 148] T (r) kinetic energy density [23, 31]
V(r) potential energy density [36, 37]
V (Ω) volume of a topological atom [34, 76, 110, 117, 126]
V volume [65, 84, 85, 148–150]
p pressure [65]
µc mobility of carriers [45, 46]
µ chemical potential [41, 75, 88, 89, 102, 109, 111, 118, 120, 148, 149] ∇ρ(r) gradient of the total electron density [21, 22, 31, 34]
∇2ρ(r) Laplacian of the total electron density [21, 23, 36, 37, 72, 75,
82, 83, 108, 109, 115, 116, 126, 148]
νi valence of the atom i [38, 39]
ω angular frequency [63]
ω rank of the critical point [32]
φ Kohn-Sham orbital [19, 23, 63]
ρ(r) total electron density [16, 18–23, 30, 31, 35, 36, 65, 67, 72, 75, 81, 82, 98, 99, 108, 109, 115, 116, 118, 123, 126, 148] ρα(r) total electron density of α electrons [31]
Symbol Description Pages σ conductivity [40, 43–46, 70, 75, 77, 88, 91, 92, 111, 113, 114, 120, 121, 130, 131, 148, 149] τ relaxation time [44, 70, 75, 77, 88, 91, 92, 103, 113, 114, 120, 121, 130, 131, 135, 148, 149]
εemf electromotive force [39]
ε electric field [40]
ϕ molecular orbital [17, 18]
% resistivity [40, 44]
ς signature of the critical point [32]
ϑ atomic orbital [26, 27]
B magnetic field in vector form [42]
E electric field in vector form [40, 42, 43] J current density in vector form [40, 43]
Ψ wave function [15–17, 24, 25]
k wave vector (k-vector) [24–26, 29, 42–44, 54, 63–67, 70, 105]
pcryst. crystal momentum [42]
q heat flow density vector [40, 43] b
F Fock operator [17, 18]
b
H Hamiltonian operator [15–17]
b
T kinetic energy operator [16, 17] b
V potential energy operator [16, 17] b softness of a bond of particular type – empirical
BVM parameter [37, 69, 96]
c global charge transfer index [35, 74, 76, 85, 117, 126, 150]
di discrepancy of sum of bond valences of ith ion from
its formal valence [39]
e elementary charge, e = 1.6021766208(98)·10−19[C] [42, 45, 46] h Planck constant, h = 6.626070040(81)·10−34[J s] = 4.135667662(25) · 10−15[eV s] [46] j current density [40] kB Boltzmann constant [41, 46, 48]
m? effective mass of carriers [29, 46] me rest mass of an electron, me = 9.10938356(11) ·
10−31[kg]
[16, 17]
n carriers concentration [45, 46]
q(Ω) partial charge of a topological atom [35, 74, 76, 82, 101, 109, 110, 117, 126, 150]
Symbol Description Pages sexpk valence of kth bond calculated from eq. 2.60 based
on experimentally measured bond length [38] stheok valence of kth bond calculated from theoretical
con-siderations based on formal ionic oxidation states [38] sij valence of a bond between ith and j th atoms [37–39]
v velocity [42, 43]
1. Introduction
The great development of human civilization in last two centuries, including population growth and industrialisation, results in tremendous increase in energy demand. According to data presented by International Energy Agency, global consumption of energy increased by ca. 60 % between 1990 and 2016 [1], and the average annual growth is about 2 %. At the same time, two thirds of all source energy input is actually wasted due to inefficiencies of conversion, as can be concluded from annually published energy flow charts by Lawrence Livermore Na-tional Laboratory, which consider energy consumption in United States [2]. The highest rate of rejected energy applies to two main segments for primary energy use – electricity generation and transportation.
Currently, most of the consumed energy originates from non-renewable sources, such as coal, natural gas or petroleum. Since the resources of these fuels are limited and their consumption still accelerates, huge increase of the role of energy production from renewable sources is required to avoid global energy crisis (apart from another important aspect of non-renewable sources consumption, i.e. greenhouse gases generation).
Sustainable energy management includes the increase of use renewable sources of energy as well as decrease of wasted energy. Nowadays, plenty of various technologies can be used to exploit waste heat, such as utilisation of heat pumps [3], to-cooling technologies [4], as well as heat-to-power generators with Organic Rankine cycle [5] and many others, including thermoelectrics [6].
The history of discovering the thermoelectric (TE) phenomena starts in XVIII century. De-spite the fact, that the effect of direct conversion of heat into electricity was primarily observed and described by Alessandro Volta (and in the literature one can find some evidences of even ear-lier observations [7]), it was named after another scientist who rediscovered this effect nearly 30 years later – Thomas Johann Seebeck. Since he observed the movement of a compass needle under the temperature difference between two joints of two wires made from different metals, he called this phenomenon thermomagnetic effect, as he did not recognised that magnetic field responsible for needle deflection was induced by an electric current in wires. The term "thermoelectricity" was later coined by Hans Christian Ørsted.
In XX century thermoelectrics started to be applied in some sophisticated niches, such as the most recognisable radioisotope TE generators for space industry needs [8] as well as nuclear powered pacemakers [9]. Since that time, some of the materials used in first applications are
still studied, but many new materials and novel concepts were developed, resulting in significant enhancement of devices efficiency, tailoring for particular working temperature range as well as environmentally-friendly non-toxic materials-based devices for sustainable technologies.
Although TE devices are still not used in large-scale applications, the development in this area is strongly dynamic in last decades, what allows one to expect a rapid increase in the importance of TE technologies on the market.
Nowadays, in silico experiment plays indispensable role in changing boundaries of knowl-edge. Quick development of theories, models and methods together with continuously increasing computing efficiency of contemporary processing units allows to conduct calculations not only agreeable with results of a real experiment, but also extend these results with additional value.
Theoretical modelling can be used to explain phenomena observed experimentally on the most primary level and shed new light on processes, the nature of which remained so far unknown. Typically low-cost in-silico experiment employed to cross-analyse wide variety of systems or conditions allows to reduce a number of real samples, syntheses, tests or trials by excluding least promising configurations. Finally, theoretical modelling enables to explore hitherto unknown research areas and predict new systems as a result.
In this thesis materials from two groups known as promising from the point of view of their thermoelectric properties – rocksalt structured chalcogenides and silicon clathrates – are studied within the density functional theory together with several additional models in order to get deeper insight into their characteristics important for materials scientists as well as chemists. In the Chapter 2 short description of theories and models used for calculations is provided. Chapter 3 contains a review of available literature about the systems analysed within this thesis to bring the current state of knowledge about the subject. The purpose for conducting a research described in this work is formulated in the Chapter 4. In the next chapter details of applied methodology are provided. Chapters 6 and 7 contain description of the results obtained from the calculations. Finally, conclusions are presented in the Chapter 8.
2. Solid state computational chemistry and principles of
thermoelectrics
Although fundamentals of quantum mechanics and quantum chemistry which determined current shape of these sciences, namely the Schrödinger’s equation and Hartree–Fock (HF) theory, were developed in 1920s and early 1930s, the huge branch of quantum chemistry devoted to solid state, so important from the point of view of materials scientists, has begun growing fast much later. This was caused in part by the fact that HF theory fails in description of solid metals due to neglecting electronic correlation, and thus results in prediction of materials treated as electron gas-like systems to be insulators [10]. As basics of solid state physics even nowadays are taught mostly in metals-focused schemes [11, 12] and historically solid state quantum theory is linked with metals electronic structure, such a drawback of HF theory resulted in temporary "separation" of molecular- and solid state-related chemistry (despite the fact, that HF theory can be actually well carried out for periodic systems [13]).
A real milestone for solid state quantum theory was the development of density functional theory (DFT), which currently serves as a starting point of nearly every theoretical analysis devoted to solid state. This thesis is no exception to this rule and thus concise description of DFT basics is given in the following section.
2.1. Calculations from first principles
Schrödinger equation [14] is known since 1926. In the early stages it has become the fun-damental equation of quantum mechanics and remains so today, although it is not the only approach to describe quantum mechanical systems [15]. Schrödinger equation universally applies to any system consisting of electrons and nuclei [10]. In its most general, time-independent form, it is given as:
b
HΨ = EΨ (2.1)
In the equation Ψ is a wave function of an analysed system,Hb is Hamilton operator and E is a total energy (which components depend on the particular definition of Hamiltonian operator) of this system and an eigenvalue of Ψ.
Hamiltonian operator is an operator corresponding to the sum of all components of energy in a system, including all kinetic and potential terms for all particles a system consists of. In the case of a single particle system it is simply the sum of operators related to kinetic and potential energies:
b
H = bT + bV (2.2)
Wave function Ψ is a function of the variables, which number is equal to the number of system’s degrees of freedom, and giving a complex number as a result. Such function intrinsically cannot be measured (it is not an observable). According to the Born rule [16] the square of wave function’s absolute value can be interpreted as a probability density of finding a particle in elementary space – electron density ρ(r) in the case of Ψ describing electron – and can be measured experimentally:
ρ(r) = |Ψ(r, t0)|2 (2.3)
2.1.1. Born–Oppenheimer approximation
Considering a system of nuclei and electrons, Hamiltonian operator in non-relativistic form is a sum of kinetic and potential energies present in a system:
b
H = bTn+ bTe+ bVne+ bVee+ bVnn (2.4)
Indices in above equation state for origin of particular operator term, i.e. total Hamiltonian operator is a sum of kinetic terms related to nuclei (n) and electrons (e) together with potential terms related to particle-particle interactions.
Born-Oppenheimer approximation [17] allows to separate motions of electrons and atomic nuclei in a system and solve Schrödinger’s equation in two subsequent steps. In the first step only electronic part of wave function is solved:
b
He(r, R)Ψ(r, R) = EeΨ(r, R) (2.5)
Variables r and R denote electronic and nuclear coordinates, respectively. Electronic Hamil-tonian operator consists of kinetic energy operator of electrons and potential energy operators of particle-particle interactions: b He = bTe+ bVne+ bVee+ bVnn (2.6) where: b Te = − X i ~2 2me ∇2r i (2.7)
b Vne= − X i X j Zie2 4π0|Ri− rj| (2.8) b Vee= X i X j>i e2 4π0|ri− rj| (2.9) b Vnn= X i X j>i ZiZje2 4π0|Ri− Rj| (2.10)
Here Zi denotes the atomic number of ith nucleus and me is a rest mass of an electron.
From eqs. (2.8) and (2.10) one can find thatHbedepends only on nuclear positions (R) and does not require information of their momenta. Because of that eq. (2.5) can be solved for chosen R. Solving electronic Schrödinger equation for a set of various R allows to obtain electronic energy eigenvalue Ee as a function of R; this function is called potential energy surface (PES).
In the second step Schrödinger equation considering nuclear motion can be solved:
h b
Tn+ Ee(R)
i
Ψ(R) = EeΨ(R) (2.11)
The transfer from general Schrödinger equation (eq. (2.1)) to its electronic form (eq. (2.6)) requires actually few mathematical operations, including transformation of Hamiltonian operator to the centre of mass system, which introduces mass-polarisation term, which is further neglected, as well as adiabatic approximation restricting the total wave function to single electronic surface. Born–Oppenheimer approximation introduces then a small error, which in majority of systems is insignificant, especially for systems with heavy nuclei [18].
The process of finding the optimum geometry of any system of nuclei and electrons (i.e. a relaxed state of a system, where no unbalanced forces acting on nuclei are present) simply utilises this approximation. It requires first solving an electronic Schrödinger equation. Then, the result is used to calculate forces acting on particular nuclei in the system. Depending on optimisation method employed, nuclei are moved a little and for such new fixed geometry equation for electrons is solved again. This process is repeated until assumed convergence of a result is reached.
2.1.2. Self-consistent field theory
The problem of analysing many-body systems (systems consisting of three or more interacting particles) is that Schrödinger equation cannot be solved analytically. HF method allows to split N-body problem to N single body problems, by introducing the approximation that N-body wave function can be represented by Slater determinant (in case of fermions, e.g. electrons) or by a permanent (for bosons) of N functions.
Hartree–Fock one-electron wave function can be found by solving the equation:
b
Here ϕi are one-electron orbitals and Fb is a Fock operator, which includes kinetic energy term of an electron, Coulomb interactions and exchange operator.
Set of these one-electron equations can be solved iteratively – based on the results of previous iteration ϕ functions are corrected and a new set of HF equations is then solved, until results become self-consistent (subsequent iterations do not change the results within defined conver-gence). This method is called self-consistent field (SCF) method. Such approach utilises also the variational principle, which states that every approximate wave function needs to have an energy that is equal (in the perfect case) or above the exact energy of analysed system.
The representation of ϕ function can be either numerical (defined on a set of points) or analytical. In practice, the former is limited to small and high symmetry systems due to low computational efficiency. Nearly all calculations express ϕ in a form of some defined basis set, i.e. set of known functions, such as Gaussian functions. In this manner every orbital is defined as a linear combination of several basis functions χ:
ϕi =
X
b
cbiχb (2.13)
Frequently, functions used in such representation are called atomic orbitals (only by conven-tion, as generally they are not solutions to atomic HF calculations) and this method is called linear combination of atomic orbitals (LCAO).
2.1.3. Hohenberg–Kohn theorems
Soon after Schrödinger published his work, Llewellyn Thomas and Enrico Fermi [19, 20] cre-ated a quantum mechanical model of the electronic structure for many-body systems. Completely separated from wave function equations, this model was formulated on the basis of electron den-sity. However, because of many imperfections and low accuracy in quantitative predictions, it was not widely applied in the form at that time.
Nevertheless, many researchers have become interested in a concept of electron density-based description of quantum mechanical models. Nearly forty years after Schrödinger equations, Ho-henberg and Kohn have formulated two theorems [21], which are now fundamental theorems of DFT. The first one puts electron density as an equivalent to wave function in the description of a quantum system in a stationary state [22]:
The external potential V (r) in a many-body system at stationary state is a unique functional of electron density: V (r) = F [ρ(r)]
With simple reductio ad absurdum proof of this theorem one can deny the possibility of existing another potential V0(r) being functional of the same electron density at stationary
state.
The second Hohenberg–Kohn theorem is an analogy to the quantum mechanical variational principle:
If F [ρ(r)] is known functional of electron density, then, to determine energy and electron density distribution in a stationary state requires only minimisation of a functional of three-dimensional function of density distribution for given potential V (r).
This means, that energy of a system in a stationary state determined using the "real" func-tional is the lowest possible energy in this system. Finding the solution for a given system one can assume, that new electron density distribution, which gives lower than previous energy from functional F [ρ(r)], is closer to the real stationary state distribution. The main problem in this approach is that the exact definition of the functional F [ρ(r)] is unknown.
2.1.4. Kohn–Sham equations
Hohenberg–Kohn theorems gave the theoretical basis to the electron density-based description of quantum mechanical models. Transfer of these theorems to a practical use was formulated by Kohn and Sham as a re-definition of Schrödinger equation. The general form of such equation for a wave function in a form of single Slater determinant created from Kohn–Sham (KS) orbitals (which are analogues to Hartree–Fock orbitals) is given as:
−~ 2 2m∇ 2+ V ef.(r) φi(r) = Eiφi(r) (2.14)
In this equation effective potential (Vef.), often called Kohn–Sham potential, is a local
poten-tial in which electrons are moving without mutual interactions. Similarly to the Hartree–Fock method, a system of N interacting electrons (described by one equation) is replaced by a set of N non-interacting electrons (described by N equations in a form of eq. (2.14)).
In analogy to Born rule (eq. (2.3)) total electron density is given as sum of squares of absolute values of all Kohn–Sham orbitals:
ρ(r) =
N
X
i
|φi(r)|2 (2.15)
The total energy of a system in given in a form of electron density functional:
E[ρ] = Ts[ρ] +
Z
dr Vext.(r)ρ(r) + EH[ρ] + Exc[ρ] (2.16)
The first component is the kinetic energy of a system, defined using KS orbitals:
Ts[ρ] = N X i=1 Z dr φ∗i(r) −~ 2 2m∇ 2 φi(r) (2.17)
EH[ρ] is Hartree energy, which stands for Coulomb electron-electron interaction, including
self-interaction, i.e. interaction of particular electron with a potential derived from its own elec-tron density. It is equal to:
EH[ρ] = e2 2 Z dr Z dr0 ρ(r)ρ(r 0) |r − r0| (2.18)
The second term of eq. (2.16) is related to external potential Vext. acting on a system,
which corresponds to electron-nuclei interactions in its simplest form. The last term is exchange-correlation energy.
KS equations require an expression of total energy of a system in a form of energies projected on particular orbitals, so that an effective potential in eq. (2.14) is given in a following form:
Vef.(r) = vext.(r) + e2
Z ρ(r0) |r − r0|dr
0+δExc[ρ]
δρ(r) (2.19)
The last term of eq. (2.19) is an exchange-correlation potential:
vxc(r) ≡
δExc[ρ]
δρ(r) (2.20)
The exact form of exchange-correlation (XC) potential (as well as corresponding Exc[ρ]) in
the Kohn–Sham equations is unknown. The problem of formulation of Exc[ρ]is still a challenge
for theoreticians developing DFT. Currently, a variety of XC functionals can be found, which are more of less accurate approximations of a real unknown functional. Some of the developed exchange-correlation functionals are supposed to be suitable for a wide range of systems, some other are dedicated to particular problems. Two following sections are devoted to give a short overview on currently most widely used functionals.
2.1.5. Exchange and correlation in density functional theory
XC energy includes correlation (C) and exchange (X) interactions between electrons. XC functional within Kohn–Sham equations consists of the following terms:
Exc[ρ] = T [ρ] − Ts[ρ] + Eee[ρ] − EH[ρ] (2.21)
Here, T [ρ] is an exact kinetic energy and Eee[ρ] is an exact electron-electron interaction,
expression of both terms are unknown.
The name of electronic correlation was given to mere Coulomb repulsion between electrons in a system. Because of the dynamical character of electrons, description of their correlated motion is actually very challenging for real systems. Several methods developed on the basis of HF method (often called as post-Hartree–Fock) are now frequently used to account for electron correlation, such as configuration interaction (CI) or Møller–Plesset perturbation theory.
The exchange between electrons is related to their spins. According to quantum mechanics laws a wave function for fermions needs to be antisymmetric to the mutual exchange of their components (electrons), what in the case of HF method intrinsically results from mathematical properties of properly formulated Slater determinant. In DFT methods this effect needs to be additionally included.
Figure 2.1: A general classification of exchange-correlation functionals based on physical terms included. Starting from the bottom, complexity of the Exc[ρ]increases together with accuracy of
reality description. "Heaven" is an ideal XC functional, perfectly describing reality. On the right side of the ladder lists of included physical terms are included: electron density ρ(r), its gradient ∇ρ(r)and Laplacian ∇2ρ(r). [26]
When two electrons of the same spin interact with each other (α−α or β−β) an exchange hole appears (at small distance between these electron) due to Pauli exclusion principle (in a system no two identical fermions can occupy the same quantum state). In contrast, if two electrons have different spins, their interaction occurs with no exchange hole.
2.1.6. Exchange-correlation functionals
In last several decades many various approaches to include exchange and correlation in multi-electron systems have emerged. A wide range of available XC functionals implemented in a variety of software for DFT calculations gives many capabilities to scientists today. At the same time, selection of optimal XC functional based on its level of theory and included terms (such as corrections to account for dispersion forces [23]), which will give reasonable results of calculations in specified time, is an inevitable step at the beginning of any in silico research.
Depending on the chosen level of theory, time of computation in ab initio methods for N-particle systems scales up to O(N7), in the case of post-Hartree-Fock perturbative methods.
In case of most DFT methods time of computation scales to index O(N3), but in some cases
scalability can be different [24]. Exponential nature of computation time increase puts limits to size of a system possible to analyse.
Classification of XC functionals is generally done based on physical terms included in a func-tional, according to Perdew et al. [25]. Because of the character of such classification (considering that increase of the degree of complexity goes together with the increase of accuracy) it can be presented in a form of Jacob’s ladder (as the biblical theme), like in a fig. 2.1.
Classification based on physical terms does not include one important aspect about various XC functionals – if given functional is empirical or not. Non-empirical functionals are created by including additional pure physical terms and constraints (real ab initio calculations are carried out using non-empirical functionals, then). In this case, the higher level of theory, the better results should be expected, but in real calculations this is not always true.
This situation changes in the case of empirical functionals. These models are, in general, parametrised according to some experimental data, to increase compliance of calculations with real, measured properties. This may provide some bias to the physical description of a model. Nevertheless, empirical functionals are useful because of their relatively good efficiency and ac-curacy, when employed to a suitable problem.
Local density approximation
The simplest approximation of a real Kohn–Sham XC functional is local density approxima-tion (LDA). This approach exploits properties of homogeneous electron gas, where dependence of XC energy to assumed gas density is well known (exchange is expressed in a simple analytical term, whereas correlation is approximated). In LDA, potential of exchange and correlation is assumed to be equal to XC potential of homogeneous gas of the same electron density, what gives the definition of a functional as follows:
ExcLDA[ρ(r)] = Z
ρ(r)Vxchom.g.[ρ(r)]dr (2.22) Here, it is assumed that in case of perfectly homogeneous electron density energy of XC should be equal to this part of energy of homogeneous electron gas. This implies that in cases of structures, where electron density distribution is close to uniform (changes in space only a little), this assumption should give good results. Metallic materials have frequently very uniformly distributed electron density over crystal space, and thus LDA approximation fits well. However, in most of other systems (such as isolated molecules, strongly covalent materials, etc.) electron density distribution significantly varies in space and this approximation is insufficient.
Generalised gradient approximation
More complex approximation of ideal XC functional are based on generalised gradient ap-proximation (GGA). The assumption here is simple: as real systems are characterised by inho-mogeneous electron density distribution, one needs to include information about variability of the electron density, what allows to get more flexible form of a functional. In this case exchange-correlation potential is based on ρ(r) as well as ∇ρ(r):
ExcGGA[ρ(r)] = Z
ρ(r)VxcGGA[ρ(r), ∇ρ(r)]dr (2.23) Some of the GGA functionals reasonably well describe properties of XC holes, for example widely used Perdew-Wang 91 (PW91) [27] and Perdew–Burke–Ernzerhof XC functional (PBE) [28]. An important variation of PBE is revised Perdew–Burke–Ernzerhof XC functional for solids
(PBEsol) [29], which was specifically designed for application in solid state calculations. Em-ploying this functional often results in well defined structure of solids, close to experimentally observed structures. This implies that electron density distribution in a structure, where inter-atomic distances are accurate, is a good starting point for further analysis of ρ(r), e.g. within topological methods.
Functionals of higher level of theory
Next step in the classification presented in the fig. 2.1 is given to meta-GGA functionals, which additionally make use of Laplacian of ρ(r). In fact, kinetic energy density derived from occupied KS orbitals includes the same physical information, as Laplacian of an electron density, so that, in meta-GGA functionals T (r) can be used instead of ∇2ρ(r):
T (r) = 1 2
X
i
|∇φi(r)|2 (2.24)
The example of meta-GGA XC functional is Tao–Perdew–Staroverov–Scuseria XC functional (TPSS).
Next level of functionals makes use of exact term for exchange:
Ex[ρ(r)] = − 1 2ρ(r) Z d3r0| P iφ∗i(r0)φi(r)|2 |r − r0| (2.25)
Such functionals are named hyper-GGA. This type of XC functional is highly popular in quantum chemical computations for molecules using localised basis sets. One of the most popular examples is hybrid Becke–three-parameter–Lee–Yang–Parr XC functional (B3LYP). It contains three empirically determined parameters adjusted to give accurate results for small molecular systems. As this functional does not satisfy condition of the uniform density limit [30], it does not give good predictions in calculations for systems characterised by high uniformity of electron density distribution, i.e. most solids, especially metallic materials.
Exchange-correlation functionals and band gap size
Band gap size is an important property of semiconducting materials from the point of view of their application, i.a. as thermoelectric materials (more details are given in section 2.5.4). This is a weak point of most of typically used functionals. Calculations using LDA or GGA functionals frequently results in underestimated gap, or even no gap [31].
Variety of methods for improving estimation of g are described in a literature. One of them
is Hubbard model [32], which introduces additional parameter for electron-electron interactions, sometimes used as a correction for popular XC functionals. The GW approximation [33] is an advanced approach to quantum chemical computations which is significantly more CPU time-consuming comparing to typical methods, however gives generally much better results, i.a. im-proved estimation of a band gap.
Relatively efficient method for better g estimation is employing modified Becke–Johnson XC
accurate results of band gap size [34]. However, as the definition of XC energy functional in this method does not exists, only XC potential, it needs to be used on the basis of results obtained from LDA or GGA computations.
2.2. Electronic structure calculations in periodic systems
Molecular systems consist of a limited number of nuclei and electrons, which need to be evaluated at a time of computations. In last decades significant increase in efficiency of computers allowed to shift upper limits of system size possible to analyse within ab initio methods, however it is still no more than a couple of thousands of atoms (with a very poor quality of used methods). Large systems, like molecules of proteins, or complex supramolecular models require using semi-empirical methods to provide reasonable results, such as using of pseudopotentials (effective potentials defined analytically or numerically as simplified form of fully ab initio potentials) or classical potentials.
Crystalline materials, however, consist of a number of interacting atoms at the level of Avo-gadro number, if one considers single grain of a crystal. Calculations for such systems using the same methods as for molecular systems are simply impossible.
When describing structure of a crystal, one needs to simplify it. Based on periodically repeat-ing orderrepeat-ing of atoms in a crystal, it is generally enough to define the smallest repeated fragment – unit cell. This way one neglects defects chaotically distributed in a crystal, but such approach allows to limit the number of atoms to countable, just by adding conditions of periodicity, which can be also exploited within ab initio methods.
2.2.1. Bloch’s theorem
In calculations for periodic systems, Bloch’s theorem is utilised to generate wave function of a crystal, consisting of crystal orbitals. In an infinite crystal, potential needs to maintain the same periodic character as crystal lattice, since subsequent unit cells are indistinguishable. So that:
V (r + R) ≡ V (r) (2.26)
Here R is a translational vector consistent with periodicity of a crystal. Periodic character of a potential puts constraints on a wave function, which need to be Bloch’s wave:
Ψ(k, r + R) = eikRu(r) (2.27) Function u(r) needs to have the same periodicity as a crystal lattice. Functions of type eikr
are called plane waves. Also, the consequence of eq. (2.27) considering case, when R = 0, it can be concluded that:
Ψ(k, r + R) = eikRΨ(k, r) (2.28)
2.2.2. Reciprocal space
The description of electronic structure of solids requires referring to wave vectors, k, which are inherently connected to reciprocal space of a crystal. When defining a vector described in a real space of a crystal, it is given in the coordinates of lattice unit vectors:
R = n1a1+ n2a2+ n3a3 (2.29)
Similarly, vector described in a reciprocal space is constructed from the reciprocal unit vectors:
K = m1g1+ m2g2+ m3g3 (2.30)
Vectors g1, g2 and g3 are defined from lattice unit vectors as:
g1 = 2π a2× a3 a1· (a2× a3) (2.31) g2 = 2π a3× a1 a2· (a3× a1) (2.32) g3 = 2π a1× a2 a3· (a1× a2) (2.33)
The denominators of eqs. (2.31) to (2.33) are essentially the volume of a unit cell.
Since crystal orbitals are dependent wave functions, it is necessary to define range of k-vectors, which is needed to fully describe electronic structure of a material. This fragment of reciprocal space is called first Brillouin zone (FBZ) and is limited to:
−π
a ≤ k ≤ π
a (2.34)
Graphical representation of a Brillouin zone if given in the fig. 2.2.
2.2.3. Band structure of solids
Wave function of electrons in a crystal needs to be described within a FBZ in order to get full information on all energy states in a crystal. Since a reciprocal space is three-dimensional, relation of energy of particular state in respect to three components of k is impossible to present in plain form of a plot. Because of that, a common practice is to choose particular set of k-points in a form of a path through FBZ according to symmetry of a reciprocal space. These plots are called band structure (BS) plots and every single curve (i.e. electronic band) on the plot represents a set of crystal orbitals.
For further consideration simple numbers instead of vectors will be used (k → k, R → R, Ψ → ψ). The simplest periodic system – one-dimensional chain of hydrogen atoms – consists of
(a) L z kz Γ kx ky W V Σ 3 X Λ Δ K M Q 1 X U X S Σ1
©bilbao crystallographic server http://www.cryst.ehu.es
(b)
Figure 2.2: Brillouin zone in a reciprocal space. (a) scheme for "cutting" first Brillouin zone from a reciprocal space: 1. Lattice node choice. 2. Placing lines to closest nodes. 3. Cutting of space with planes perpendicular to these lines and equally distanced to both connected nodes. [35]; (b) first Brilloin zone for F d¯3m real space (so that reciprocal space group is Im¯3m) with denoted special symmetry points and lines. Based on such diagram one can choose a set of k-points to analyse electronic structure of a system. The most characteristic point in FBZ is reciprocal space origin, k = (0, 0, 0), always denoted as Γ. [36]
nuclei equally spaced by the distance a and electrons, where each one was originated from a single 1s atomic orbital. The nth atomic orbital, ϑn, originates from the atom at position R = n · a.
The Bloch’s sum representing a single band over all orbitals is:
ψ(k, R) =X n eikRϑn= X n eiknaϑn (2.35)
Value of k can be limited according to eq. (2.34). For k = 0, the exponential part for each orbital is equal to 1, so that Bloch’s sum will be following:
ψ(0, na) =X
n
ϑn= ϑ1+ ϑ2+ ϑ3· · · (2.36)
The same result will give any k = 2πl/a, l ∈ Z. At the edge of FBZ k = π/a, value of exponential part is 1 for even n and −1 for odd n:
ψ(π/a, na) =X
n
eiπnϑn= −ϑ1+ ϑ2− ϑ3+ ϑ4− ϑ5· · · (2.37)
For k = π/2a and k = −π/2a:
ψ(π/2a, na) =X
n
inϑn= iϑ1− ϑ2− iϑ3+ ϑ4+ iϑ5· · · (2.38)
ψ(−π/2a, na) =X
n
Figure 2.3: Band structure of one-dimensional hydrogen chain described in text, presented for range of k ∈ (0, π/a). Chains of 1s atomic orbitals are presented on the right. Colours red/blue are given to show difference in sign as component of Bloch’s sum.
According to Kramers theorem energy states of at k and −k are the same, so that ψn(π/2a, na)and ψn(−π/2a, na)have the same energy [10]. In order to give them more chemical
meaning one can make linear combination, either by addition, or subtraction:
ψ+≡ ψ(π/2a, na) + ψ(−π/2a, na) ∼ −ϑ2+ ϑ4− ϑ6+ ϑ8· · · (2.40)
ψ− ≡ ψ(π/2a, na) − ψ(−π/2a, na) ∼ ϑ1− ϑ3+ ϑ5− ϑ7· · · (2.41)
Considerations given above are conveniently presented in the fig. 2.3.
Band structure analysis
Band structures typically consist of numerous different bands, which are the result of over-lapping various atomic orbitals with each other. Example calculated band structures for several materials are presented in the fig. 2.4. The rate of atomic orbitals overlapping affects band disper-sion, i.e. range of energy spanned by this band – typically, the stronger overlapping, the greater dispersion. Because of that, different types of materials show different shape of band structure. For strongly separated states between entities (orbitals which do not take part in bonding) dis-persion of band will be small. It is typically found in deep core states and in case of weakly bounded molecular crystals (fig. 2.4(d)). Bands representing bonding states in materials with strong bonding, such as covalent materials, as well as strongly delocalised states in metals, are generally spanning wide energy range (figs. 2.4(b) and 2.4(c)).
The band gap size (if there is any) can be simply read from BS. It is the difference between energy of the lowest state at conduction band (CB) (which is, actually, a set of bands above Fermi level, F) and the highest state at valence band (VB) (analogically, set of bands below
(a) (b)
(c) (d)
Figure 2.4: Example band structures for several types of materials: (a) GaAs; (b) diamond; (c) bcc α-iron (ferrite); (d) crystalline urea. Red horizontal line represents Fermi level – level of highest occupied state. The values of energies on y-axes are relative to Fermi level.
g= min(CB)− max(VB) (2.42)
Since electrical conduction is realised by carriers occupying states close to Fermi level, analysis of top valence bands and bottom conduction bands can give important information about these carriers, i.a. allows to estimate their effective mass, m?. This is the apparent mass of the charge
carrier moving in the periodic potential; carrier "behaves" as if it had a different mass than when it is moving in vacuum. Typically, effective mass is given in relative units compared to particle’s rest mass. The simplest parabolic approximation of band peaks at maximum of VB or minimum of CB (assuming isotropic dispersion of bands around these points) is following:
(k) = 0+~ 2k2
2m? (2.43)
More general case, which includes anisotropy of properties is three-dimensional parabolic dispersion relation: (k) = 0+ ~ 2 2m? x (kx− k0,x)2+ ~ 2 2m? y (ky− k0,y)2+ ~ 2 2m? z (kz− k0,z)2 (2.44)
Conductivity effective mass, which can be used in calculations of carriers mobil-ity/conductivity, is a harmonic mean of m?
x, m?y and m?z: m?cond.= 3 1 m? x + 1 m? y + 1 m? z −1 (2.45)
2.2.4. Modelling of point defects in solids
As it was stated at the beginning of the section 2.2, calculations in periodic systems need to be reduced to particular tractable unit cell. If one introduces a defect in a system, it is then repeated periodically in every subsequent unit cell. In case of complex systems, where unit cell is large, this generally should not bias the results of calculations of related properties. However, for small unit cell distance between repeated defect may be small and thus introduce significant error to calculations due to mutual interactions. In such case it is necessary to define a supercell, i.e.larger part of a periodic system constructed from several neighbouring unit cells, for example twice larger in every direction (2 × 2 × 2). Such approach greatly increases requirements of computations, but allows to analyse more complex cases. Modelling of charged defects typically requires larger supercells in order to reduce long-range Coulomb interaction, but even some neutral defects may introduce significant interaction due to structure distortion [37].
Also, when modelling defects this way one needs to take into account how composition changes. It is extremely hard, or even impossible, to model in conventional way e.g. very small amounts of dopant in a system, because it requires defining large supercells. Because of that frequently ab initio modelling of point defects is carried out for exaggerated dopant amounts [38, 39].
When analysing systems with multiple various defects it becomes necessary to consider all unique configurations for positioning them in a structure, since they may interact with each other. This greatly increases a number of required calculations, but also allows to extract important information about defects mutual influence.
Besides the conventional method for point defects modelling, some other techniques were developed which allow to analyse an influence of introducing particular dopant. One of them is virtual crystal approximation (VCA), which is a simple method for modelling doping by defining a virtual atom which represents interpolated properties of the actual components in a material. In some cases this approach gives comparable results to calculations of ordered supercells [40].
Another example is coherent potential approximation (CPA), which it the method for mod-elling physical properties of substitutionally random systems by configurational averaging. In this technique an affective medium is defined for including the influence of dopant in a form of coherent potential [41]. Because of that employing CPA in calculations enables to model also very small amounts of dopants, beyond the limits of typical ab initio calculations [42].
2.3. Quantum theory of atoms in molecules
From the obvious reasons, electrons within quantum mechanical description of molecules and crystals are indistinguishable. Molecular or crystal orbitals are just mathematical functions defined within a whole space of a system, and thus cannot be directly ascribed, in general, to particular atoms forming this system [43]. This creates a problem how to define an atom as a separate subsystem. A solution to this problem is important from the point of view of every chemist since any chemical debate on matter and related processes uses the concepts of an atom and a bond, such as simple prediction of a molecule geometry employing still very popular valence shell electron pair repulsion (VSEPR) model [44].
These two, and other fundamental concepts in chemistry originated in early years of the XXthcentury, before the development of the quantum mechanics. The inconsistencies in classical
definitions concerning modern knowledge created the necessity to develop appropriate reliable mathematical models.
In the last few decades efforts of many researchers have resulted in a variety of theoretical methods to describe concepts of atoms and bonds they create. As an example, large branch of DFT has been developed for over 40 years called "conceptual DFT" which directly aims to give sharp definitions of many chemical concepts in use, such as electronegativity [45]. In order to extract the idea of a chemical bond, a sequence of natural localised orbital sets binding atomic orbitals and molecular orbitals was described, where the set of natural bond orbitals (NBOs) is the most popular member of this sequence. Another model, energy decomposition analysis (EDA), provides a detailed quantitative description of a bond in terms of several energy components [46]. Great importance in the manner of re-defining chemical concepts needs to be given to methods of topological analysis based on electron density (ρ(r)), which most popular and widely used are quantum theory of atoms in molecules (QTAiM) [47] and electron localisation function (ELF)
[48]. The latter defines a measure of the electron localisation based on kinetic energy density (T (r)) and electron spin density (the difference between total electron density of one spin and total electron density of the second, ρα(r) − ρβ(r)). This method allows to show such elements
in a molecule as core and valence electrons of particular atoms, lone pairs and covalent bonds, and can be treated as a useful tool to visualise ideas of VSEPR theory. Both, QTAiM ale ELF, are components of so-called quantum chemical topology (QCT), which encompasses all methods sharing the same idea of topological analysis.
QTAiM is the model introduced and developed by Richard Bader [47], which gives precise definition of an atom in a system like molecule or crystal and relates variety of properties con-cerning atoms and bonds in this system with electron density and its derivatives. Despite it was not the first theory which was based on the analysis of the real space of a system (this idea was introduced by Raymond Daudel and was called "loge theory" [49]), it became one of the most important models in the field of computational chemistry, because it served very detailed, yet simple and elegant mathematical description with applicability limited only by computation efficiency.
2.3.1. Principles of Bader’s approach
In general, the topology of the electron density is governed by the attractive forces of the atomic nuclei, what results in clear local maxima of the electron density at position of every nuclei; it is the basic feature of every system. In Bader’s approach, topological analysis is based on so-called topological critical points (CPs). They are defined as points with non-zero value of electron density (ρ(r) > 0) and density gradient (∇ρ(r)) being zero vector:
∇ρ(r) = i∂ρ ∂x+ j ∂ρ ∂y + k ∂ρ ∂z = 0 (2.46)
It is important for the right side of eq. (2.46) to be zero vector, as it implies that all partial derivatives of gradient operator need to be equal to 0, not just their simple sum. Local maximum of the electron density at position of an atomic nucleus is one of the types of CPs and is called a nuclear critical point (NCP). In real calculations the approximation related to neglect of a nucleus size results in sharp peak in a point of nucleus position, so that, in fact, a derivative of electron density does not exists at this point. In this manner position of nucleus is not a CP in mathematical sense of eq. (2.46). However, maxima at nuclei points exhibit topological aspect of CP, hence the convention to mark them as NCPs.
In order to determine if particular critical point represents local maximum, minimum or a saddle point of an electron density, one needs to consider second derivatives of ρ(r) in respect to all directions in space, i.e. elements of the Hessian matrix at CP denoted as rcp:
H(rcp) = ∂2ρ ∂x2 ∂2ρ ∂x ∂y ∂2ρ ∂x ∂z ∂2ρ ∂y ∂x ∂2ρ ∂y2 ∂2ρ ∂y ∂z ∂2ρ ∂z ∂x ∂2ρ ∂z ∂y ∂2ρ ∂z2 (2.47)
This matrix is symmetric and real, therefore can be diagonalised. The rotation matrix to transfer the coordinate system (r(x, y, z)− > r(x0, y0, z0)) in order to obtain diagonalised Hessian
is a unitary matrix U, constructed using three calculated eigenvalues from equations Hui= λiui
(i = 1, 2, 3). A transformation U−1HU = Λresults in Hessian in diagonalised form, given as:
Λ = ∂2ρ ∂x02 0 0 0 ∂ 2ρ ∂y02 0 0 0 ∂ 2ρ ∂z02 r0=r cp = λ1 0 0 0 λ2 0 0 0 λ3 (2.48)
where: λ1, λ2 and λ3 are the curvatures of the electron density in respect to x0, y0 and z0.
It is important to notice that trace of the Hessian is independent to coordinate system rotations and its value is the Laplacian of the electron density, ∇2ρ(r):
∇2ρ(r) = ∂ 2ρ(r) ∂x2 + ∂2ρ(r) ∂y2 + ∂2ρ(r) ∂z2 = λ1+ λ2+ λ3 (2.49)
2.3.2. Topological critical points
Classification of CPs is based on values of two parameters, which are ascribed depending on the components of the Laplacian. First parameter is rank (ω), which is a number of non-zero values of λ at this point. Critical points with ω < 3 are mathematically unstable and thus sensitive to small distortions of electron density related to e.g. changes of the positions of atomic nuclei (point disappears of bifurcates). Presence of such point in a system indicates instability of current electron density topology, and can be noticed at saddle points of the chemical reactions, but in equilibrium state only ω = 3 points are observed.
The signature (ς) is a sum of density curvature signs, as defined by the following equation:
ς =
3
X
1
sgn(λi) (2.50)
Using such classification four types of stable (with rank ω = 3) critical points are: – (3, −3) – local maximum of the electron density;
– (3, −1) – maximum of the electron density in two directions and a minimum in a third one, when a coordinate system is rotated to give diagonalised Hessian (like in eq. (2.48));
Figure 2.5: Graphical representation of critical points distribution in a fcc structure of silicon. Green spheres stand for BCPs, orange – RCPs and blue – CCPs. Paths created by small yellow spheres represent bond paths between adjacent silicon atoms.
– (3, 1) – analogically, maximum of the electron density in one direction and a minimum in two other directions;
– (3, 3) – local minimum of the electron density.
Each of these types of critical points represents different element in chemical structure. NCPs mentioned already state for positions of the nuclei and their classification is (3, −3). Bond critical points (BCPs) classified as (3, −1) are placed in space between adjacent atoms and represent bonds. The set of points of local maximum of electron density between two bonded nuclei con-necting their NCPs and including BCP is called a bond path. BCP is a point of minimum at a bond path, and also a point of bond path crossing with zero flux surface between bonded atoms. The third type of CP is ring critical point (RCP) (3, 1), which typically is surrounded by several BCPs and appears inside rings. Point representing local minimum of electron density is called cage critical point (CCP). Example of critical points distribution in a periodic structure is shown in the fig. 2.5.
The number of critical points of each type is strictly correlated. This relationship, depending on the system type, is called Morse sum (in the case of periodic systems) or Poincaré-Hopf relationship (in the case of molecules), and is defined as:
Figure 2.6: Atomic basins in PCl3 molecule. Atomic basin of chlorine is coloured green, surface
of the phosphorus basins are orange. Black curves connecting adjacent NCPs are bond paths.
nN CP − nBCP + nRCP − nCCP = 0 (periodic systems) 1 (isolated molecules) (2.51)
This relationship has to be always satisfied; if the result obtained from calculation does not fulfil this rule, one can be sure this calculation needs to be corrected. It is also accidentally possible to match this relationship by missing e.g. both a RCP and BCP, however probability of such situation is quite small and, in practice, when eq. (2.51) is satisfied it is taken as a proof of calculations consistency.
2.3.3. Topological atoms
Each atomic nucleus in an analysed system is a local maximum of the total electron density. As a natural consequence of this fact, Bader’s model divides space of a system into subregions, Ω, which represents topological atoms. Each subregion has one point of maximum and is bounded by surface of zero flux S(Ω) of the gradient of electron density field. It can be written in the following equation:
∇ρ(r) · n(r) = 0 (2.52)
where n(r) stands for a normal vector to a S(Ω) at point r. In the fig. 2.6 example of a molecule partitioned into atomic basins is shown.
A topological atom is defined as a particular nucleus with belonging atomic basin. In some cases it is possible to observe local maxima of electron density at position not assigned to any nucleus. Such points are denoted as non-nuclear attractor and topologically exhibit the same features as ordinary NCPs, including own basin. As a consequence, they are called pseudoatoms. They are especially important in characteristics of metallic bonding [50].
Analysis of the properties related to topological atoms is an important part of Bader’s model. The most basic example is the atomic volume, V (Ω). In QTAiM, a system partitioned into atomic subregions exhibits atomic additivity, i.e. expected value of particular property of a whole system
equals the sum of individual atomic contributions. This also means that every point in space of a system belongs to a particular topological atom (or, eventually, a non-nuclear attractor).
Integration of total electron density over a whole space of an atomic basin gives the electron population of this atom:
N (Ω) = Z
Ω
ρ(r)dr (2.53)
Comparison the electron population with atomic number (ZΩ) gives us the information about
partial charge of this atom:
q(Ω) = ZΩ− N (Ω) (2.54)
Analysis of the values of partial charges of atoms in a system shows transfer of the electron density between topological atoms related to broadly defined electronegativity of the elements. In comparison to other methods of atomic partial charges approximation, results obtained within QTAiM are close to formal oxidation states for most of ionic systems [51].
A measure of averaged global charge transfer, which can be compared between different types of systems is an index defined as [52]:
c = 1 N N X Ω=1 q(Ω) OS(Ω) (2.55)
OS(Ω)is a formal oxidation state of atom Ω. The value of parameter c gives us the information about the rate of ionicity, where c ≈ 1 for ideal ionic material, and is equal to zero in homoatomic covalent systems by definition.
2.3.4. Chemical bond and its properties
In QTAiM characteristic of chemical bonds in a system is focused on the analysis of total electron density field at bond critical points. BCP is a point, where "space" of a chemical bond is represented the most accurately – it is placed on a bond path and at a zero flux surface which separates bonded atoms.
The value of ρ(r) at BCP can be used as a specific measure of a strength of a chemical bond. In the most basic definition the electron density is bound with a bond order BO in the following equation:
BO = eA(ρ(r)−B) (2.56)
The parameters A and B are fitted empirically to particular bond type (depend on types of bonded elements).
In the most simple interpretation, large values of ρ(r) represent strong bonds (e.g. covalent), small electron density values at BCP represent weak closed-shell bonds (e.g. van der Waals, hydrogen bonding).
Table 2.1: Values of ρ(r) and ∇2ρ(r)at BCPs as well as curvatures λ
iin various types of chemical
bonds (all parameters are given in atomic units (a.u.)). Compound ρ(r) ∇2ρ(r) λ 1, λ2 λ3 NaCl 0.011 0.052 -0.008 0.068 Cdc 0.239 -0.673 -0.436 0.198 Sidc 0.083 -0.128 -0.071 0.014 Nabcc 0.005 0.000 0.000 0.001 Arfcc 0.003 0.014 -0.002 0.019
Evaluation of the Laplacian of the electron density at BCP is very important as well. From the definition of Laplacian its value shows us if the electron density is accumulated at this point (i.e. negative curvatures of the density dominates in sum, so that ∇2ρ(r) < 0), or the BCP is just
more of the minimum (∇2ρ(r) > 0). The former case represents all bonds where electrons are
shared between bonded atoms (covalent bonds), the latter is observed in closed-shell interactions (for example ionic).
In some cases, evaluation of a bond character based on ρ(r) and ∇2ρ(r)is not enough to give
unambiguous conclusions and requires extended comparison or employing more sophisticated methods. One of the drawbacks of QTAiM theory which can be found in a literature is its low sensitivity on localisation/delocalisation of a bond, related to metallic character of this bond. ELF model can be employed in such case. Methods developed especially to analyse localisation of an electron in metallic and covalent-metallic materials are even more appropriate, e.g. electron localizability index (ELI) [53].
In the table 2.1 values of topological parameters at BCP for several different types of chemical bonds are presented. As can be seen, a strong covalent C–C bond is characterised by high value of the electron density as well as significant negative value of Laplacian. In the metallic system (Na), the electron density distribution is flat around BCP (values of all λ are close to zero). In ionic NaCl Laplacian of the electron density has clear positive value. In a typical ionic-covalent material, β-cristobalite, ρ(r) value is high, as well as ∇2ρ(r) 0.
Together with total electron density distribution, the densities of energy of particular type at bond critical point can give us important information about the stability of the bond. The virial theorem in its local formulation in a stationary state binds together the Laplacian of ρ(r) with potential energy density V(r) and gradient kinetic energy density G(r):
~2 4m
∇2ρ(r) = 2G(r) + V(r) (2.57)
Cremer and Kraka [54] showed benefits while evaluating the simple sum of both energy densities, which is the total electronic energy density He(r):
