Currently, DFT can treat up to 100 atoms in routine applications, sometimes even more, and has been successfully applied to molecular dynamics simulations up to several picoseconds. In order to extend the applicability of the simulations to even more complex system such as large biochromophores, solid state physics, or molecules surrounded by solvate molecules, approximations to DFT such as the density functional tight binding (DFTB) method have been developed.[56–59] DFTB has been shown to provide a quite accurate description of ground state properties such as molecular geometries, vibrational frequencies, and reaction energies comparable in accuracy to full DFT.[58,61] It has been successfully applied to a wide range of problems in the fields of biomolecules, surfaces, interfaces, and point and extended defects in solid-state systems.[59,97,98]
The DFTB method is derived from DFT by choosing a reference densityρ0as a super-position of neutral atomic densitiesρ0 =PAρ0A and by expanding the DFT exchange-correlation energy up to second order:[58]
E = Here,φi are the KS orbitals,ni are their occupation numbers, ˆH0is the KS-Hamiltonian evaluated at the reference density, EXC and VXC are the exchange-correlation energy and potential, andEii are the core-core repulsions. Eq. 1.9 serves as a starting point for further approximations leading to the tight binding version of DFT.
The second term on the right hand side of Eq. 1.9, which describes the energy contri-bution due to the density fluctuation, can be decomposed in atom centered monopole
contributions:
where the charge fluctuations ∆qA on atomA are estimated from the Mulliken charge analysis. γAB is the so-calledγ-functional defined as:
γAB =Z Z 1 Here, FA denotes the normalized spherical density distribution located on A, which means that the angular deformation of the charge density in second order is neglected.
In the short-range limit|r−r0| →0,E2nd describes the electron-electron interaction on atomA. In this case,γABcan be approximated as a Hubbard-type interaction depending only on the Hubbard parameter (also known as chemical hardness)UAleading to:
E2nd≈ 1
2UA∆q2A. (1.12)
The first term on the right hand side of Eq. 1.9 involves the summation over Kohn-Sham orbitals, which are expanded as linear combination of a minimal basis set in DFTB.
Employing the basis functions bν, the first order Hamiltonian terms in Eq. 1.9 can be expressed as an eigenvalue equation: with the KS molecular orbital coefficients ciµ and the overlap matrix elements Sµν = hbµ|bνi. The Hamiltonian matrix elements have the form:
Hµν0 =DbµHˆ0bνE. (1.14) The diagonal elements thus correspond to atomic KS eigenvalues and the non-diagonal Hamiltonian matrix elements are calculated in a two-center approximation:
Hµν0 =Dbµ
Tˆ+Vef fhρ0A+ρ0Bibν
E, (1.15)
where bµ and bν are centered on atomsA and B, respectively, andVef f is the effective KS potential.
The four terms in the second line of Eq. 1.9 depend only on the neutral atomic densities
and inter-atomic distances and are therefore collected in a repulsive potential Erep:
whereErep can be approximated as a sum of short-range two body potentialsUAB: Erep = 1
2 X
AB
UAB(RAB). (1.17)
In practice, UAB is usually fitted from the difference of the total DFT energy and the electronic part of the DFTB energy with respect to the bond length RAB of an atom pair for an adequate set of reference systems. The determination of the pair potentials exhibits an effort ofN2 for N sorts of atoms.
Employing the above described definitions and approximations, the total DFTB energy can be summarized as: Applying the variational principle leads to the KS eigenvalue problem:
X
ν
(Hµν−iSµν)ciν = 0 (1.19) with the KS orbital energiesi, and the Hamilton matrix elements Hµν given as:
Hµν =Hµν0 +1 2Sµν
X
C
(γAC+γBC) ∆qC. (1.20) Notice that the Hµν are calculated only once for all possible combinations of elements with a DFT functional for a dense grid of two atomic distances and are tabulated after that. The Sµν have to be calculated from the atomic orbitals, since the basis set is not orthogonal, but this is also only performed once. Subsequently, a tabulated form is used.
Since no integral evaluation is necessary during the calculations, the remaining compu-tational effort for the determination of the total DFTB energy is the iterative solution of the eigenvalue problem in Eq. 1.20. Thus, the computational costs are dramatically reduced compared to DFT. Therefore, DFTB allows the calculation of systems up to several hundreds of atoms.
1.3.1 Limitations of DFTB
The molecular geometries optimized by DFTB are comparable to those obtained from DFT, while the vibrational properties are not sufficiently accurate, especially if DFTB is
not used in the self-consistent charge (SCC) version. The inaccuracy for the vibrational frequencies has been improved by applying a special parametrization of the repulsive potential which employs experimental data.[99] A further improvement of energies and frequencies seems to be possible by optimizing the strategies for the parametrization of Erep, e.g. by a more extensive use of experimental data or by employing genetic algo-rithms for the fitting procedure as suggested by Knaup et al.[100]
In current DFTB, the Hubbard parameters (cf. Eq. 1.12) are assumed to be constant, thus neglecting their dependence on atomic charge. However, this dependence is impor-tant for the correct description of deprotonation energies or the total energies of ions.
This deficiency might be overcome by introducing the derivative of the chemical hard-ness with respect to the atomic charge, as indicated by results of Elstner et al.[101]Other routes to improve the accuracy of DFTB are to advance the schemes for the evaluation of the atomic charges ∆qA or to further optimize the γ-functional.
Since DFTB is an approximation based on GGA functionals, it also shares all of the shortcomings of current DFT-GGA functionals (cf. Section 1.1.1), such as the problem of over-polarizability in extended conjugated systems and the problem of van der Waals interactions.[101] The latter problem has been addressed in an “ad hoc” way by adding an empirical dispersion correction to the total energy.[102]