research area has been dispersion corrections to DFT.(65, 66) While no single method
“just works” in any chemical environment, some variant of DFT has been successfully applied in nearly all areas of chemistry and chemical physics.(40) For a more in depth introduction, see e.g. the second half of the book by Koch and Holthausen.(57)
3.6 Basis sets
Wavefunction based methods, such as the HF theory, do need to express this many-body wavefunction in some basis. In Kohn-Sham DFT, the electron density is generated from orbital densities. Instead of simply optimizing a free form function on a spatial grid, i.e. using the position basis, the orbitals are built from linear combinations of specially adapted functions.
ϕi(r) =X
µ
Cµiχµi(r) (3.38)
These basis functions χ may be derived from atomic orbitals, obtained e.g. by solving the Schr¨odinger equation for the hydrogen atom. Hence, a so called Slater Type Orbital (STO) has the form
χ(r, θ, φ) =Arn−1e−ζrYlm(θ, φ) (3.39) The radial part of the solutions for the hydrogen atom contains an exponentially decaying part as well as a polynomial.
3.6.1 Gaussian type orbitals
The exponential part is inconvenient from a numerical point of view, as integrals e.g of the type in eq. 3.13 are not easily evaluated. Therefore, typically Gaussian type functions (GTFs), also abbreviated as GTOs (Gaussian type orbitals) are chosen for this purpose, as integrals over products of differently centered Gaussians may be easily evaluated using the Gaussian product rule. A single Gaussian does not accurately reproduce the sharp cusp at the center of an exponential, neither does it have the correct asymptotic decay. However, using linear combinations of Gaussians one may approximate atomic orbitals. Linear combinations of Gaussian functions with fixed
coefficients may be used as one atomic orbital type basis function. This is referred to as a contracted GTO (CGTO). In this case
χ(r, θ, φ) =Ylm(θ, φ)rlX
i
αiAi(l, α)e−ζir2 (3.40) Where the coefficients ci and the exponentsζi predetermined in order to optimally fit atomic orbitals. Instead of spherical coordinates, the GTOs may also be expressed in cartesian coordinates. Atomic orbital type basis functions, such as GTOs have localized centers in real space.
3.6.2 Plane waves
Another basis especially suited for periodic systems are plane waves. These are by nature periodic and localized in reciprocal space.
χG(r) =AeiG·r (3.41)
The normalization factor A is given by √1
Ω with Ω as the volume of the unit cell.
In a three dimensional lattice, G contains the reciprocal lattice vectors. The peri-odicity of the plane waves is especially advantageous in the domain of solid state physics, where crystal periodicity enforces periodic wavefunctions. According to the Bloch theorem(67), the wavefunctions of a particle in a periodic potential will have the form of
ψnk(r) =e[ik·r]unk(r) (3.42) Where u has the same periodicity as the potential, e.g. the crystal unit cell. The different wavevectors k are quantum numbers with discrete energy levels k. They correspond to a wavefunction phase factor with a different periodicity than the unit cell itself, e.g. several multiples of a unit cell. Finding allkin the first Brillouin zone of the crystal will give the band structure. In the case of a metallic system these conduction bands have nonzero occupation, hence an electron can move freely over several unit cells. When treating insulators, as is the case in this thesis, it is sufficient to consider the electronic 0 K ground state atk= 0. This is called the Γ-point approximation. In the following sections, thekdependency will be omitted for the sake of simplicity.
3.6 Basis sets
The plane waves (PW) form a complete basis - the momentum basis for a particular periodic system. As they are not localized in real but in reciprocal space, they avoid some problems of atom centered basis sets (see below). The Kohn-Sham orbitals are expressed in PW as follows
ϕi(r) = 1
√Ω X
G
Ci,Ge[iG·r] (3.43)
As in practice the basis set expansion is cut off at some point, the basis is incomplete which leads to certain systematic errors. When using GTOs, errors arise from the position dependence of the basis functions, leading to effects dependent on the relative positions of nuclei to each other. One example for this is the basis set superposition error (BSSE). Another are the Pulay forces. The use of plane waves avoids these problems, as the basis is dependent on the unit cell G and not the atomic positions.
However, the energy may instead unphysically depend on the position of a molecule in the periodic box because of the finite grid resolution. When using PW, the basis set cutoff is usually given in terms of the kinetic energy. This kinetic energy corresponds to the energy of the free particle whose wavefunction isχG, e.g. 280 Ry. As the reciprocal unit cellGis three dimensional, this energy cutoffEcis expressed in terms of a sphere
1
2|G|2 ≤Ec (3.44)
Another natural advantage of plane waves is the ease of computing the Hartree potential (i.e. mean field classical electrostatic potential defined in eq. 3.13). In real space, the Hartree potentialνH is written as
νH(r) =
Z ρ(r0)
|r−r0|dr0 (3.45)
In absence of nuclear charge density, in reciprocal space the equation becomes νH(G) = 4πρ(G)
G2 (3.46)
and accordingly the classical electrostatic energy is
J[ρ] = 2πΩ X
G<GC
ρ∗(G)ρ(G)
G2 (3.47)
Where the Fourier component of the potential for G = 0 is not included. It cor-responds to an average, which will result only in a shift. In the case of a charged periodic cell, the potential diverges. In such a situation, a neutralizing charge has to be introduced, e.g. included in the G = 0 background. The external potential due to the nuclei can also be computed in reciprocal space in a similar manner, by adding the nuclear charge density to the density. In practice short-range interactions are often treated differently, in a local representation. The density ρ(r) can be accumulated in real space by summing the plane waves on a real space grid.
ρ(r) = 2 Ω
n/2
X
i
X
G
Ci,G∗ Ci,GeiGr (3.48) In order to obtain ρ(G) from ρ(r), a Fast Fourier Transform (FFT) may be used, yieldingνH(G) from the real space density with a scaling ofO(NlogN). Plane waves are by nature orthogonal to each other and their kinetic energy expectation value is trivial to compute. Another advantage is the very good representation of diffuse electron densities. However, the electron density near the atomic centers varies rapidly, so that a very high number of plane waves would be necessary to represent it. This problem is especially important for elements with high electronegativity and even worse for the KS-orbitals of core electrons. This is why pseudopotentials are used in most practical plane wave based computations - they are explained in the next section.