## Abstract

This paper reviews the basic principles of the density-functional tight-binding (DFTB) method, which is based on density-functional theory as formulated by Hohenberg, Kohn and Sham (KS-DFT). DFTB consists of a series of models that are derived from a Taylor series expansion of the KS-DFT total energy. In the lowest order (DFTB1), densities and potentials are written as superpositions of atomic densities and potentials. The Kohn–Sham orbitals are then expanded to a set of localized atom-centred functions, which are obtained for spherical symmetric spin-unpolarized neutral atoms self-consistently. The whole Hamilton and overlap matrices contain one- and two-centre contributions only. Therefore, they can be calculated and tabulated in advance as functions of the distance between atomic pairs. The second contributions to DFTB1, the DFT double counting terms, are summarized together with nuclear repulsion energy terms and can be rewritten as the sum of pairwise repulsive terms. The second-order (DFTB2) and third-order (DFTB3) terms in the energy expansion correspond to a self-consistent representation, where the deviation of the ground-state density from the reference density is represented by charge monopoles only. This leads to a computationally efficient representation in terms of atomic charges (Mulliken), chemical hardness (Hubbard) parameters and scaled Coulomb laws. Therefore, no additional adjustable parameters enter the DFTB2 and DFTB3 formalism. The handling of parameters, the efficiency, the performance and extensions of DFTB are briefly discussed.

## 1. Introduction

Density-functional theory (DFT) within the orbital formulation by Kohn & Sham (KS) [1] is a very powerful method for the calculation of physical and chemical properties of molecules, clusters and condensed matter systems. It is a very successful complement to the wave function-based methods, which build upon Hartree–Fock (HF) theory to include electron correlation.

Although efficient DFT-KS computer codes and the access to increasingly powerful computers allow one to study systems with remarkable size and complexity, there is demand for even faster but approximate methods to reach much larger system sizes, extended time scales in molecular dynamics simulations and to treat a multitude of structures of one molecular system. Classical force field methods (molecular mechanics, MM), on the other hand, are several orders of magnitude faster than DFT(HF), allowing one to treat hundreds of thousands to millions of atoms and to follow their dynamics beyond the nanosecond time scale. Such classical force fields contain a large number of empirical parameters and can be very accurate for the system they have been parametrized for, but may suffer from a limited transferability, from the unavailability of accurate force field parameters for a number of elements, and the missing quantum effects in such types of simulations.

Semi-empirical (SE) quantum chemical methods are located in the middle between DFT(HF) and MM methods, being roughly three orders of magnitude faster than DFT(HF) and slower than MM, containing a number of element-specific parameters, more than DFT methods, which are often parametrized nowadays as well, but significantly less than that in MM. These parameters require a reasonable effort in their determination and are sufficiently transferable, but most importantly, the SE methods can describe nearly all quantum effects accessible by DFT or HF and post-HF methods with few exceptions. One developmental line started from HF by neglecting and parametrizing the interaction integrals and resulted in the very successful neglect of diatomic differential overlap (NDDO) type of approaches, the most well known called MNDO, AM1 and PM3. Research along this line has been very active over the years, and several improved methods have been proposed, such as PM6 [2], PDDG-PM3 [3] or OMx [4], to name a few.

The empirical tight-binding (TB) [5] methods represent the counterpart to the SE quantum-chemical methods in the condensed phase with similar limitations. It was shown in the 1980s that a ‘TB’ scheme, i.e. the application of limited atomic orbital (AO) basis-set representation, can be used as the basis for a nearly parameter-free approximate treatment within the KS-DFT [6]. These early ideas have been further developed during the past 20 years, resulting in a density-functional theory-based tight-binding (DFTB) method. The DFTB method combines the computational efficiency of the SE ‘TB’ methods; however its parametrization is less empirical and more straightforward to conduct because it is basically related to a number of DFT calculations. In contrast to the NDDO-type methods, which start from HF and introduce correlation effects solely through parametrization, electron correlation in DFTB is included from the beginning owing to its descent from DFT. On the other hand, this strategy imports also all the DFT failures resulting from the use of approximate exchange correlation functionals.

The DFTB method has been applied to study large molecules (e.g. biomolecules), clusters, nanostructures and condensed-matter systems with a wide range of elements. The purpose of this paper is to describe the main ideas of DFTB and its potential, illustrated by a few selected applications. Finally, its conceptual limitations as well as extensions are discussed.

## 2. Theoretical approach

DFTB can be derived from a Taylor series expansion of the KS density functional total energy [1] around a properly chosen reference density *ρ*(** r**). Instead of finding the electron density

*ρ*(

**) that minimizes the energy a reference density**

*r**ρ*

^{0}is assumed which is perturbed by some density fluctuation,

*ρ*(

**)=**

*r**ρ*

^{0}(

**)+**

*r**δρ*(

**). The exchange–correlation energy functional is then expanded in a Taylor series up to the third order and the total energy can be written as 2.1Depending on the inclusion of terms from this expansion, different models appear, which are built successively on top of each other. Over the years, several DFTB models have been implemented, starting from the first-order non-self-consistent DFTB1 [7,8] (originally called DFTB), the second-order DFTB2 (originally called SCC-DFTB) [9] and the recent extension to the third-order DFTB3 [10–13].**

*r*### (a) DFTB1

The DFTB1 method takes the first two contributions of equation (2.1), *E*^{0}[*ρ*_{0}] and *E*^{1}[*ρ*_{0},*δρ*], into account.

DFTB is based on a linear combination of atomic orbital (LCAO) ansatz of the KS orbitals
2.2The AOs are obtained from DFT calculations of the corresponding atoms. The basis is restricted to the valence shell of the atoms. Such restriction requires an orthogonalization to the core orbitals of the atoms, which can be achieved by a Schmidt orthogonalization
2.3where |*ϕ*_{μ}) is the valence AO *μ* at atom *a* and is a core orbital at atom *b*, as obtained from the corresponding atomic calculations. This orthogonalization is only necessary for the case *a*≠*b*, because the AOs of the atomic DFT calculation represent an orthonormal set of functions. The LCAO ansatz leads to a general eigenvalue problem of the form
2.4with the Hamiltonian matrix elements and the overlap matrix elements *S*_{μν}. Considering the orthogonality construction (2.3) and writing the effective KS potential *V* _{eff} as a superposition of atomic-like potentials (), one can write
2.5where *ϵ*_{κc} is the energy of a core state *κ* at centre *c*. Together with its core correction term the potential can be interpreted as a pseudo-potential. The pseudo-potential appears only in the three-centre terms of the Hamiltonian matrix elements (*a*≠*b*≠*c*) and in the ‘diagonal terms’ *ν*,*μ*∈*a* with *c*≠*a*, whereas the ‘full’ potential appears in all the other terms
2.6The neglect of the three-centre terms and pseudo-potential contributions leads to the two-centre approximation for the Hamilton matrix elements
2.7The approximations lead to the same structure of the secular equations as in (non-orthogonal) TB schemes or the extended Hückel method; however all matrix elements are calculated within DFT, because the AO basis set *ϕ*_{μ} is computed by solving the DFT-KS equations for atoms. Pure AOs would be too diffuse for a minimal AO basis set; therefore, the atomic KS equations are usually solved applying an additional (harmonic) potential to the atomic KS equations, which leads to ‘compressed’ AOs or optimized AOs, as introduced by Eschrig [14]
2.8

The confinement radius is taken to be roughly two times the covalent radius of the atom [8]. The potentials *V*^{a}_{eff}=*V* _{eff}(*ρ*_{a}) contain the electron density of the neutral atoms *a*. This density is usually also determined from atomic KS equations, in the same way as the AO basis set. However, a slightly different confinement radius can be chosen for these initial densities *ρ*_{a}. Therefore, two parameters per element have to be determined at this stage, the confinement radius for the AO basis, *r*_{0}, and the confinement radius for the initial density, .

Now, with AO basis and initial density determined, the KS equations can be solved (equation (2.4)) leading to the energy (*n*_{i}: occupation number of KS orbital *i*)
2.9This is the electronic energy of the DFTB method, i.e. the sum of the occupied KS energies. To get an expression for the total energy, the *E*^{0} term has to be approximated. *E*^{0}[*ρ*_{0}] consists of the DFT ‘double counting’ contributions in the first line of equation (2.1) and depends only on the reference density *ρ*_{0}, i.e. the superposition of the neutral atomic densities. This means that this term is not dependent on the specific chemical environment, i.e. it can be determined for a ‘reference system’ and then be applied to other molecules in different chemical environments. This is the key to transferability of the parameters. In DFTB, this term is approximated by a sum of pair potentials called repulsive energy term [15]
2.10which is either determined by comparison with DFT calculations [8] or fitted to empirical data [16]. The total energy then reads
2.11Forces are readily computed as derivatives of the total energy, which allows very efficient geometry optimization or molecular dynamics simulations for molecular systems. This results, on the one hand, from the minimal basis representation and on the other hand from the efficient handling of the matrix elements. The *H*_{μν} and *S*_{μν} are computed once and stored, i.e. they do not have to be computed during the runtime of the DFTB program and they are read in from tables. Further, DFTB1 is a non-self-consistent TB method, i.e. the KS equations are solved only once. This makes DFTB1 about 5–10 times faster than DFTB2 and DFTB3, where a self-consistent solution of the eigenvalue problem has to be computed. DFTB1 is suitable for systems where the charge transfer between the atoms is small, typically homonuclear systems or systems with atoms of similar electronegatitvity. Therefore, DFTB1 is very well suited, for example, for the description of hydrocarbons, because the higher order terms are small here. On the other hand, DFTB1 can also treat systems where a complete charge transfer between the atoms occurs, as e.g. in NaCl [10]. For systems exhibiting a delicate charge balance, higher order terms have to be considered.

### (b) DFTB2 and DFTB3

DFTB2 and DFTB3 approximate the *E*^{2} and *E*^{3} terms in equation (2.1) further. First, the density fluctuations are written as a superposition of atomic contributions
2.12and second, the atomic-like density fluctuations are expanded in a multipole expansion, however, only keeping the monopole term
2.13Practically, the second-order integral (third line in equation (2.1)) is evaluated assuming an exponentially decaying charge density
2.14and neglecting the term ∂^{2}*E*_{xc}/∂*ρ*(*r*)∂*ρ*(*r*′) for the moment, i.e. only the Hartree integral is analytically evaluated. This leads to an analytical function *γ*_{ab}, and the second-order term finally reads (for *a*≠*b*)
2.15The Hartree term therefore describes the interaction of the charge density fluctuations *δρ*_{a} and *δρ*_{b}, which reduces to a Coulomb interaction of partial charges Δ*q*_{a} and Δ*q*_{b} for large distances, i.e. *γ*_{ab} approaches 1/*R*_{ab} for large distances.

For short distances, one has to deal with the electron–electron interaction within one atom, and here the exchange–correlation contributions become important. In this case, *E*^{2} describes the electron–electron interaction within one atom, and *E*^{2} becomes proportional to and *γ*_{ab} describes the effective on-site electron–electron interaction. The exponent *τ*_{a} is chosen such that the on-site value of the *γ*-function properly describes the atomic chemical hardness (or alternatively the Hubbard parameter as calculated from DFT) and thus implicitly takes into account the exchange–correlation contribution to the second-order term. The chemical hardness of an atom, or the Hubbard parameter, is the second derivative of the energy with respect to the charge density, i.e. exactly the third term in equation (2.1). Then, the Hartree energy can be evaluated analytically, and after introducing and approximating for the *E*_{xc} contribution this term boils down to an effective interaction between atomic charges, computed using the Mulliken scheme.

For *R*_{ab}=0, we find and the function *γ*_{aa} is simply the second derivative of the total energy with respect to the charge of an atom, i.e. it equals the so-called Hubbard parameter *U*_{a}, which is twice the chemical hardness of an atom [9]. From the functional form of *γ*_{ab}, one finds that the Slater exponent of an atom *a* is related to the Hubbard parameter as
2.16i.e. the width of the atomic charge density is inversely proportional to its chemical hardness. This relation is intuitive in that more diffuse atoms (or anions) have a smaller chemical hardness. A problem, however, arises because the function *γ*_{ab} encodes a constant proportionality, which turns out not to be valid across the periodic table. In fact, at least every row of the periodic table should be treated slightly differently. The deviation is most pronounced between hydrogen and the first-row elements; therefore a modified has been proposed [17,10]. This modified function can be applied within DFTB2, but it is not the ‘default’ option. Within DFTB3, it has become default and therefore is a key ingredient of the DFTB3 methodology [11–13].

For *E*^{3}, the same approximations are introduced as for *E*^{2}. Originally, only the diagonal terms were included [17,10,11]. The third-order terms describe the change of the chemical hardness of an atom with its charge state [10], i.e. in the third order, a new parameter is introduced, the chemical hardness derivative . This parameter can be computed from DFT or fitted in order to improve the performance of the model. A function *Γ*_{ab} results as the derivative of the γ-function with respect to charge by introducing the Hubbard derivative parameter. In that way, the third-order terms within SCC-DFTB can be rather seen as a robust way to introduce the charge dependence capturing some deficiencies of problematic approximations within the second-order formalism, namely, the small size of the pseudo-AO basis as well as the very simplified density fluctuation scheme.

With all these approximations, the SCC-DFTB total energy in the third order is given by
2.17The derivative of this expression with respect to the molecular orbital coefficients leads to the corresponding KS equations
2.18and
2.19where *S*_{μν} is the overlap matrix. The Hamilton matrix elements depend on the Mulliken charges *q*_{a} (Δ*q*_{a}=*q*_{a}−*Z*_{a}) which in turn depend on the molecular orbital coefficients *c*_{μi}, and thus these equations have to be solved self-consistently.

### (c) Parameters and efficiency

To parametrize the various DFTB models, one first has to determine four parameters per atom type. These are the confinement radii for the AOs *ϕ*_{μ}, which is called *r*_{0}, the confinement radii for the atomic densities , which is called , the atomic Hubbard parameters *U*_{a} and its derivative , which both can be computed from DFT. The repulsive potentials are two-body contributions, therefore they are much more involved and some effort has to be invested in the parametrization process. Typically, on determining for every pair of atoms the total energy *E*^{tot} versus bond distance *R* for specific bonds using a high level DFT or *ab initio* method determines (*E*^{1}+*E*^{2}+*E*^{3}) also versus bond distance, and then computes *E*_{rep} as
2.20The determination of *r*_{0} and is an empirical procedure and can be quite involved [16,13], while *U*_{a} and can be easily computed in principle. For the modified function *γ*^{h} one additional parameter appears, which is fitted to reproduce the water dimer binding energy [13].

The substantial advantage of using DFTB2 is its time/performance efficiency. Table 1 presents benchmark calculations for the CPU-time of a single-point energy calculation.^{1}

DFTB2 is at least 250 times faster than RI-PBE and more than 1000 times faster than B3LYP. This acceleration is mainly owing to two factors: first the usage of a minimal basis set within SCC-DFTB and second the tabulation and neglect of integrals. For the water cluster (H_{2}O)_{48}, for example, *N*=288 basis function is needed for a minimal basis set and *N*=864 basis function for the 6-31G(d) basis set. The time-limiting step for obtaining the energy with all methods discussed here is a matrix-diagonalization, which scales with *N*^{3}. Thus, an acceleration just from using the minimal basis of a factor of 27 is achieved. The remaining factor is owing to the tabulation and neglect of integrals; in this example this factor is roughly 10 and 40 for comparison with RI-PBE and B3LYP, respectively.

### (d) Extensions

Within DFT, the correct description of the weak van der Waals interactions, the London dispersion interactions, are not included in the DFTB approach as formulated above. However, the van der Waals interactions may be included in the DFTB approach *a posteriori*, which means that the van der Waals contributions are calculated afterwards and the energy *E*_{disp} is added to the total energy. The first implementation [20,21] showed a very promising performance, stimulating the combination with empirical dispersion also on the full DFT level [22–25]. Recently, Grimme has parametrized his D3 [25] correction for DFTB3, leading to an excellent performance of DFTB3-D3 for a large set of hydrogen- and van der Waals-bonded molecules.

As is well known, the popular DFT functionals fall short in a proper handling of strong electron correlation effects in molecules or solids. Practical ways to cure the problems related with electron localization are the LDA+U approach and the application of self-interaction correction schemes. Both have been implemented into the DFTB methodology [26]. The proper description of the correct balance between charge delocalization and polarization, for example in charge transfer complexes, is also a challenge for DFT. Rapacioli *et al.* [27] adapted recently a configuration interaction method, based on constrained DFT calculations, into the DFTB approach which allows one to investigate charge resonances in molecular complexes and to describe the proper dissociation behaviour.

Electronic excitations can be described using time-dependent density-functional response theory. The corresponding approach within the DFTB scheme has been put forward by Niehaus *et al.* [28]. The time-dependent density-functional approach was also adapted for the description of non-adiabatic processes [29] and later on combined with Tully's stochastic surface hopping algorithm [30].

For heavy atoms, a scalar relativistic treatment is used [31] and the spin–orbit coupling can also be considered [32].

In combination with the non-equilibrium Green's function technique, the DFTB approach was also extended for the calculation of electron transport properties on the molecular scale [33].

Although the DFTB approach is about two to three orders of magnitude faster than the DFT one, for large systems (more than about 1000 atoms) the scaling behaviour in the diagonalization methods for the solution of the general eigenvalue problems becomes notable. Therefore, the so-called order-N scaling schemes [34,35] or the Car–Parrinello treatment [36] have been implemented within DFTB. For calculations of even larger systems, a combination (QM/MM) of the DFTB method, as a quantum mechanical (QM) method, with empirical force fields (MM) has also been realized [37–39]. This scheme has further been extended to include also a continuum electrostatics environment in the DFTB/MM-GSBP scheme [40].

## 3. Performance and tests

DFTB is a set of parametrized computational models which require a careful parametrization before application. This has a number of implications. First, DFTB models are based on DFT; therefore these models also inherit most of the well-known DFT shortcomings. As the Hamilton matrix is computed using DFT-GGA (with PBE functional), the problem of self-interaction is also found in DFTB. This results in a number of problems, showing up in a problematic description of the KS energy spectrum, the overpolarizability problem, etc. The problematic description of van der Waals complexes was recognized early and a correction has been proposed (see above), which has later on also been adopted for full self-consistent DFT calculations [22,23]. Second, on top of DFT, the above-described approximations are introduced. This results in a limited flexibility of the computational scheme, leading to a restricted transferability of the parameters. It has been shown that for certain properties, the accuracy of DFT can be achieved. This holds, for example, for geometries and molecular geometries, or for geometries and vibrational frequencies. However, no parametrization could be derived which approaches DFT accuracy for both energies and vibrational frequencies simultaneously [16,13]. Therefore, different parametrizations have been proposed.

DFTB requires two sorts of parameters. First, there are the electronic parameters concerning the wave function and density confinement, which enter the *E*^{1} term, and the Hubbard parameters and its derivatives, which are usually computed from DFT and are therefore easy to determine. The confinement radii are a little more involved, because they can be chosen to optimize molecular properties (e.g. [16,13]). However, these parameters allow only a fine-tuning, i.e. they cannot change the performance of the models drastically. Second, there are the repulsive potential parameters constituting *E*^{0}. These are fitted to theoretical or experimental data and this fitting process is crucial for DFTB accuracy. However, only bonds are affected (no bond angles, dihedrals, etc.), i.e. the bond energies, bond distances and stretch vibrational frequencies can be modified.

DFTB has not been parametrized for the whole periodic table (for available parameters, see http://www.dftb.org). Basically, there are parameters available for organic systems and for inorganic systems, which are not all compatible to date. For materials science applications, there exist parameter sets also for a number of element combinations (parameter set matsci—for references, see http://www.dftb.org), which have been applied to several inorganic systems, based on the DFTB1 approximation (e.g. [41,42] and references therein). DFTB2 has been particulariy parametrized for O, N, C, H [9], S [43] and Mg [44], a parameter set called mio. These ONCH parameters have been tested quite intensely with respect to their performance for atomization and reaction energies, molecular geometries, vibrational frequencies, ionization potentials and electron affinities, dipole moments and for non-covalent interactions [45–47]. DFTB2 performs very well for reaction energies [45], while for heats of formation other SE methods, in particular the newly developed OMx and PDDG-PM3 methods, turned out to be slightly superior.

Concerning geometries, the DFTB2 performance is close to that of full DFT-GGA, while for vibrational geoemetries larger deviations with respect to DFT have been reported. A special parametrization for vibrational frequencies turned out to be very successful leading also to a very accurate description for these type of properties [48]. A deeper analysis showed that there is indeed an optimization conflict: it is not possible within DFTB2 (and also DFTB3) to achieve simultaneously an accuracy comparable to that of DFT. Therefore, special parametrizations for vibrational frequencies have been proposed for DFTB2 [16] and DFTB3 [13].

Ionization potentials and electron affinities for small molecules are not very accurate using small basis set methods. In contrast to the NDO-type SE methods, these properties do not enter the DFTB parametrization process. Therefore, theses values may be less accurate and a careful testing is required before the application to a new class of molecules. Dipole moments in DFTB are computed using a Mulliken population analysis, leading to an underestimation of these properties. Using parametrized charges similar to those supplied by the CM3 scheme, a huge improvement can be found [49,50]. However, the derivatives of charges and polarizabilities with respect to nuclear coordinates, as needed in order to compute IR and Raman intensities, cannot be improved by such simple extension as CM3 charges or the chemical potential equilibration extension [50,51].

The main shortcomings of DFTB2 are its overbinding, its performance for hydrogen-bonded complexes and charged species, which show up, for example, in the proton affinity values. Here, the extension to the third order and the modified function has led to a major improvement. In the new DFTB3 parametrization [13], the overbinding has been removed leading to an excellent performance for reaction energies and a clear improvement for non-bonded interactions. Note that DFTB3/3OB is even better than DFT-GGA when only a small, double zeta-type, basis set is applied [13]. In particular hydrogen bonding energies, proton affinities and proton transfer barriers are very well described.

## 4. Conclusion

The DFTB models allow one now to treat a broad range of systems (molecules, clusters and solids) and their computational efficiency enables one to consider systems with up to several thousand atoms. Applications have been described including elements all over the periodic table; however a systematic parametrization across the periodic table is still missing.

With the realization of the third-order expansion in the total energy (DFTB3) an excellent performance with respect to accuracy and computational efficiency for organic and biomolecules could be reached being about two to three orders faster than DFT-GGA methods with small double zeta-type basis sets, depending on the implementation and system studied. In many cases, DFTB reaches the DFT accuracy, however not for all properties within one parametrization. The computational efficiency comes at a certain price, which has been described in this review. The derivation from DFT-GGA makes DFTB inherit the well-known DFT problems, which can be remedied only partially. In addition to this, the use of a minimal basis set, the limited flexibility owing to the monopole approximation at DFTB2/DFTB3 level and the neglect of three-centre integrals mark approximations which lead to a deterioration of the accuracy when compared with full DFT calculations. Unfortunately, there is really no free lunch.

## Footnotes

One contribution of 8 to a Theme Issue ‘Density functional theory across chemistry, physics and biology’.

↵1 All calculations were carried out on a single processor of a standard desktop PC. For SCC-DFTB, the DFTB+ code (http://www.dftb.org) was used. The DFT values were obtained using the TURBOMOLE program package [18]. For the PBE calculations, the resolution of the identity (RI) integral evaluation has been used [19]. As basis set for the DFT methods, we chose 6-31G(d) which is a rather small basis set for practical use.

- © 2014 The Author(s) Published by the Royal Society. All rights reserved.