Numerous systems and materials have been successfully investigated within the DFTB ap-proach [89, 170–176]. But not for every case the chosen input density is satisfactory. For systems with considerably charged atoms the error introduced by the assumed neutral den-sity becomes appreciable. Such problematic systems include heteronuclear molecules and polar semiconductors. A systematic extension of the method is accomplished by including up to the second term of the energy functional (3.5), that is, the second-order in the density difference
∆ρ. This leads to a self-consistent scheme with the iterative update of the electron density until convergence is reached. Such a extension is called self-consistent-charge DFTB (SCC-DFTB).
In the following we will consider a spin-unrestricted scheme [177]. This generalization implies to write the second-order contribution to the energy as
E(2) = 1 2
X
στ
Z Z
∆ρσ(r)fhxcστ[ρ0](r,r0) ∆ρτ(r0)drdr0, (3.18) where the sum runs over the spin variables σ =↑,↓, and ∆ρσ = ρσ −ρ0,σ. The input density ρ0 = [ρ0,↑, ρ0,↓] is chosen spin-unpolarized (ρ0,↑(r)−ρ0,↓(r) = 0, ∀r). The exchange-correlation component of fhxc now depends on the spin-dependent electron density,
fxcστ[ρ](r,r0) = δ2Exc[ρ]
δρσ(r)δρτ(r0). (3.19)
For the succeeding formulation, the Greek indices µ, ν, κ, λ will denote atomic orbitals (AO).
Let us also abbreviate a general two-point integral over a kernelg(r,r0) in the following form:
(f|g|h) = Z Z
f(r)g(r,r0)h(r0)drdr0. (3.20)
For atomic orbital productsf(r) =φµ(r)φν(r) andh(r0) = φκ(r0)φλ(r0) the shorthand (µν|g|κλ) will be also used. Expanding the KS orbitals into the AO set, ψsσ = P
µcσµsφµ, E(2) can be written as
E(2) = 1 2
X
στ
X
µνκλ
∆Pµνσ (µν|fhxcστ[ρ0]|κλ) ∆Pκλτ , (3.21) where ∆Pµνσ =Pµνσ −Pµν0,σ denotes the AO density matrix of the difference density ∆ρσ(r), with
Pµνσ =X
st
cσµscσνtPstσ, Pµν0,σ =X
st
cσµscσνtPst0,σ. (3.22) Here Pstσ = hΨ|aˆ†tσˆasσ|Ψi are the MO density matrix elements, Ψ being the ground state KS determinant. This leads to Pstσ = nsσδst for the orthogonal KS functions ψsσ. The term Pst0,σ designates the MO density matrix of the reference system.
By applying the variational principle to the modified energy functional, E =E[ρ0] +E(2), we obtain the secular equation (3.3) with a modified DFTB Hamiltonian,Hµν =Hµν0 +Hµν1 , where
Hµνσ1 =X
τ
X
κλ
(µν|fhxcστ[ρ0]|κλ) ∆Pκλτ . (3.23) Thus, second-order corrections due to charge fluctuations leads to a KS Hamiltonian which includes up to first-order corrections. The Hamiltonian is now self-consistent and so the secular equations are solved iteratively.
3.3.1 Mulliken and monopole approximations
So far the SCC formalism has been described exactly. Next, two important approximations will be applied for the evaluation of the multicenter integrals appearing in (3.21) and (3.23), which will simplify considerably the scheme. First, the Mulliken approximation is applied. This amounts to set φµφν ≈ 12Sµν(|φµ|2 +|φν|2), using the known AO overlap integrals Sµν. Thus, the general four-center integrals are written in terms of two-center integrals,
(µν|fhxcστ[ρ0]|κλ)≈ 1
4SµνSκλ X
α∈{µ,ν}
X
β∈{κ,λ}
(αα|fhxcστ[ρ0]|ββ). (3.24)
By substituting Eq. (3.24) into (3.21), the second order contribution to the energy reads, E(2) = 1
2 X
στ
X
µν
∆ ˜Pµµσ (µµ|fhxcστ [ρ0]|νν) ∆ ˜Pνντ , (3.25)
where ˜Pµνσ are the elements of the dual density matrix defined as P˜σ = 1
2(Pσ·S+S·Pσ). (3.26)
The concept of dual representation of the density matrix was introduced by Han and coworkers in their implementation of the LDA+U method [178]. The authors show that for a nonorthog-onal orbital basis this matrix satisfies exactly the sum rule, P
σTr( ˜Pσ) =Ne, where Ne is the
total number of electrons. They also pointed out, that the use of the dual density matrix is consistent with the Mulliken population analysis. Thus, for example, the trace of the atom-A block of ˜Pσ returns the Mulliken population for that atom:
qσA=X
µ∈A
P˜σ =X
s
nsσX
µ∈A
X
ν
cσµscσνsSµν. (3.27)
A further simplification to Eq. (3.25) is obtained by spherical averaging over AO products. Such monopole approximation ensures that the final total energy expression is invariant with respect to arbitrary rotations of the molecular frame. To this end, the functions FAl are introduced:
FAl(|r−RA|) = 1 2l+ 1
m=l
X
m=−l
|φµ(r−RA)|2, (3.28) where l and m denote the angular momentum and magnetic quantum number of AO φµ, centered on atom A. We then have for µ={Alm},ν ={Bl0m0}:
(µµ|fhxcστ[ρ0]|νν)≈(FAl|fhxcστ[ρ0]|FBl0) = ΓστAl,Bl0 (3.29) introducing the shorthand notation Γ. Thus, the second-order energy can be expressed as
E(2) = 1 2
X
στ
X
AB,ll0
∆qσAlΓστAl,Bl0∆qBlτ 0, (3.30) where qσAl is the Mulliken population for atomA, angular momentum l and spin σ.
In practice, rather than working with spin-dependent charge fluctuations, it is preferred to evaluate the net charge and spin populations during the self-consistent process. This is achieved after transformation from the set {ρ↑, ρ↓} to the total density ρ = ρ↑+ρ↓ and magnetization m = ρ↑−ρ↓. By this change of variables, the XC kernel can be split and one arrives at (see Appendix A)
ΓστAl,Bl0 =γAl,Bl0 +δσδτδABWAl,l0, (3.31) where δσ = 2δσ↑−1 and the parameters
γAl,Bl0 = (FAl|fhxc[ρ0]|FBl0) (3.32) WAl,l0 = FAl
δ2Exc[ρ, m]
δm(r)δm(r0) ρ0,0
FAl0
!
. (3.33)
The constants WAl,l0 depend only on the XC kernel, but not on the long-range Coulomb in-teraction. Moreover, as the reference density ρ0 is built from neutral spin-unpolarized atomic densities and the XC functional is an even functional inm, there are no integrals that involve mixed derivatives of the XC energy with respect to both the density and magnetization. The parameters γAl,Bl0 and WAl,l0 are known in DFTB as the γ-functional and spin coupling con-stants, respectively [179]. The latter are treated as strictly on-site parameters, whereas γAl,Bl0 is calculated for every atom pair using an interpolation formula, that depends on the distance RAB between the atomsA and B and the atomic Hubbard-like parametersγAl,Al and γBl0,Bl0:
γAl,Bl0 = 1
RAB −S(RAB, γAl,Al, γBl0,Bl0). (3.34)
S is a short-range function ensuring the correct convergence at RAB = 0 and the correct R−1 behavior at large interatomic distances.
Traditionally, the Hubbard parameters are not computed directly from the integrals (3.32), but from total energy derivatives according to γAl,Al = δ2E/δn2, where E denotes the DFT total energy of atom A and n refers to the occupation of the shell with angular momentum l. The derivative is evaluated numerically by full self-consistent field calculations at perturbed occupations. Due to orbital relaxation, Hubbard parameters obtained in this way are roughly 10-20 % smaller than the ones from a direct integral evaluation. It has been pointed out that use of a low-quality basis set leads to overestimated Hubbard parameters as the missing basis functions would otherwise screen its value [180]. Thus, the approximation used for determining γAl,Al in DFTB partly compensates for the error introduced by the employed minimal basis.
In part because of this error compensation, DFTB often returns better results than DFT itself when employing a minimal basis set.
By substituting Eq. (3.31) into (3.30), the second-order energy term finally reads E(2) = 1
2 X
AB,ll0
∆qAlγAl,Bl0∆qBl0 +1 2
X
A,ll0
∆mAlWAl,l0∆mBl0, (3.35) whereqAl =qAl↑ +qAl↓ andmAl =qAl↑ −qAl↓ are the Mulliken net charge and spin population of shell l on atom A. Likewise, by applying the Mulliken (3.24) and monopole (3.29) approximations to Eq. (3.23), one arrives at
Hµνσ1 = 1
2SµνX
τ
X
γ
[(µµ|fhxcστ [ρ0]|γγ) + (νν|fhxcστ[ρ0]|γγ)] ∆ ˜Pγγτ
= 1
2SµνX
τ
X
Cl00
ΓστAl,Cl00+ ΓστBl0,Cl00
∆qτCl00
= 1
2SµνX
Cl00
(γAl,Cl00+γBl0,Cl00) ∆qCl00 +1
2δσSµν
X
l00
(WAl,l00∆mAl00+WBl0,l00∆mBl00). (3.36) It should be noted that atomic blocks of the Hamiltonian matrix are diagonal due to the orthogonality of the AO centered at the same atom. In the SCC scheme, the repulsive potentials are determined as a function of the distance by taking the difference between DFT cohesive energies and the corresponding SCC-DFTB electronic energies, P
iniεi +E(2), for suitable reference systems. SCC-DFTB has been applied to large biological systems, atomic clusters and solids with great success [181–184]. Nowadays, the use of non-self-consistent DFTB has been diminished in favor of its SCC variant, for which the term DFTB will stand in the following.