Density-functional tight-binding: basic concepts and applications to molecules and clusters

ABSTRACT

The scope of this article is to present an overview of the Density Functional based Tight Binding (DFTB) method and its applications. The paper introduces the basics of DFTB and its standard formulation up to second order. It also addresses methodological developments such as third order expansion, inclusion of non-covalent interactions, schemes to solve the self-interaction error, implementation of long-range short-range separation, treatment of excited states via the time-dependent DFTB scheme, inclusion of DFTB in hybrid high-level/low level schemes (DFT/DFTB or DFTB/MM), fragment decomposition of large systems, large scale potential energy landscape exploration with molecular dynamics in ground or excited states, non-adiabatic dynamics. A number of applications are reviewed, focusing on -(i)- the variety of systems that have been studied such as small molecules, large molecules and biomolecules, bare or functionalized clusters, supported or embedded systems, and -(ii)- properties and processes, such as vibrational spectroscopy, collisions, fragmentation, thermodynamics or non-adiabatic dynamics. Finally outlines and perspectives are given.

GRAPHICAL ABSTRACT

KEYWORDS:

• DFTB
• electronic structure
• simulation
• molecules
• clusters

PACS:

• 31.15.E–
• 31.15.Ct
• 36.40.–c
• 82.00.00
• 34.00.00

1. Introduction

Since the demonstration by Hohenberg and Kohn [1] of the theoretical grounding of the Density Functional Theory (DFT) [2–4], stating that the energy of any electronic system is a universal functional of the density $\rho$ and the proposal of the Kohn-Sham scheme [5] to find the density, DFT has proved ubiquitous in the theoretical description of electronic system properties of atoms, molecules and condensed matter [6,7]. It has become a choice tool for atomic-scale simulations in Chemistry and Material Science [6–8]. In the Kohn-Sham formulation, the energy of the actual many interacting electrons system is shown to be equivalent to that of a fictitious system of independent electrons within an effective potential involving the interaction with the nuclei (and possibly external potentials) complemented by the electron-electron Coulomb interaction and the exchange-correlation functional $E_{xc}[\rho ]$

(1) $E[\rho ]=\sum _k{n}_k⟨{\phi }_k\left|-\frac{1}{2}\mathrm{\Delta }\right|{\phi }_k⟩+V_{\mathrm{e}\mathrm{x}\mathrm{t}}[\rho ]+\frac{1}{2}\int \frac{\rho (\mathbf{r})\rho (\mathbf{r}\prime )}{\left|\mathbf{r}\prime -\mathbf{r}\right|}{\mathrm{d}}^3\mathbf{r}{\mathrm{d}}^3{\mathbf{r}}^{\prime }+E_{\mathrm{x}\mathrm{c}}[\rho ]+\frac{1}{2}\sum _{a,b}\frac{Z_a{Z}_b}{\left|{\mathbf{R}}_a-{\mathbf{R}}_b\right|}$(1)

The first term is the kinetic energy of independent electrons in orbitals ${\phi }_k$ weighted by their occupation numbers. V_ext is the functional contribution associated with the external potential v_ext. Applying the variational theorem, the resolution is obtained in terms of the mean-field type Kohn-Sham (KS) equation

(2) $\left(-\frac{1}{2}\mathrm{\Delta }+v_{ext}(\mathbf{r})+\int \frac{\rho (\mathbf{r}\prime )}{|\mathbf{r}\prime -\mathbf{r}|}{\mathrm{d}}^3\mathbf{r}\prime +\frac{\delta E_{xc}}{\delta \rho (\mathbf{r})}\right){\phi }_k={ϵ}_k{\phi }_k$(2)

The left hand side of the above equation is the Kohn-Sham operator $H^{KS}=\frac{\delta E}{\delta \rho }$ consisting of the sum of the kinetic contribution and the Kohn-Sham potential $v_{KS}$

(3) $v_{KS}[\rho ]=v_{ext}(\mathbf{r})+\int \frac{\rho (\mathbf{r}\prime )}{|\mathbf{r}\prime -\mathbf{r}|}{\mathrm{d}}^3\mathbf{r}\prime +\frac{\delta E_{xc}}{\delta \rho (\mathbf{r})}$(3)

The density (normalized to the number of electrons) is obtained from the individual orbitals

(4) $\rho (\mathbf{r})=\sum _k{n}_k|{\phi }_k(\mathbf{r}){|}^2$(4)

The Kohn-Sham operator depends on the orbitals via the density and must hence be solved self-consistently. While the Kohn-Sham equation is mathematically very similar to the Hartree-Fock equation, a major difference lies in the fact that it formally incorporates the electron-electron correlation. On the opposite, the Hartree-Fock energy must be complemented by a wavefunction type many-body correlation contribution based on multi-configurational schemes with a generally unfavorable dependence to the number of electrons. Conversely to many-body wavefunctions which are functions of coordinates in space $R^{3N}$, the electronic density is only a function of variables in $R^3$. Hence, the resolution of the KS equation is much simpler and computationally much more efficient than Configuration Interaction type schemes, which explains the success of DFT. Using linear scaling algorithms and High Performance Computing systems, DFT is now able to deal with a few thousands of atoms and a few tens of thousands of electrons at least for a single geometry. Of course, the main theoretical handicap of DFT is that the exchange-correlation functional remains unknown. This brings various drawbacks in many applications of DFT such as the self-interaction error (SIE) [4,9–13], and consequent inherent failures like improper description of the charge localization in extended compounds, ill-behaved dissociation or an incorrect energy derivative with the number of electrons. The account of dispersion forces is also problematic in standard DFT functionals. This situation has led to the proposal of a forest of functionals, some of them taking advantage of theoretical grounding, other empirically determined over reference training sets. This has sometimes questioned the practice of DFT as a first principle theory. Many progresses are currently done to design improved functionals, in particular based around the concept of long-range correction (LC) through a short-range/long-range separation [4,14–17] and its account through double hybrid functionals [17]. Correction of SIE and improvements of functionals are also major challenges in the representation of excited states via the time-dependent version of DFT (TD-DFT), in particular to properly describe Rydberg states or charge transfer excitations.

Despite the favorable computational adaptation of DFT and dedicated progress to achieve linear scaling, there is always a need from the computational point of view for even more efficient techniques. This is the case if one aims at modelling larger systems in the nanoscale domain for instance or running Molecular Dynamics (MD) or Monte Carlo (MC) simulations for medium size systems with the scope of reaching statistical convergence, which requires calculations of energies and energy gradients that must be repeated up to ${10}^6$-${10}^8$ times or even more. The development of approximate schemes, still treating electrons quantum-mechanically, has always been a challenge since the early years of quantum chemistry. There have been essentially two ways for designing such schemes. One is offered by most of the approximate single electron descriptions, which start with very simple elements and can be further complexified in a bottom-up strategy. The second one, more recent and efficient, tends to be theoretically derived in a top-down approximation scheme, from well established mean-field theories, formerly Hartree-Fock and now DFT. It is in this last scheme that the Density Functional based Tight Binding (DFTB) formalism [18–20] has been developed over the two-three decades, now described in a number of review [20–24] or introductory [25] articles. The position of DFTB among other simulation methods in terms of size and simulated time scales is shown in Figure 1. The scope of the present article is (i) to provide an overview of the principles and advances of DFTB in the domain of electronic structure and molecular simulation and (ii) to illustrate applications to molecules, clusters and nanoparticles.

Figure 1. Size and time scales of simulation methods used in chemistry and material science

Section 2 introduces the basic formalism and approximations of DFTB. Section 3 describes developments and extensions such as description of non-covalent forces, improvement of electrostatics, inclusion of DFTB in hybrid methods or determination of electronic excited states. The use of DFTB in large scale simulations (global optimization, dynamics in ground and excited states or thermodynamics) is also commented. After reporting the accuracy of DFTB on small molecules, section 4 overviews applications to more involved classes of systems such as biomolecules, bare or functionalized clusters and nanoparticles, or supported/embedded systems. Note that the number of articles within the DFTB framework is now too large to allow for a fully exhaustive account in the present review article. Hence, the application sections should only be considered as an attempt to provide representative DFTB applications to various fields of chemistry and molecular physics. Finally, outlines and perspectives are given in the last section. Throughout the paper, we will in general use $a$, $b$, $c$, $d$ to label atoms, greek letters $\mu$, $\nu$, $\lambda$, $\tau$ … to label atomic orbitals, $i$, $j$, $k$, $l$ … for molecular orbitals, and capital letters $A$, $B$, $C$ … to address fragment systems. $\mathbf{R}$, ${\mathbf{R}}_a$, and $\mathbf{r}$ will label global nuclei coordinates, nuclei coordinates of atom $a$ and electronic coordinates, respectively.

2. The density-functional based tight-binding approach: basic concepts

2.1. A brief overview of tight-binding theories

Prior to describe the principles of the DFTB method in details, we provide in this subsection a brief general framework for Tight-Binding (TB) theories. Simplified quantum methods for electronic structure rely on several general approximations. A first one concerns the restriction of the Hamiltonian to a subclass of electrons directly involved in the electronic properties of interest. Consideration of the valence electrons only is also related to the physics and chemistry underlying frozen cores and pseudopotential schemes in ab initio calculations. In general, model valence Hamiltonians are defined in linear combination of atomic orbital (LCAO)-type basis sets, so-called minimal in the sense that each valence orbital $\mu$ of atom $a$ is defined by a single atomic function ${\varphi }_{a\mu }$. This is a basic assumption of early quantum semi-empirical methods, as featured by the Hückel [26] or extended-Hückel Hamiltonians [27–30] of quantum chemistry or the tight-binding equivalent in solid state [31–33] and surface physics [34,35] corresponding to one-electron pictures. Restriction to the valence space is also the basis of semi-empirical, multi- or mono-configurational approximations of quantum chemistry such as CNDO [36], MNDO [37], AM1 [38] and PM3 [39]. It remains the basis of the modern tight-binding versions [21,40]. In all these schemes, the basis set is implicit and the Hamiltonian is defined in the matrix form. Transferability and flexibility are accounted for by the dependence of the matrix elements upon geometry [41].

A generic electronic TB Hamiltonian is defined by its matrix elements

(5) $H_{a\mu ,b\nu }=⟨{\varphi }_{a\mu }|\tilde{H}|{\varphi }_{b\nu }⟩$(5)

expressed in the minimal LCAO representation. The diagonal elements have the meaning of effective single-electron atomic energy levels associated with the valence shell atomic orbitals, possibly screened by an effective potential $\tilde{V}$ not necessarily explicited:

(6) $H_{a\mu ,a\mu }=⟨{\varphi }_{a\mu }|-\frac{\mathrm{\Delta }}{2}+\tilde{V}|{\varphi }_{a\mu }⟩={\epsilon }_{a\mu }$(6)

while the interatomic off-diagonal elements between orbitals of different atoms ($a\ne b$), called hopping integrals, describe electron delocalization primarily induced by the (screened) kinetic energy operator

(7) $H_{a\mu ,b\nu }=⟨{\varphi }_{a\mu }|-\frac{\mathrm{\Delta }}{2}+\tilde{V}|{\varphi }_{b\nu }⟩$(7)

The on-site off diagonal elements are generally zero.

In a LCAO non-orthogonal basis, the tight-binding eigenvalue problem is solved for the orbitals ${\phi }_k$ and energies ${\epsilon }_k$

(8) ${\phi }_k=\sum _{a\mu }{c}_{a\mu }^k{\varphi }_{a\mu }$(8)

via the set of secular equations

(9) $\sum _{b\nu }(H_{a\mu ,b\nu }-{\epsilon }_k{S}_{a\mu ,b\nu })c_{b\nu }^k=0,\mathrm{\forall }a\mu$(9)

If the atomic basis functions are supposed to be orthogonal, thus giving rise to the orthogonal tight-binding scheme, the one-electron levels and orbitals are simply obtained by diagonalizing the effective Hamiltonian matrix.

Labelling ${\rho }_{a\mu ,b\nu }$ the one-particle density matrix elements by the elements of the one-particle density matrix $\hat{\rho }$, the sum of the valence electrons energies is

(10) $\sum _k{n}_k{\epsilon }_k=\sum _k\sum _{a\mu ,b\nu }{n}_k{c}_{a\mu }^{k\ast }{c}_{b\nu }^k{H}_{a\mu ,b\nu }=\sum _{a\mu ,b\nu }{\rho }_{a\mu ,b\nu }{H}_{a\mu ,b\nu }$(10)

Finally, the total TB energy can be cast under the very general form, consistent with DFT:

(11) $E[\rho ]=V_{\mathrm{r}\mathrm{e}\mathrm{p}}(\mathbf{R})+\sum _k{n}_k{\epsilon }_k+G[\rho ]$(11)

where $V_{\mathrm{r}\mathrm{e}\mathrm{p}}(\mathbf{R})$ essentially describes the short-range repulsion of the ionic cores, the sum of the single electron energies defines the band energy and the functional contribution of the density $G[\rho ]$ provides an account of all residual contributions, namely the exchange and correlation energies (in particular the dispersion contribution) that are not included in the effective band contribution, as well as the double-counting corrections (the most important being the double counting of Coulomb terms when relevant).

In the simplest version with no electrostatics and no self-consistency included, $\tilde{V}$ is supposed to account for electron screening. In the case of ionic or iono-covalent systems or systems with significant charge fluctuations, interactions between on-site charges can be taken into account, either perturbatively [42–44] or self-consistently [21,40,45–51]. Tight-binding methods may also be considered according to the origin of their parametrization: either semi-empirical tight-binding, where simple functional forms are used for the matrix elements fitted to reproduce ab initio or experimental data, or ab initio tight-binding, where the formalism, functions and inputs are fully derived from first principles references [47].

2.2. From DFT to DFTB

The basic idea of DFTB consists in an expansion of the density $\rho (r)={\rho }_0(r)+\delta \rho (r)$ around a reference density ${\rho }_0(r)$

$\begin{array}{rl}& E[\rho (r)]=E[{\rho }_0(\mathbf{r})]+\int {\left.\frac{\delta E[\rho (\mathbf{r})]}{\delta \rho (\mathbf{r})}\right|}_{{\rho }_0}\delta \rho (\mathbf{r})\\ & +\frac{1}{2}\int \int {\left.\frac{{\delta }^{2}E[\rho (\mathbf{r})]}{\delta \rho (\mathbf{r})\delta \rho ({\mathbf{r}}^{\prime })}\right|}_{{\rho }_0}\delta \rho (\mathbf{r})\delta \rho ({\mathbf{r}}^{\prime })...\\ & +\frac{1}{p!}\int \int ...\int {\left.\frac{\delta E[\rho ({\mathbf{r}}^{\prime })]}{\delta \rho (\mathbf{r})\delta \rho ({\mathbf{r}}^{\prime })...\delta \rho ({\mathbf{r}}^{(\mathbf{p})})}\right|}_{{\rho }_0}\delta \rho (\mathbf{r})\delta \rho ({\mathbf{r}}^{\prime \prime })...\delta \rho ({\mathbf{r}}^{(\mathbf{p})})+...(12)\end{array}$

In current DFTB schemes, the superposed density of the atoms (isolated or in a confined potential) is taken as starting point ${\rho }_0$. Collecting the terms which depend on ${\rho }_0$ only in a so-called repulsive energy contribution, one has

(13) $E^{rep}=E({\rho }_0)-\int {\left.\frac{\delta E}{\delta \rho }\right|}_{{\rho }_0}{\rho }_0(\mathbf{r})d\mathbf{r}$(13)

Using the expression of the Kohn-Sham operator, the terms depending on $\rho$ only provide the so-called band-energy, which was the basis of the initial version of DFTB or DFTB1 [18] (including the above repulsion energy).

(14) $E^{(1)}=\int {\left.\frac{\delta E}{\delta \rho }\right|}_{{\rho }_0}\rho (\mathbf{r})d\mathbf{r}$(14)

The second order dependence upon density fluctuation of the Coulomb and of the exchange-correlation energy only appears in the second order term, namely

(15) $E^{(2)}=\frac{1}{2}\int \int \left(\frac{1}{|\mathbf{r}-\mathbf{r}\prime |}+{\left.\frac{{\delta }^2{E}_{xc}}{\delta \rho (\mathbf{r})\delta \rho (\mathbf{r}\prime )}\right|}_{{\rho }_0}\right)\delta \rho (\mathbf{r})\delta \rho (\mathbf{r}\prime )d\mathbf{r}d\mathbf{r}\prime$(15)

This provides the second order or DFTB2 expansion, namely

(16) $E^{DFTB2}=E^{rep}+\sum _i{n}_i⟨{\phi }_i|H_0^{KS}|{\phi }_i⟩>+\frac{1}{2}\int \int \left(\frac{1}{|\mathbf{r}-\mathbf{r}\prime |}+{\left.\frac{{\delta }^2{E}_{xc}}{\delta \rho (\mathbf{r})\delta \rho (\mathbf{r}\prime )}\right|}_{{\rho }_0}\right)\delta \rho (\mathbf{r})\delta \rho (\mathbf{r}\prime )d\mathbf{r}d\mathbf{r}\prime$(16)

which is the most widely spread DFTB scheme, also called self-consistent charge DFTB (SCC-DFTB) [19,20]. The next step consists in expressing the molecular orbitals as linear combinations of atomic orbitals, consequently defining the matrix elements of the Kohn-Sham operator for the reference density

(17) $H_{a\mu ,b\nu }^0=⟨{\varphi }_{a\mu }|H_0^{KS}|{\varphi }_{b\nu }⟩$(17)

Another approximation consists in replacing the 3D continuous electronic density by a set of discretized atomic electron populations. Assuming a nonpolar expansion of the density fluctuation $\delta \rho (r)$ over the atomic centers

(18) $\delta \rho (\mathbf{r})=\sum _a\mathrm{\Delta }{q}_a{F}_0(\mathbf{r}-{\mathbf{R}}_a)$(18)

the electrostatic situation is described by atomic charges fluctuations $\mathrm{\Delta }{q}_a$ with respect to the atomic neutral references. In the standard versions of DFTB, Mulliken’s charges are used [52]. One should note here that atomic charges are not observables and their definition is arbitrary (see below section 3.1).

One can then introduce the two-electron integrals ${\gamma }_{ab}$ as

(19) ${\gamma }_{ab}=\int \int \left(\frac{1}{|\mathbf{r}-\mathbf{r}\prime |}+{\left.\frac{{\delta }^2{E}_{xc}}{\delta \rho (\mathbf{r})\delta \rho (\mathbf{r}\prime )}\right|}_{{\rho }_0}\right)F_0(\mathbf{r}-{\mathbf{R}}_a)F_0(\mathbf{r}\prime -{\mathbf{R}}_b)d\mathbf{r}d\mathbf{r}\prime$(19)

and the total DFTB2 energy reads

(20) $E=E^{rep}+\sum _i{n}_i\sum _{a\mu ,b\nu }{H}_{a\mu ,b\nu }^0{c}_{a\mu }^i{c}_{b\nu }^i+\frac{1}{2}\sum _{a,b}{\gamma }_{ab}\mathrm{\Delta }{q}_a\mathrm{\Delta }{q}_b$(20)

The next approximation consists in retaining the two-center contributions only in the matrix elements. These terms are then estimated making use of the superposition of pair reference atomic densities ${\rho }_0={\rho }_0^a+{\rho }_0^b$. The second order expression for the KS operator is thus

(21) $H_{a\mu ,b\nu }^{DFTB2}=H_{a\mu ,b\nu }^0+H_{a\mu ,b\nu }^1=H_{a\mu ,b\nu }^0+\frac{1}{2}{S}_{a\mu ,b\nu }\sum _{c\ne a,b}\mathrm{\Delta }{q}_c({\gamma }_{ac}+{\gamma }_{bc})$(21)

Also the repulsive contribution $E^{rep}$ is usually taken as a sum of pair potentials

(22) $E^{rep}=\sum _{a,b}{u}_{ab}^{rep}(|{\mathbf{R}}_a-{\mathbf{R}}_b|)$(22)

Finally, the last standard approximation is to consider minimal valence sets only (although auxiliary bases [53] and extended basis sets [54] have been considered also), namely $s$ set for H and He, $s,p$ set for the second and third row elements, $s,p,d$ set for transition elements and $s,p,d,f$ for rare earths.

The expansion of DFTB was carried out up to third order (DFTB3) by Elstner and co-workers [55]

(23) $E^{(3)}=\frac{1}{6}\int \int \int {\left.\frac{{\delta }^3{E}^{xc}[\rho ]}{\delta \rho (\mathbf{r})\delta \rho (\mathbf{r}\prime )\delta \rho (\mathbf{r}\prime \prime )}\right|}_{{\rho }_0}\delta \rho (\mathbf{r})\delta \rho (\mathbf{r}\prime )\delta \rho (\mathbf{r}\prime \prime )d\mathbf{r}d\mathbf{r}\prime d\mathbf{r}\prime \prime$(23)

$=\frac{1}{6}\int \int \int {\left.\frac{\delta }{\delta \rho ({\mathbf{r}}^{\prime \prime })}\frac{{\delta }^2{E}^{xc}[\rho ]}{\delta \rho (\mathbf{r})\delta \rho ({\mathbf{r}}^{\prime })}\right|}_{{\rho }_0}\delta \rho (\mathbf{r})\delta \rho ({\mathbf{r}}^{\prime })\delta \rho ({\mathbf{r}}^{\prime \prime })d\mathbf{r}d{\mathbf{r}}^{\prime }d{\mathbf{r}}^{\prime \prime }$

Submitting the third order terms to the DFTB approximations (retaining only two body terms) yields the following expression

(24) $E^{(3)}=\frac{1}{6}\sum _{abc}\mathrm{\Delta }{q}_a\mathrm{\Delta }{q}_b\mathrm{\Delta }{q}_c\frac{d{\gamma }_{ab}}{d{q}_c}=\frac{1}{6}\sum _{ab}(\mathrm{\Delta }{q}_a^2\mathrm{\Delta }{q}_b{\mathrm{\Gamma }}_{ab}+\mathrm{\Delta }{q}_a\mathrm{\Delta }{q}_b^2{\mathrm{\Gamma }}_{ba})$(24)

$\mathrm{w}\mathrm{i}\mathrm{t}\mathrm{h}{\mathrm{\Gamma }}_{ab}={\left.\frac{d{\gamma }_{ab}}{d{q}_b}\right|}_{q_a^0}$

The matrix elements of the KS operator are

(25) $H_{a\mu ,b\nu }^{DFTB3}=H_{a\mu ,b\nu }^{DFTB2}+S_{a\mu ,b\nu }\sum _{c\ne a,b}\left[\frac{1}{3}(\mathrm{\Delta }{q}_a{\mathrm{\Gamma }}_{ac}+\mathrm{\Delta }{q}_b{\gamma }_{bc})+\frac{1}{6}\mathrm{\Delta }{q}_c({\mathrm{\Gamma }}_{ca}+{\mathrm{\Gamma }}_{cb})\right]$(25)

This introduces a dependence on the atomic charges via an integral that explicitly depends itself on the other atomic charges. Combined with a modification of the $\gamma$ matrix, DFTB3 was shown [55] to provide an additional flexibility and, in particular, better proton affinities for systems involving C, H, O, N, P and other elements important for chemistry in gas phase or in solvents and, in particular, water. In contrast, DFTB3 only brings a poor improvement of the reaction barriers for proton transfer [55].

2.3. Parametrization issues

The parametrization of the matrix elements $H_{a\mu ,b\nu }^{KS}$ is achieved from DFT calculation. One starts from atomic calculations to determine the atomic KS orbitals ${\varphi }_{a\mu }$ and eigenvalues ${ϵ}_{a\mu }$

(26) $H_{a\mu ,a\mu }^{KS}={ϵ}_{a\mu }$(26)

In principle, the above atomic orbitals could provide the LCAO basis to span the DFTB Hamiltonian. These atomic orbitals are actually constrained by the addition of a confinement potential to the Kohn-Sham atomic operator under the form

(27) $v^{con}={\left(\frac{r}{r_0}\right)}^m$(27)

This confinement potential may yield better transferability. The resolution of the KS equation in the presence of this potential thus defines confined atomic orbitals ${\tilde{\varphi }}_{a\mu }$ which will be taken as the actual DFTB/LCAO basis set.

The overlap integrals S${ }_{a\mu ,b\nu }$ and the off-diagonal elements of the Hamiltonian are determined from the equivalent DFT matrix elements of the atom pairs over the above frozen atomic basis, along the inter-atomic distance $R=|{\mathbf{R}}_a-{\mathbf{R}}_b|$

(28) $H_{a\mu ,b\nu }^0=⟨{\tilde{\varphi }}_{a\mu }|H_0^{KS}|{\tilde{\varphi }}_{b\nu }⟩S_{a\mu ,b\nu }=⟨{\tilde{\varphi }}_{a\mu }|{\tilde{\varphi }}_{b\nu }⟩$(28)

The on-site second order contributions ${\gamma }_{aa}$ are identified with the atom Hubbard parameters $U_a$ and taken as the difference between the first ionization potential (IP) and the electron affinity (EA) of atom $a$

(29) ${\gamma }_{aa}=U_a=IP(a)-EA(a)$(29)

The two-center integrals ${\gamma }_{ab}$ (b$\ne a$) could in principle be calculated numerically from the exact expression provided the atomic charges and the expansion functions are known. In practice, they are expressed via an analytical damped Coulomb formula.

(30) ${\gamma }_{ab}=\frac{1}{R_{ab}}-f(U_a,U_b,R_{ab})$(30)

depending on the on-site integrals $U$_a and $U$_b.

The parametrization of the repulsive term is certainly the most delicate. The initial and somewhat consistent recipe should determine this term as the difference between the purely electronic DFTB contribution to the interaction energy $\mathrm{\Delta }{E}_{ab}^{DFTB(elec)}$ and the total DFT interaction energy $\mathrm{\Delta }{E}_{ab}^{DFT}$ of a given pair of atoms

(31) $u_{ab}^{rep}(R_{ab})=\mathrm{\Delta }{E}_{ab}^{DFT}(R_{ab})-\mathrm{\Delta }{E}_{ab}^{DFTB(elec)}(R_{ab})$(31)

Let us mention a number of attempts to improve the transferability of the parametrization beyond this basic recipe. For instance, constraints on the confinement potentials of the atomic orbitals have been used to optimize bulk electronic band spectra of all elements throughout the periodic table [56,57]. Also several authors have developed automatized algorithms [58–61] to optimize the repulsive terms in multiproperty fits to various ensembles of observables such as molecular binding energies, equilibrium geometries, bulk data band structure, elastic constants or to develop parameters dedicated to specific chemical environment [62]. Some authors also reported on-the-fly parametrization mapping the DFTB parameters on the DFT data during global optimization simulations [63]. Recently, a new scheme has been pioneered with the use of machine-learning algorithms to develop optimized parametrizations [64].

The parameters, most of the time tabulated pointwise, are finally interpolated via spline functions or polynomials. The main parameter sets available are the mio set [20], the matsci set [65], the 3ob set (adapted to DFTB3) [66], the pbc set [67] (adapted to periodic calculations) and that of Wahiduzzaman et al. [56] for the electronic matrix elements throughout the periodic table. Note that there is a dependency between the electronic version of DFTB and the repulsive potentials. In the following, if not specified, DFTB will be used as a generic name referring either to DFTB1, DFTB2 or DFTB3.

3. Extensions of density-functional tight-binding

3.1. Non-covalent interactions

Due to its formulation in minimal basis sets and considering the present quality of the DFT functionals from which it is parametrized, DFTB tends to underestimate or even almost ignore non-covalent contributions to the energy. This includes in particular the polarization energy and the London dispersion energy. In low dimensional systems, such as 1D or 2D systems for instance, whereas the calculation of longitudinal polarizabilities can benefit of the presence of neighboring bases (mediated by the hopping integrals), the calculation of perpendicular polarizabilities may be considerably hindered due to the atomic point charge definition used in the second order term and the absence of basis sets in the orthogonal direction. In addition, the description of electrostatic fluctuations in weakly bound systems may be poorly described via the Mulliken charges. Improvement of electric dipole polarizabilities and polarization energies in the framework of DFTB2 [68] and DFTB3 [53,69] was proposed within the so-called Chemical Potential Equalization (CPE) scheme. The principle is based on an expansion of the energy as a response to the field in the vicinity of the field-less DFTB density

(32) $\begin{array}{rl}& \mathrm{\Delta }{E}^{CPE}=\int {\left[\frac{\delta E}{\delta \rho (\mathbf{r})}\right]}_{{\rho }_{DFTB}}\delta {\rho }^{CPE}(\mathbf{r})d\mathbf{r}\\ & +\frac{1}{2}\int \int {\left[\frac{{\delta }^{2}E}{\delta \rho (\mathbf{r})\delta (\mathbf{r}\prime )}\right]}_{{\rho }_{DFTB}}\delta {\rho }^{CPE}(\mathbf{r})\delta {\rho }^{CPE}(\mathbf{r}\prime )d\mathbf{r}d\mathbf{r}\prime \end{array}$(32)

The response density is itself expanded over $p$-type atomic-centered Gaussian functions

(33) $\delta {\rho }^{CPE}(\mathbf{r})=\sum _j{d}_j{g}_j(\mathbf{r})$(33)

Within the DFTB approximation of charge densities by discrete atomic charges, the minimization of the CPE energy is made via the resolution of a system of linear equations, from which the $d_j$ coefficients are determined. The CPE implementation yields a modification of the Hamiltonian matrix $H_{a\mu ,bv}$

(34) $\mathrm{\Delta }{H}_{a\mu ,b\nu }^{CPE}=\frac{1}{2}{S}_{a\mu ,b\nu }\left(\frac{\mathrm{\partial }\mathrm{\Delta }{E}^{CPE}}{\mathrm{\partial }\mathrm{\Delta }{q}_a}+\frac{\mathrm{\partial }\mathrm{\Delta }{E}^{CPE}}{\mathrm{\partial }\mathrm{\Delta }{q}_b}\right)$(34)

The DFTB3/CPE response was shown to improve intermolecular interactions involving charged and highly polarizable molecules [69].

An alternative scheme for improving polarization can be formulated in analogy with the effective core polarization operators in ab initio treatment. It consists in adding phenomenological atomic contributions to the DFTB energy

(35) $E^{pol}=\sum _a-\frac{1}{2}{\alpha }_a{\left[\sum _{b\ne a}f(R_{ab})\frac{\mathrm{\Delta }{q}_b{\mathbf{R}}_{ab}}{R_{ab}^3}\right]}^2$(35)

This expression accounts for the polarization of atom a due to the resultant electric field created by all other atomic charges. ${\alpha }_a$ is the polarizability (or possibly an effective polarizability) of atom $a$, $f(R_{ab})$ a cut-off function to prevent short distance divergence. $E^{pol}$ can be incorporated in the SCC convergence. It does not require any extra basis but may yield some overestimation of polarization contributions since the atomic polarizability correction is isotropic and may be, at least partially, superfluous (case of longitudinal polarizabilities for instance). Note however that it can be extremely helpful to properly describe MM atoms as polarizable centers in the case of combination of DFTB with MM force fields, for example in the treatment of cryogenic matrices [70].

Continuous theoretical efforts are made to derive DFT functionals describing the London dispersion [71–78]. A more phenomenological approach used in a number of applications [79–84] and systemized by Grimme et al. [85] consists in adding to the total energy specific pair-additive dispersion contributions with $1/R_{ab}^6$, $1/R_{ab}^8$ … long-range behaviour. This empirical approach was first applied for DFTB by Elsner et al. [86]. As, in standard DFTB, the dispersion energy is almost completely absent, due to the reduced basis and the functionals used for parametrization, very little double counting of the dispersion energy is expected. As for polarization, a damping cut-off is necessary to avoid attractive divergence at short distance. The form of the cut-off is strongly related to the parametrization of the repulsive potential [86–88].

As an example, the benzene dimer, unstable at the DFTB2 level, becomes stable when dispersion interactions are added [88]. Benchmarks of intermolecular interactions have been done by Christensen et al. [69] combining DFTB3, CPE and the D3 form of Grimme’s dispersion [89,90]

(36) $E^{disp}=-\sum _{a<b}\sum _{k=3,4}{s}_{2k}\frac{C_{2k}^{ab}}{R_{ab}^{2k}+{[f_{ab}(R_{ab})]}^{2k}}$(36)

with $C_{2k}^{ab}$ the 2k-order dispersion coefficient for the atom pair ab, $s_{2k}$ a scaling factor and $f_{ab}$ a damping function.

Finally, the energy can still be improved by modifying the Coulomb interaction. In its formulation, DFTB makes use of Mulliken definition of atomic point charges to define second and third order terms responsible for the long-range Coulomb interaction between charges fluctuations. This difference with DFT, where the Coulomb interaction is calculated from explicit 3D electronic densities, can be problematic in the case of non-covalently bonded systems, due to a delicate balance between different small contributions in the interaction energy. Among the other definitions of the atomic charges (Bader [91], Löwdin [92], …), the Class IV – Charge Model 3 (CM3) developed by the group of Truhlar [93], easy to implement within the DFTB scheme, corrects the Mulliken charges to take into account a more relevant bond polarization

(37) $\mathrm{\Delta }{q}_a^{CM3}=\mathrm{\Delta }{q}_a^{Mull}+\sum _{b\ne a}^{atoms}[D_{ab}{B}_{ab}+C_{ab}{B}_{ab}^2]$(37)

where $B_{ab}$ is the Mayer’s bond order [94] along bond $ab$ and $D_{ab}$, $C_{ab}$ are empirical parameters. The use of CM3 charges instead of Mulliken charges, first introduced in DFTB as an a posteriori correction of molecular dipole values to compute IR spectra [95], was also shown to improve the long-range Coulomb interactions when used instead of Mulliken charges in DFTB equations [88]. An alternative definition of charges for DFTB was further proposed recently [96]. Let us finally mention that it was also proposed to introduce additional multipoles in the DFTB scheme to describe systems interacting with an electric field [97].

3.2. Spin-polarized DFTB

DFTB was initially formulated within the restricted scheme, corresponding to closed shells in which pairs of electrons $\alpha$ and $\beta$ share the same spatial orbital. DFTB has also been formulated within spin-polarized (unrestricted DFTB) versions [98,99] with possibly different energies ${ϵ}_{i\sigma }$ and orbitals ${\phi }_{i\sigma }$ for different values of the spin-projection $\sigma$. Kohler et al. [98,99] published an atomic shell-resolved formulation. The spin-polarization (magnetization) density $m(\mathbf{r})={\rho }_{\alpha }(\mathbf{r})-{\rho }_{\beta }(\mathbf{r})$ is discretized over the atomic centers and shell-resolved, defining atomic spin-polarization differences $m_{al}=n_{al\alpha }-n_{al\beta }$ ($n_{al\alpha }$ is the electron population with $\alpha$ spin in shell $l$ of atom $a$). Consistently, the charge populations $q_{al}$ and the on-site electron-electron integrals $U_{al}$ become shell-dependent as well as the two-center integrals ${\gamma }_{al,bl\prime }$ which are functions of the $U_{al}$ parameters. The spin-polarized DFTB energy (SDFTB) at second order reads

(38) $E^{DFTB2,spin-pol}=E_{rep}+\sum _{i,\sigma }<{\phi }_{i\sigma }|H^0|{\phi }_{i\sigma }>+\frac{1}{2}\sum _{al,bl\prime }{\gamma }_{al,bl\prime }\mathrm{\Delta }{q}_{al}\mathrm{\Delta }{q}_{bl\prime }+\frac{1}{2}\sum _{all\prime }{m}_{al}{m}_{al\prime }{W}_{all\prime }$(38)

where the $W_{all\prime }$ parameters are shell-dependent atomic constants which, similarly to the Hubbard constants, can be derived from Janak’s theorem [100].

(39) $W_{all\prime }=\frac{1}{2}\left(\frac{\mathrm{\partial }{ϵ}_{al\alpha }}{\mathrm{\partial }{n}_{al\prime \alpha }}-\frac{\mathrm{\partial }{ϵ}_{al\alpha }}{\mathrm{\partial }{n}_{al\prime \beta }}\right)$(39)

The above SDFTB energy corresponds to the electron spin-dependent operator

(40) $\begin{array}{c}{H}_{a\mu ,b\nu ,\sigma }^{DFTB2,spin-pol}=H_{a\mu ,b\nu }^0+\frac{1}{2}{S}_{a\mu ,b\nu }\sum _{c,l\prime \prime \in c}\mathrm{\Delta }{q}_{cl\prime \prime }({\gamma }_{a{l}_{\mu },cl\prime \prime }+{\gamma }_{b{l}_v^{\prime },cl\prime \prime })\\ +{\delta }_{\sigma }\frac{1}{2}{S}_{a\mu ,b\nu }\left(\sum _{l^{\prime }\in a}{m}_{al\prime }{W}_{a{l}_{\mu }l\prime }+\sum _{l^{\prime }\in b}{m}_{b{l}^{\prime }}{W}_{b{l}_{\nu }{l}^{\prime }}\right)\end{array}$(40)

where index $l_{\mu }$ indicates the shell associated with orbital $\mu$ on a given atom.

Note that Melix et al. [101] use a version resolved to atoms only where the spin-polarized DFTB energy is

(41) $E^{DFTB2,spin-pol}=E_{rep}+\sum _{i,\sigma }<{\phi }_{i\sigma }|H^0|{\phi }_{i\sigma }>+\frac{1}{2}\sum _{ab}{\gamma }_{ab}\mathrm{\Delta }{q}_a\mathrm{\Delta }{q}_b+\frac{1}{2}\sum m_a^2{W}_a$(41)

where $W_a$ is now a single atomic constant related to the population derivative of the highest occupied atomic orbital and $m_a=n_{a\alpha }-n_{a\beta }$ the difference between total populations with $\alpha$ and $\beta$ spins on atom $a$.

3.3. Self-interaction correction schemes

Most of standard DFT functionals undergo self-interaction error (SIE) which stems from the fact that the self-exchange contributions in the functionals do not cancel the self-Coulomb contribution. In its original formulation, DFTB meets the same problem. The SIE is responsible for several major errors of standard DFT (and LDA in particular), namely (i) the deviation of the asymptotic potential from $-\frac{1}{r}$ which induces electron overdelocalization, (ii) an underestimation of the HOMO-LUMO gap and (iii) the non linearity and derivative continuity of the energy dependence of the system upon the number of electrons [102].

Several schemes have been proposed to cure the SIE of standard DFT, involving full self-interaction corrections [103], the GW formalism [104], or using hybrid functionals including a part of Hartree-Fock exchange [105]. Other schemes to correct LDA calculations consist in adding corrections $\mathrm{\Delta }{E}^{SIC}$ calculated within the Hubbard model and on-site electron-electron effective interactions $U_a$. This has yielded the LDA+U schemes which have also been declined using $l$-resolved electron-electron screened interactions $U_{al}-J_{al}$ [106]. An alternative so-called pseudo-SIC scheme [107–109] consists in expressing the corrections via the projections of the KS orbitals onto atomic states concerned with the highly correlated shells ($d$ and/or $f$ electrons). Houharine et al. [110] transposed those LDA+U and pseudo-SIC corrections within the spin-polarized DFTB formalism. For example, the pseudo-SIC correction reads

(42) $\mathrm{\Delta }{E}^{pSIC}=-\alpha \sum _a\sum _{l\in a}\frac{(U_{al}-J_{al})}{2}\sum _{\sigma }\sum _{a\mu ,a\nu \in l}(n_{a\mu ,a\nu }^{\sigma }{)}^2$(42)

where $U_{al}-J_{al}$ is taken from atomic DFT calculations and $n_{a\mu ,b\nu }^{\sigma }$ is a matrix generalization of the basis functions Mulliken atomic occupation numbers for a given shell $l$ and a given spin projection $\sigma$. $\alpha$ is here an empirical scaling parameter. Analogous expressions were given for the LDA+U schemes either in the fully localized (FLL) or in the mean-field (AMF) limits. All these corrections rely on the fact that the largest contribution to the SIE is that corresponding to electrons in localized shells. Those contributions to the energy may bring significant improvement. For instance they allow for a gap opening in the strongly correlated antiferromagnetic phase II of bulk NiO, even though the gap remains underestimated. Conversey the corrected magnetic moments show magnitudes comparable with the experimental ones. Further corrections, based on the trace of the idempotent expression $\hat{\rho }-\hat{\rho }S\hat{\rho }$ were proposed to tackle the derivative continuities of the energy as a function of the electron number. Test calculations over several aromatic molecules with CuS substitutive contacts show that such corrections strongly increase the HOMO-LUMO gap which becomes quite consistent with its thermodynamic charge definition $E(N+1)-2E(N)+E(N-1)$.

Another extension of DFTB in relationship with the SIE problem concerns specific classes of systems such as cationic molecular clusters which consist of well identified subsystems. In such cases, delocalization can be strongly overestimated in DFTB as in standard DFT. The single electron picture may also present incorrect dissociation and, since it equally distributes the charge on the separated subsystems (case of two identical subsystems), it may induce spurious Coulomb repulsion at intermediate and long distance separation [10]. Those drawbacks can be circumvented when combining DFTB with Configuration Interaction within a valence bond framework, namely describing the global system via a multiconfigurational wavefunction expanded on charge-localized configurations:

(43) ${\mathrm{\Psi }}_0^{+}=\sum _A{C}_A{a}_A^{HOMO}{\mathrm{\Psi }}_0=\sum _A{C}_A{\mathrm{\Psi }}_A^{+}$(43)

where ${\mathrm{\Psi }}_0$ is the wavefunction of the neutral cluster and $a_A^{HOMO}$ the electron annihilation operator of the HOMO on fragment $A$. The CI problem is then restricted to a secular equation in the charge localized basis

(44) $\sum _B(H_{AB}^{CI}-E_0{S}_{AB}^{CI})C_B=0$(44)

where $H^{CI}$ and $S^{CI}$ are the Hamiltonian and overlap matrices respectively in the charge-localized configurations basis {${\mathrm{\Psi }}_A^{+}$}. The dimension of the CI matrix is only the number of fragments. In this approach, the diagonal terms of the Hamiltonian represent the energies of fragment-localized charge configurations, while the non dynamical correlation arising from the charge resonance and determining the extension of charge fluctuation is mediated by the hopping integrals in the CI resolution. Note that this valence bond CI formulation is well suited to investigate hole transfer through extended system since it provides a naturally quasi-diabatic framework where the hole dynamics is promoted by the hopping integrals [111,112].

A similar partitioning scheme was the principle of the DFTB coarse-grained based approach developed by Elstner et al. [113–116] to study charge transfer in DNA. In this approach, the MOs are calculated independently for each fragment (the fragment orbital approach [117,118]). The diagonal elements are estimated from DFTB2 single particle energies and the hopping term between two fragments is calculated as

(45) $H_{AB}^{CI}=<{\phi }_{HOMO}^A|H^0|{\phi }_{HOMO}^B>$(45)

where ${\phi }_{HOMO}^A$ is the HOMO of the charged fragment $A$ in configuration ${\mathrm{\Psi }}_A^{+}$ and $H^0$ the DFTB1 Kohn-Sham Hamiltonian. The charge mobility in DNA could be described by non-adiabatic MD in a mean field approach with a refined version of this coarse grain model [119].

An alternative scheme in a similar philosophy is that adapted from the constrained-DFT scheme [120–122], in which the orbitals of the charge-localized configurations ${\mathrm{\Psi }}_A^{+}$ are calculated variationally within the DFTB scheme, minimizing a Lagrangian with respect to the orbitals ${\phi }_i^A$ with constraint of charge localization on a given fragment $A$

(46) $\begin{array}{rl}& {\mathcal{L}}^A=E(\{{\phi }_i^A\})-\sum _{ij}{ϵ}_{ij}^A\left(<{\phi }_i^A|{\phi }_j^A>-{\delta }_{ij}\right)\\ & +{\lambda }^A\left(\sum _i<{\phi }_i^A|P^A|{\phi }_j^A>-N_A\right)\end{array}$(46)

where $E(\{{\phi }_i^A\})$ is the DFT energy and the second term ensures the MO orthonormality constraint. The last term is the expression of the charge localization constraint, with ${\lambda }^A$ a Lagrange parameter, $P^A$ a projector of the density on the fragment carrying the charge and $N_A$ the number of electrons fixing the charge localization on fragment $A$. Following Wu and Van Voohris [123,124] the hopping integrals can be computed from the different charge-localized MO coefficients and the Lagrange constraints parameters. The combination of this approach within the DFTB approximations gives the DFTB-CI method [125,126]. This approach differs from the previous coarse-grained one [113] in the sense that each charge localized configuration is calculated self-consistently, thus including relaxation and polarization of the neutral fragments by the charged one. From the computational point of view, the Lagrangian optimization has to be repeated for each fragment, which is more time-consuming than simple DFTB.

3.4. Long-range corrected DFTB

The long-range corrected DFT scheme (LC-DFT) has also been quite fruitful in curing DFT deficiencies. It is based on a range separation of the electron-electron Coulomb interaction. The short-range part is treated via a DFT exchange-correlation functional while the long-range contribution can receive a better treatment, for instance via exact Hartree Fock exchange, contributing to cancellation of the SIE. LC-DFT achievements are obviously more general since they also address issues of long-range correlation either via a higher level correlation functional or even via combinations with Wavefunction type calculations [4,14–16] in order to deal with the dynamical and non-dynamical contributions to electronic correlation. The long-range corrected DFTB scheme (LC-DFTB) was formulated by Lutsker et al. [127] using a Yukawa long-range/short-range type separation of the Coulomb operator

(47) $\frac{1}{r_{12}}=\frac{exp(-\omega r_{12})}{r_{12}}+\frac{1-exp(-\omega r_{12})}{r_{12}}$(47)

This scheme depends on a separation range parameter $\omega$. Using the specific DFTB approximations, the Hamiltonian can be cast as

(48) $H_{a\mu ,b\nu }=H_{a\mu ,b\nu }^0+\frac{1}{4}\sum _{c\lambda ,d\tau }{\rho }_{c\lambda ,d\tau }{S}_{a\mu ,b\nu }{S}_{c\lambda ,d\tau }({\gamma }_{ac}+{\gamma }_{ad}+{\gamma }_{bc}+{\gamma }_{bd})$(48)

(49) $-\frac{1}{8}\sum _{c\lambda ,d\tau }{\rho }_{c\lambda ,d\tau }{S}_{a\mu ,c\lambda }{S}_{d\tau ,b\nu }({\gamma }_{ab}^{lr}+{\gamma }_{cd}^{lr}+{\gamma }_{ad}^{lr}+{\gamma }_{db}^{lr})$(49)

where ${\gamma }_{ab}$ is the two center second order integral calculated with the full Coulomb potential while ${\gamma }_{ab}^{lr}$ is calculated with the long-range part only.

Lutsker et al. [127] benchmarked applications with DFTB parameters extracted from LC-DFT calculations involving the LDA exchange functional and the local PBE form of correlation for a set of organic molecules. They showed that, similarly to LC-DFT schemes, LC-DFTB largely cures the delocalization problem attributed to SIE. As a consequence, a number of properties of the systems are significantly improved, such as the energy of the frontier orbitals, and consequently the estimations of the ionization potentials based on the HOMO energies, the HOMO-LUMO gap, or electrical properties (longitudinal polarizabilities of polyacenes). The LC-DFTB also significantly improves the density of states with respect to photoelectron spectroscopy data. The ordering of the orbitals in delicate cases can still turn out to be incorrect and electron affinities still in default, either due to inherent DFTB approximations (minimal basis set, retain of two-center integrals only) or to the PBE-based parametrization. The improvement of excited electronic states with the LC-DFTB correction is discussed in section 3.6.

3.5. DFTB in hybrid and QM-MM methods

DFTB has also been involved in schemes were the most active atoms/molecules are treated via a higher level quantum-mechanical (QM) scheme while the largest part of the system (large molecule or solvent) is treated at a lower level of approximation, generally via molecular mechanics (MM) potentials or force fields (FF). It should be noted that DFTB, involving two-center approximations, atom-based charges and two-atom repulsive interactions, is very well suited for combination with force fields. The inclusion of point charges in force field is quite straightforward since DFTB is itself based on point charges for the QM atoms. Thus there have been adaptations of DFTB (QM method) within various MM packages such as CHARMM [128], AMBER [129] or GROMACS [130].

Another type of QM-MM combination was adapted to investigate the dynamics of molecules or clusters in a cryogenic environment, namely rare gas inert matrices. This scheme relies on the definition of (possibly) anisotropic two-body interactions between the active atoms and the rare gas atoms added to the DFTB-KS operator in the AO basis, the description of inert atoms interaction (Rg-Rg) via a pair potential, and the inclusion of the polarization response of the Rg atoms. Inclusion of the latter can be handled via atomic polarization operators (see Equation 35)(35) $E^{pol}=\sum _a-\frac{1}{2}{\alpha }_a{\left[\sum _{b\ne a}f(R_{ab})\frac{\mathrm{\Delta }{q}_b{\mathbf{R}}_{ab}}{R_{ab}^3}\right]}^2$(35) which can be finalized adding the following contributions to the initial electrostatic/exchange correlation contributions to the DFTB2 $\gamma$ matrix

(50) ${\gamma }_{ab}^{pol}=-\sum _{c\in Rg}{\alpha }_c{f}_{ac}(R_{ac})f_{bc}(R_{bc})\frac{{\mathbf{R}}_{ab}{\mathbf{R}}_{bc}}{R_{ab}^6{R}_{bc}^6}$(50)

where the rare gas atomic polarization ${\alpha }_c$ and the cut-off functions $f_{ac}(R)$ between active atoms and Rg inert atoms are introduced [70]. Such scheme proved able to describe the influence of the matrix on the structures of molecular complexes such as water clusters in interaction with polycyclic aromatic hydrocarbons [131]. Another combination has also been explored combining DFTB as the low level description with DFT as the high level method [132].

Finally, let us mention that environmental effects can also be taken into account through a polarizable continuum model (for both ground and excited states) [133].

3.6. Excited states and time-dependent DFTB

In the framework of Density Functional Theory, the access to excited states is given by the electronic response, based on the time-dependent Kohn-Sham equation

(51) $i\frac{\mathrm{\partial }{\phi }_j(\mathbf{r},t)}{\mathrm{\partial }t}=\left(-\frac{1}{2}\mathrm{\Delta }+v_{KS}[\rho (\mathbf{r},t)]\right){\phi }_j(\mathbf{r},t)$(51)

The linear response TD-DFTB was originally developed by Niehaus et al. [134] as a DFTB analogue of the linear response TD-DFT [135,136]. Excitation energies are given as the eigenvalues ${\mathrm{\Omega }}_J$ of the following matrix equation

(52) $\left(\begin{array}{cc}\mathbf{A}& \mathbf{B}\\ \mathbf{B}& \mathbf{A}\end{array}\right)\left(\begin{array}{c}\mathbf{X}\\ \mathbf{Y}\end{array}\right)={\mathrm{\Omega }}_J\left(\begin{array}{cc}\mathbf{I}& 0\\ 0& -\mathbf{I}\end{array}\right)\left(\begin{array}{c}\mathbf{X}\\ \mathbf{Y}\end{array}\right)$(52)

where $\mathbf{I}$ is the identity matrix, $\mathbf{A}$ and $\mathbf{B}$ are matrices with the following elements

(53) $A_{ik,jl}=({ϵ}_j-{ϵ}_k){\delta }_{ij}{\delta }_{kl}+2K_{ik,jl}$(53)

(54) $B_{ik,jl}=2K_{ik,jl}$(54)

where indices $i,j$ and $k,l$ label occupied and virtual orbitals respectively, with energies ${ϵ}_i$, ${ϵ}_j$ and ${ϵ}_k$, ${ϵ}_l$. The coupling matrices $\mathbf{K}$, depending on the spin configuration, are determined within the DFTB scheme [134] using the Mulliken approximation to compute transition dipoles. The first application of the linear response TD-DFTB was reported in ref [134]. Absorption spectra were computed for neutral polyacenes ranging in size from naphthalene to heptacene and compared with experimental as well as TD-DFT data. Vibrationally resolved UV/Vis spectra of various aromatic and polar molecules were calculated using TD-DFTB excitation energies and analytical gradients in ref [137]. The results of TD-DFTB were found in a very good agreement with the TD-DFT calculations using local functionals.

Several extensions were developed in the framework of the linear response TD-DFTB. Spin-unrestricted TD-DFTB [138,139] has been implemented in order to study absorption spectra of open-shell systems. Conventional TD-DFTB fails to properly describe PES for charge transfer states. TD-DFTB was combined with LC-DFTB [140–142] to benefit from the range separation improvement for excited states that, in particular, leads to the recovering of a correct -$1/r$ behaviour of the potential. Also, incorporation of intra-atomic exchange integrals [139,143] was shown to improve the transitions energies both towards triplet and singlet TD-DFTB states. Calculation of spin-orbit coupling was interfaced by Gao et al. [144] for TD-DFT approaches, including TD-DFTB. From a computational efficiency point of view, intensity-selected TD-DFTB has been introduced by Rüger et al. [145], delivering similar accuracy as the linear response TD-DFTB, but at a lower computational cost. More details about the TD-DFTB method as well as some other examples of applications can be found in the review paper of T. A. Niehaus [146].

Further improvements were done in order to derive intermolecular excitonic transfer couplings according to the Förster mechanism, implying a formulation of the interaction integral between the transition dipoles of the interacting molecules $A$ and $B$

(55) $J_{AB}^m=\int \int \frac{<{\mathrm{\Psi }}_A^0|\hat{\rho }(\mathbf{r})|{\mathrm{\Psi }}_A^m><{\mathrm{\Psi }}_B^m|\hat{\rho }(\mathbf{r}\prime )|{\mathrm{\Psi }}_B^0>}{|\mathbf{r}-\mathbf{r}\prime |}d\mathbf{r}d\mathbf{r}\prime$(55)

where ${\mathrm{\Psi }}_A^m$ is the intramolecular excited state on $A$ correlated with the exciton band. Within the DFTB formalism this integral becomes [147,148]

(56) $J_{AB}^m=\sum _{a\in A}\sum _{b\in B}{Q}_a^m{\gamma }_{ab}{Q}_b^m$(56)

where quantities $Q_a^m$ are atomic many-body transition charges determined within the TD-DFTB scheme.

Another extension has also been opened for charged molecular clusters in the framework of the DFTB-CI scheme (see above). Initially developed to investigate the ground state, it also delivers excited states as higher roots of the CI matrix. The formalism has been extended in order to provide a better description of the ionic excited states considering in the basis of charged localized configurations, not only the removal of an electron from the HOMO of the charged fragment, but also electron removal from sub-HOMO occupied orbitals ${\phi }_i^A$, yielding a more general wavefunction [149]

(57) ${\mathrm{\Psi }}_0^{+}=\sum _{A,i\in occ}{c}_{Ai}{a}_{Ai}^{†}{a}_{HOMO}{\mathrm{\Psi }}_A^{+}=\sum _{A,i\in occ}{c}_{Ai}{\mathrm{\Psi }}_{Ai}^{+}$(57)

This improvement $vs$ the simple initial scheme restricted to the HOMO orbital becomes important for clusters or stacks of large molecules, presenting a small orbital separation below the HOMO. Moreover, it allows to incorporate not only the excited states of the charge transfer band, but also those correlated with local excitations on the fragment ions, and their coupling. This scheme has been applied to ionic clusters of polyaromatic hydrocarbon molecules and shown to yield satisfactory excited states potential energy surface in the full geometry range up to intermolecular dissociation [149].

3.7. Global exploration of the energy landscape and dynamics

Global exploration of the potential energy surface (PES) or energy landscape is now standard either using Monte Carlo (MC) or molecular dynamics (MD) evolution schemes. While MC only requires the knowledge of the total DFTB energy, the energy gradient is needed in MD. In the widely used DFTB2 approximation, the expression of the gradient is

(58) $\begin{array}{rl}& {\mathbf{F}}_a=-\sum _b\frac{\mathrm{\partial }{u}_{ab}^{rep}}{\mathrm{\partial }{\mathbf{R}}_a}-\sum _i\sum _{a\mu ,b\nu }{n}_i{c}_{a\mu }^i{c}_{b\nu }^i\left(\frac{\mathrm{\partial }{H}_{a\mu ,b\nu }^0}{\mathrm{\partial }{\mathbf{R}}_a}-\left({ϵ}_i-\frac{H_{a\mu ,b\nu }^1}{S_{a\mu ,b\nu }}\right)\frac{\mathrm{\partial }{S}_{a\mu ,b\nu }}{\mathrm{\partial }{\mathbf{R}}_a}\right)\\ & -\mathrm{\Delta }{q}_a\sum _b\frac{{\gamma }_{ab}}{\mathrm{\partial }{\mathbf{R}}_a}\mathrm{\Delta }{q}_b\end{array}$(58)

Note that ground state PES gradients are also available in various extended versions of DFTB such as DFTB3 [55], spin-polarized DFTB [99], CI-DFTB [150] or when LDA+U or pSIC-corrections are included [110].

In large systems like extended and/or flexible molecules, atomic or molecular clusters, structural intuition is delicate, due to the large number of degrees of freedom. Finding the most stable structure (global minimum) and possibly secondary metastable minima might become a challenging task [151] and requires global optimization (GO) schemes with no a priori knowledge of the final structure. A variety of them have been coupled with DFTB and often require the computation of millions of single point energies and possibly gradients for various geometries. A first family of GO schemes rely on genetic algorithms [152] often used to search for atomic cluster structures [63,153–157]. Simulated annealing [158] as well as basin-hopping schemes [159,160] have also often been used either in their standard form [161–163] or improved versions like the modified basin hopping [164,165] or the Tsinghua global minimum algorithms [166]. Other approaches rely on the exploration of the complex potential energy surface (PES) with either MC or MD simulations, which are combined with regular local optimization of the visited geometries as done for ammonium/water clusters [167]. Reaching the bottom of the lowest energy PES basin requires low temperature exploration, but, in such case, the system might be trapped in local minima with vanishing possibility to overcome barriers. An alternative consists in running several simulations at different temperatures [168] and to allow for replica exchange (RE) between the latter following a Boltzmann criterion leading to Parallel Tempering (PT) schemes for MC [169] or MD [170,171]. In the context of DFTB, Parallel-Tempering schemes have appeared quite powerful in finding local minima for atomic and molecular clusters [172–175].

Obviously, MD is also be used to follow the dynamical aspects of the system, for instance to simulate a reaction, collision and/or fragmentation (see section 4.7). A Car-Parrinello version of DFTB molecular dynamics was also implemented [176] as well as biased dynamics schemes like metadynamics [177–179]. Thermodynamical quantities can also be calculated. For instance, DFTB has been combined with the multiple histogram method of Labastie and Whetten [180] to derive the entropy and the heat capacity curves of finite clusters and complexes [181].

IR spectra can be determined in the harmonic approximation, calculating the eigenmodes of the mass-weighted Hessian matrix. However, MD allows to go beyond the harmonic approximation, integrating the IR absorption spectra at finite temperature on-the-fly via the Fourier transform of the autocorrelation function of the electric dipole $\mathit{\mu }$ along the trajectories [182]

(59) $I(\omega )\propto {\omega }^2{\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}\mathrm{d}t⟨\mu (0).\mu (t)⟩e{}^{i\omega t}$(59)

where $<>$ indicates a statistical average to minimize spurious correlations. Let us mention that anharmonic effect can also be obtained from a posteriori treatment of cubic and quartic derivatives of the PES [183,184]. However, the quartic constant can only be obtained at the DFT level for small systems, whereas their computation at the DFTB level could allow for the application of such approaches to larger molecules [185,186].

Finally, recent advances concern the dynamics of excited states. In order to propagate the classical trajectory on a given excited PES, the TD-DFTB excited states energy gradients were developed. The derivation relies on the so-called Z-vector method, which was initially introduced by Furche and Ahlrichs [187,188] to compute analytical forces for the TD-DFT excited states. The procedure was further used to derive TD-DFTB gradients by Heringer et al. [189,190] and led to the final expression published in ref [137].

Non-Adiabatic Molecular Dynamics (NAMD) coupling electronic and nuclear motions has been implemented in the framework of mixed approaches within a DFTB/TD-DFTB quantum description of the electrons and classical nuclei.

Mostly two directions have been followed. In the first approach, the electronic motion is described by the explicit propagation of electronic wavepackets mediated by the time-dependent Schrödinger equation (or equivalent schemes) while the nuclei are propagated in a mean time-dependent electronic potential $E[\rho (\mathbf{R},\mathbf{r},t)]$ (Ehrenfest-like propagation of nuclei).

(60) $M_a\frac{d^2{\mathbf{R}}_{\mathbf{a}}}{d{t}^2}=-{\mathrm{\nabla }}_{a}E[\rho (\mathbf{R}\mathbf{(}\mathbf{t}\mathbf{)},\mathbf{r},t)]$(60)

where $\rho (\mathbf{R}\mathbf{(}\mathbf{t}\mathbf{)},\mathbf{r},t)$ is now the time-dependent electronic density corresponding to molecular orbitals ${\phi }_i(\mathbf{R}\mathbf{(}\mathbf{t}\mathbf{)},\mathbf{r},t)$ obeying the time-dependent DFTB equation. A version of mean potential non-adiabatic dynamics with DFTB was first derived by Niehaus et al. via a variational treatment of the equation of motion (EOM) and the definition of a Lagrangian from which the time-dependant equations can be derived [191]. Mean potential NAMD schemes can also be derived based on the Liouville-von Neumann equation

(61) $i\frac{\mathrm{\partial }\hat{\rho }(\mathbf{R}\mathbf{(}\mathbf{t}\mathbf{)},\mathbf{r},t)}{\mathrm{\partial }t}=[H_{KS},\hat{\rho }(\mathbf{R}\mathbf{(}\mathbf{t}\mathbf{)},\mathbf{r},t)]$(61)

One may cite the NAMD scheme derived by Jakowski [192] and other developments made in the context of electronic transport [193,194].

The second approach relies on the Tully’s Trajectory Surface Hopping (TSH) scheme [195,196]. Here, the motion is propagated on the adiabatic PES of the TD-DFTB excited states, with probabilities to hop between states ${\mathrm{\Psi }}_m$ and ${\mathrm{\Psi }}_n$ determined by the non-adiabatic couplings

(62) $<{\mathrm{\Psi }}_m|\frac{\mathrm{\partial }}{\mathrm{\partial }Q}|{\mathrm{\Psi }}_n>$(62)

along some relevant coordinate $Q$ (possibly a generalized coordinate along the trajectory).

The first article describing methodological as well as development aspects of TSH (in the fewest-switches or FSSH version) coupled with TD-DFTB electronic structure calculation was published by Mitrić et al. [197]. DFTB, as a density functional method, is not initially designed to use wavefunctions to compute properties. Nevertheless, the most common practice is to use the excited state wavefunctions associated with the single excitation configuration interaction (CIS) approximation spanning the TD-DFTB excited states to determine the non-adiabatic couplings presented above [197–201]. This can be achieved through the calculation of the overlap of the CIS electronic wavefunctions between nuclear time steps $t$ and $t+\mathrm{\Delta }t$. This procedure is described within the framework of TD-DFTB by Humeniuk and Mitrić [200]. Several implementations of FSSH are available within various open-source DFTB codes, such as DFTBaby [200], DFTB+ coupled with the NewtonX or PYXAID packages [201,202] and DeMonNano [203].

4. Applications

4.1. Small molecules

Small and medium size molecules can be treated safely via DFT or wavefunction methods. Nevertheless, determination of their ground state properties (structure, energetics, dipole moments, binding energies, vibrational spectra, proton affinities, hydrogen bonds, proton transfer barriers) provides benchmarks for checking the accuracy of DFTB $vs$ DFT, wavefunction calculations (MP2, MP4, Coupled-Cluster or multi-reference CI) or experimental data. Moreover, generic small molecules are often building blocks of larger and/or new systems for which one may expect some transferability. Finally, since reference data are available they also allow to evaluate the various DFTB improvements including the parametrization issues.

In the early DFTB2 versions, the average performances for a set of small organic molecules [204] were found to be 0.017 Å for bond lengths, 2 degrees for bond angles, 5 kcal/mol for dissociation energies and relative errors in the range 6–7 percent on harmonic vibrational frequencies. Recent studies focused on the barrier heights and energetics of reactions with organic molecules [205,206]. The description of the isomers (epimers) of glucose at the DFTB level has also been compared with DFT and wavefunction results: the agreement between structural parameters was shown to be good except when hydrogen bonds are present [207]. The goal was to study large carbohydrate networks which would be out of reach with DFT approaches. Very systematic benchmarks were produced recently to assess the accuracy of the DFTB3 and LC-DFTB2 methods [24,208] covering reference molecule sets. So far, the DFTB3 level appears as the DFTB reference, including benchmarks of proton affinities and hydrogen bonding in organic and biological molecules [209]. Systematic benchmarks of DFTB3 (with the corresponding OB3 parametrization and possibly completed by the addition of the D3 dispersion), LC-DFTB2 with re-optimized parametrization (named OB2${ }^{0.3}$) and DFTB2 (with the $mio$ parameters set), have been recently performed in particular for about 70 neutral closed shells molecules containing C, H, N, and O including the G2/97 set [210]. Structural reference data originate from DFTB-B3LYP/6-31+G(d,p) calculations while the G3B3 data [211] are the reference for energetics. All DFTB methods perform quite well for geometries. The mean absolute deviations (MAD) of bond lenghts $vs$ B3LYP calculations are around 0.01 Å with DFTB2, and around 0.005 Å with DFTB3 and LC-DFTB2, deviations for bond angles are in the range 0.6–0.7 degrees for all methods, while deviations for dihedral angles are within 2–3 degrees. Atomization energies have been compared with the reference data G3B3/MP2 [211]. A net improvement is observed for DFTB3 and LC-DFTB2 methods with a mean deviation of 5–6 kcal/mol $vs$ 20 kcal/mol for the initial DFTB2 scheme. Conversely, all DFTB methods provide deviations of reaction energies in the range 8–10 kcal/mol. Frequencies of selected stretch vibrations show a much better accuracy with DFTB3 and LC-DFTB2 methods, with MADs of 35 and 42 cm${ }^{-1}$ respectively, than with DFTB2 (MAD of 156 cm${ }^{-1}$). Comparisons against the experimental molecular data of the Jorgensen set [212] for geometries and energetics show deviations with the same order of magnitude as above, while the dipole moments deviations are in the range 0.3–0.4 Debye, whatever the DFTB level. Other benchmarks have been done for molecule subsets (closed shells including C, H, O and N atoms only) of the GMTKN0 database [213–216] dedicated to main group thermochemistry and non-covalent interactions of small molecules and even proteins. Errors on a set of reaction energies obtained with DFTB2 and LC-DFTB2 are in the range 0.5–14 kcal/mol, while DFTB3 performs slightly better. Hydrogen binding energies show mean deviations of 3.5 kcal/mol with DFTB3 and 5–6 kcal for LC-DFTB2. Deviations for proton affinities of acidic oxygen (nitrogen) species are 3.7 (17.4) kcal/mol with DFTB2, 3.7 kcal/mol (6.9 kcal/mol and 2.9 kcal/mol, respectively, with modified NH parameters) for DFTB3 and around 8.5 kcal/mol with LC-DFTB2, while proton transfer barriers are in the range 2–3 kcal/mol with DFTB2 and LC-DFTB2 instead of 1 kcal/mol for DFTB3. Finally non-covalent interactions in molecular complexes corresponding to the S66 set [217] were benchmarked against the CCSDT/CBS limit, showing a deviation of 0.82 kcal/mol and around 2.3 kcal/mol for LC-DFTB2 with dispersion.

Other families of molecules outside the above sets have been investigated. Geometries and relative energies were determined for organometallic complexes, the electronic structure of which may be delicate to describe [218–221]. Investigating a series of organometallic complexes with SDFTB2, Zheng et al. [220] estimated an average accuracy of 0.1 Å for bond lengths, 10 degrees for bond angles, finding significant average errors on dissociation energies (25–50 kcal/mol) and on transition energies between spin isomers (10–40 kcal/mol). More recently, it was shown on the example of zinc and manganese complexes [219], that the DFTB3 level (here with $l$-dependent Hubbard integrals) strongly reduces the mean errors down to 0.03 Å for the bond lengths and 2–5 kcal/mol for the energetics, referencing to B3LYP and even G3B3/MP2 data, the largest errors corresponding to interactions of the metal ions with highly charged or polarizable ligands.

One can also cite the specific case of halogens. Kubar et al. [222] benchmarked SDFTB2 parametrization against the experimental CCCBDB database [223] for a series of halogen-containing organic molecules and found absolute errors of 0.045 Å for bond lengths, below 3.6 degrees for bond angles, and 26 and 16  $\mathrm{c}{\mathrm{m}}^{-1}$  for stretching and bending modes respectively. Conversely, reaction energies could present significant errors, in the range 3–30 kcal/mol depending on the type of rearrangement. Kubillus et al. [224] benchmarked DFTB3+D3(X) results against the specific $X40$ halogen database of Rezac et al. [225] and showed that, depending on the halogen atom, DFTB3+D3 could provide mean absolute errors smaller than 0.035 Å and 3 degrees for bond lengths and bond angles respectively, and 25–45 cm${ }^{-1}$ for vibrational frequencies with a larger error ($\approx$ 108 cm${ }^{-1}$) for bromine. Atomization energy errors were found in the range 5–17 kcal/mol, significantly large, however somewhat better than PBE/def2-sv results for Cl and F.

The performance of DFTB regarding the computation of ionization potentials and electron affinities has also been evaluated. Darghouth et al. achieved DFTB3 calculations on a set of small and medium size organic molecules with potential photovoltaic interest [226], comparing with experimental data. Determination of total energy differences ($\mathrm{\Delta }SCF$) gave deviations within $\pm 0.75$ eV and $\pm 0.49$ eV for IPs and EAs respectively, while somewhat better results were obtained when using Koopmans’s theorem, namely IPs and EAs errors within $\pm 0.45$ eV and $\pm$ 0.33 $\mathrm{e}\mathrm{V}$, respectively.

Let us finally mention the case of pure individual carbon clusters, for which the electronic structure, the relative energies and vibrational spectra have been investigated [227–230]. Such systems have sustained a lot of interest due to their relevance in the astrophysical, atmospherical and nanomaterial domains. One can cite for instance the important case of buckminsterfullerene C${ }_{60}$ which has been detected in space. The interest of DFTB for small and medium size molecules is that its efficiency allows the description of large populations of isomers. For instance, an automatic search of benzene isomers has led to the identification of 7000 isomers and 26229 transition structures [231]. DFTB was also used to explore hundreds of thousands of carbon clusters isomers containing 24 to 60 carbon atoms, allowing a classification into structural families and a statistical characterization of their spectroscopic properties (see Figure 2) [232].

Figure 2. Distributions of the samples of C${ }_{60}$ isomers based on specific order parameters. Left: two-dimensional distribution as a function of the $sp2$ hybridization fraction and asphericity parameter $\beta$. The boxes classify the four structural families, cage, planar polycyclic aromatic (PPA), pretzel, and branched. Right: one-dimensional distributions as a function of isomer energy (top), number of 6-member cycles (medium), and gyration radius (bottom). Reproduced from [232]

4.2. Large molecules and biomolecules

One of the main goals behind the development of DFTB was the possibility of modeling systems much larger than those accessible in DFT, while maintaining an electronic scale description of the systems studied. In this framework, many studies have focused on the modeling of nucleic acids and proteins [86,233]. In the case of nucleic acids, most DFTB studies are concerned with the interaction of DNA fragments with different systems. Examples include investigations of the interaction between small DNA fragments and anticancer drugs [234–236], and also between a DNA basis and a carbon nanotube [237]. Charge transport through a short DNA oligomer has also been investigated [238]. It should be noted that some authors have reported that the DFTB2+D method fails to adequately describe deoxyribose and ribose sugar ring pucker [239,240]. In the case of enzymes, studies involving DFTB mainly concern reaction mechanisms carried out using the QM/MM method, with DFTB making it possible to include in the QM reactive zone much more reactive groups than the DFT/MM calculations [24]. The implementation of the DFTB method in codes widely used in hybrid DFT/MM calculations has considerably facilitated access to this method for such hybrid studies. Very different enzymatic mechanisms have been explored, such as proton transfer reactions or proton storage [241,242], histone methylation [243], C-terminal residue cleavage [244], amide hydrolysis [245], glycosylation/deglycosylation [246,247], inactivation of a new tuberculosis target [248], hydrolysis of organophosphorus [51] or proton-coupled electron transfer reactions [249]. One can also cite DFTB studies aimed at investigating substrate promiscuity [250], ion binding and transport by membrane proteins [251], proton distribution over multiple binding sites of a membrane protein [252] or evaluating the pKa of protein residues [253]. The efficiency of the DFTB/MM method also allows the comparison of catalytic pathways [254,255] and the contribution to protein design [256]. Note that it has been reported that, although the DFTB2 method is accurate with regard to protein structure, it sometimes differs from more precise calculations with regard to the electronic states on which it converges [132]. Even the DFTB3 level does not allow a good evaluation of vertical transition energies in the case of the Red Fluorescent protein [257]. Some studies focus on other biologically relevant systems, such as drug [258] or plasma species [259]. To further reduce the computational cost of such biochemical processes studies, different research groups are working at coupling DFTB with linear scaling methods, such as the fragment molecular orbital (FMO) one [260,261].

4.3 Clusters and nanoparticles

DFTB has been used to investigate various clusters including sodium [262], ceria [295], cadmium sulfides [233,264], boron [166], silver and gold [155,157,165,172,173,267–272], ZnO [273], molybdenum disulfide [274], iron [154,275] or nanodiamond [276,277]. In addition to the necessary work dedicated to specific DFTB parametrization for these systems [155,156,172,173,268–270,278], a number of studies have been devoted to their structural characterisation [63,153,154,157,161,165,268,278]. Figure 3 illustrates examples of investigated structures for silver cluster Ag${ }_{561}$ [172]. An interesting question is the evolution with size of the competition between ordered and disordered structures [157,165,173,272]. For instance, global exploration performed at the DFTB level followed by local optimization at the DFT level, suggested that Au${ }_{55}$ presents cavities [173] (recently confirmed by two other DFT studies [63,279]), and showed that the amorphous forms of Au${ }_{147}$ are expected to be more stable than the regular icosahedral ones, or at least very competitive with this latter at low temperature [272] (see Figure 4). Shi et al. evidenced the presence of a core/shell structuration in Pt-Ru alloys [155].

Let us also mention the original approach based on machine learning to correlate the structure/morphology of silver NPs (with diameters up to 4.9 nm) and their electron transfer properties [280]. The magnetic properties of clusters have also been investigated evidencing strong changes with the number of atoms for small iron clusters [154,275].

Figure 3. Icosahedral and cuboctahedral structures of Ag${ }_{561}$

Figure 4. Core-shell like organization of the lowest energy Au${ }_{147}$ isomer (left): surface atoms (middle) and core atoms (right) only. Adapted from reference [272] with the permission of AIP publishing

In addition to atomic clusters, molecular clusters have also been investigated within the DFTB framework. This requires to go beyond simple second order DFTB for a proper treatment of intermolecular interactions including various corrections as describe in section 3. The characterisation of the most stable structures for water clusters provide a picture of the isomer excitation spectra strongly depending on the number of molecules. The ordering found for those isomers with DFTB turn out to be essentially correct. For instance in the DFTB studies of Simon et al., the most stable water octamer is a cube, the next isomer lying 20 kJ.mol${ }^{-1}$ above, whereas the most stable hexamer is a prism followed by 4 other isomers within 9 kJ/mol [281,282]. Interestingly, this structural size dependence induces different thermodynamic behaviors with higher melting temperatures for the octamer than for the hexamer (180K vs 80K [181]). The evolution of water clusters IR spectra with temperature was also investigated [282].

Understanding the interactions between water clusters and molecules is of prime interest as it can be regarded as a step towards the understanding of solvation. Besides, the interaction of water clusters with carbonaceous particles, and in particular polycyclic aromatic hydrocarbons (PAHs), has sustained a lot of interest lately due to their relevance in both atmospherical science and astrochemistry. The PES of water clusters in interaction with planar PAHs was explored with MD [281,282] and PTMD [181] simulations. The lowest energy structures of PAH-(H${ }_2$O)${ }_n$ clusters were determined for planar PAHs [281–283]. Figure 5 reports the lowest energy structures of corannulene (non planar PAH) in interaction with small water clusters C${ }_{20}$H${ }_{10}$-(H${ }_2$O)${ }_n$ (n = 1–8) obtained after PTMD simulations using a similar GO procedure as for C${ }_{16}$H${ }_{10}$-(H${ }_2$O)${ }_n$ clusters [96]. The interaction of the water clusters with the concave face of corannulene is the most energetically favorable, as previously shown for a single water molecule [284]. Interestingly, the water trimer tends to linearize, this is due to its interaction with the edge hydrogens, and such an effect is due to the finite-size of the systems [281,283]. Finite-temperature conformational dynamics of water clusters adsorbed on PAH were also studied [281,282] as well as the influence of PAH adsorption on the IR spectra of water clusters [281–283] and on their thermodynamic properties (heat capacities) [181].

Figure 5. Most stable structures of C${ }_{20}$H${ }_{10}$-(H${ }_2$O)${ }_n$ (n = 1–8) obtained after PTMD/DFTB and local DFTB optimization following the procedure detailed in ref [96]

Water clusters containing impurities, such as ammonium [167] or hydroxyde group [162] have also been considered within DFTB. New isomers were reported in the case of sulfate containing clusters $({\mathrm{H}}_2\mathrm{O})\mathrm{S}{\mathrm{O}}_4^{2-}$ and $({\mathrm{H}}_2\mathrm{O}){\mathrm{H}}_2\mathrm{S}{\mathrm{O}}_4^{2-}$ [174]. The application of DFTB to model protonated water clusters was first reported by Goyal et al. [285]. Korchagina et al. [175] showed that the cluster $({\mathrm{H}}_2\mathrm{O}{)}_{21}{\mathrm{H}}^{+}$ is particularly stable, in agreement with reference calculations [286,287], and present a specific behavior of the heat capacity curves also observed experimentally. The main differences between the IR spectra of pure and protonated water clusters have also been studied [288].

When molecular clusters are singly ionized, alternative DFTB-CI schemes (see section 3.3) may be considered to properly describe the charge and excitation resonance over the different units. Its combination with global exploration schemes allowed to identify the most stable structures of cationic pyrene (Py) clusters, showing that the charge is delocalised over a dimer or trimer core [150], and to compute their electronic spectra [149]. This model was further used to interpret various experiments concerned with thermal evaporation of $\mathrm{P}{\mathrm{y}}_{\mathrm{n}}^{+}$ clusters [289], photodissociation of $\mathrm{P}{\mathrm{y}}_2^{+}$ [290], combined photoionisation and dissociation of $\mathrm{P}{\mathrm{y}}_2$ [291] and the determination of $\mathrm{P}{\mathrm{y}}_{\mathrm{n}}$ ionisation potentials (see Figure 6) [292].

Figure 6. Left: DFTB most stable structure of the cationic pyrene heptamer with all-atom relaxation [150]. Right: experimental and computed ionization potentials for pyrene clusters. Adapted from reference  [292] (https://doi.org/10/1021/acs.jpclett.7b01546, further permissions should be directed to the ACS)

Finally, let us note that the ability of DFTB for describing ionic clusters (clusters of ion pairs) and nanoparticles has recently been reported [293,294].

4.4. Functionalized clusters

The accuracy of the DFTB approach to model bare metal systems, inorganic particles [295–297] as well as organic molecules [205,209] combined with the transferability of the DFTB potentiel over different chemical systems, makes it a valuable tool to describe functionalized clusters and hybrid organic-inorganic systems. Hence, over the last 15 years, this strength of the DFTB approach has led to a number of studies devoted to functionalized clusters.

A number of them focused on metal particles, in particular gold and silver. In the case of gold, the study of thiolates has been of utmost importance as they are often used to stabilise gold particles. In this context, attachment of thiolates on gold clusters were first studied at the DFTB level by Mäkinen et al. [298]. The authors first validated the DFTB approach against experimental and DFT data for three systems: Au${ }_{25}$(SMe)${ }_{18}^{-}$, Au${ }_{102}$(SMe)${ }_{44}$ and Au${ }_{144}$(SMe)${ }_{60}$ and on Au${ }_{102}$(p-MBA)${ }_{44}$ (p-MBA = para-mercaptobenzoic acid). Then, they demonstrated its ability to accurately describe the low-energy structures of Au_m(SMe)_n species as well as qualitatively describe their electronic structure. A similar study was latter conducted by Fihey et al. who developed a new set of DFTB parameters for Au-X (X = Au, H, C, S, N, O) elements in order to better describe the interaction of thiolates and other molecules with gold particles [269]. Those parameters were validated by considering two species: Au${ }_3$SCH${ }_3$ and Au${ }_{25}$SCH${ }_3$ for which structural, energetic and electronic properties were calculated and compared to DFT results. Castro et al. also applied the DFTB approach to describe amino-acids grafted on gold clusters [299]. As for thiolathes, DFTB leads to geometries and adsorption energies that are in good agreement with DFT results, which allowed the authors to study the electron-acceptor and electron-donor character of several amino-acids grafted to gold clusters. In the case of silver, an elegant application of DFTB was conducted by Douglas-Gallardo et al. who tried to rationalize the impact of two adsorbates, water and 1,4-benzoquinone, on the surface plasmon resonance (SPR) band of silver particles of various sizes [300]. This study was a continuation of a previous work devoted to bare icosahedral silver nanoparticles undergoing strong laser pulses [301]. The characteristic of this SPR band, in particular excitation energy and line width, are key in the application of plasmonic particles. However, experiments can have difficulties in probing such properties as they strongly depend on size [302,303], morphology [302–304] and chemical environment [302] of the particles. Combining real-time excited-state dynamics and DFTB, Douglas-Gallardo and co-workers were able to draw a linear relationship between the surface plasmon excitation energy and the inverse cube root of the cluster number of atoms as well as the impact of the adsorbate molecule by studying five different cluster sizes: Ag${ }_{55}$, Ag${ }_{147}$, Ag${ }_{309}$, Ag${ }_{561}$ and Ag${ }_{923}$. In a similar spirit, using real-time excited-state dynamics and DFTB, part of these authors also studied the impact of oxidation on the plasmonic properties of aluminum nanoclusters [305]. To do so, they first simulate the optical absorption spectra of five bare icosahedral aluminum nanostructures: Al${ }_{55}$, Al${ }_{147}$, Al${ }_{309}$, Al${ }_{561}$ and Al${ }_{923}$. Then, focusing on Al${ }_{561}$, MD simulation were performed to describe the structure of Al${ }_{561}$ at different stage of oxidation, from which absorption spectra were re-evaluated. The resulting SPR band displays a red-shift, a broadening and a decrease in intensity that get stronger as oxidation state increases. This was shown to result from the presence of oxygen and not from the symmetry loss.

Figure 7. (a) Schematic representation of the atomic structure of a naphthalenediol-TiO${ }_2$ complex. Superimposed are the corresponding HOMO (red) and LUMO (blue). (b) Time-dependent population of the HOMO and higher-energy orbitals for naphthalenediol-TiO${ }_2$ subject to a continuous laser-type perturbation. Naphthalenediol-TiO${ }_2$ undergoes a direct injection mechanism where population exchange occurs between the HOMO and a manifold of high-energy orbitals. Adapted with permission from [306]. Copyright (2012) American Chemical Society

DFTB has also been applied to model the behavior of dyes grafted on inorganic particles, mainly TiO${ }_2$, under light excitation to understand charge injection mechanisms in dye-sensitized solar cells (DSSC) [306–308]. Indeed, in Grätzel-type solar cells, photoexcitation of the grafted dyes leads to the injection of electrons into the conduction band of the semiconductor. Understanding this mechanism is thus a key step in developing more efficient DSSC. To provide an atomistic-scale description of this process, electron photoinjection was described at the DFTB level for various dyes: alizarin, coumarin C343, derivatives of aniline, naphthalenediol [306], catechol, cresol [307] on a TiO${ }_2$ cluster and 4-nitrophenyl-acetylacetonate and coumarin 343 on a polyoxotitanate particle [308]. Note that Fuertes et al. also studied at the DFTB level the optical properties of bare TiO${ }_2$ particles [266]. These various studies allowed to understand the different steps of the electron transfer from the dye to the inorganic particles for both type I and type II mechanisms and the influence of the excitation wavelength. As a representative example, Figure 7 shows how the electronic structure of a naphthalenediol-TiO${ }_2$ system evolves when subject to a laser-type perturbation. The population exchange between the HOMO and an excited state of the dye followed by an electron transfer to the conduction band of the semiconductor is characteristic of a indirect injection mechanism as opposed to a direct mechanism where the exchange directly occurs from the HOMO of the dye to the semiconductor conduction band [306].

4.5. Supported or embedded systems

DFTB has been widely used to study the adsorption of organic molecules on oxide surfaces. First of all, the adsorption of small molecules such as CO${ }_2$ or NH${ }_3$ on ZnO was studied and the results were found to be in agreement with both DFT and experimental data [296]. Subsequently, the adsorption modes of larger molecules were studied. The grafting of a zwitterionic amino acid (glycine) on a germinal hydroxylated silica surface showed a domination of the adsorption through the carboxylic acid group vs the NH${ }_3^{+}$ one in an explicit water environment [309]. The effect of water has been investigated in the case of the adsorption on TiO${ }_2$ of a serine molecule, an amino acid slightly larger than glycine. It was found that the presence of water weakens the O-Ti bonds and H-bonds existing between the -COO${ }^{-}$/-OH groups and the surface [310]. The effect of the grafting of an organic molecule on the surface gap has also been explored and was found to be negligible in the case of an acetic acid molecule adsorbed either on a crystalline oxide surface (anatase (101), rutile (110) and (B)-TiO${ }_2$ (001)) or on an amorphous one ((a)-TiO${ }_2$) [311]. More recently, the development of new DFTB parameters has also made it possible to study the adsorption of organic molecules on metal surfaces. One can for example mention a study of the adsorption of a corrosion inhibitor (chalcone derivative) on a Fe(110) surface in which the $\pi$ molecular orbitals were found to play a major role in the adsorption phenomenon [312]. DFTB was also developed in order to study adsorption of organic molecules on carbon surfaces, for example transition metal complexes (porphyrin and porphycene) on graphene [313] or small molecules (H${ }_2$O, CH${ }_4$, NH${ }_3$) on defective carbon nanotubes which were all found to physisorb on the nanotubes, except NH${ }_3$ which also chemisorbs [314]. Optical properties of natural pigments (flavonols) adsorbed on boron nitride nanotubes were also analyzed using DFTB (Figure 8) [315]. Some DFTB surface adsorption studies have also given rise to reactivity studies, for example water splitting on anatase (001) [316] or H${ }_2$ dissociation on plutonium [317].

Figure 8. The near site view of the flavonols outside (A) and inside (B) the BNNT surface. The closest contact distance is also shown. Reproduced from [315] with permission of John Wiley and Sons

The adsorption of a PAH on a water ice surface and its influence on the PAH properties are relevant topics for interstellar chemistry. In dense molecular clouds, PAHs are likely to condense on grains covered by H${ }_2$O rich ice mantles with exposure to ionizing radiation [318], and a rich heterogeneous photochemistry on interstellar grains is expected to occur [319]. This motivated experimental studies where PAHs in an icy environment are irradiated with UV-photons leading to the following statements; -(i)- the interaction with the ice leads to a decrease of the ionisation energy of the PAH by 1.5 to 2 eV [320,321] and (ii)- the photo-initiated reactions of PAHs with water on the ice surface [322,323], even at low energy, could be ion-mediated [324]. In this context, Michoulier et al. [96] determined the effect of ice on the ionization energies (IEs) of PAHs using DFTB and constrained DFTB schemes [96] for a series of PAHs from naphthalene (C${ }_{10}$H${ }_8$) to ovalene (C${ }_{32}$H${ }_{14}$) on different types of ices, crystalline (hexagonal Ih and cubic Ic) and amorphous (low density amorphous LDA). They also observed a correlation between the presence (resp. absence) of dangling OH (dOH) bonds interacting with the PAH and the increase (resp. decrease) of the PAH ionisation energy [96]. The conclusion is that the small magnitude of the IE variation, that is at most 0.8 eV for amorphous ice (the experimental type of ice) cannot account for the experimental results. Actually, the electron ejected from the PAH could be transfered to the water ice or recombine with impurities such as the OH radicals. A future theoretical challenge will be to treat such an electron transfer process.

Furthermore, in the astrophysical context, the IR signature of the adsorption of PAH on water ice is an issue of paramount relevance with the imminent launch of the James Webb Space Telescope, which will aim at providing high resolution IR spectra from various regions of the interstellar medium. Therefore, diagnostics for the presence of PAHs condensed on water ice need to be established beforehand. Using the efficiency of DFTB, combinations of harmonic IR spectra of several PAH-amorphous ice systems possessing various PAH-surface interacting structures was computed. The shifts of the dOH bond frequencies induced by the adsorption of the PAH were found to range from −70 to −85 cm${ }^{-1}$ depending on the PAH, in good agreement with experimental results [325]. Further details about the description of water based systems with DFTB can be found in a previous review [326].

Figure 9. Structure of a water dimer interacting with coronene within an argon rare gas matrix subpiece treated via a DFTB-MM scheme [335]

Beyond the adsorption of single molecules, the DFTB method, due to its low computational cost, also allows for the study of extended monolayers. In this framework, the impact of an organic molecule layer on the tunneling current was studied in the case of a PTCDA (3,4,9,9,10-perylene tetracarboxylic dianhydride) monolayer on a 2 × 1 S-passivated GaAs (100) surface. The presence of the layer was found to reduce by one order of magnitude the current with respect to the free surface, in agreement with experimental data [327]. Monolayers (OH, HS and S) were also added in a DFTB study of the sulfidization-amine flotation mechanism of smithsonite in order to model the hydration effect of water and the sulfidization effect on the ZnCO${ }_3$ (101) surface [328]. The structural study of a water monolayer on oxide surfaces has led to several DFTB studies, for example on ZnO [329], on a TiO${ }_2$ anatase surface [330] or on an alumina surface on which it has been found that water dissociates rapidly, leading to an -OH group coverage of about 4.2 groups/nm${ }^2$ [331]. Finally, one can also find a DFTB study of graphene formation on a surface-molten copper surface. In the latter, the authors explains the high quality of a graphene layer grown on Cu by the fact that the high mobility and rapid diffusion of surface Cu atoms induce defect-healing during graphene growth [332].

The deposition of clusters on surfaces has also led to a few studies at the DFTB level. Structural and energetic changes were reported when potassium clusters up to 20 atoms adsorb on a potassium surface K(110) or K(100) [333], the interaction energy being found to dominate the structural reorganization one. MgO supported Au islands were also studied [334,335]. In these islands, the inner atoms were found to remain neutral while the perimeter ones were found to be negatively charged. The specific role played by the peripheral atoms during adsorption and reaction processes was attributed to this charge accumulation coupled with a high density of state.

Finally, structural properties and IR spectroscopy of carbonaceous molecules, water molecules and complexes embedded in cryogenic argon matrix were investigated via the DFTB-MM model described in section 3.5[70,131,335]. The structuration of a water dimer/coronene complex within the argon matrix is illustrated in Figure 9. Fine effects such as the modification of the energetic order of the (H${ }_2$O)${ }_6$ isomers with respect to the gas phase were shown. Besides, MD simulations using the DFTB-MM model allowed to show the influence of (even low) temperature (10 K) on the IR spectrum of a single water molecule embedded in the Ar matrix: red shifts and broadening experimentally observed with respect to the gas phase could be interpreted [335].

4.6. Vibrational spectroscopy

Determining theoretical vibrational spectra of large systems is an important issue as such spectra are among the most popular diagnostics for the presence of species in laboratory experiments, in the earth atmosphere or in space. The determination of vibrational spectra requires the description of charge fluctuation. The use of DFTB2 (possibly with extensions) or DFTB3 thus appears as a convenient approach to compute the vibrational spectra of large molecular systems or clusters as well as the anharmonic effects due to the PES on the spectra.

IR or Raman vibrational spectra can be modeled in the double harmonic approximation. The normal modes are obtained by diagonalizing the full weighted hessian matrix while intensities are obtained by evaluating the variations of the dipole moments (IR) [336] or the changes of the molecular polarizability tensor (Raman) induced by the normal mode oscillations [337]. Vibrational spectra at the DFTB2 level were benchmarked on small molecules with respect to hybrid DFT methods in particular [338], showing that the approach could be used to compute the vibrational spectra of large organic molecules. For instance, the structures of the isomers of oxidized graphene nanoflakes were differentiated by their IR spectra and a correlation was established between stability and IR data [339]. The IR spectra of carbon clusters, either individual structures [227–230,340], or populations (families of isomers) of astrophysical interest [232] have also been determined.

When the internal energy increases or/and when systems exhibit a floppy behaviour, as for instance molecular clusters or systems of biological interest, anharmonic effects due to the shape of the PES are likely to become non negligible. Anharmonic effects on vibrational spectra can be obtained from on-the-fly MD computing the time correlation function of the dipole moment (IR) or of the polarizability (Raman) [341]. The DFTB approach is convenient because long simulations are possible and convergence of spectra in terms of positions and intensities can be reached in reasonable computational time (ns scale) for systems of several tens of atoms [177]. This approach allows to describe the expected redshift of the modes (when no coupling occurs). The example of the out-of-plane CH mode (${\gamma }_{CH}$) of PAHs is quite illustrative. A linear fit of the shift of the latter mode as a function of the internal energy (kinetic temperature) yields the anharmonicity coefficient, the value of which determined at the DFTB level was comparable to the experimental one [342]. This approach was applied to complexes of astrophysical relevance such as SiPAH and FePAH, for which increasing the energy leads to an enhanced motion of the atom (Si, Fe) on the PAH surface [221,342]. In the case of Si, this leads to a merging of the ${\gamma }_{CH}$ modes, that are resonant at two different energies at low temperature and thus induce a deviation from linearity of the function ${\nu }_{{\gamma }_{CH}}(T)$ [342]. Using the same approach, it was found that the influence of the coordination of water clusters on PAHs led to a modification of the anharmonicity of the ${\gamma }_{CH}$ mode, and that this could be a fingerprint of the edge-coordination of the water cluster on the PAH [283]. In the case of a water molecule (described at the DFTB level) surrounded by a rare gas matrix (described with a force field FF), it was shown from MD//DFTB/FF simulations that at low temperature ($\sim$10 K), the water molecule rotates inside the matrix (in agreement with experimental results at low concentration of water), and that leads to red shifts and broadening of the water stretching modes [335] (for a review, see ref [326]).

4.7. Reactivity and fragmentation

The efficiency of DFTB allows for dynamical reactivity studies that can be achieved either through MD/DFTB simulations or through biased molecular dynamics techniques [343] such as umbrella sampling [344] and metadynamics [345]. Statistical convergence on averaged properties can be reached taking into account explicitly the electronic structure for quite large systems. We can cite for instance the unimolecular reactivity of isolated molecular systems in the gas phase such as the isomerisation [177] and dissociation at high energy [346–348] of PAH radical cations. MD/DFTB simulations provide insights into statistical dissociation branching ratios and pathways. The competition between isomerisation and dissociation was shown (see as an illustrative example some isomers and cationic fragments structures of cationic perylene [C${ }_{20}$H${ }_{12}{]}^{+}$ in Figure 10). Comparison with experimental results reporting collision induced dissociation of PAHs [348] or competition between hydrogenation and dissociation of PAHs [347] gave satisfactory results and allowed to cross-benchmark the approaches.

Figure 10. Snapshots retrieved from MD/DFTB simulations describing the evolution of cationic perylene [C${ }_{20}$H${ }_{12}{]}^{+}$ at high energy ($\sim$24–26 eV of internal energy): the formation of a fulvene-type isomer was observed, as well as losses of H, H${ }_2$ and ${\mathrm{C}}_2{\mathrm{H}}_2$, the expected statistical dissociation pathways for PAH radical cations [346]

The low energy conformational dynamics of water clusters, isolated and adsorbed onto a molecular PAH was addressed [181,281,282]. Bimolecular reactions were also investigated via collision dynamics simulations, for instance the collision of H with CO adsorbed on water clusters [349] or the hydrogen uptake of carbon fullerene cages and boron doped heterofullerene [350]. Finally, MD/DFTB simulations at high temperature in simulation chambers were performed to study the growth of carbonaceous systems: formation of large carbonaceous species with various structural orders formed from mixtures of benzene varying the H/C ratio [351–353], growth of carbon nanotubes, possibly catalyzed by a metallic clusters (iron [354,355]), on a SiC surface [356], or formation of metallofullerenes [357].

4.8. Thermodynamics

Some studies have been concerned with the evolution of structural properties with temperature, as well as the determination of the heat capacities of clusters, taking advantage of Parallel Tempering strategies. For instance Choi et al. [163]. simulated the caloric curve of the water octamer. Note that although the qualitative evolution is expected to be well reproduced, one should keep in mind that the simulation results may depend on the type of DFTB and parametrization used [163,358]. A subsequent work was published by Oliveira et al. who redetermined the caloric curves of the water hexamer and heptamer [181]. They also investigated in details the microscopic nature of the phase transition at melting, fingerprinting in particular the evolution of the isomer populations. They furthermore investigated the effect of depositing water clusters on a graphite type substrate modeled as a coronene molecule. Other DFTB thermodynamical studies were concerned with metallic systems and in particular silver and gold clusters. The effect of charge on the doubly magic (electronically closed shell and geometrically a symmetric pyramid) cluster Au${ }_{20}$ was investigated [359] as well as the the correlation between the isomer spectra features and the nature of the solid-to-liquid transition [360], from the comparison between the caloric curves of structurally ordered systems (Au${ }_{20}$, Ag${ }_{55}$) and those of disordered cases (Ag${ }_{20}$, Au${ }_{55}$).

4.9. Dynamics in excited states

Figure 11. Example of FSSH molecular dynamics simulation for neutral polyacenes [203]. Population dynamics averaged over 63 trajectories following excitation to the brightest excited $S_{10}$ state in pentacene (left panel) and hexacene (right panel). Adapted by permission of the PCCP owner societies

The TD-DFTB method was successfully used to study the charge migration in the caffeine molecule induced by an ionizing XUV pulse [361]. In addition to the simulation of exciton dynamics in molecular clusters [200,362–364] reported in section 4.3, the FSSH scheme for non-adiabatic dynamics has been used to simulate excimer formation in the pyrene dimer [365] or relaxation of excited fluorene oligomers [200]. Relaxation dynamics enhanced by transition density analysis has been investigated by Stojanović et al. for two cycloparaphenylene molecules (labelled [8]CPP and [10]CPP) in ref [201]). Other authors have studied the intraband electron and hole relaxation as well as nonradiative electron-hole recombination in a CdSe quantum dot and the (10,5) semiconducting carbon nanotube [202]. The version of FSSH coupled to TD-DFTB in the DeMonNano code was used to investigate the relaxation mechanisms in neutral polyacenes (see Figure 11) ranging in size from naphthalene to heptacene, showing an alternation in decay times of the brightest singlet state with the number of aromatic cycles. More details about the implementation as well as discussion about the observed size effect can be found in ref [203].

Electronic excited states of molecular clusters have also been investigated via DFTB-based schemes. The excitation energy transfer in molecular aggregates has been described through a Frenkel Hamiltonian whose parameters are computed from TD-DFTB [147,148,366,367]. The combination of non-adiabatic dynamics with long-range corrected DFTB [200] has been used to simulate the dynamical evolution of excitons in clusters of tetracene [362] and perylene diimides [363]. The dynamical coupling between local and charge transfer excitons in pentacene clusters was also investigated [364].

Another promising application of DFTB for large metal NPs concerns plasmonics [300,301,368]. For instance, the sub-picosecond breathing-like radial oscillations following a laser pulse excitation have been evidenced for silver NPs up to 309 atoms [301].

5. Outlines and perspectives

The Density Functional based Tight Binding Theory is now more than 25 years old. With respect to many other usual Tight Binding theories, it displays several advantages. One is that it is based on a formal expansion of the energy as a function of the density. Thus, it can be expanded and improved by considering significant terms at higher orders of the expansions, which provides a theoretical basis for upgrade. Being derived from DFT, DFTB exhibits the drawbacks inherent to the former, such as being practically a mean field theory since the exact exchange-correlation functional remains unknown, or suffering from self-interaction errors. In the same time, it has also benefited from many methodological developments adapted from DFT, such as the long-range/short- range separation scheme or the time-dependent version which provides access to excited states, visible/UV spectra and non-adiabatic dynamics. Important initial weaknesses, such as poor treatment of non-covalent interactions, have been cured through various complementary schemes.

DFTB has been now implemented in several packages such as DFTB+ [369], DeMonNano [370], ADF [371], Amber [372], Gromacs [373], Gaussian [374], DFTBaby [200], CP2K [375] where various functionalities are available. Parameters are now available for a large set of elements, even though the problem of the determination and transferability of the repulsive form must still undergo further progress. Of course, many applications have been made for standard atoms C, H, O, N, P, Si, etc. for which the transferability of various DFTB parameter sets has been tested, possibly combined with various versions of DFTB. For other elements, for instance transition metals or even heavier elements, transferability is still to be fully assessed. Machine learning might be useful to finalize the parametrization work [64].

In the domain of DFTB-MM methodologies, combination of DFTB with polarizable force-fields for liquids, and in particular water, would certainly yield a desirable advance for molecules in liquid phase, and even for chemistry with ice. Multi-spatial shell treatments (the active system and a near shell of water molecules treated explicitly with DFTB, the other ones addressed via accurate polarizable force fields) may also improve the study of reactivity in cases where the solvent is likely to participate in the process.

With the development of TD-DFTB and related formalisms, photochemistry and electron transfer processes become feasible for quite large systems. In the field of excited states, an obvious lack concerns Rydberg states which cannot be reached in DFTB, based on valence orbitals only. It could be interesting to include diffuse basis functions that would make at least the low Rydberg states available. Also, DFTB is based on a LCAO expansion and is thus a theory for bound states. As in LCAO-based methods, the continuum is only poorly represented by a discrete set of virtual orbitals, even worse with DFTB. Development of matrix coupling to the continuum could make it able to describe molecular physics processes involving unbound electrons (ionization, electronic attachment).

Figure 12. Box of 350000 water molecules treated via a DFTB cluster division algorithms [376]. Adapted with permission from (J. Chem. Theory Comput. 2014, 106, 2344–2354). Copyright (2014) American Chemical Society

Many important processes involving light atoms require a quantum description of nuclei motion, for instance flexible molecules, reactions associated with proton transfer, water dynamics and ice dynamics. Implementation of quantum dynamics of nuclei via the Path Integral Molecular Dynamics (PIMD) with DFTB electronic structure was reported recently [377]. PIMD yields a system of replicas which multiplies the actual number of degrees of freedom by a factor between 8 and 32, depending on the target accuracy. Development of PIMD within the DFTB framework for highly parallel computing architectures should make nuclei quantum dynamics affordable even for rather large and complex systems in gas phase.

Despite the fact that the present paper is essentially devoted to finite systems, it is important to mention that DFTB in various distributions is implemented in periodic version to address crystals and condensed matter. DFTB offers the possibility to achieve bulk matter simulations using large unit cells (above 10${ }^3$ atoms). This can be of primary importance for investigating the dynamics of default propagation in pure metal and alloys at the microscopic scale. The detailed interaction of atoms, molecules or clusters with surfaces can also be investigated via DFTB within the periodic framework. Deposition of clusters on surfaces may drastically change their structural, spectroscopic, chemical or thermodynamical properties. Such studies also lead to the conception of nanodevices including nanostructuration, nanowires, nanotransport [194]. A neighbouring topic is the collision of atoms or molecules with metal surfaces which may exhibit quite complex electron-surface dynamical coupling involving phonons, plamons and hole-electron pairs excitation. Such complex physics can be addressed by DFTB considering explicitly all the atoms of the active systems and of the surface slabs. Methodological developments can also be thought by combining classical phenomenological description accounting for electron-pair excitation and DFTB via a dissipative dynamics in the ground state [378].

Finally, a word can be said about computational efficiency. Standard DFTB2 is 10${ }^2$ to 10${ }^3$ times faster than even local functionals, and even more if compared with higher-level functionals such as hybrid, double hybrid or LC-corrected functionals. Algorithmic schemes achieving linear scaling with the number of atoms in solving the DFTB Hamiltonian [21,376,379–382] such as the Divide and Conquer techniques [21,381,382] or cluster type algorithms [376] have now proved the feasibility of calculations on extremely large systems up to one million atoms at least for covalent or intermolecular complexes (see Figure 12: a box of 350000 water molecules), even though one should mention that the case of metals remains more delicate due to electronic delocalization. Even if large scale dynamical simulations on such huge systems are not yet practicable, DFTB certainly stands as a promising method to address simulations of systems with up to 10000 atoms on the next generation of High-Performance Computing architectures, which would be quite helpful for theoretical investigation of properties and processes involved in the chemistry and physics of large molecular systems, possibly biomolecules, or in nanoparticle physics.

Correction Statement

This article has been republished with minor changes. These changes do not impact the academic content of the article.

Acknowledgments

The European Union (EU) and Horizon 2020 are acknowledged for funding awarded under the Marie Skłodowska-Curie action to the EUROPAH consortium, grant number 722346 for E.P and TCCM project ITN-EJD-642294 for C. M-A. The authors thank the CALMIP HPC center for continuous support. L. D. acknowledges ERC for support under grant ERC-2013- Syg-610256-NANOCOSMOS and EUR grant NanoX n 080876 in the framework of the Programme des Investissements d’Avenir. The authors would like to thank the colleagues who have permitted reproduction of pictures from their original articles.

Disclosure statement

No potential conflict of interest was reported by the authors.