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

Fernand Spiegelman, Nathalie Tarrat, Jérôme Cuny, Léo Dontot, Evgeny Posenitskiy, Carles Martí, Aude Simon, Mathias Rapacioli

► To cite this version:

Fernand Spiegelman, Nathalie Tarrat, Jérôme Cuny, Léo Dontot, Evgeny Posenitskiy, et al.. Density-functional tight-binding: basic concepts and applications to molecules and clusters. Advances in Physics: X, 2020, 5 (1), pp.1710252. 10.1080/23746149.2019.1710252. hal-04737151v2

HAL Id: hal-04737151

https://hal.science/hal-04737151v2

Submitted on 15 Oct 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Density-functional tight-binding: basic concepts and applications to molecules and clusters

Fernand Spiegelman, Nathalie Tarrat, Jérôme Cuny, Leo Dontot, Evgeny Posenitskiy, Carles Martí, Aude Simon & Mathias Rapacioli

**To cite this article:** Fernand Spiegelman, Nathalie Tarrat, Jérôme Cuny, Leo Dontot, Evgeny Posenitskiy, Carles Martí, Aude Simon & Mathias Rapacioli (2020) Density-functional tight-binding: basic concepts and applications to molecules and clusters, Advances in Physics: X, 5:1, 1710252, DOI: 10.1080/23746149.2019.1710252

**To link to this article:** https://doi.org/10.1080/23746149.2019.1710252

<figure><img src="image_2.png" /></figure>

© 2020 The Author(s). Published by Informa UK Limited, trading as Taylor & Francis Group.

<figure><img src="image_3.png" /></figure>

Published online: 18 Feb 2020.

<figure><img src="image_4.png" /></figure>

Submit your article to this journal

<figure><img src="image_5.png" /></figure>

Article views: 18858

<figure><img src="image_6.png" /></figure>

View related articles

<figure><img src="image_7.png" /></figure>

View Crossmark data

CrossMark

<figure><img src="image_8.png" /></figure>

Citing articles: 45 View citing articles

<figure><img src="image_9.png" /></figure>
<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /></figure>

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

Fernand Spiegelman IDª, Nathalie Tarrat IDᵇ, Jérôme Cuny IDª, Leo Dontotª, Evgeny Posenitskiy IDᶜ, Carles Martí IDª,ᵃ, Aude Simon IDª and Mathias Rapacioli IDª

ªLaboratoire de Chimie et Physique Quantiques LCPQ/IRSAMC, UMR5626, Université de Toulouse (UPS) and CNRS, Toulouse, France; ᵇCEMES, Université de Toulouse (UPS), CNRS, UPR 8011, Toulouse, Toulouse, France; ᶜLaboratoire Collisions Agrégats et Réactivité LCAR/IRSAMC, UMR5589, Université de Toulouse (UPS) and CNRS, Toulouse, France; ᵈLaboratoire de Chimie, UMR5182, Ecole Normale Supérieure de Lyon, Université de Lyon and CNRS, Lyon, France

## 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.

<figure><img src="image_12.png" /></figure>

## ARTICLE HISTORY

Received 17 October 2019
Accepted 19 December 2019

## KEYWORDS

DFTB; electronic structure;
simulation; molecules;
clusters

## PACS

31.15.E-; 31.15.Ct; 36.40.-c;
82.00.00; 34.00.00
<figure><img src="image_13.png" /></figure>

# 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]$

$$E[\rho] = \sum_k n_k \left\langle \varphi_k \left| -\frac{1}{2} \Delta \left| \varphi_k \right\rangle \right\rangle + V_{\text{ext}}[\rho] + \frac{1}{2} \int \frac{\rho(\mathbf{r})\rho(\mathbf{r}')}{|\mathbf{r}' - \mathbf{r}|} d^3\mathbf{r} d^3\mathbf{r}' + E_{xc}[\rho] \\ + \frac{1}{2} \sum_{a,b} \frac{Z_a Z_b}{|\mathbf{R}_a - \mathbf{R}_b|} \qquad (1)$$

The first term is the kinetic energy of independent electrons in orbitals $\varphi_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

$$\left( -\frac{1}{2}\Delta + v_{\text{ext}}(\mathbf{r}) + \int \frac{\rho(\mathbf{r}')}{|\mathbf{r}' - \mathbf{r}|} d^3\mathbf{r}' + \frac{\delta E_{xc}}{\delta\rho(\mathbf{r})} \right) \varphi_k = \epsilon_k \varphi_k \qquad (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}$

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

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

$$\rho(\mathbf{r}) = \sum_k n_k |\varphi_k(\mathbf{r})|^2 \qquad (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
<figure><img src="image_14.png" /></figure>

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.
<figure><img src="image_15.png" /></figure>

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.

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

<figure><img src="image_16.png" /><figcaption>Figure 1. Size and time scales of simulation methods used in chemistry and material science.</figcaption></figure>
<figure><img src="image_17.png" /></figure>

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 μ, ν, λ, τ ... to label atomic orbitals, i, j, k, l ... for molecular orbitals, and capital letters A, B, C ... to address fragment systems. R, Ra, and r will label global nuclei coordinates, nuclei coordinates of atom a and electro- nic 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 μ of atom a is defined by a single atomic function $\phi_{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

$$H_{a\mu,bv} = \langle \phi_{a\mu} | \tilde{H} | \phi_{bv} \rangle \quad (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 explicit:
<figure><img src="image_18.png" /></figure>

$$H_{a\mu,a\mu} = \langle \phi_{a\mu} | -\frac{\Delta}{2} + \tilde{V} | \phi_{a\mu} \rangle = \epsilon_{a\mu} \quad (6)$$

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

$$H_{a\mu,b\nu} = \langle \phi_{a\mu} | -\frac{\Delta}{2} + \tilde{V} | \phi_{b\nu} \rangle \quad (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$

$$\phi_k = \sum_{a\mu} c_{a\mu}^k \phi_{a\mu} \quad (8)$$

via the set of secular equations

$$\sum_{b\nu} (H_{a\mu,b\nu} - \epsilon_k S_{a\mu,b\nu}) c_{b\nu}^k = 0, \forall a\mu \quad (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

$$\sum_k n_k \epsilon_k = \sum_k \sum_{a\mu,b\nu} n_k c_{a\mu}^{k*} c_{b\nu}^k H_{a\mu,b\nu} = \sum_{a\mu,b\nu} \rho_{a\mu,b\nu} H_{a\mu,b\nu} \quad (10)$$

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

$$E[\rho] = V_{\text{rep}}(\mathbf{R}) + \sum_k n_k \epsilon_k + G[\rho] \quad (11)$$

where $V_{\text{rep}}(\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
<figure><img src="image_19.png" /></figure>

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 $ρ(r) = ρ_0(r) + δρ(r)$ around a reference density $ρ_0(r)$

$$ 
\begin{aligned}
E[ρ(r)] = E[ρ_0(\mathbf{r})] &+ \left. \frac{\delta E[ρ(\mathbf{r})]}{\delta ρ(\mathbf{r})} \right|_{ρ_0} δρ(\mathbf{r}) \\
&+ \frac{1}{2} \left. \int \int \frac{\delta^2 E[ρ(\mathbf{r})]}{\delta ρ(\mathbf{r}) δρ(\mathbf{r}')}
\right|_{ρ_0} δρ(\mathbf{r}) δρ(\mathbf{r}') \dots \\
&+ \frac{1}{ρ!} \int \dots \int \left. \frac{\delta E[ρ(\mathbf{r}')]}{\delta ρ(\mathbf{r}) δρ(\mathbf{r}') \dots δρ(\mathbf{r}^{(p)})}
\right|_{ρ_0} δρ(\mathbf{r}) δρ(\mathbf{r}'') \dots δρ(\mathbf{r}^{(p)}) + \dots
\end{aligned} 
\quad (12) $$

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

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

Using the expression of the Kohn-Sham operator, the terms depending on $ρ$ 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).

$$ E^{(1)} = \left. \frac{\delta E}{\delta \rho} \right|_{\rho_0} \rho(\mathbf{r}) d\mathbf{r} \quad (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

$$ E^{(2)} = \frac{1}{2} \iint \left( \frac{1}{|\mathbf{r} - \mathbf{r}'|} + \frac{\delta^2 E_{xc}}{\delta \rho(\mathbf{r}) \delta \rho(\mathbf{r}'')} \right) \delta \rho(\mathbf{r}) \delta \rho(\mathbf{r}') d\mathbf{r} d\mathbf{r}' \quad (15) $$

This provides the second order or DFTB2 expansion, namely
<figure><img src="image_20.png" /></figure>

$$E^{\text{DFTB2}} = E^{\text{rep}} + \sum_{i} n_{i} \langle \varphi_{i} | H_{0}^{\text{KS}} | \varphi_{i} \rangle + \frac{1}{2} \iint \left( \frac{1}{|\mathbf{r} - \mathbf{r}'|} + \frac{\delta^{2} E_{xc}}{\delta\rho(\mathbf{r})\delta\rho(\mathbf{r}')} \bigg|_{\rho_{0}} \right) \delta\rho(\mathbf{r})\delta\rho(\mathbf{r}') d\mathbf{r} d\mathbf{r}' \quad (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

$$H_{a\mu,b\nu}^{0} = \langle \phi_{a\mu} | H_{0}^{KS} | \phi_{b\nu} \rangle \quad (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

$$\delta\rho(\mathbf{r}) = \sum_{a} \Delta q_{a} F_{0}(\mathbf{r} - \mathbf{R}_{a}) \quad (18)$$

the electrostatic situation is described by atomic charges fluctuations $\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

$$\gamma_{ab} = \iint \left( \frac{1}{|\mathbf{r} - \mathbf{r}'|} + \frac{\delta^2 E_{xc}}{\delta\rho(\mathbf{r})\delta\rho(\mathbf{r}')} \bigg|_{\rho_0} \right) F_0(\mathbf{r} - \mathbf{R}_a) F_0(\mathbf{r}' - \mathbf{R}_b) d\mathbf{r} d\mathbf{r}' \quad (19)$$

and the total DFTB2 energy reads

$$E = E^{\text{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} \Delta q_{a} \Delta q_{b} \quad (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

$$H_{a\mu,b\nu}^{\text{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 \neq a,b} \Delta q_c (\gamma_{ac} + \gamma_{bc}) \quad (21)$$
<figure><img src="image_21.png" /></figure>

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

$$E^{rep} = \sum_{a,b} u_{ab}^{rep} (|\mathbf{R}_a - \mathbf{R}_b|) \qquad (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]

$$E^{(3)} = \frac{1}{6} \iint \int \frac{\delta^3 E^{xc}[\rho]}{\delta\rho(\mathbf{r})\delta\rho(\mathbf{r}')\delta\rho(\mathbf{r}'')} \Big|_{\rho_0} \delta\rho(\mathbf{r})\delta\rho(\mathbf{r}')\delta\rho(\mathbf{r}'') d\mathbf{r} d\mathbf{r}' d\mathbf{r}'' \quad (23) \\ = \frac{1}{6} \iint \int \frac{\delta}{\delta\rho(\mathbf{r}'')} \frac{\delta^2 E^{xc}[\rho]}{\delta\rho(\mathbf{r})\delta\rho(\mathbf{r}')} \Big|_{\rho_0} \delta\rho(\mathbf{r})\delta\rho(\mathbf{r}')\delta\rho(\mathbf{r}'') d\mathbf{r} d\mathbf{r}' d\mathbf{r}''$$

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

$$E^{(3)} = \frac{1}{6} \sum_{abc} \Delta q_a \Delta q_b \Delta q_c \frac{d\gamma_{ab}}{dq_c} = \frac{1}{6} \sum_{ab} (\Delta q_a^2 \Delta q_b \Gamma_{ab} + \Delta q_a \Delta q_b^2 \Gamma_{ba}) \quad (24)$$

with $\Gamma_{ab} = \frac{d\gamma_{ab}}{dq_b} \Big|_{q_a^0}$

The matrix elements of the KS operator are

$$H_{a\mu,bv}^{DFTB3} = H_{a\mu,bv}^{DFTB2} + S_{a\mu,bv} \sum_{c \neq a,b} \left[ \frac{1}{3} (\Delta q_a \Gamma_{ac} + \Delta q_b \gamma_{bc}) + \frac{1}{6} \Delta q_c (\Gamma_{ca} + \Gamma_{cb}) \right] \quad (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].
<figure><img src="image_22.png" /></figure>

## 2.3. Parametrization issues

The parametrization of the matrix elements $H_{a\mu,bv}^{KS}$ is achieved from DFT calculation. One starts from atomic calculations to determine the atomic KS orbitals $\phi_{a\mu}$ and eigenvalues $\epsilon_{a\mu}$

$$H_{a\mu,a\mu}^{KS} = \epsilon_{a\mu} \qquad (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

$$v^{\text{con}} = \left(\frac{r}{r_0}\right)^m \qquad (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{\phi}_{a\mu}$ which will be taken as the actual DFTB/LCAO basis set.

The overlap integrals $S_{a\mu,bv}$ 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|$

$$H_{a\mu,bv}^{0} = \langle \tilde{\phi}_{a\mu} | H_{0}^{KS} | \tilde{\phi}_{bv} \rangle \qquad S_{a\mu,bv} = \langle \tilde{\phi}_{a\mu} | \tilde{\phi}_{bv} \rangle \qquad (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

$$\gamma_{aa} = U_a = IP(a) - EA(a) \qquad (29)$$

The two-center integrals $\gamma_{ab}$ ($b \neq 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.

$$\gamma_{ab} = \frac{1}{R_{ab}} - f(U_a, U_b, R_{ab}) \qquad (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 $\Delta E_{ab}^{\text{DFTB(elec)}}$ and the total DFT interaction energy $\Delta E_{ab}^{\text{DFT}}$ of a given pair of atoms
<figure><img src="image_23.png" /></figure>

$$u_{ab}^{rep}(R_{ab}) = \Delta E_{ab}^{DFT}(R_{ab}) - \Delta E_{ab}^{DFTB(elec)}(R_{ab}) \quad (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 multiproduct 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
<figure><img src="image_24.png" /></figure>

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

$$ \Delta E^{\text{CPE}} = \int \left[ \frac{\delta E}{\delta \rho(\mathbf{r})} \right]_{\rho_{\text{DFTB}}} \delta \rho^{\text{CPE}}(\mathbf{r}) d\mathbf{r} + \frac{1}{2} \iint \left[ \frac{\delta^2 E}{\delta \rho(\mathbf{r}) \delta(\mathbf{r}')} \right]_{\rho_{\text{DFTB}}} \delta \rho^{\text{CPE}}(\mathbf{r}) \delta \rho^{\text{CPE}}(\mathbf{r}') d\mathbf{r} d\mathbf{r}' \quad (32) $$

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

$$ \delta \rho^{\text{CPE}}(\mathbf{r}) = \sum_{j} d_{j} g_{j}(\mathbf{r}) \quad (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}$

$$ \Delta H_{a\mu,bv}^{\text{CPE}} = \frac{1}{2} S_{a\mu,bv} \left( \frac{\partial \Delta E^{\text{CPE}}}{\partial \Delta q_a} + \frac{\partial \Delta E^{\text{CPE}}}{\partial \Delta q_b} \right) \quad (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

$$ E^{\text{pol}} = \sum_{a} -\frac{1}{2} \alpha_{a} \left[ \sum_{b \neq a} f(R_{ab}) \frac{\Delta q_{b} \mathbf{R}_{ab}}{R_{ab}^{3}} \right]^{2} \quad (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].
<figure><img src="image_25.png" /></figure>

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]

$$E^{disp} = -\sum_{a<b} \sum_{k=3,4} s_{2k} \frac{C_{2k}^{ab}}{R_{ab}^{2k} + [f_{ab}(R_{ab})]^{2k}} \quad (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

$$\Delta q_a^{CM3} = \Delta q_a^{Mull} + \sum_{b \neq a}^{atoms} [D_{ab}B_{ab} + C_{ab}B_{ab}^2] \quad (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
<figure><img src="image_26.png" /></figure>

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 α and β share the same spatial orbital. DFTB has also been formulated within spin-polarized (unrestricted DFTB) versions [98,99] with possibly different energies $\epsilon_{i\sigma}$ and orbitals $\varphi_{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'}$ which are functions of the $U_{al}$ parameters. The spin-polarized DFTB energy (SDFTB) at second order reads

$$E^{\text{DFTB2,spin-pol}} = E_{\text{rep}} + \sum_{i,\sigma} \langle \varphi_{i\sigma} | H^0 | \varphi_{i\sigma} \rangle + \frac{1}{2} \sum_{al,bl'} \gamma_{al,bl'} \Delta q_{al} \Delta q_{bl'} \\ + \frac{1}{2} \sum_{all'} m_{al} m_{all'} W_{all'} \quad (38)$$

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

$$W_{all'} = \frac{1}{2} \left( \frac{\partial \epsilon_{al\alpha}}{\partial n_{al'\alpha}} - \frac{\partial \epsilon_{al\beta}}{\partial n_{al'l\beta}} \right) \quad (39)$$

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

$$H_{a\mu,b\nu,\sigma}^{\text{DFTB2,spin-pol}} = H_{a\mu,b\nu}^{0} + \frac{1}{2} S_{a\mu,b\nu} \sum_{c,l'' \in c} \Delta q_{cl''} (\gamma_{a\mu,cl''} + \gamma_{bl'',cl''}) \\ + \delta_{\sigma} \frac{1}{2} S_{a\mu,b\nu} \left( \sum_{l' \in a} m_{al'} W_{al\mu l'} + \sum_{l' \in b} m_{bl'} W_{bl\nu l'} \right) \quad (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
<figure><img src="image_27.png" /></figure>

$$E^{\text{DFTB2,spin-pol}} = E_{\text{rep}} + \sum_{i,\sigma} \langle \varphi_{i\sigma} | H^0 | \varphi_{i\sigma} \rangle + \frac{1}{2} \sum_{ab} \gamma_{ab} \Delta q_a \Delta q_b \\ + \frac{1}{2} \sum m_a^2 W_a \qquad (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_{aa} - n_{ab}$ 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 $\Delta E^{\text{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

$$\Delta E^{\text{pSIC}} = -\alpha \sum_a \sum_{l \in a} \frac{(U_{al} - J_{al})}{2} \sum_{\sigma} \sum_{a\mu,av \in l} (n_{a\mu,av}^{\sigma})^2 \qquad (42)$$

where $U_{al} - J_{al}$ is taken from atomic DFT calculations and $n_{a\mu,bv}^{\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
<figure><img src="image_28.png" /></figure>

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}\hat{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:

$$ \Psi_0^+ = \sum_A C_A a_A^{\text{HOMO}} \Psi_0 = \sum_A C_A \Psi_A^+ \quad (43) $$

where $\Psi_0$ is the wavefunction of the neutral cluster and $a_A^{\text{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

$$ \sum_B (H_{AB}^{\text{CI}} - E_0 S_{AB}^{\text{CI}}) C_B = 0 \quad (44) $$

where $H_{AB}^{\text{CI}}$ and $S_{AB}^{\text{CI}}$ are the Hamiltonian and overlap matrices respectively in the charge-localized configurations basis {$\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

$$H_{AB}^{CI} = \langle \varphi_{HOMO}^{A} | H^{0} | \varphi_{HOMO}^{B} \rangle \quad (45)$$

where $\varphi_{HOMO}^{A}$ is the HOMO of the charged fragment A in configuration $\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 $\Psi_{A}^{+}$ are calculated variationally within the DFTB scheme, minimizing a Lagrangian with respect to the orbitals $\varphi_{i}^{A}$ with constraint of charge localization on a given fragment A

$$\mathcal{L}^A = E(\{\varphi_i^A\}) - \sum_{ij} \epsilon_{ij}^A (\langle \varphi_i^A | \varphi_j^A \rangle - \delta_{ij}) \\ + \lambda^A \left( \sum_i \langle \varphi_i^A | P^A | \varphi_j^A \rangle - N_A \right) \quad (46)$$

where $E(\{\varphi_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.
<figure><img src="image_29.png" /></figure>

### 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

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

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

$$ 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}) \quad (48) $$

$$ -\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}) \quad (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 where 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) which can be finalized adding the following contributions to the initial electrostatic/exchange correlation contributions to the DFTB2 $\gamma$ matrix

$$ \gamma_{ab}^{pol} = -\sum_{c \in R_g} \alpha_c f_{ac}(R_{ac}) f_{bc}(R_{bc}) \frac{\mathbf{R}_{ab} \mathbf{R}_{bc}}{R_{ab}^6 R_{bc}^6} \quad (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].
<figure><img src="image_30.png" /></figure>

### 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

$$i \frac{\partial \varphi_j(\mathbf{r}, t)}{\partial t} = \left( -\frac{1}{2} \Delta + v_{KS} [\rho(\mathbf{r}, t)] \right) \varphi_j(\mathbf{r}, t) \quad (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 $\Omega_j$ of the following matrix equation

$$\begin{pmatrix} \mathbf{A} & \mathbf{B} \\ \mathbf{B} & \mathbf{A} \end{pmatrix} \begin{pmatrix} \mathbf{X} \\ \mathbf{Y} \end{pmatrix} = \Omega_j \begin{pmatrix} \mathbf{I} & \mathbf{0} \\ \mathbf{0} & -\mathbf{I} \end{pmatrix} \begin{pmatrix} \mathbf{X} \\ \mathbf{Y} \end{pmatrix} \quad (52)$$

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

$$A_{ik,jl} = (\epsilon_j - \epsilon_k)\delta_{ij}\delta_{kl} + 2K_{ik,jl} \quad (53)$$

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

where indices $i, j$ and $k, l$ label occupied and virtual orbitals respectively, with energies $\epsilon_i, \epsilon_j$ and $\epsilon_k, \epsilon_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

$$J_{AB}^{m} = \iint \frac{\langle \Psi_{A}^{0} | \hat{\rho}(\mathbf{r}) | \Psi_{A}^{m} \rangle - \langle \Psi_{B}^{m} | \hat{\rho}(\mathbf{r}') | \Psi_{B}^{0} \rangle}{|\mathbf{r} - \mathbf{r}'|} d\mathbf{r} d\mathbf{r}' \quad (55)$$

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

$$J_{AB}^{m} = \sum_{a \in A} \sum_{b \in B} Q_{a}^{m} \gamma_{ab} Q_{b}^{m} \quad (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]

$$\Psi_{0}^{+} = \sum_{A, i \in occ} c_{Ai} a_{Ai}^{\dagger} a_{HOMO} \Psi_{A}^{+} = \sum_{A, i \in occ} c_{Ai} \Psi_{Ai}^{+} \quad (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].
<figure><img src="image_31.png" /></figure>

### 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

$$ F_a = -\sum_b \frac{\partial u_{ab}^{rep}}{\partial \mathbf{R}_a} - \sum_i \sum_{a\mu,bv} n_i c_{a\mu}^i c_{bv}^i \left( \frac{\partial H_{a\mu,bv}^0}{\partial \mathbf{R}_a} - \left( \epsilon_i - \frac{H_{a\mu,bv}^1}{S_{a\mu,bv}} \right) \frac{\partial S_{a\mu,bv}}{\partial \mathbf{R}_a} \right) \\ - \Delta q_a \sum_b \frac{\gamma_{ab}}{\partial \mathbf{R}_a} \Delta q_b \qquad (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 $\mu$ along the trajectories [182]

$$I(\omega) \propto \omega^2 \int_{-\infty}^{+\infty} dt \langle \mu(0) \cdot \mu(t) \rangle e^{i\omega t} \quad (59)$$

where $\langle \rangle$ 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).

$$M_a \frac{d^2 \mathbf{R}_a}{dt^2} = -\nabla_a E[\rho(\mathbf{R}(t), \mathbf{r}, t)] \quad (60)$$
<figure><img src="image_32.png" /></figure>

where $\rho(\mathbf{R}(t), \mathbf{r}, t)$ is now the time-dependent electronic density corresponding to molecular orbitals $\phi_i(\mathbf{R}(t), \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-dependent equations can be derived [191]. Mean potential NAMD schemes can also be derived based on the Liouville-von Neumann equation

$$i \frac{\partial \hat{\rho}(\mathbf{R}(t), \mathbf{r}, t)}{\partial t} = [H_{KS}, \hat{\rho}(\mathbf{R}(t), \mathbf{r}, t)] \quad (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 $\Psi_m$ and $\Psi_n$ determined by the non-adiabatic couplings

$$\langle \Psi_m | \frac{\partial}{\partial Q} | \Psi_n \rangle \quad (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 + \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
<figure><img src="image_33.png" /></figure>

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 OB20.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 lengths 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⁻¹ respectively, than with DFTB2 (MAD of 156 cm⁻¹). 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.
<figure><img src="image_34.png" /></figure>

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 $cm^{-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 (≈ 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.
<figure><img src="image_35.png" /></figure>

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 ($\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 eV, 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><img src="image_36.png" /><figcaption>**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].</figcaption></figure>
<figure><img src="image_37.png" /></figure>

## 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
<figure><img src="image_38.png" /></figure>

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 structure in Pt-Ru alloys [155].

<figure><img src="image_39.png" /><figcaption>**Figure 3.** Icosahedral and cuboctahedral structures of $Ag_{561}$.</figcaption></figure>

<figure><img src="image_40.png" /><figcaption>**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.</figcaption></figure>
<figure><img src="image_41.png" /></figure>

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].

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 described 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⁻¹ 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₂O)ₙ 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₂₀H₁₀-(H₂O)ₙ (n = 1–8) obtained after PTMD simulations using a similar GO procedure as for C₁₆H₁₀-(H₂O)ₙ 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].

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 (H₂O)SO₄²⁻ and (H₂O)H₂SO₄²⁻ [174]. The application of DFTB to
<figure><img src="image_42.png" /></figure>

<figure><img src="image_43.png" /><figcaption>**Figure 5.** Most stable structures of C₂₀H₁₀-(H₂O)ₙ (n = 1–8) obtained after PTMD/DFTB and local DFTB optimization following the procedure detailed in ref [96].</figcaption></figure>

model protonated water clusters was first reported by Goyal et al. [285]. Korchagina et al. [175] showed that the cluster (H₂O)₂₁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 Pyₙ⁺ clusters [289], photodissociation of Py₂⁺ [290], combined photoionisation and dissociation of Py₂ [291] and the determination of Pyₙ ionisation potentials (see Figure 6) [292].

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].
<figure><img src="image_44.png" /></figure>

<figure><img src="image_45.png" /><figcaption>**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).</figcaption></figure>

### 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 potential 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 = paramercaptobenzoic 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_3SCH_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 thiolates, DFTB
<figure><img src="image_46.png" /></figure>

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 adsor-bates, 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 under-going 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<sub>55</sub>, Ag<sub>147</sub>, Ag<sub>309</sub>, Ag<sub>561</sub> and Ag<sub>923</sub>. 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<sub>55</sub>, Al<sub>147</sub>, Al<sub>309</sub>, Al<sub>561</sub> and Al<sub>923</sub>. Then, focusing on Al<sub>561</sub>, MD simulation were performed to describe the structure of Al<sub>561</sub> 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.

DFTB has also been applied to model the behavior of dyes grafted on inorganic particles, mainly TiO<sub>2</sub>, 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<sub>2</sub> 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<sub>2</sub> 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
<figure><img src="image_47.png" /></figure>

wavelength. As a representative example, Figure 7 shows how the electronic structure of a naphthalenediol-TiO₂ 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 an 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].

<figure><img src="image_48.png" /></figure>

**Figure 7.** (a) Schematic representation of the atomic structure of a naphthalenediol-TiO₂ complex. Superimposed are the corresponding HOMO (red) and LUMO (blue). (b) Time-dependent population of the HOMO and higher-energy orbitals for naphthalenediol-TiO₂ subject to a continuous laser-type perturbation. Naphthalenediol-TiO₂ 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.
<figure><img src="image_49.png" /></figure>

### 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₂ or NH₃ 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₂ 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₂ (001)) or on an amorphous one ((a)-TiO₂) [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 π 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₂O, CH₄, NH₃) on defective carbon nanotubes which were all found to physisorb on the nanotubes, except NH₃ 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₂ dissociation on plutonium [317].

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₂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
<figure><img src="image_50.png" /></figure>

<figure><img src="image_51.png" /><figcaption>**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.</figcaption></figure>

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
<figure><img src="image_52.png" /></figure>

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⁻¹ 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].

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₃ (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₂ 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² [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
<figure><img src="image_53.png" /></figure>

<figure><img src="image_54.png" /><figcaption>**Figure 9.** Structure of a water dimer interacting with coronene within an argon rare gas matrix subpiece treated via a DFTB-MM scheme [335].</figcaption></figure>

the argon matrix is illustrated in Figure 9. Fine effects such as the modification of the energetic order of the $(H_2O)_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
<figure><img src="image_55.png" /></figure>

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 ($γ_{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 $γ_{CH}$ modes, that are resonant at two different energies at low temperature and thus induce a deviation from linearity of the function $ν_{γ_{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 $γ_{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 (~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
<figure><img src="image_56.png" /></figure>

[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₂₀H₁₂]⁺ 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.

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

<figure><img src="image_57.png" /><figcaption>**Figure 10.** Snapshots retrieved from MD/DFTB simulations describing the evolution of cationic perylene [C₂₀H₁₂]⁺ at high energy (~24–26 eV of internal energy): the formation of a fulvene-type isomer was observed, as well as losses of H, H₂ and C₂H₂, the expected statistical dissociation pathways for PAH radical cations [346].</figcaption></figure>
<figure><img src="image_58.png" /></figure>

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 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

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
<figure><img src="image_59.png" /></figure>

<figure><img src="image_60.png" /><figcaption>**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.</figcaption></figure>

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
<figure><img src="image_61.png" /></figure>

practically a mean field theory since the exact exchange-correlation func-
tional 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 weak-
nesses, 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, photochem-
istry 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 con-
tinuum 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).

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
<figure><img src="image_62.png" /></figure>

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.
<figure><img src="image_63.png" /><figcaption>**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.</figcaption></figure>

**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.

**ORCID**

Fernand Spiegelman http://orcid.org/0000-0001-9412-2866

Nathalie Tarrat http://orcid.org/0000-0002-0909-0001

Jérôme Cuny http://orcid.org/0000-0002-7882-9156

Evgeny Posenitskiy http://orcid.org/0000-0002-1623-0594

Carles Martí http://orcid.org/0000-0003-3341-320X

Aude Simon http://orcid.org/0000-0002-2315-9490

Mathias Rapacioli http://orcid.org/0000-0003-2394-6694
<figure><img src="image_64.png" /></figure>

References

[1] Hohenberg P, Kohn W. Inhomogeneous electron gas. Phys Rev. 1964;136:B864-B871.

[2] Parr RG, Yang W. Density-functional theory of atoms and molecules. Oxford, NY: Oxford University Press; 1989.

[3] Gross EKU, Ullrich CA, Gossman U. Density functional theory. Gross EKU, Dreizler RM, ed. New York: Plenum Press; 1995.

[4] Savin A. Recent developments and applications of modern density functional theory. Seminario J. Amsterdam: Elsevier; 1996. p. 327-357.

[5] Kohn W, Sham LJ. Self-consistent equations including exchange and correlation effects. Phys Rev. 1965;140:A1133-A1138.

[6] Becke AD. Perspective: fifty years of density-functional theory in chemical physics. J Chem Phys. 2014;140:18A301.

[7] Mardirossian N, Head-Gordon M. Thirty years of density functional theory in computational chemistry: an overview and extensive assessment of 200 density functionals. Mol Phys. 2017;115:2315-2372.

[8] Neugebauer J, Hickel T. Density functional theory in materials science. Wiley Interdiscip Rev Comput Mol Sci. 2013;3:438-448.

[9] Bally T, Sastry GN. Incorrect dissociation behavior of radical ions in density func-tional calculations. J Phys Chem A. 1997;101:7923-7925.

[10] Gräfenstein J, Kraka E, Cremer D. The impact of the self-interaction error on the density functional theory description of dissociating radical cations: ionic and cova-lent dissociation limits. J Chem Phys. 2004;120:524-539.

[11] Gräfenstein J, Kraka E, Cremer D. Effect of the self-interaction error for three-electron bonds: on the development of new exchange-correlation functionals. Phys Chem Chem Phys. 2004;6:1096-1112.

[12] Gräfenstein J, Cremer D. The self-interaction error and the description of non-dynamic electron correlation in density functional theory. Theory Chem Acc. 2009;123:171-182.

[13] Bao JL, Gagliardi L, Truhlar DG. Self-interaction error in density functional theory: an appraisal. J Phys Chem Lett. 2018;9:2353-2358.

[14] Leininger T, Stoll H, Werner H-J, et al. Combining long-range configuration inter-action with short-range density functionals. Chem Phys Lett. 1997;275:151-160.

[15] Goll E, Werner H-J, Stoll H. A short-range gradient-corrected density functional in long-range coupled-cluster calculations for rare gas dimers. Phys Chem Chem Phys. 2005;7:3917-3923.

[16] Goll E, Werner HJ, Stoll H. Short-range density functionals in combination with local long-range ab initio methods: application to non-bonded complexes. Chem Phys. 2008;346:257.

[17] Toulouse J, Colonna F, Savin A. Long-range-short-range separation of the electron-electron interaction in density-functional theory. Phys Rev A. 2004;70:062505.

[18] Porezag D, Frauenheim T, Köhler T, et al. Construction of tight-binding-like poten-tials on the basis of density-functional theory: application to carbon. Phys Rev B. 1995;51:12947-12957.

[19] Seifert G, Porezag D, Frauenheim T. Calculations of molecules, clusters, and solids with a simplified LCAO-DFT-LDA scheme. Int J Quantum Chem. 1996;58:185-192.

[20] Elstner M, Porezag D, Jungnickel G, et al. Self-consistent-charge density-functional tight-binding method for simulations of complex materials properties. Phys Rev B. 1998;58:7260-7268.
<figure><img src="image_65.png" /></figure>

[21] Frauenheim T, Seifert G, Elstner M, et al. Atomistic simulations of complex materials: ground-state and excited-state properties. J Phys Cond Mat. 2002;14:3015.

[22] Elstner M, Seifert G. Density functional tight binding. Philos Trans R Soc A Math Phys Eng Sci. 2014;372:20120483.

[23] Seifert G, Joswig J-O. Density-functional tight binding-an approximate density-functional theory method. Wiley Interdiscip Rev Comput Mol Sci. 2012;2:456-465.

[24] Cui Q, Elstner M. Density functional tight binding: values of semi-empirical methods in an ab initio era. Phys Chem Chem Phys. 2014;16:14368-14377.

[25] Koskinen P, Makinen V. Density-functional tight-binding for beginners. Comput Mater Sci. 2009;47:237-253.

[26] Hückel EH. Quantentheoretische beiträge zum benzolproblem. Z Phys. 1931;70:204-286.

[27] Wolfsberg M, Helmholz L. The spectra and electronic structure of the tetrahedral ions $MnO_4^-$, $CrO_4^{2-}$, and $ClO_4^-$. J Chem Phys. 1952;20:837-843.

[28] Hoffmann R. An extended Hückel theory. I. Hydrocarbons. J Chem Phys. 1963;39:1397-16.

[29] Hoffmann R. Extended Hückel theory. III. Compounds of boron and nitrogen. J Chem Phys. 1964;40:2474-2477.

[30] Hoffmann R. Extended hückel theory-v. Tetrahedron. 1966;22:521-538.

[31] Chadi DJ. (110) surface atomic structures of covalent and ionic semiconductors. Phys Rev B. 1979;19:2074-2082.

[32] Harrison WA. In: Electronic structure and the properties of solids. Freeman WH, editor. Dover Editions; 1980.

[33] Ziman JM. The physics of metals. London, NY: Cambridge University Press; 1969.

[34] Desjonquères MC, Spanjaard D. Concepts in surface physics, Springer series in surface sciences. Berlin: Springer Verlag; 1993.

[35] Noguera C. Clean oxide surfaces: a theoretical review. Chem Phys Sol Surf. 2001;9:35-86.

[36] Pople JA, Santry DP, Segal GA. Approximate self-consistent molecular orbital theory. I. Invariant procedures. J Chem Phys. 1965;43:S129-S135.

[37] Dewar MJS, Thiel W. Ground states of molecules. 38. The MNDO method. Approximations and parameters. J Am Chem Soc. 1977;99:4899-4907.

[38] Dewar MJS, Zoebisch EG, Healy EF, et al. Development and use of quantum mechanical molecular models. 76. AM1: a new general purpose quantum mechanical molecular model. J Am Chem Soc. 1985;107:3902-3909.

[39] Stewart JJP. Optimization of parameters for semiempirical methods I. Method. J Comp Chem. 1989;10:209-220.

[40] Goringe CM, Bowler DR, Hernández E. Tight-binding modelling of materials. Rep Prog Phys. 1997;60:1447.

[41] Slater JC, Koster GF. Simplified LCAO method for the periodic potential problem. Phys Rev. 1954;94:1498-1524.

[42] Majewski JA, Vogl P. Simple model for structural properties and crystal stability of sp-bonded solids. Phys Rev B. 1987;35:9666-9682.

[43] Kohyama M, Yamamoto R, Ebata Y, et al. Atomic forces in the self-consistent tight-binding model. Phys Status Solidi B. 1989;152:533-541.

[44] Liu F. Self-consistent tight-binding method. Phys Rev B. 1995;52:10677-10680.

[45] Sankey OF, Niklewski DJ. Ab initio multicenter tight-binding model for molecular-dynamics simulations and other applications in covalent systems. Phys Rev B. 1989;40:3979-3995.
<figure><img src="image_66.png" /></figure>

[46] Tang MS, Wang CZ, Chan CT, et al. Environment-dependent tight-binding potential model. Phys Rev B. 1996;53:979-982.

[47] Horsfield AP, Bratkovsky AM. Ab initio tight binding. J Phys Condens Mater. 2000;12:R1.

[48] Khakimov ZM, Tereshchuk PL, Sulaymanov NT, et al. Nonconventional tight-binding method for the calculation of the total energy and spectroscopic energies of atomic clusters: transferable parameters for silicon. Phys Rev B. 2005;72:115335.

[49] Leahy C, Yu M, Jayanthi CS, et al. Coherent treatment of the self-consistency and the environment-dependency in a semi-empirical Hamiltonian: applications to bulk silicon, silicon surfaces, and silicon clusters. Phys Rev B. 2006;74:155408-155413.

[50] Areshkin DA, Shenderova OA, Schall JD, et al. A self-consistent tight binding model for hydrocarbon systems: application to quantum transport simulation. J Phys Cond Mat. 2004;16:6851-6866.

[51] Zlobin A, Mokrushina Y, Terekhov S, et al. QM/MM description of newly selected catalytic bioscavengers against organophosphorus compounds revealed reactivation stimulus mediated by histidine residue in the acyl-binding loop. Front Pharmacol. 2018;9:834.

[52] Mulliken RS. Electronic population analysis on LCAO-MO molecular wave functions. I. J Chem Phys. 1955;23:1833-1840.

[53] Kaminski S, Giese TJ, Gaus M, et al. Extended polarization in third-order SCC-DFTB from chemical-potential equalization. J Phys Chem A. 2012;116:9131-9141.

[54] Srinivasan SG, Goldman N, Tamblyn I, et al. A density functional tight binding model with an extended basis set and three-body repulsion for hydrogen under extreme thermodynamic conditions. J Phys Chem A. 2014;118:5520-5528.

[55] Yang Y, Yu H, York D, et al. Extension of the self-consistent-charge density-functional tight-binding method: third-order expansion of the density functional theory total energy and introduction of a modified effective coulomb interaction. J Phys Chem A. 2007;111:10861-10873.

[56] Wahiduzzaman M, Oliveira AF, Philipsen P, et al. DFTB parameters for the periodic table: part 1, electronic structure. J Chem Theor Comp. 2013;9:4006-4017.

[57] Oliveira AF, Philipsen P, Heine T. DFTB parameters for the periodic table, part 2: energies and energy gradients from hydrogen to calcium. J Chem Theor Comp. 2015;11:5209-5218.

[58] Gaus M, Chou C-P, Witek H, et al. Automatized parametrization of SCC-DFTB repulsive potentials: application to hydrocarbons. J Phys Chem A. 2009;113:11866-11881.

[59] Chou C-P, Nishimura Y, Fan C-C, et al. Automatized parameterization of DFTB using particle swarm optimization. J Chem Theor Comp. 2016;12:53-64.

[60] Bodrog Z, Aradi B, Frauenheim T. Automated repulsive parametrization for the DFTB method. J Chem Theor Comp. 2011;7:2654-2664.

[61] Lourenço MP, da Silva MC, Oliveira AF, et al. FASP: a framework for automation of Slater-koster file parameterization. Theor Chem Acc. 2016;135:250.

[62] Doemer M, Liberatore E, Knaup JM, et al. In situ parameterisation of SCC-DFTB repulsive potentials by iterative Boltzmann inversion. Mol Phys. 2013;111:3595-3607.

[63] Van den Bossche M. DFTB-assisted global structure optimization of 13- and 55-atom late transition metal clusters. J Phys Chem A. 2019;123:3038-3045.

[64] Kranz JJ, Kubillus M, Ramakrishnan R, et al. Generalized density-functional tight-binding repulsive potentials from unsupervised machine learning. J Chem Theor Comp. 2018;14:2341-2352.
<figure><img src="image_67.png" /></figure>

[65] Frenzel J, Oliveira AF, Jardillier N, et al. Semi-relativistic, self-consistent charge Slater-Koster tables for density-functional based tight-binding (DFTB) for materials science simulations. 2004–2009

[66] Gaus M, Goez A, Elstner M. Parametrization and benchmark of DFTB3 for organic molecules. J Chem Theor Comp. 2013;9:338–354.

[67] Köhler C, Hajnal Z, Deák P, et al. Theoretical investigation of carbon defects and diffusion in α-quartz. Phys Rev B. 2001;64:085333.

[68] Giese TJ, York DM. Density-functional expansion methods: grand challenges. Theor Chem Acc. 2012;131:1145.

[69] Christensen AS, Elstner M, Cui Q. Improving intermolecular interactions in DFTB3 using extended polarization from chemical-potential equalization. J Chem Phys. 2015;143:084123.

[70] Iftner C, Simon A, Korchagina K, et al. A density functional tight binding/force field approach to the interaction of molecules with rare gas clusters: application to $(C_6H_6)^{t/0}Ar_n$ clusters. J Chem Phys. 2014;140:034301.

[71] Zhao Y, Truhlar D. The M06 suite of density functionals for main group thermochemistry, thermochemical kinetics, noncovalent interactions, excited states, and transition elements: two new functionals and systematic testing of four M06-class functionals and 12 other functionals. Theor Chem Acc. 2008;120:215–241.

[72] Dion M, Rydberg H, Schröder E, et al. Van der Waals density functional for general geometries. Phys Rev Lett. 2004;92:246401–246404.

[73] Sato T, Tsuneda T, Hirao K. Long-range corrected density functional study on weakly bound systems: balanced descriptions of various types of molecular interactions. J Chem Phys. 2007;126:234114.

[74] Langreth DC, Dion M, Rydberg H, et al. Van der Waals density functional theory with applications. Int J Quant Chem. 2005;101:599–610.

[75] Chakarova-Kack SD, Schroder E, Lundqvist BI, et al. Application of van der Waals density functional to an extended system: adsorption of benzene and naphthalene on graphite. Phys Rev Lett. 2006;96:146107.

[76] Thonhauser T, Cooper VR, Li S, et al. Van der Waals density functional: self-consistent potential and the nature of the van der Waals bond. Phys Rev B. 2007;76:125112.

[77] von Lilienfeld OA, Tavernelli I, Rothlisberger U, et al. Optimization of effective atom centered potentials for London dispersion forces in density functional theory. Phys Rev Lett. 2004;93:153004.

[78] Gräfenstein J, Cremer D. An efficient algorithm for the density-functional theory treatment of dispersion interactions. J Chem Phys. 2009;130:124105–124116.

[79] Lewis JP, Sankey OF. Geometry and energetics of DNA basepairs and triplets from first principles quantum molecular relaxations. Biophys J. 1995;69:1068–1076.

[80] Meijer EJ, Sprik M. A density-functional study of the intermolecular interactions of benzene. J Chem Phys. 1996;105:8684–8689.

[81] Gianturco FA, Paesani F, Laranjeira MF, et al. Intermolecular forces from density functional theory. III. A multiproerty analysis for the Ar(1S)-CO(1Σ) interaction. J Chem Phys. 1999;110:7832–7845.

[82] Wu Q, Yang W. Empirical correction to density functional theory for van der Waals interactions. J Chem Phys. 2002;116:515–524.

[83] Zimmerli U, Parrinello M, Koumoutsakos P. Dispersion corrections to density functionals for water aromatic interactions. J Chem Phys. 2004;120:2693–2699.
<figure><img src="image_68.png" /></figure>

[84] Goursot A, Mineva T, Kevorkyants R, et al. Interaction between n-alkane chains: applicability of the empirically corrected density functional theory for van der Waals complexes. J Chem Theor Comput. 2007;3:755-763.

[85] Grimme S. Accurate description of van der Waals complexes by density functional theory including empirical corrections. J Comput Chem. 2004;25:1463-1473.

[86] Elstner M. The SCC-DFTB method and its application to biological systems. Theory Chem Acc. 2006;116:316-325.

[87] Zhechkov L, Heine T, Patchovskii S, et al. An efficient a posteriori treatment for dispersion interaction in density-functional-based tight binding. J Chem Theor Comput. 2005;1:841-847.

[88] Rapacioli M, Spiegelman F, Talbi D, et al. Correction for dispersion and Coulombic interactions in molecular clusters with density functional derived methods: application to polycyclic aromatic hydrocarbon clusters. J Chem Phys. 2009;130:244304-244310.

[89] Grimme S. Semiempirical GGA-type density functional constructed with a long-range dispersion correction. J Comput Chem. 2006;27:1787-1799.

[90] Grimme S, Antony J, Ehrlich S, et al. A consistent and accurate ab initio parametrization of density functional dispersion correction (DFT-D) for the 94 elements H-Pu. J Chem Phys. 2010;132:154104.

[91] Bader RF. Atoms in molecules: a quantum theory. Oxford: Oxford University Press; 1990.

[92] Lowdin P-O. On the non-orthogonality problem connected with the use of atomic wave functions in the theory of molecules and crystals. J Chem Phys. 1950;18:365-375.

[93] Li J, Zhu T, Cramer C, et al. New class IV charge model for extracting accurate partial charges from wave functions. J Phys Chem A. 1998;102:1820-1831.

[94] Mayer I. Charge, bond order and valence in the AB initio SCF theory. Chem Phys Lett. 1983;97:270-274.

[95] Kalinowski J, Lesyng B, Thompson J, et al. Class IV charge model for the self-consistent charge density-functional tight-binding method. J Phys Chem A. 2004;108:2545-2549.

[96] Michoulier E, Ben Amor N, Rapacioli M, et al. Theoretical determination of adsorption and ionisation energies of polycyclic aromatic hydrocarbons on water ice. Phys Chem Chem Phys. 2018;20:11941-11953.

[97] Wu Y, Ilie A, Crampin S. Self-consistent charge and dipole density functional tight binding method and application to carbon-based systems. Comput Matter Sci. 2017;134:206-213.

[98] Köhler C, Seifert G, Gerstmann U, et al. Approximate density-functional calculations of spin densities in large molecular systems and complex solids. Phys Chem Chem Phys. 2001;3:5109-5114.

[99] Köhler C, Frauenheim T, Hourahine B, et al. Treatment of collinear and noncollinear electron spin within an approximate density functional based method. J Phys Chem A. 2007;111:5622-5629.

[100] Janak JF. Proof that $δe/δn = ε$ in density-functional theory. Phys Rev B. 1978;18:7165-7168.

[101] Melix P, Oliveira A, Rüger R, et al. Spin polarization in SCC-DFTB. Theory Chem Acc. 2016;135:232.

[102] Mori-Sanchez P, Cohen AJ, Yang W. Many-electron self-interaction error in approximate density functionals. J Chem Phys. 2006;125:201102.
<figure><img src="image_69.png" /></figure>

[103] Perdew JP, Zunger A. Self-interaction correction to density-functional approximations for many-electron systems. Phys Rev B. 1981;23:5048-5079.

[104] Aryasetiawan F, Gunnarsson O. The GW method. Rep Prog Phys. 1998;61:237-312.

[105] Moreira IDPR, Illas F, Martin RL. Effect of Fock exchange on the electronic structure and magnetic coupling in NiO. Phys Rev B. 2002;65:155102.

[106] Anisimov VI, Aryasetiawan F, Lichtenstein AI. First-principles calculations of the electronic structure and spectra of strongly correlated systems: the LDA+ U method. J Phys. 1997;9:767-808.

[107] Vogel D, Krüger P, Pollmann J. Self-interaction and relaxation-corrected pseudopotentials for II-VI semiconductors. Phys Rev B. 1996;54:5495-5511.

[108] Vogel D, Krüger P, Pollmann J. Ab initio electronic structure of silver halides calculated with self-interaction and relaxation-corrected pseudopotentials. Phys Rev B. 1998;58:3865-3869.

[109] Filippetti A, Spaldin NA. Self-interaction-corrected pseudopotential scheme for magnetic and strongly-correlated systems. Phys Rev B. 2003;67:125109.

[110] Hourahine B, Sanna S, Aradi B, et al. Self-interaction and strong correlation in DFTB. J Phys Chem A. 2007;111:5671-5677.

[111] Van Voorhis T, Kowalczyk T, Kaduk B, et al. The diabatic picture of electron transfer, reaction barriers, and molecular dynamics. Annu Rev Phys Chem. 2010;61:149-170.

[112] de la Lande A, Salahub DR. Derivation of interpretative models for long range electron transfer from constrained density functional theory. J Mol Struct. 2010;943:115-120.

[113] Kubař T, Woiczikowski PB, Cuniberti G, et al. Efficient calculation of charge-transfer matrix elements for hole transfer in DNA. J Phys Chem B. 2008;112:7937-7947.

[114] Kubař T, Elstner M. What governs the charge transfer in DNA? The role of DNA conformation and environment. J Phys Chem B. 2008;112:8788-8798.

[115] Kubař T, Kleinekathöfer U, Elstner M. Solvent fluctuations drive the hole transfer in DNA: a mixed quantum-classical study. J Phys Chem B. 2009;113:13107-13117.

[116] Woiczikowski PB, Kubař T, Gutiérrez R, et al. Combined density functional theory and Landauer approach for hole transfer in DNA along classical molecular dynamics trajectories. J Chem Phys. 2009;130:215104-215114.

[117] Endres RG, Cox DL, Singh RRP. Colloquium: the quest for high-conductance DNA. Rev Mod Phys. 2004;76:195-214.

[118] Senthilkumar K, Grozema FC, Guerra CF, et al. Absolute rates of hole transfer in DNA. J Am Chem Soc. 2005;127:14894-14903.

[119] Kubař T, Elstner M. Coarse-grained time-dependent density functional simulation of charge transfer in complex systems: application to hole transfer in DNA. J Phys Chem B. 2010;114:11221-11240.

[120] Wu Q, Van Voorhis T. Direct optimization method to study constrained systems within density-functional theory. Phys Rev A. 2005;72:024502-4.

[121] Wu Q, Van Voorhis T. Direct calculation of electron transfer parameters through constrained density functional theory. J Phys Chem A. 2006;110:9212-9218.

[122] Wu Q, Van Voorhis T. Constrained density functional theory and its application in long-range electron transfer. J Chem Theor Comput. 2006;2:765-774.

[123] Wu Q, Van Voorhis T. Extracting electron transfer coupling elements from constrained density functional theory. J Chem Phys. 2006;125:164105-164109.

[124] Wu Q, Cheng C-L, Van Voorhis T. Configuration interaction based on constrained density functional theory: a multireference method. J Chem Phys. 2007;127:164119.

[125] Rapacioli M, Spiegelman F. Modelling singly ionized coronene clusters. Eur Phys J D. 2009;52:55-58.
<figure><img src="image_70.png" /></figure>

[126] Rapacioli M, Spiegelman F, Scemama A, et al. Modeling charge resonance in cationic molecular clusters: combining DFT-tight binding with configuration interaction. J Chem Theor Comput. 2011;7:44-55.

[127] Lutsker V, Aradi B, Niehaus TA. Implementation and benchmark of a long-range corrected functional in the density functional based tight-binding method. J Chem Phys. 2015;143:184107.

[128] Cui Q, Elstner M, Kaxiras E, et al. A QM/MM implementation of the self-consistent charge density functional tight binding (SCC-DFTB) method. J Phys Chem B. 2001;105:569-585.

[129] Seabra GDM, Walker RC, Elstner M, et al. Implementation of the SCC-DFTB method for hybrid QM/MM simulations within the Amber molecular dynamics package. J Phys Chem A. 2007;111:5655-5664.

[130] Kubar T, Welke K, Groenhof G. New QM/MM implementation of the DFTB3 method in the gromacs package. J Comput Chem. 2015;36:1978-1989.

[131] Simon A, Noble JA, Rouaut G, et al. Formation of coronene:water complexes: FTIR study in argon matrices and theoretical characterisation. Phys Chem Chem Phys. 2017;19:8516-8529.

[132] Lundberg M, Sasakura Y, Zheng G, et al. Case studies of ONIOM(DFT:DFTB) and ONIOM(DFT:DFTB:MM) for enzymes and enzyme mimics. J Chem Theor Comp. 2010;6:1413-1427.

[133] Nishimoto Y. DFTB/PCM applied to ground and excited state potential energy surfaces. J Phys Chem A. 2016;120:771-784.

[134] Niehaus TA, Suhai S, Della Sala F, et al. Tight-binding approach to time-dependent density-functional response theory. Phys Rev B. 2001;63:085108.

[135] Casida ME. In: Delano P Chong, editor. Time-dependent density functional response theory for molecules in recent advances in density functional methods. World Scientific; 1995. p. 155-192.

[136] Casida ME. Time-dependent density-functional theory for molecules and molecular solids. J Mol Struct. 2009;914:3-18.

[137] Rüger R, Niehaus T, van Lenthe E, et al. Vibrationally resolved UV/Vis spectroscopy with time-dependent density functional based tight binding. J Chem Phys. 2016;145:184102.

[138] Trani F, Scalmani G, Zheng G, et al. Time-dependent density functional tight binding: new formulation and benchmark of excited states. J Chem Theor Comp. 2011;7:3304-3313.

[139] Domnguez A, Aradi B, Frauenheim T, et al. Extensions of the time-dependent density functional based tight-binding approach. J Chem Theor Comp. 2013;9:4901-4914.

[140] Humeniuk A, Mitrić R. Long-range correction for tight-binding TD-DFT. J Chem Phys. 2015;143:134120.

[141] Kranz JJ, Elstner M, Aradi B, et al. Time-dependent extension of the long-range corrected density functional based tight-binding method. J Chem Theor Comput. 2017;13:1737-1747.

[142] Chen R-X, Aquino AJA, Sue AC-H, et al. Characterization of charge transfer in excited states of extended clusters of π-stacked donor and acceptor complexes in lock-arm supramolecular ordering. J Phys Chem A. 2019;123:4532-4542.

[143] Freitag J, Domínguez A, Niehaus TA, et al. Nitrogen (II) oxide charge transfer complexes on TiO₂: a new source for visible-light activity. J Phys Chem C. 2015;119:4488-4501.

[144] Gao X, Bai S, Fazzi D, et al. Evaluation of spin-orbit couplings with linear-response time-dependent density functional methods. J Chem Theor Comp. 2017;13:515-524.
<figure><img src="image_71.png" /></figure>

[145] Rüger R, van Lenthe E, Lu Y, et al. Efficient calculation of electronic absorption spectra by means of intensity-selected time-dependent density functional tight binding. J Chem Theor Comp. 2015;11:157-167.

[146] Niehaus T. Approximate time-dependent density functional theory. J Mol Struct Theory Chem. 2009;914:38-49.

[147] Plötz P-A, Niehaus T, Kühn O. A new efficient method for calculation of Frenkel exciton parameters in molecular aggregates. J Chem Phys. 2014;140:174101.

[148] Plötz P-A, Megow J, Niehaus T, et al. Spectral densities for Frenkel exciton dynamics in molecular crystals: a TD-DFTB approach. J Chem Phys. 2017;146:084112.

[149] Dontot L, Suaud N, Rapacioli M, et al. An extended DFTB-CI model for charge-transfer excited states in cationic molecular clusters: model studies versus ab initio calculations in small PAH clusters. Phys Chem Chem Phys. 2016;18:3545-3557.

[150] Dontot L, Spiegelman F, Rapacioli M. Structures and energetics of neutral and cationic pyrene clusters. J Phys Chem A. 2019;123:9531-9543.

[151] Wales DJ, Miller MA, Walsh TR. Archetypal energy landscapes. Nature. 1998;394:758-760.

[152] Johnston RL. Evolving better nanoparticles: genetic algorithms for optimising cluster geometries. Dalton Trans. 2003;22:4193-4207.

[153] Bobadova-Parvanova P, Jackson KA, Srinivas S, et al. Scanning the potential energy surface of iron clusters: a novel search strategy. J Chem Phys. 2002;116:3576-3587.

[154] Köhler C, Seifert G, Frauenheim T. Density functional based calculations for Fe<sub>n</sub> (n<32). Chem Phys. 2005;309:23-31.

[155] Shi H, Koskinen P, Ramasubramaniam A. Self-consistent charge density-functional tight-binding parametrization for Pt-Ru alloys. J Phys Chem A. 2017;121:2497-2502.

[156] Wu L, Zhang L, Qi Y. Structures and electronic properties of a Si<sub>55</sub> cluster on DFTB calculations. Mater Sci Forum. 2015;815:49-53.

[157] Dong Y, Springborg M. Unbiased determination of structural and electronic properties of gold clusters with up to 58 atoms. J Phys Chem C. 2007;111:12528-12535.

[158] Selli D, Fazio G, Di Valentin C. Modelling realistic TiO<sub>2</sub> nanospheres: a benchmark study of SCC-DFTB against hybrid DFT. J Chem Phys. 2017;147:164701.

[159] Li Z, Scheraga HA. Monte Carlo-minimization approach to the multiple-minima problem in protein folding. Proc Natl Acad Sci. 1987;84:6611-6615.

[160] Wales DJ, Doye JPK. Global optimization by basin-hopping and the lowest energy structures of Lennard-Jones clusters containing up to 110 atoms. J Phys Chem A. 1997;101:5111-5116.

[161] Aktürk A, Sebetci A. BH-DFTB/DFT calculations for iron clusters. AIP Adv. 2016;6:055103.

[162] Choi TH, Liang R, Maupin CM, et al. Application of the SCC-DFTB method to hydroxide water clusters and aqueous hydroxide solutions. J Phys Chem B. 2013;117:5165-5179.

[163] Choi TH. Simulation of the (H<sub>2</sub>O)<sub>8</sub> cluster with the SCC-DFTB electronic structure method. Chem Phys Lett. 2012;543:45-49.

[164] Yen TW, Lai SK. Use of density functional theory method to calculate structures of neutral carbon clusters C<sub>n</sub> (3≤n≤24) and study their variability of structural forms. J Chem Phys. 2015;142:084313.

[165] Yen T-W, Lim T-L, Yoon T-L, et al. Studying the varied shapes of gold clusters by an elegant optimization algorithm that hybridizes the density functional tight-binding theory and the density functional theory. Comput Phys Commun. 2017;220:143-149.
<figure><img src="image_72.png" /></figure>

[166] Chen X, Zhao Y-F, Wang L-S, et al. Recent progresses of global minimum searches of nanoclusters with a constrained Basin-Hopping algorithm in the TGMin program. Comput Theor Chem. 2017;1107:57-65.

[167] Zhao Y, Mautner M, Gonzalez C. Ionic hydrogen-bond networks and ion solvation. 1. An efficient Monte Carlo/Quantum mechanical method for structural search and energy computations: ammonium/water. J Phys Chem A. 2009;113:2967-2974.

[168] Earl DJ, Deem MW. Parallel tempering: theory, applications, and new perspectives. Phys Chem Chem Phys. 2005;7:3910-3916.

[169] Swendsen RH, Wang J-S. Replica Monte Carlo simulation of spin-glasses. Phys Rev Lett. 1986;57:2607-2609.

[170] Sugita Y, Okamoto Y. Replica-exchange molecular dynamics method for protein folding. Chem Phys Lett. 1999;314:141-151.

[171] Sugita Y, Okamoto Y. Replica-exchange multicanonical algorithm and multicanonical replica-exchange method for simulating systems with rough energy landscape. Chem Phys Lett. 2000;329:261-270.

[172] Oliveira LFL, Tarrat N, Cuny J, et al. Benchmarking density functional based tight-binding for silver and gold materials: from small clusters to bulk. J Phys Chem A. 2016;120:8469-8483.

[173] Tarrat N, Rapacioli M, Cuny J, et al. Global optimization of neutral and charged 20- and 55-atom silver and gold clusters at the DFTB level. Comput Theor Chem. 2017;1107:102-114.

[174] Korchagina KA, Simon A, Rapacioli M, et al. Structural characterization of sulfur-containing water clusters using a density-functional based tight-binding approach. J Phys Chem A. 2016;120:9089-9100.

[175] Korchagina K, Simon A, Rapacioli M, et al. Theoretical investigation of the solid-liquid phase transition in protonated water clusters. Phys Chem Chem Phys. 2017;19:27288-27298.

[176] Rapacioli M, Barthel R, Heine T, et al. Car-Parrinello treatment for an approximate density-functional theory method. J Chem Phys. 2007;126:124103.

[177] Rapacioli M, Simon A, Marshall CCM, et al. Cationic methylene-pyrene isomers and isomerization pathways: finite temperature theoretical studies. J Phys Chem A. 2015;119:12845-12854.

[178] Cuny J, Korchagina K, Menakbi C, et al. Metadynamics combined with auxiliary density functional and density functional tight-binding methods: alanine dipeptide as a case study. J Mol Model. 2017;23:72.

[179] Mitchell I, Aradi B, Page AJ. Density functional tight binding-based free energy simulations in the DFTB+ program. J Comp Chem. 2018;39:2452-2458.

[180] Labastie P, Whetten RL. Statistical thermodynamics of the cluster solid-liquid transition. Phys Rev Lett. 1990;65:1567-1570.

[181] Oliveira LFL, Cuny J, Moriniere M, et al. Phase changes of the water hexamer and octamer in the gas phase and adsorbed on polycyclic aromatic hydrocarbons. Phys Chem Chem Phys. 2015;17:17079-17089.

[182] McQuarrie DA. In: Harper, Row, editor. Statistical mechanics, Harper and Row. New York: University Science Books; 1976.

[183] Candian A, Mackie CJ. Anharmonic interstellar PAH molecules. Int J Quantum Chem. 2017;117:146-150.

[184] Mulas G, Falvo C, Cassam-Chennai P, et al. Anharmonic vibrational spectroscopy of polycyclic aromatic hydrocarbons (PAHs). J Chem Phys. 2018;149:144102.
<figure><img src="image_73.png" /></figure>

[185] Gamboa A, Rapacioli M, Spiegelman F. Automatic differentiation of the energy within self-consistent tight-binding methods. J Chem Theory Comput. 2013;9:3900-3907.

[186] Witek HA, Irle S, Morokuma K. Analytical second-order geometrical derivatives of energy for the self-consistent-charge density-functional tight-binding method. J Chem Phys. 2004;121:5163-5170.

[187] Furche F, Ahlrichs R. Adiabatic time-dependent density functional methods for excited state properties. J Chem Phys. 2002;117:7433-7447.

[188] Furche F, Ahlrichs R. Erratum: "Time-dependent density functional methods for excited state properties" [J. Chem. Phys. 117, 7433 (2002)]. J Chem Phys. 2004;121:12772-12773.

[189] Heringer D, Niehaus TA, Wanko M, et al. Analytical excited state forces for the time-dependent density-functional tight-binding method. J Comp Chem. 2007;28:2589-2601.

[190] Heringer D, Niehaus TA, Wanko M, et al. Analytical excited state forces for the time-dependent density-functional tight-binding method [J. Comp. Chem. 28, 2589]. J Comp Chem. 2011;33:593-593.

[191] Niehaus T, Heringer D, Torralva B, et al. Importance of electronic self-consistency in the TDDFT based treatment of nonadiabatic molecular dynamics. Eur Phys J D. 2005;35:467-477.

[192] Jakowski J, Morokuma K. Liouville-von Neumann molecular dynamics. J Chem Phys. 2009;130:224106.

[193] Wang Y, Yam C-Y, Frauenheim T, et al. An efficient method for quantum transport simulations in the time domain. Chem Phys. 2011;391:69-77.

[194] Oppenländer C, Korff B, Frauenheim T, et al. Atomistic modeling of dynamical quantum transport. Phys Stat Solid B. 2013;250:2349-2354.

[195] Tully JC. Molecular dynamics with electronic transitions. J Chem Phys. 1990;93:1061-1071.

[196] Hammes-Schiffer S, Tully JC. Proton transfer in solution: molecular dynamics with quantum transitions. J Chem Phys. 1994;101:4657-4667.

[197] Mitrić R, Werner U, Wohlgemuth M, et al. Nonadiabatic dynamics within time-dependent density functional tight binding method. J Phys Chem A. 2009;113:12700-12705.

[198] Tapavicza E, Tavernelli I, Rothlisberger U. Trajectory surface hopping within linear response time-dependent density-functional theory. Phys Rev Lett. 2007;98:023001.

[199] Tapavicza E, Tavernelli I, Rothlisberger U, et al. Mixed time-dependent density-functional theory/classical trajectory surface hopping study of oxirane photochemistry. J Chem Phys. 2008;129:124108.

[200] Humeniuk A, Mitrić R. DFTBaby: a software package for non-adiabatic molecular dynamics simulations based on long-range corrected tight-binding TD-DFT(B). Comput Phys Commun. 2017;221:174-202.

[201] Stojanović L, Aziz SG, Hilal RH, et al. Nonadiabatic dynamics of cycloparaphenylenes with TD-DFTB surface hopping. J Chem Theor Comput. 2017;13:5846-5860.

[202] Pal S, Trivedi DJ, Akimov AV, et al. Nonadiabatic molecular dynamics for thousand atom systems: a tight-binding approach toward PYXAID. J Chem Theor Comput. 2016;12:1436-1448.

[203] Posenitskiy E, Rapacioli M, Lepetit B, et al. Non-adiabatic molecular dynamics investigation of the size dependence of the electronic relaxation in polyacenes. Phys Chem Chem Phys. 2019;21:12139-12149.
<figure><img src="image_74.png" /></figure>

[204] Krüger T, Elstner M, Schiffels P, et al. Validation of the density-functional based tight-binding approximation method for the calculation of reaction energies and other data. J Chem Phys. 2005;122:114110.

[205] Gruden M, Andjeklović L, Jissy AK, et al. Benchmarking density functional tight binding models for barrier heights and reaction energetics of organic molecules. J Comput Chem. 2017;38:2171-2185.

[206] Perriot R, Negre CFA, McGrane SD, et al. Density functional tight binding calculations for the simulation of shocked nitromethane with LATTE-LAMMPS. AIP Conf Proc. 2018;1979:050014.

[207] Lee KH, Schnupf U, Sumpter BG, et al. Performance of density-functional tight-binding in comparison to Ab initio and first-principles methods for isomer geometries and energies of glucose epimers in vacuo and solution. ACS Omega. 2018;3:16899-16915.

[208] Vuong VQ, Akkarapattiakal Kuriappan J, Kubillus M, et al. Parametrization and benchmark of long-range corrected DFTB2 for organic molecules. J Chem Theory Comput. 2018;14:115-125.

[209] Gaus M, Cui Q, Elstner M. Density functional tight binding: application to organic and biological molecules. Wiley Interdiscip Rev Comput Mol Sci. 2014;4:49-61.

[210] Curtiss LA, Raghavachari K, Redfern PC, et al. Assessment of Gaussian-2 and density functional theories for the computation of enthalpies of formation. J Chem Phys. 1997;106:1063-1079.

[211] Baboul AG, Curtiss LA, Redfern PC, et al. Gaussian-3 theory using density functional geometries and zero-point energies. J Chem Phys. 1999;110:7650-7657.

[212] Sattelmeyer KW, Tubert-Brohman I, Jorgensen WL. NO-MNDO: reintroduction of the overlap matrix into MNDO. J Chem Theory Comput. 2006;2:413-419.

[213] Goerigk L, Grimme S. A general database for main group thermochemistry, kinetics, and noncovalent interactions- assessment of common and reparameterized (meta-) GGA density functionals. J Chem Theory Comput. 2010;6:107-126.

[214] Goerigk L, Grimme S. Efficient and accurate double-hybrid-meta-GGA density functionals—evaluation with the extended GMTKN30 database for general main group thermochemistry, kinetics, and noncovalent interactions. J Chem Theory Comput. 2011;7:291-309.

[215] Manna D, Kesharwani MK, Sylvetsky N, et al. Conventional and explicitly correlated ab initio benchmark study on water clusters: revision of the BEGDB and WATER27 data sets. J Chem Theory Comput. 2017;13:3136-3152.

[216] Hostas J, Rezac J. Accurate DFT-D3 calculations in a small basis set. J Chem Theory Comput. 2017;13:3575-3585.

[217] Rezac J, Riley KE, Hobza P. S66: a well-balanced database of benchmark interaction energies relevant to biomolecular structures. J Chem Theory Comput. 2011;7:2427-2438.

[218] Gaus M, Jin H, Demapan D, et al. DFTB3 parametrization for copper: the importance of orbital angular momentum dependence of hubbard parameters. J Chem Theor Comp. 2015;11:4205-4219.

[219] Lu X, Gaus M, Elstner M, et al. Parametrization of DFTB3/3OB for magnesium and zinc for chemical and biological applications. J Phys Chem B. 2015;119:1062-1082.

[220] Zheng G, Witek HA, Bobadova-Parvanova P, et al. Parameter calibration of transition-metal elements for the spin-polarized self-consistent-charge density-functional tight-binding (DFTB) method: Sc, Ti, Fe, Co, and Ni. J Chem Theor Comput. 2007;3:1349-1367.
<figure><img src="image_75.png" /></figure>

[221] Simon A, Rapacioli M, Lanza M, et al. Molecular dynamics simulations on [FePAH] +π-complexes of astrophysical interest: anharmonic infrared spectroscopy. Phys Chem Chem Phys. 2011;13:3359-3374.

[222] Kubar T, Bodrog Z, Gaus M, et al. Parametrization of the SCC-DFTB method for halogens. J Chem Theor Comp. 2013;9:2939-2949.

[223] NIST Standard Reference Database. 2019

[224] Kubillus M, Kubar T, Gaus M, et al. Parameterization of the DFTB3 method for Br, Ca, Cl, F, I, K, and Na in organic and biological systems. J Chem Theor Comp. 2015;11:332-342.

[225] Rezac J, Riley KE, Hobza P. Benchmark calculations of noncovalent interactions of halogenated molecules. J Chem Theory Comput. 2012;8:4285-4292.

[226] Darghouth A, Casida M, Taouali W, et al. Assessment of density-functional tight-binding ionization potentials and electron affinities of molecules of interest for organic solar cells against first-principles GW calculations. Computation. 2015;3:616656.

[227] Zheng G, Irle S, Morokuma K. Performance of the DFTB method in comparison to DFT and semiempirical methods for geometries and energies of C₂₀–C₈₆ fullerene isomers. Chem Phys Lett. 2005;412:210-216.

[228] Witek HA, Irle S, Zheng G, et al. Modeling carbon nanostructures with the self-consistent charge density-functional tight-binding method: vibrational spectra and electronic structure of C₂₈, C₆₀, and C₇₀. J Chem Phys. 2006;125:214706.

[229] Malolepsza E, Witek HA, Irle S. Comparison of geometric, electronic, and vibrational properties for isomers of small fullerenes C₂₀–C₃₆. J Phys Chem A. 2007;111:6649-6657.

[230] Malolepsza E, Lee Y-P, Witek HA, et al. Comparison of geometric, electronic, and vibrational properties for all pentagon/hexagon-bearing isomers of fullerenes C₃₈, C₄₀, and C₄₂. Int J Quant Chem. 2009;109:1999-2011.

[231] Tokoyama H, Yamakado H, Maeda S, et al. Exploration of isomers of benzene by GRRM/SCC-DFTB. Chem Lett. 2014;43:702-704.

[232] Dubosq C, Falvo C, Calvo F, et al. Mapping the structural diversity of C60 carbon clusters and their infrared spectra. Astron Astrophys. 2019;625:L11.

[233] Elstner M, Frauenheim T, Kaxiras E, et al. A self-consistent charge density-functional based tight-binding scheme for large biomolecules. Phys Stat Solidi (B). 2000;217:357-376.

[234] Eynollahi S, Riahi S, Ganjali MR, et al. Density functional theory studies on the geometry and electronic properties of Mitomycin C, DNA Base Pairs and their complex. Int J Electrochem Sci. 2010;5:12.

[235] Riahi S, Eynollahi S, Ganjali MR. Computational studies on effects of MDMA as an anticancer drug on DNA. Chem Biol Drug Des. 2010;76:425-432.

[236] Riahi S, Eynollahi S, Ganjali MR. Interaction of emodin with DNA bases: a density functional theory. J Theory Comput Chem. 2010;09:875-888.

[237] Chehelamirani M, da Silva M, Salahub D. Ultrasonically treated liquid interfaces for progress in cleaning and separation processes. Phys Chem Chem Phys. 2016;18:21-46.

[238] Gutiérrez R, Caetano R, Woiczikowski PB, et al. Structural fluctuations and quantum transport through DNA molecular wires: a combined molecular dynamics and model Hamiltonian approach. New J Phys. 2010;12:023022.

[239] Huang M, Giese TJ, Lee T-S, et al. Improvement of DNA and RNA sugar pucker profiles from semiempirical quantum methods. J Chem Theor Comput. 2014;10:1538-1545.
<figure><img src="image_76.png" /></figure>

[240] Sattelle BM, Almond A. Less is more when simulating unsulfated glycosaminoglycan 3D-structure: comparison of GLYCAM06/TIP3P, PM3-CARB1/TIP3P, and SCC-DFTB-D/TIP3P predictions with experiment. J Comput Chem. 2010;31:2932-2947.

[241] Elstner M, Frauenheim T, Suhai S. An approximate DFT method for QM/MM simulations of biological structures and processes. J Mol Struct. 2003;632:29-41.

[242] Phatak P, Ghosh N, Yu H, et al. Amino acids with an intermolecular proton bond as proton storage site in bacteriorhodopsin. Proc Natl Acad Sci. 2008;105:19672-19677.

[243] Guo H-B, Guo H. Mechanism of histone methylation catalyzed by protein lysine methyltransferase SET7/9 and origin of product specificity. Proc Natl Acad Sci. 2007;104:8797-8802.

[244] Xu D, Guo H. Quantum mechanical/molecular mechanical and density functional theory studies of a prototypical zinc peptidase (Carboxypeptidase A) suggest a general acid-general base mechanism. J Am Chem Soc. 2009;131:9780-9788.

[245] Lodola A, Branduardi D, De Vivo M, et al. A catalytic mechanism for cysteine N-terminal nucleophile hydrolases, as revealed by free energy simulations. PLoS One. 2012;7:e32397.

[246] Lian P, Guo H-B, Smith JC, et al. Catalytic mechanism and origin of high activity of cellulase TmCel12A at high temperature: a quantum mechanical/molecular mechanical study. Cellulose. 2013;21:937-949.

[247] Rogers IL, Naidoo KJ. Producing DFT/MM enzyme reaction trajectories from SCC-DFTB/MM driving forces to probe the underlying electronics of a glycosyltransferase reaction. J Comp Chem. 2017;38:1789-1798.

[248] Silva JRA, Govender T, Maguire GEM, et al. Simulating the inhibition reaction of Mycobacterium tuberculosis L, D-transpeptidase 2 by carbapenems. Chem Commun. 2015;51:12560-12562.

[249] Gillet N, Elstner M, Kubar T. Coupled-perturbed DFTB-QM/MM metadynamics: application to proton-coupled electron transfer. J Chem Phys. 2018;149:072328.

[250] Yao J, Guo H, Chaiprasongsuk M, et al. Substrate-assisted catalysis in the reaction catalyzed by salicylic acid binding protein 2 (SABP2), a potential mechanism of substrate discrimination for some promiscuous enzymes. Biochemistry. 2015;54:5366-5375.

[251] Zhekova HR, Ngo V, Chagas da Silva M, et al. Selective ion binding and transport by membrane proteins - a computational perspective. Coord Chem Rev. 2017;345:108-136.

[252] Wolf S, Freier E, Gerwert K. A delocalized proton-binding site within a membrane protein. Biophys J. 2014;107:174-184.

[253] Riccardi D, Schaefer P, Cui Q. pKacalculations in solution and proteins with QM/MM free energy perturbation simulations: a quantitative test of QM/MM protocols. J Phys Chem B. 2005;109:17715-17733.

[254] Guo H, Wlodawer A, Guo H. A general acid-base mechanism for the stabilization of a tetrahedral adduct in a serine-carboxyl peptidase: a computational study. J Am Chem Soc. 2005;127:15662-15663.

[255] Xiao K, Wang X, Yu H. Nutritional and lipidomics biomarkers of docosahexaenoic acid-based multivitamin therapy in pediatric NASH. Sci Rep. 2019;9:2045-2322.

[256] Wang X, Li R, Cui W, et al. QM/MM free energy simulations of an efficient gluten hydrolase (Kuma030) implicate for a reactant-state based protein-design strategy for general acid/base catalysis. Sci Rep. 2018;8:7042.

[257] Kulakova AM, Khrenova MG, Nemukhin AV. Simulation of spectra of red fluorescent protein mutants. Moscow Univ Chem Bull. 2018;73:212-215.
<figure><img src="image_77.png" /></figure>

[258] Bende A, Ioan T. Nitrogen substituted phenothiazine derivatives: modelling of molecular self-assembling. Int J Mol Sci. 2011;12:3102-3116.

[259] Bogaerts A, Khosravian N, der Paal JV, et al. Multi-level molecular modelling for plasma medicine. J Phys D Appl Phys. 2015;49:054002.

[260] Nishimoto Y, Fedorov DG. The fragment molecular orbital method combined with density-functional tight-binding and the polarizable continuum model. Phys Chem Chem Phys. 2016;18:22047-22061.

[261] Vuong VQ, Nishimoto Y, Fedorov DG, et al. The fragment molecular orbital method based on long-range corrected density-functional tight-binding. J Chem Theor Comput. 2019;15:3008-3020.

[262] Tevekeliyska V, Dong Y, Springborg M, et al. Structural and energetic properties of sodium clusters. Eur Phys J D. 2007;43:19-22.

[263] Frenzel J. Structural, electronic and optical properties of cadmium sulfide nanoparticles. [phd]. Dresden: Technische Universität Dresden; 2006.

[264] Joswig J-O, Springborg M, Seifert G. Structural and electronic properties of cadmium sulfide clusters. J Phys Chem B. 2000;104:2617-2622.

[265] Popov IS, Enyashin AN, Rempel AA. Size dependent content of structural vacancies within TiO nanoparticles: quantum-chemical DFTB study. Superlattices Microstruct. 2018;113:459-465.

[266] Fuertes VC, Negre CFA, Oviedo MB, et al. A theoretical study of the optical properties of nanostructured TiO₂. J Phys Condens Matter. 2013;25:115304.

[267] Cuny J, Tarrat N, Spiegelman F, et al. Stochastic models of multi-channel particulate transport with blockage. J Phys Condens Matter. 2018;30:303001.

[268] Koskinen P, Häkkinen H, Seifert G, et al. Density-functional based tight-binding study of small gold clusters. New J Phys. 2006;8:9.

[269] Fihey A, Hettich C, Touzeau J, et al. SCC-DFTB parameters for simulating hybrid gold-thiolates compounds. J Comput Chem. 2015;36:2075-2087.

[270] Van den Bossche M, Grönbeck H, Hammer B. Tight-binding approximation-enhanced global optimization. J Chem Theor Comput. 2018;14:2797-2807.

[271] Koskinen P, Korhonen T. Plenty of motion at the bottom: atomically thin liquid gold membrane. Nanoscale. 2015;7:10140-10145.

[272] Tarrat N, Rapacioli M, Spiegelman F. Au₁₄₇ nanoparticles: Ordered or amorphous? J Chem Phys. 2018;148:204308.

[273] Minh N, Tuoc V, Lien L. Density functional based tight binding study on wurtzite ZnO prismatic nanoparticles. Commun Phys. 2011;21:235.

[274] Enyashin A, Gemming S, Seifert G. Nanosized allotropes of molybdenum disulfide. Eur Phys J Spec Top. 2007;149:103-125.

[275] Köhler C, Frauenheim T. Magnetic moment of iron clusters with 109, 110, 111, and 147 atoms. J Comp Theor Nanosci. 2007;4:264-269.

[276] Yuan T, Larsson K. Theoretical study of size effects on surface chemical properties for nanoscale diamond particles. J Phys Chem C. 2014;118:26061-26069.

[277] Li W, Irle S, Witek HA. Convergence in the evolution of nanodiamond raman spectra with particle size: a theoretical investigation. ACS Nano. 2010;4:4475-4486.

[278] Rincon L, Hasmy A, Marquez M, et al. A perturbatively corrected tight-binding method with hybridization: application to gold nanoparticles. Chem Phys Lett. 2011;503:171-175.

[279] Schebarchov D, Baletto F, Wales DJ. Structure, thermodynamics, and rearrangement mechanisms in gold clusters—insights from the energy landscapes framework. Nanoscale. 2018;10:2004-2016.
<figure><img src="image_78.png" /></figure>

[280] Sun B, Fernandez M, Barnard AS. Machine learning for silver nanoparticle electron transfer property prediction. J Chem Inf Model. 2017;57:2413-2423.

[281] Simon A, Spiegelman F. Water clusters adsorbed on polycyclic aromatic hydrocarbons: energetics and conformational dynamics. J Chem Phys. 2013;138:194309.

[282] Simon A, Spiegelman F. Conformational dynamics and finite-temperature infrared spectra of the water octamer adsorbed on coronene. Comput Theor Chem. 2013;1021:54-61.

[283] Simon A, Rapacioli M, Mascetti J, et al. Vibrational spectroscopy and molecular dynamics of water monomers and dimers adsorbed on polycyclic aromatic hydrocarbons. Phys Chem Chem Phys. 2012;14:6771-6786.

[284] Perez C, Steber AL, Rijs AM, et al. Corannulene and its complex with water: a tiny cup of water. Phys Chem Chem Phys. 2017;19:14214-14223.

[285] Goyal P, Elstner M, Cui Q. Application of the SCC-DFTB method to neutral and protonated water clusters and bulk water. J Phys Chem B. 2011;115:6790-6805.

[286] Wu -C-C, Lin C-K, Chang H-C, et al. Protonated clathrate cages enclosing neutral water molecules: H⁺(H₂O)₂₁ and H⁺(H₂O)₂₈. J Chem Phys. 2005;122:074315.

[287] Iyengar SS, Petersen MK, Day TJF, et al. The properties of ion-water clusters. I. The protonated 21-water cluster. J Chem Phys. 2005;123:084309.

[288] Yu H, Cui Q. The vibrational spectra of protonated water clusters: a benchmark for self-consistent-charge density-functional tight binding. J Chem Phys. 2007;127:234504.

[289] Zamith S, Ji M-C, LHermite J-M, et al. Thermal evaporation of pyrene clusters. J Chem Phys. 2019;151:194303.

[290] Bernard J, Al-Mogeeth A, Dontot L, et al. in prep.

[291] Garcia G, Nahon L, Dontot L, et al. in prep.

[292] Joblin C, Dontot L, Garcia GA, et al. Size effect in the ionization energy of PAH clusters. J Phys Chem Lett. 2017;8:3697-3702.

[293] Addicoat MA, Fukuoka S, Page AJ, et al. Stochastic structure determination for conformationally flexible heterogenous molecular clusters: application to ionic liquids. J Comp Chem. 2013;34:2591-2600.

[294] Addicoat MA, Stefanovic R, Webber GB, et al. Assessment of the density functional tight binding method for protic ionic liquids. J Chem Theor Comput. 2014;10:4633-4643.

[295] Kullgren J, Wolf MJ, Hermansson K, et al. Self-consistent-charge density-functional tight-binding (SCC-DFTB) parameters for ceria in 0D to 3D. J Phys Chem C. 2017;121:4593-4607.

[296] Moreira NH, Dolgonos G, Aradi B, et al. Toward an accurate density-functional tight-binding description of zinc-containing compounds. J Chem Theor Comput. 2009;5:605-614.

[297] Li W, Kotsis K, Manzhos S. Comparative density functional theory and density functional tight binding study of arginine and arginine-rich cell penetrating peptide TAT adsorption on anatase TiO 2. Phys Chem Chem Phys. 2016;18:19902-19917.

[298] Mäkinen V, Koskinen P, Häkkinen H. Modeling thiolate-protected gold clusters with density-functional tight-binding. Eur Phys J D. 2013;67:38.

[299] Dominguez-Castro A, Hernández D, Guzmán F. Insights into the interactions of biomolecules with small gold clusters: a theoretical study from a DFTB perspective. Theory Chem Acc. 2017;136:84.

[300] Douglas-Gallardo OA, Berdakin M, Sánchez CG. Atomistic insights into chemical interface damping of surface plasmon excitations in silver nanoclusters. J Phys Chem C. 2016;120:24389-24399.
<figure><img src="image_79.png" /></figure>

[301] Bonafe FP, Aradi B, Guan M, et al. Plasmon-driven sub-picosecond breathing of metal nanoparticles. Nanoscale. 2017;9:12391-12397.

[302] Kelly KL, Coronado E, Zhao LL, et al. The optical properties of metal nanoparticles: the influence of size, shape, and dielectric environment. J Phys Chem B. 2003;107:668-677.

[303] Eustis S, El-Sayed MA. Why gold nanoparticles are more precious than pretty gold: noble metal surface plasmon resonance and its enhancement of the radiative and nonradiative properties of nanocrystals of different shapes. Chem Soc Rev. 2006;35:209-217.

[304] Millstone JE, Hurst SJ, Métraux GS, et al. Colloidal gold and silver triangular nanoprisms. Small. 2009;5:646-664.

[305] Douglas-Gallardo OA, Soldano GJ, Mariscal MM, et al. Effects of oxidation on the plasmonic properties of aluminum nanoclusters. Nanoscale. 2017;9:17471-17480.

[306] Oviedo MB, Zarate X, Negre CFA, et al. Quantum dynamical simulations as a tool for predicting photoinjection mechanisms in dye-sensitized TiO₂ solar cells. J Phys Chem Lett. 2012;3:2548-2555.

[307] Negre CFA, Fuertes VC, Oviedo MB, et al. Quantum dynamics of light-induced charge injection in a model dye-nanoparticle complex. J Phys Chem C. 2012;116:14748-14753.

[308] Negre CFA, Young KJ, Oviedo MB, et al. Photoelectrochemical hole injection revealed in polyoxotitanate nanocrystals functionalized with organic adsorbates. J Am Chem Soc. 2014;136:16420-16429.

[309] Zhao YL, Köppen S, Frauenheim T. An SCC-DFTB/MD study of the adsorption of zwitterionic glycine on a geminal hydroxylated silica surface in an explicit water environment. J Phys Chem C. 2011;115:9615-9621.

[310] Zhao Y-L, Wang C-H, Zhai Y, et al. Selective adsorption of L-serine functional groups on the anatase TiO₂ (101) surface in benthic microbial fuel cells. Phys Chem Chem Phys. 2014;16:20806-20817.

[311] Manzhos S, Giorgi G, Yamashita K. A density functional tight binding study of acetic acid adsorption on crystalline and amorphous surfaces of titania. Molecules. 2015;20:3371-3388.

[312] Guo L, Qi C, Zheng X, et al. Toward understanding the adsorption mechanism of large size organic corrosion inhibitors on an Fe(110) surface using the DFTB method. RSC Adv. 2017;7:29042-29050.

[313] Kanematsu Y, Gohara K, Yamada H, et al. Applicability of density functional tight binding method with dispersion correction to investigate the adsorption of porphyrin/porphycene metal complexes on graphene. Chem Lett. 2017;46:51-52.

[314] Li D, Wang F, Zhang Z, et al. The nature of small molecules adsorbed on defective carbon nanotubes. R Soc Open Sci. 2019;6:190727.

[315] Fan G, Zhu S, Xu H. Density-functional theory study of the interaction mechanism and optical properties of flavonols on the boron nitride nanotubes. Int J Quantum Chem. 2018;118:e25514.

[316] Dolgonos G, Aradi B, Moreira NH, et al. An improved self-consistent-charge density-functional tight-binding (SCC-DFTB) set of parameters for simulation of bulk and molecular systems involving titanium. J Chem Theory Comput. 2010;6:266-278.

[317] Goldman N, Aradi B, Lindsey RK, et al. Development of a multicenter density functional tight binding model for plutonium surface hydriding. J Chem Theory Comput. 2018;14:2652-2660.
<figure><img src="image_80.png" /></figure>

[318] Bernstein MP, Sandford SA, Mattioda AL, et al. Near-and mid-infrared laboratory spectra of PAH cations in solid H₂O. Astrophys J. 2007;664:1264.

[319] Osberg KI. Photochemistry and astrochemistry: photochemical pathways to interstellar complex organic molecules. Chem Rev. 2016;116:9631-9663.

[320] Gudipati MS, Allamandola LJ. Polycyclic aromatic hydrocarbon ionization energy lowering in water ices. Astrophys J. 2004;615:L177-L180.

[321] Woon DE, Park JY. Photoionization of benzene and small polycyclic aromatic hydrocarbons in ultraviolet-processed astrophysical ices: a computational study. Astrophys J. 2004;607:342-345.

[322] Guennoun Z, Aupetit C, Mascetti J. Photochemistry of coronene with water at 10 K: first tentative identification by infrared spectroscopy of oxygen containing coronene products. Phys Chem Chem Phys. 2011;13:7340-7347.

[323] Guennoun Z, Aupetit C, Mascetti J. Photochemistry of pyrene with water at low temperature: study of atmospherical and astrochemical interest. J Phys Chem A. 2011;115:1844-1852.

[324] Bouwman J, Paardekooper DM, Cuppen HM, et al. Real-time optical spectroscopy of vacuum ultraviolet irradiated pyrene: H₂O interstellar ice. Astrophys J. 2009;700:56.

[325] Michoulier E, Toubin C, Simon A, Mascetti J, Aupetit C, Noble JA. Perturbation of the Surface of Amorphous Solid Water by the Adsorption of Polycyclic Aromatic Hydrocarbons. J Phys Chem C. 2019, in press. doi:10.1021/acs.jpcc.9b09499

[326] Simon A, Rapacioli M, Michoulier E, et al. Contribution of the density-functional-based tight-binding scheme to the description of water clusters: methods, applications and extension to bulk systems. Mol Sim. 2019;45:249-268.

[327] Szücs B, Hajnal Z, Scholz R, et al. Theoretical study of the adsorption of a PTCDA monolayer on S-passivated GaAs(100). Appl Surf Sci. 2004;234:173-177.

[328] Chen Y, Liu M, Chen J, et al. A density functional based tight binding (DFTB+) study on the sulfidization-amine flotation mechanism of smithsonite. Appl Surf Sci. 2018;458:454-463.

[329] Hellström M, Jorner K, Bryngelsson M, et al. An SCC-DFTB repulsive potential for various ZnO polymorphs and the ZnO-water system. J Phys Chem C. 2013;117:17004-17015.

[330] Selli D, Fazio G, Seifert G, et al. Water multilayers on TiO₂ (101) anatase surface: assessment of a DFTB-based method. J Chem Theory Comput. 2017;13:3862-3873.

[331] Prasetyo N, Hofer TS. Adsorption and dissociation of water molecules at the α-Al₂O₃ (0001) surface: a 2-dimensional hybrid self-consistent charge density functional based tight-binding/molecular mechanics molecular dynamics (2D SCC-DFTB/MM MD) simulation study. Comput Matter Sci. 2019;164:195-204.

[332] Li H-B, Page AJ, Hettich C, et al. Graphene nucleation on a surface-molten copper catalyst: quantum chemical molecular dynamics simulations. Chem Sci. 2014;5:3493-3500.

[333] Abdalla S, Springborg M, Dong Y. Isolated and deposited potassium clusters: energetic and structural properties. Surf Sci. 2013;608:255-264.

[334] Stiehler C, Pan Y, Schneider W-D, et al. Electron quantization in arbitrarily shaped gold islands on MgO thin films. Phys Rev B. 2013;88:115415.

[335] Simon A, Iftner C, Mascetti J, et al. Water clusters in an argon matrix: infrared spectra from molecular dynamics simulations with a self-consistent charge density functional-based tight binding/force-field potential. J Phys Chem A. 2015;119:2449-2467.
<figure><img src="image_81.png" /></figure>

[336] Witek H, Morokuma K, Stradomska A. Modeling vibrational spectra using the self-consistent charge density-functional tight-binding method II: infrared spectra. J Theory Comput Chem. 2005;4:639-655.

[337] Witek H, Morokuma K, Stradomska A. Modeling vibrational spectra using the self-consistent charge density-functional tight-binding method. I. Raman spectra. J Chem Phys. 2004;121:5171-5178.

[338] Witek H, Morokuma K. Systematic study of vibrational frequencies calculated with the self-consistent charge density functional tight-binding method. J Comput Chem. 2004;25:1858-1864.

[339] Page AJ, Chou C-P, Pham BQ, et al. Quantum chemical investigation of epoxide and ether groups in graphene oxide and their vibrational spectra. Phys Chem Chem Phys. 2013;15:3725-3735.

[340] Dos Santos HF, De Souza LA, De Almeida WB, et al. Structure, stability, and infrared spectrum of capped carbon cones: a DFTB study. J Phys Chem C. 2014;118:24761-24768.

[341] Kaminski S, Gaus M, Phatak P, et al. Vibrational raman spectra from the self-consistent charge density functional tight binding method via classical time-correlation functions. J Chem Theor Comp. 2010;6:1240-1255.

[342] Joalland B, Rapacioli M, Simon A, et al. Molecular dynamics simulations of anharmonic infrared spectra of $[sipah]⁺$ $[SiPah]π$-complexes. J Phys Chem A. 2010;114:5846-5854.

[343] Pietrucci F. Strategies for the exploration of free energy landscapes: unity in diversity and challenges ahead. Rev Phys. 2017;2:32-45.

[344] Torrie GM, Valleau JP. Nonphysical sampling distributions in Monte Carlo free-energy estimation: umbrella sampling. J Comput Phys. 1977;23:187-199.

[345] Laio A, Parrinello M. Escaping free-energy minima. Proc Natl Acad Sci. 2002;99:12562-12566.

[346] Simon A, Rapacioli M, Rouaut G, et al. Dissociation of polycyclic aromatic hydrocarbons: molecular dynamics studies. Phil Trans R Soc A. 2017;375:20160195.

[347] Rapacioli M, Cazaux S, Foley N, et al. Atomic hydrogen interactions with gas-phase coronene cations: hydrogenation versus fragmentation. Phys Chem Chem Phys. 2018;20:22427-22438.

[348] Simon A, Champeaux JP, Rapacioli M, et al. Dissociation of polycyclic aromatic hydrocarbons at high energy: MD/DFTB simulations versus collision experiments. Theory Chem Acc. 2018;137:106.

[349] Korchagina KA, Spiegelman F, Cuny J. Molecular dynamics study of the collision-induced reaction of H with CO on small water clusters. J Phys Chem A. 2017;121:9485-9494.

[350] Javier Dominguez-Gutierrez F, Krstic PS, Irle S, et al. Low-energy hydrogen uptake by small-cage $C_n$ and $C_{n-1}B$ fullerenes. Carbon. 2018;134:189-198.

[351] Qian H-J, van Duin ACT, Morokuma K, et al. Reactive molecular dynamics simulation of fullerene combustion synthesis: reaxFF vs DFTB potentials. J Chem Theor Comput. 2011;7:2040-2048.

[352] Saha B, Shindo S, Irle S, et al. Quantum chemical molecular dynamics simulations of dynamic fullerene self-assembly in benzene combustion. ACS Nano. 2009;3:2241-2257.

[353] Saha B, Irle S, Morokuma K. Formation mechanism of polycyclic aromatic hydrocarbons in benzene combustion: quantum chemical molecular dynamics simulations. J Chem Phys. 2010;132:224303.
<figure><img src="image_82.png" /></figure>

[354] Ohta Y, Okamoto Y, Irle S, et al. Rapid growth of a single-walled carbon nanotube on an iron cluster: density-functional tight-binding molecular dynamics simulations. ACS Nano. 2008;2:1437-1444.

[355] Zheng G, Irle S, Morokuma K. Fe/C interactions during SWNT growth with C₂ feedstock molecules: a quantum chemical molecular dynamics study. J Nanosci Nanotechnol. 2006;6:1259-1270.

[356] Wang Z, Irle S, Zheng G, et al. Carbon nanotubes grow on the C face of SiC (0001) during sublimation decomposition: quantum chemical molecular dynamics simulations. J Phys Chem C. 2007;111:12960-12972.

[357] Deng Q, Heine T, Irle S, et al. Self-assembly of endohedral metallofullerenes: a decisive role of cooling gas and metal-carbon bonding. Nanoscale. 2016;8:3796-3808.

[358] Adeagbo WA, Entel P. Determination of melting of water clusters using density functional theory. Phase Transitions. 2004;77:63-79.

[359] Rapacioli M, Tarrat N, Spiegelman F. Melting of the Au₂₀ gold cluster: does charge matter? J Phys Chem A. 2018;122:4092-4098.

[360] Rapacioli M, Spiegelman F, Tarrat N. Evidencing the relationship between isomer spectra and melting: the 20- and 55-atom silver and gold cluster cases. Phys Chem Chem Phys. 2019;21:24857-24866.

[361] Niehaus TA, Meziane M, Lepine F, et al. Pulse shape and molecular orientation determine the attosecond charge migration in Caffeine. Eur Phys J B. 2018;91:152.

[362] Titov E, Humeniuk A, Mitrić R. Exciton localization in excited-state dynamics of a tetracene trimer: a surface hopping LC-TDDFTB study. Phys Chem Chem Phys. 2018;20:25995-26007.

[363] Medrano CR, Oviedo MB, Sánchez CG. Photoinduced charge-transfer dynamics simulations in noncovalently bonded molecular aggregates. Phys Chem Chem Phys. 2016;18:14840-14849.

[364] Darghouth AAMHM, Correa GC, Juillard S, et al. Davydov-type excitonic effects on the absorption spectra of parallel-stacked and herringbone aggregates of pentacene: time-dependent density-functional theory and time-dependent density-functional tight binding. J Chem Phys. 2018;149:134111.

[365] Hoche J, Schmitt H-C, Humeniuk A, et al. The mechanism of excimer formation: an experimental and theoretical study on the pyrene dimer. Phys Chem Chem Phys. 2017;19:25002-25015.

[366] Plötz P-A, Megow J, Niehaus T, et al. All-DFTB approach to the parametrization of the system-bath hamiltonian describing exciton-vibrational dynamics of molecular assemblies. J Chem Theory Comput. 2018;14:5001-5010.

[367] Plötz PA, Polyutov SP, Ivanov SD, et al. Biphasic aggregation of a perylene bisimide dye identified by exciton-vibrational spectra. Phys Chem Chem Phys. 2016;18:25110-25119.

[368] Negre CFA, Perassi EM, Coronado EA, et al. Quantum dynamical simulations of local field enhancement in metal nanoparticles. J Phys Condens Matter. 2013;25:125304.

[369] Aradi B, Hourahine B, Frauenheim T. DFTB+, a sparse matrix-based implementation of the DFTB method. J Phys Chem A. 2007;111:5678-5684.

[370] Heine T, Rapacioli M, Patchkovskii S, et al. deMon-nano experiment. 2009. Available from: http://demon-nano.ups-tlse.fr/

[371] Te Velde G, Bickelhaupt FM, Baerends EJ, et al. Chemistry with ADF. J Comput Chem. 2001;22:931-967.

[372] Salomon-Ferrer R, Case DA, Walker RC. An overview of the Amber biomolecular simulation package. Wiley Interdiscip Rev Comput Mol Sci. 2013;3:198-210.
<figure><img src="image_83.png" /></figure>

[373] Berendsen H, van der Spoel D, van Drunen R. GROMACS: a message-passing parallel molecular dynamics implementation. Comput Phys Commun. 1995;91:43-56.

[374] Frisch MJ, Trucks GW, Schlegel HB, et al. Gaussian16 Revision C.01. Wallingford CT: Gaussian Inc.; 2016.

[375] Hutter J, Iannuzzi M, Schiffmann F, et al. Wiley Interdiscip. Rev Comput Mol Sci. 2014;4:15-25.

[376] Scemama A, Renon N, Rapacioli M. A sparse self-consistent field algorithm and its parallel implementation: application to density-functional-based tight binding. J Chem Theory Comput. 2014;10:2344-2354.

[377] Kosugi K, Nakano H, Sato H. SCC-DFTB-PIMD method to evaluate a multidimensional quantum free-energy surface for a proton-transfer reaction. J Chem Theory Comput. 2019;15:4965-4973.

[378] Head-Gordon M, Gordon M, Tully JC. Molecular dynamics with electronic frictions. J Chem Phys. 1995;103:10137-10145.

[379] Giese TJ, Chen H, Dissanayake T, et al. A variational linear-scaling framework to build practical, efficient next-generation orbital-based quantum force fields. J Chem Theory Comput. 2013;9:1417-1427.

[380] Zhang H, Zapol P, Dixon DA, et al. Shift-and-invert parallel spectral transformation eigensolver: massively parallel performance for density-functional based tight-binding. J Comput Chem. 2015;37:448-459.

[381] Nishizawa H, Nishimura Y, Kobayashi M, et al. Three pillars for achieving quantum mechanical molecular dynamics simulations of huge systems: divide-and-conquer, density-functional tight-binding, and massively parallel computation. J Comput Chem. 2016;37:1983-1992.

[382] Nishimura Y, Nakai H. Dcdftbmd: divide-and-conquer density functional tight-binding program for huge-system quantum mechanical molecular dynamics simulations. J Comput Chem. 2019;40:1538-1549.