# Automatic Differentiation of the Energy within Self-consistent Tight-Binding Methods

Antonio Gamboa, Mathias Rapacioli, Fernand Spiegelman

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

► To cite this version:

Antonio Gamboa, Mathias Rapacioli, Fernand Spiegelman. Automatic Differentiation of the Energy within Self-consistent Tight-Binding Methods. Journal of Chemical Theory and Computation, 2013, 9 (9), pp.3900-3907. 10.1021/ct400214b. hal-00868882

HAL Id: hal-00868882

https://hal.science/hal-00868882

Submitted on 3 Jun 2021

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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Automatic Differentiation of the Energy within Self-consistent Tight-Binding Methods

Antonio Gamboa,\* Mathias Rapacioli, and Fernand Spiegelman

Université de Toulouse, UPS, LCPQ (Laboratoire de Chimie et Physique Quantiques), IRSAMC, 118 Route de Narbonne, F-31062 Toulouse, France

CNRS, LCPQ (Laboratoire de Chimie et Physique Quantiques), IRSAMC, F-31062 Toulouse, France

ABSTRACT: We present and implement the calculation of analytical n-order geometric derivatives of the energy obtained within the framework of the density functional based tight binding approach. The use of automatic differentiation techniques allows a unique implementation for the calculation of derivatives up to any order providing that the computational facilities are sufficient. As first applications, the derivatives are used to build an analytical potential energy surface around the optimized geometry of acetylene. We also discuss the relevant anharmonic contributions that have to be considered when building such an analytical potential energy surface for acetylene, ethylene, ethane, benzene, and naphtalene.

# 1. INTRODUCTION

High-order response properties of molecular systems subjected to perturbations have received a growing interest in the past few years (for a recent review, see ref 1). For instance, the response to an electric field allows the description of nonlinear spectroscopy,² through calculation of hyperpolarizabilities constants. Computing high-order energy derivatives with respect to geometric displacements can be used in optimization schemes or to build an analytical local potential energy surface (PES) with a Taylor expansion. Such a surface can then be used with approaches demanding fast energy computation to access quantum nuclear effects, like wave packet dynamics, quantum vibrational Monte Carlo, or calculation of infrared (IR) spectra including quantum effects with many-body mean-field or multiconfigurational methods for quantum vibrational perturbation theory.³⁻⁵

In the past decades, density functional theory (DFT) combined with the development of computational facilities has opened the route to the study of large systems. In many ab initio code using DFT, the third-order derivatives can be obtained as well as some partial fourth-order ones. In other calculations, harmonic constants may be derived with high level methods such as CCSD(T) and corrections with DFT third and fourth partial order derivatives are achieved to provide analytical surfaces.⁶

Tight-binding (TB) methods are much faster approaches that conserve the quantum description with explicit use of molecular orbitals (MO). Developed as the early quantum semiempirical methods (like the extended Hückel Hamiltonians⁷), these approaches, which evolved to more complex formulations including for instance self-consistent processes, can be derived from ab initio methods (TB-HF from Hartree-

Fock⁸ or DFTB from DFT⁹⁻¹²). In all TB approaches, the solution of the electronic problem is achieved by the diagonalization (possibly iteratively) of the Hamiltonian matrix (eventually expressed in a nonorthogonal basis). The computational efficiency of TB approaches relies on the use of a minimal atomic basis set and the fact that the Hamiltonian and overlap matrices elements are obtained from two-body parametrized functions or interpolated from precalculated points.

Among these TB approaches, density functional based tight-binding (DFTB⁹⁻¹²) is a computationally efficient DFT scheme in which all integrals are parametrized from DFT calculations. The second-order DFTB, also called self-consistent charge (SCC¹¹) DFTB, is a refinement of this approach providing a more accurate description of moderately charged systems requiring a self-consistent diagonalization scheme. This approach has been widely used over the past decades to extract properties of systems for which DFT turns out to be too time-consuming in fields ranging from solid-state structures¹³ to biomolecules.¹⁴

Analytical second-order SCC-DFTB energy derivatives with respect to geometry have been formulated by Witek et al.¹⁵ and used to extract physical properties, in particular IR spectra in the harmonic approximation.¹⁶,¹⁷ In this work, we describe how higher-order geometric derivatives of the energy can be computed analytically within the SCC-DFTB framework. The use of automatic differentiation enables a unique implementation that allows the calculation of these derivatives up to any
<figure><img src="image_3.png" /><figcaption>Figure 1. Hamiltonian and overlap matrix elements as functions of the distance for C-H interaction obtained from the analytical fit eq 6 and some of points obtained from the DFT integrals used to fit these parameters.</figcaption></figure>

order, provided that the CPU time and storage capabilities are sufficient.

Although the approach is applied here to SCC-DFTB, it is general for any self-consistent (or not) TB scheme. It can also be applied to achieve high-order derivations with respect to any perturbation (for instance electric or magnetic field) once the TB matrices elements derivatives with respect to this perturbation are known.

The paper is organized as follows. In the next section, after a brief presentation of the DFTB approach, we present the formalism allowing the calculation of any order properties derivatives within this framework as well as an analytic fit of DFTB parameters. We then benchmark the approach by comparing numerical and analytical calculation of derivatives. In the last section, the derivatives are used to build an analytical PES. We discuss the n-order derivative terms, that are quantitatively relevant on some examples.

## 2. METHOD

2.1. SCC-DFTB. In this section, we briefly review the basics of the SCC-DFTB method.⁹⁻¹² It is derived from the Kohn-Sham (KS) DFT through an expansion of the DFT energy functional with respect to the electronic density around a reference density up to second-order. All three center integrals are neglected. With these approximations, the SCC-DFTB is given by

$$E^{\text{SCC-DFTB}} = \sum_{\alpha<\beta}^{\text{atoms}} E_{\alpha\beta}^{\text{rep}} + \sum_{i} n_i \langle \psi_i | \hat{H}^0 | \psi_i \rangle + \frac{1}{2} \sum_{\alpha,\beta}^{\text{atoms}} \gamma_{\alpha\beta} \Delta q_{\alpha} \Delta q_{\beta} \quad (1)$$

where $\hat{H}^0$ is the KS operator at the reference density, $n_i$ is the electronic occupation of the KS MO $\Psi_i$, and $E_{\alpha\beta}^{\text{rep}}$ is a repulsive potential between atoms $\alpha$ and $\beta$. The last term is the second-order contribution to the energy expressed through atomic pair

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

contributions involving a $\gamma$ function of interatomic distances and atomic charges fluctuations $\Delta q_{\alpha}$.

The molecular orbitals are developed on a minimal atomic orbital (AO) basis set {$\phi_{\mu}$}

$$\psi_i = \sum_{\mu} c_{i\mu} \phi_{\mu} \quad (2)$$

leading to the matrix formulation

$$E^{\text{SCC-DFTB}} = \sum_{\alpha<\beta} E_{\alpha\beta}^{\text{rep}} + \sum_{i}^{\text{occ}} \sum_{\mu\nu} n_i c_{i\mu} c_{i\nu} H_{\mu\nu}^{0} + \frac{1}{2} \sum_{\alpha\beta} \gamma_{\alpha\beta} \Delta q_{\alpha} \Delta q_{\beta} \quad (3)$$

where $H^0$ and S are the KS and overlap matrices expressed in the AO basis. Atomic charge fluctuations are calculated with the Mulliken analysis approach $\Delta q_{\alpha} = q_{\alpha} - q_{\alpha}^0$ where $q_{\alpha}^0$ is the number of valence electrons in the neutral free atom and

$$q_{\alpha} = \sum_{i} n_{i} \sum_{\mu \in \alpha} \sum_{\nu} c_{i\mu} c_{i\nu} S_{\mu\nu}$$

As the atomic charge fluctuations depend on the MO coefficients, the energy minimization is performed through self-consistent cycles by solving the eigenvalue equation:

$$\sum_{\nu} c_{i\nu} (H_{\mu\nu} - \epsilon_{i} S_{\mu\nu}) = 0 \quad (4)$$

with

$$H_{\mu \in \alpha, \nu \in \beta} = H_{\mu\nu}^{0} + H_{\mu\nu}^{1} = H_{\mu\nu}^{0} + \frac{1}{2} S_{\mu\nu} \sum_{\xi} (\gamma_{\alpha\xi} + \gamma_{\xi\beta}) \Delta q_{\xi} \quad (5)$$

The matrix elements of $H^0$, S, and $E_{\alpha\beta}^{\text{rep}}$ are interpolated from two-body DFT calculations.

Additional terms can be added to account for London dispersion ($E^{\text{disp}}$) forces as a sum over atomic pairs.¹⁸⁻²⁰ In
recent years, order terms with respect to the density have been introduced.²¹

**2.2. Parameter Fitting.** In the next section, *n*-order response properties resulting from a geometric perturbation ζ (geometrical change, electric field, etc.) will be derived analytically assuming the knowledge of the DFTB parameter derivatives with respect to this perturbation: $\partial H_{\mu\nu}^0 / \partial \zeta$, $\partial S_{\mu\nu} / \partial \zeta$, $\partial E_{\alpha\beta}^{\text{exp}} / \partial \zeta$, $\partial \gamma_{\alpha\beta} / \partial \zeta$.

DFTB matrices elements are computed at the DFT level for selected interatomic distances. Clearly, one cannot use low-order splines for high-order derivations and global analytic fitting of the DFTB parameters is relevant. The analytical fitting has been discussed by several authors either with Chebyshev polynomials²² or standard polynomials.²³ One has to be aware of possible artifacts and transferability problems for very small matrices elements and out of the DFT tabulated ranges at long interatomic distances. Long-range artifacts can be in principle corrected via cutoff functions. The essential scope of the present paper is to present the analytical derivatives of the DFTB energy. This algebra is valid for any kind of derivable fitting functions. In the present work, we have expressed the interatomic elements of the Hamiltonian matrix at the reference density $H^0$ and the overlap matrix S as a polynomial function of the interatomic distance through a polynomial expansion:

$$ \Theta_k(R) = \sum_{n=-5}^{15} a_{nk} R^n \quad (6) $$

The subindex *k* specifies the different types of interacting orbitals. For instance, in the case of second-row elements, the possible orbital interactions are s-s, s-p, p-s, (p-p)σ and (p-p)π. As an example, we show in Figure 1 the results of the fitting process for the Hamiltonian and overlap matrices elements corresponding to s-s and s-p atomic basis orbitals for interaction of C-H atoms. We have used high-order polynomials in order to ensure the fit accuracy better than 9 × 10⁻⁵ Hartrees for all the DFT tabulated points. This is important to compute high-order derivatives. The *n*-derivatives of the $H^0$ and S matrices are obtained by calculating the derivatives of the products of this function with the rotation matrices corresponding to various angular momenta.²⁴ The fact that the DFTB method only takes in account two-body interactions simplifies the calculations because all the derivatives that involve coordinates of three or more different atoms are zero.

**2.3. Analytical Derivatives.** We first consider the calculation of *n*-order derivatives with respect to a single coordinate in the standard second-order SCC-DFTB. We follow here the approach of Masmoudi et al.²⁵ and extend it to the self-consistent scheme. In order to obtain the *n*-derivatives of the MO coefficients $c_{\mu\nu}$ and the eigenvalues $\epsilon_{\mu}$, we start calculating the *n*-derivative of eq 4:

$$ \left[ \sum_{\nu} c_{\mu\nu} (H_{\mu\nu} - \epsilon_{\mu} S_{\mu\nu}) \right]^n = 0 \quad (7) $$

The superscript stands for the *n*-order derivative $\partial^n / \partial\zeta^n$. We suppose now that we have solved eq 7 for all orders $k < n$. So, in this equation, the only unknown terms are the highest-order ones $c_{\mu\nu}^n$ and $\epsilon_{\mu}^n$.

In order to calculate this terms, the derivatives of the coefficients $c_{\mu\nu}^n$ are expressed as a linear combination of the MO coefficients $c_{ij\mu}$:

$$ c_{\mu\nu}^n = \sum_m c_{m\nu} u_{mi} \quad (8) $$

Introducing this expansion in eq 7, using the formula of differentiation of a product

$$ (AB)^n = \sum_{r=0}^{n} \binom{n}{r} A^r B^{n-r} \quad (9) $$

where $\binom{n}{r}$ is the binomial coefficient $n!/(r!(n-r)!)$ and multiplying by $c_{ij\mu}$ coefficients and summing over $\mu$ leads to the following expression:

$$ u_{ji}(\epsilon_i - \epsilon_j) = \sum_{\mu\nu} \sum_{r=0}^{n-1} \binom{n}{r} c_{j\mu} c_{i\nu}^r (H_{\mu\nu} - \epsilon_i S_{\mu\nu})^{(n-r)} \quad (10) $$

If we extract the term corresponding to $r=0$ from the sum and we apply again eq 9 to expand the product of $(\epsilon_i S_{\mu\nu})^{(n-r)}$, we can separate eq 10 in two cases. For $i \neq j$, we have

$$ u_{ji}(\epsilon_i - \epsilon_j) = \sum_{\mu\nu} c_{j\mu} c_{i\nu} \left( H_{\mu\nu}^n - \sum_{k=1}^n \binom{n}{k} \epsilon_i^{n-k} S_{\mu\nu}^k \right) + \sum_{\mu\nu} \sum_{r=1}^{n-1} \binom{n}{r} c_{j\mu} c_{i\nu}^r (H_{\mu\nu} - \epsilon_i S_{\mu\nu})^{(n-r)} \quad (11) $$

and for $i=j$,

$$ \epsilon_i^n = \sum_{\mu\nu} c_{i\mu} c_{i\nu} \left( H_{\mu\nu}^n - \sum_{k=1}^n \binom{n}{k} \epsilon_i^{n-k} S_{\mu\nu}^k \right) + \sum_{\mu\nu} \sum_{r=1}^{n-1} \binom{n}{r} c_{i\mu} c_{i\nu}^r (H_{\mu\nu} - \epsilon_i S_{\mu\nu})^{(n-r)} \quad (12) $$

where we are making use of the fact that $\sum_{\mu\nu} c_{i\mu} c_{j\nu} S_{\mu\nu} = \delta_{ij}$.

These two sets of equations can be solved independently. For the non-SCC case, all the terms in the right-hand side (rhs) of eq 11 are known, as $H_{\mu\nu} = H_{\mu\nu}^0$. Therefore, this equation gives the solution for each $u_{if}$ for $i \neq j$. On the other hand, in the SCC scheme, the Hamiltonian matrix H depends on $\Delta q^n$ and therefore on $u_{if}$. The contribution of the *n*-order derivatives of the Hamiltonian involving $u_{if}$ are identified and transferred to the left-hand side (lhs) of the equation, leading to the linear system:

$$ \sum_{lm} A_{ijlm} u_{lm} = R'(c^{k<n}, \epsilon^{k<n}, H^{0k \le n}, S^{k \le n}) \quad (13) $$

At first order, this system is known as the coupled perturbed equations (see, for instance, ref 26). The full expressions of $A_{ijlm}$ and $R'$ are given in Appendix A. We notice that, although this system has to be solved for each order of the $c_{i\mu}$ derivatives, the matrix $A_{ijlm}$ will allways remain the same and only a single inversion of this matrix is necessary for the calculations of all derivatives.

To get the diagonal terms $u_{ij\mu}$ we use the normalization condition.

$$ \left( \sum_{\mu\nu} c_{i\mu} S_{\mu\nu} c_{i\nu} \right)^{n_1 \dots n_N} = 0 \quad (14) $$

that leads to
$$u_{it} = -\frac{1}{2} \sum_{\mu\nu} \left[ \sum_{r=0}^{n-1} \sum_{p=0}^{r} \binom{n}{r} \binom{r}{p} c_{i\mu}^{p} c_{i\nu}^{r} - p S_{\mu\nu}^{n-r} + \sum_{r=1}^{n-1} \binom{n}{r} c_{i\mu}^{r} c_{i\nu}^{n-r} S_{\mu\nu} \right] \quad (15)$$

Once the MO coefficients are determined, eq 12 can be solved to obtain the eigenvalues derivatives $\epsilon_i^n$.

In the general case, namely when we want to calculate a derivative according to more than a single coordinate, eq 7 is replaced by

$$[\sum_{\nu} c_{i\nu}(H_{\mu\nu} - \epsilon_i S_{\mu\nu})]^{n_1, \dots, n_N} = 0 \quad (16)$$

where the general superscript $n_1, \dots, n_N$ stands for the following derivatives:

$$\frac{\partial^{n_1+\dots+n_N}}{\partial \xi_1^{n_1} \dots \partial \xi_N^{n_N}}$$

Equation 10 is replaced by

$$u_{ji}(\epsilon_i - \epsilon_j) = \sum_{\mu\nu} \sum_{r_1=0}^{n_1-1} \sum_{r_2=0}^{n_2-1} \dots \sum_{r_N=0}^{n_N-1} \binom{n_1}{r_1} \binom{n_N}{r_N} c_{j\mu} c_{i\nu}^{r_1, \dots, r_N} \times (H_{\mu\nu} - \epsilon_i S_{\mu\nu})^{(n_1-r_1, \dots, n_N-r_N)} \quad (17)$$

The procedure for solving the latter is exactly the same as the derivation according to a single variable. Besides, in the SCC case, the matrix $A_{ijlm}$ that needs to be inverted to solve the linear problem equivalent to eq 13, is exactly the same as the one for the case of derivation with respect to a single variable. This property is very favorable to save computational efforts.

Another fact that simplifies the calculation of the derivatives is the Wigner's $2n + 1$ rule,²⁷ according to which, the knowledge of the *n*-derivative of the coefficients is sufficient to calculate the $(2n+1)$-order derivative for the energy.

The applications in the next section section are carried on using second-order SCC-DFTB. Nevertheless, for sake of completeness, we show that analytical derivatives can also be carried on in the case of the third-order expansion of DFTB with respect to the density (DFTB3 approach²¹). The algebra of the derivatives is given in Appendix B.

### 3. BENCHMARKS

The algorithm presented in the previous section has been implemented in a preliminary version of the deMonNano²⁸ code, a version of the deMon suite of programs devoted to large systems calculations with the DFTB approach. To benchmark the algorithm and its implementation, *n*-order derivatives have been calculated with the analytical approach presented previously and through finite differences. Both numerical and analytical derivatives calculations have been performed with the analytical parameter fitting eq 6. The choice of the geometric displacement step, taken here as 0.01 bohr, results from a compromise between the fact that this value should be as small as possible to respect the differentiation definition but not too small to avoid numerical errors. We are calculating this derivatives in two steps. In the first one, we make an evaluation of the energy in a grid given by all the displacements: $(\pm\delta x, \pm2\delta x, \dots, \pm n\delta x)$, where *n* is the maximum order we want to calculate. The energy is then stored for all

points of this grid. In the second step, we use this stored values to calculate the derivatives. This way of calculating numerical derivatives makes it fast, but it has the drawback of requiring a very big amount of memory as we increase the size of the system.

As can be seen in Table 1, the relative error remains smaller than 0.1% up to third-order derivatives and increases up to 1–

<table><caption>**Table 1. Maximum Relative Deviation between Numerical and Analytical *n*-Order (*n* = 2, 3, 4) Derivatives for Acetylene, Ethylene, Ethane, and Benzene Moleculesª**</caption><thead><tr><th>molecule</th><th>second-order</th><th>third-order</th><th>fourth-order</th></tr></thead><tbody><tr><td>acetylene</td><td>9.8419 &times; 10<sup>-4</sup></td><td>1.0348 &times; 10<sup>-4</sup></td><td>6.8994 &times; 10<sup>-2</sup></td></tr><tr><td>ethylene</td><td>1.0037 &times; 10<sup>-4</sup></td><td>9.3118 &times; 10<sup>-4</sup></td><td>1.1079 &times; 10<sup>-3</sup></td></tr><tr><td>ethane</td><td>1.2842 &times; 10<sup>-4</sup></td><td>6.2790 &times; 10<sup>-4</sup></td><td>3.1483 &times; 10<sup>-2</sup></td></tr><tr><td>benzene</td><td>1.9457 &times; 10<sup>-4</sup></td><td>6.8088 &times; 10<sup>-4</sup></td><td>1.2265 &times; 10<sup>-2</sup></td></tr></tbody></table>

ªAs the relative errors have no meaning for very small derivatives, we only considered derivatives that have an absolute value of at least 1% of the biggest derivative.

7% for fourth-order derivative calculations. We also mentioned that the analytical differentiation is much faster than performing finite difference calculations as can be seen in Figure 2. This is

<figure><img src="image_5.png" /><figcaption>**Figure 2.** Time scale (seconds) comparison between calculation of analytical (dashed lines) and numerical (solid lines) as a function of the basis size (i.e., the dimension of DFTB matrices).</figcaption></figure>

true, even when we are calculating the numerical derivatives with the efficient scheme previously described. Nevertheless, the numerical fourth-order derivatives for naphthalene were not calculated because the memory and time demands were too big. This fact indicates that, if we want to calculate derivatives up to fourth-order for big systems, one would need to do it without storing the values of the energy in a grid, and the time differences between analytical and numerical derivatives would be even larger than the ones showed in Figure 2.

### 4. APPLICATIONS

In this section, we present selected applications of *n*-order derivative calculations. We first present on a simple system (C₂H₂) how going up to fourth-order is mandatory to have a reasonably good description of the PES within a Taylor expansion. Finally, we discuss the distributions of derivative absolute values, depending on the number of involved normal modes for a variety of systems, namely acetylene, ethylene, ethane, benzene, and naphthalene.
<figure><img src="image_6.png" /><figcaption>**Figure 3.** PES of the acetylene molecule as function of C–C and C–H bond lengths obtained from single point calculations (left). In the upper panel, we show, from left to right, the PES obtained from Taylor’s expansion using analytical derivatives up to second-, third-, and fourth-order. In the lower panel, we plot the absolute value of the difference between the full PES and the *n*-order approximation.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 4.** Absolute values of third- and fourth-order derivatives of the energy with respect to frequency weighted modes coordinates, $\tilde{N}_{\alpha}$. The values of the derivatives are normalized with respect to the maximum value of the derivative in each case. Derivatives with respect to 1, 2, 3, and 4 different coordinates are separated by dashed lines.</figcaption></figure>

**4.1. Building an Analytical Potential Energy Surface of Acetylene.** As a first application of high-order derivative calculations, we focus on analytical PES reduced to two possible geometric deformations, namely the C–C and C–H bond length. Figure 3 present the corresponding potential energy subsurface, obtained either from single point DFTB calculations either from Taylor expansion of the PES using *n*-order ($n = 2–4$) derivatives. The first plot shows that, around the minimum energy structure, iso-energetic contours look like a slightly deformed ellipse. Obviously, the harmonic approximation can not go beyond the ellipse shape description of the PES. This is not the case with the third-order expansion, the non-symmetric shape of the PES being recovered. However, two fictitious local metastable minima appear at intermediate distances from the most stable configuration. Such minima could have strong artificial effects for instance if some physical properties would be calculated by molecular dynamics on this approximated surface. These errors are washed by including fourth-order terms in the expansion, the local metastable minima being pushed much further from the most stable configuration.

**4.2. Relevant Terms in *n*-Order Derivatives.** Calculating high-order response properties in a DFT framework can be extremely computationally costly for systems even with a relatively small number of degrees of freedom. A solution often used is to calculate derivatives according to a selected number of modes. For instance, the calculation of fourth-order derivatives is often restricted to the dependence of one or two different normal modes $N_i$. Taking advantage of the computational cost of DFTB and analytical high-order response implementation, it is possible to calculate all derivatives up to fourth-order for several molecular systems and to identify the significant terms. Figure 4 represents, for acetylene, ethylene, ethane, benzene, and naphtalene, the third- and fourth-order energy derivatives with respect to frequency scaled normal mode coordinates: $\tilde{N}_{\alpha} = N_{\alpha}(\omega_{\alpha})^{1/2}$, $N_{\alpha}$ being the normal mode resulting from the diagonalization of the mass weighted Hessian matrix.

For the third-order derivatives of the five studied systems, some derivatives with respect to two or three different normal modes are of the same order of magnitude (or even larger) as the ones with respect to one single normal mode.
For fourth-order derivatives, the situation is different between the three systems. In some systems, as acetylene and ethane, the derivatives according to one or two normal modes are the most important. However in the cases of ethylene, benzene, and naphtalene, derivatives with respect to 4 different coordinates are the same order of magnitude as the previous ones.

Such examples suggests that when only some values of derivatives are used to compute physical properties from a fitted surface, it is more consistent to calculate all derivatives and to use a threshold to eliminate some derivatives than to select it according to the number of involved normal modes. In a DFT framework, calculating all derivatives would have a computationally expensive cost. In such cases, DFTB could be used as a precalculation to select a limited number of derivatives to be calculated at the DFT level.

## 5. CONCLUSION

In this paper, we have derived general equation to calculate the *n*-order response of the energy of a system with respect to a geometric perturbation, calculated within SCC-DFTB approach. The implementation of these equations in the deMonNano code makes use of automatic differentiation which enable the calculation of derivatives up to any order provided the fact that memory and computational time is available. Due to storage limitation, only derivatives up to fourth-order have been calculated in this work.

The implementation has been benchmarked by comparing the results with derivatives computed by the finite differences method. The analytical calculation is much faster than the latter approach. Such derivatives can be used to obtain an approximated PES of the system through a Taylor expansion up to a given order.

We notice that the equations reported in this paper are general and can be used to calculate the response with respect to any perturbation as long as the DFTB parameter derivatives with respect to this perturbation are known. This opens the route for instance to calculation of hyperpolarisabilities constants. The scope of this work was to express analytical *n*-order derivatives of DFTB properties to avoid its numerical calculations, without addressing the quality of DFTB PES. One should however keep in mind that DFTB is an approximated DFT scheme and that the quality of the *n*-order derivatives of any property relies on the DFTB accuracy itself, which can always be improved by refining parameters, introducing higher-order terms<sup>29</sup> or some specific corrections for instance to treat long-range dispersion and Coulombic interactions.<sup>18-20</sup>

As an illustration, we have discussed how including third- and fourth-order terms in the building of an analytical PES based on a Taylor expansion improves the local description of the acetylene PES. The role of high-order terms has been shown to have a significant effect. Finally, we have discussed the relative values of the derivatives for third- and fourth-order energy derivatives with respect to 1, ..., 4 normal mode deformations for acetylene, ethylene, ethane, benzene, and naphtalene. We showed that fourth-order derivatives involving two normal modes are for some systems less important than the same order derivatives involving 3 or 4 normal modes. As calculating derivatives with respect to all coordinates might be prohibitive at the DFT level and as the selection of the significant terms can hardly be done on the number of normal modes involved, we suggest that a fast DFTB high-order derivative calculation

could be performed in order to select the DFT derivatives to be calculated, when this level of description is required.

## APPENDIX A: FULL EXPRESSION OF THE TERMS $A_{ijlm}$ AND $R'$ IN EQUATION 13

In this appendix, we give the steps to obtain eq 13 from eq 11 for the SCC case. The left part of eq 13 contains all the terms that involve $u_{lm}$ in eq 11. The hamiltonian $H_{\mu\nu}$ is the sum of a non-self-consistent part $H_{\mu\nu}^0$ and the SCC component, $H_{\mu\nu}^1$ given by eq 5. The *n*-derivative of $H^1$ is given by

$$ 
\begin{aligned}
[H_{\mu,\nu}^1]^n &= \left[ \frac{1}{2} S_{\mu\nu} \sum_{\xi} (\gamma_{\alpha\xi} + \gamma_{\xi\beta}) \Delta q_{\xi} \right]^n \\
&= \frac{1}{2} \sum_{\xi} \sum_{r=1}^{n} \binom{n}{r} [S_{\mu\nu} \sum_{\xi} (\gamma_{\alpha\xi} + \gamma_{\xi\beta})]^{\prime} \Delta q_{\xi}^{n-r} \\
&\quad + \frac{1}{2} S_{\mu\nu} \sum_{\xi} (\gamma_{\alpha\xi} + \gamma_{\xi\beta}) \Delta q_{\xi}^{n}
\end{aligned} 
\quad (18) $$

and

$$ \Delta q_{\xi}^{n} = g_{\xi}^{n} + \sum_{l}^{\text{occ}} n_{l} \sum_{\omega \in \xi} \sum_{\chi} S_{\omega\chi}(c_{l\omega}c_{l\chi}^{n} + c_{l\chi}c_{l\omega}^{n}) \quad (19) $$

where

$$ g_{\xi}^{n} = \sum_{l}^{\text{occ}} n_{l} \sum_{\omega \in \xi} \sum_{\chi} \left[ \sum_{r=1}^{n} \binom{n}{r} S_{\omega\chi}^{r} (c_{l\omega} c_{l\chi})^{n-r} \right. \\ \left. + S_{\omega\chi} \sum_{r=1}^{n-1} \binom{n}{r} c_{l\omega}^{r} c_{l\chi}^{n-r} \right] \quad (20) $$

The second term in the expression for $\Delta q_{\xi}^n$ is the only one in which appear the *n*-order derivatives of the coefficients $c_{ij\omega}$ while $g_{\xi}$ contains only lower-order terms. The matrix $A_{ijlm}$ is obtained after expressing the derivatives $c_{ij\mu}^n$ in terms of the coefficients $u_{if}$ using eq 8 and grouping all these terms in the rhs of eq 11. The final expression is the following:

$$ 
\begin{aligned}
A_{ijlm} = (\varepsilon_l - \varepsilon_m) \delta_{il} \delta_{jm} & - \frac{1}{2} \sum_{\mu\nu} c_{ij\mu} c_{ij\nu} S_{\mu\nu} \sum_{\xi} (\gamma_{\alpha\xi} + \gamma_{\xi\beta}) \\ & \times n_m \sum_{\omega \in \xi} \sum_{\chi} S_{\omega\chi} (c_{l\chi} c_{m\omega} + c_{l\omega} c_{m\chi})
\end{aligned} 
\quad (21) $$

We notice that the expression of this matrix is equivalent to the one used to obtain the MO coefficient derivatives in the context of the DFTB-CI<sup>30</sup>

The rhs of eq 13 contains all the other terms and can be expressed as

$$ R'(c^{k<n}, \varepsilon^{k<n}, H^{0k \le n}, S^{k \le n}) = R_{ij}^1 + R_{ij}^2 + R_{ij}^3 + R_{ij}^4 \quad (22) $$

where

$$ R_{ij}^{1} = \sum_{\mu\nu} \sum_{r=1}^{n-1} \binom{n}{r} c_{j\mu} c_{i\nu}^{r} (H_{\mu\nu} - \varepsilon_{i} S_{\mu\nu})^{n-r} \quad (23) $$

$$ R_{ij}^{2} = - \sum_{\mu\nu} c_{j\mu} c_{i\nu} \sum_{r=1}^{n-1} \binom{n}{r} \varepsilon_{i}^{n-r} S_{\mu\nu}^{r} \quad (24) $$
$$R_{ij}^3 = \frac{1}{2} \sum_{\mu\nu} c_{\mu\nu} c_{\nu\mu} \sum_{\xi} \sum_{r=1}^{n} \binom{n}{r} [S_{\mu\nu} \sum_{\xi} (\gamma_{\alpha\xi} + \gamma_{\xi\beta})]^r \Delta q_{\xi}^{n-r} \quad (25)$$

and

$$R_{ij}^4 = \frac{1}{2} \sum_{\mu\nu} c_{\mu\nu} c_{\nu\mu} S_{\mu\nu} \sum_{\xi} (\gamma_{\alpha\xi} + \gamma_{\xi\beta}) g_{\xi} \quad (26)$$

## APPENDIX B: FULL EXPRESSION OF THE TERMS $A_{ijLM}$ AND $R'$ IN THE DFTB3 SCHEME

In the third-order extension of the SCC-DFTB method, the Hamiltonian can be expressed as

$$H_{\mu \in \alpha, \nu \in \beta} = H_{\mu\nu}^{2nd} + \sum_{\xi} \Delta q_{\xi} L_{\xi\mu\nu} \quad (27)$$

where $H_{\mu\nu}^{2nd}$ is the second-order Hamiltonian given by equation eq 5, and

$$L_{\xi\mu \in \alpha\nu \in \beta} = \frac{1}{3} (\Delta q_{\alpha} \Gamma_{\alpha\xi} + \Delta q_{\beta} \Gamma_{\beta\xi}) + \frac{\Delta q_{\xi}}{6} (\Gamma_{\xi\alpha} + \Gamma_{\xi\beta}) \quad (28)$$

where

$$\Gamma_{\xi\alpha} = \left. \frac{\partial \gamma_{\xi\alpha}}{\partial q_{\xi}} \right|_{q_{\xi}^0} \quad (29)$$

Applying the same procedure used in the second-order case, we arrive to the correspondent expressions for the left and right parts in this extended scheme. For the matrix $A_{ijlm}$, we have

$$A_{ijlm} = A_{ijlm}^{2nd} - \sum_{\mu\nu} c_{\mu\nu} c_{\nu\mu} S_{\mu\nu} n_m \times \left\{ \sum_{\xi} (L_{\xi\mu\nu} + M_{\xi\mu\nu}) \right. \\ \left. \sum_{\omega \in \xi} \sum_{\chi} S_{\omega\chi} (c_{l\chi} c_{m\omega} + c_{l\omega} c_{m\chi}) + \sum_{\xi} \frac{\Delta q_{\xi}}{3} \Gamma_{\alpha\xi} \sum_{\omega \in \alpha} \sum_{\chi} S_{\omega\chi} (c_{l\chi} c_{m\omega} + c_{l\omega} c_{m\chi}) + \sum_{\xi} \frac{\Delta q_{\xi}}{3} \Gamma_{\beta\xi} \sum_{\omega \in \beta} \sum_{\chi} S_{\omega\chi} (c_{l\chi} c_{m\omega} + c_{l\omega} c_{m\chi}) \right\} \quad (30)$$

where $A_{ijlm}^{2nd}$ is given for eq 21 and

$$M_{\xi\mu\nu} = \frac{\Delta q_{\xi}}{6} (\Gamma_{\xi\alpha} + \Gamma_{\xi\beta}) \quad (31)$$

The expression for the right part $R'$ is given by

$$R'(c^{k<n}, \epsilon^{k<n}, H^{0k\le n}, S^{k\le n}) = R_{ij}^{2nd} + T_{ij}^{1} + T_{ij}^{2} + T_{ij}^{3} + T_{ij}^{4} + T_{ij}^{5} \quad (32)$$

where the second-order part is given by eq 22 and the rest of the terms are

$$T_{ij}^{1} = \sum_{\mu\nu} \sum_{r=1}^{n-1} \binom{n}{r} c_{\mu\nu} c_{\nu\mu}^{r} \left( \sum_{\xi} \Delta q_{\xi} L_{\xi\mu\nu} \right)^{n-r} \quad (33)$$

$$T_{ij}^{2} = \sum_{\mu\nu} c_{\mu\nu} c_{\nu\mu} \sum_{\xi} \sum_{r=1}^{n} \binom{n}{r} S_{\mu\nu}^{r} (\Delta q_{\xi} L_{\xi\mu\nu})^{n-r} \quad (34)$$

$$T_{ij}^{3} = \sum_{\mu\nu} c_{\mu\nu} c_{\nu\mu} \sum_{\xi} \sum_{r=1}^{n-1} \binom{n}{r} S_{\mu\nu} \Delta q_{\xi}^{r} L_{\xi\mu\nu}^{n-r} \quad (35)$$

being

$$T_{ij}^{4} = \sum_{\mu\nu} c_{\mu\nu} c_{\nu\mu} \sum_{\xi} S_{\mu\nu} \Delta q_{\xi} T_{\xi\mu\nu}^{4} \quad (36)$$

where

$$T_{\xi\mu\nu}^{4} = \sum_{r=0}^{n-1} \binom{n}{r} \left( \frac{\Delta q_{\alpha}^{r}}{3} \Gamma_{\alpha\xi}^{n-r} + \frac{\Delta q_{\beta}^{r}}{3} \Gamma_{\beta\xi}^{n-r} + \frac{\Delta q_{\xi}^{r}}{6} (\Gamma_{\xi\alpha} + \Gamma_{\xi\beta})^{n-r} \right) \quad (37)$$

and

$$T_{ij}^{5} = \sum_{\mu\nu} c_{\mu\nu} c_{\nu\mu} S_{\mu\nu} \sum_{\xi} \left\{ g_{\xi}^{n} \left( L_{\xi\mu\nu} + \frac{\Delta q_{\xi}}{6} (\Gamma_{\xi\alpha} + \Gamma_{\xi\beta}) \right) + \frac{\Delta q_{\xi}}{3} (g_{\alpha}^{n} \Gamma_{\alpha\xi} + g_{\beta}^{n} \Gamma_{\beta\xi}) \right\} \quad (38)$$

## ASSOCIATED CONTENT

### Supporting Information

Values of the coefficients $a_{nk}$ in eq 6. This material is available free of charge via the Internet at http://pubs.acs.org/.

## AUTHOR INFORMATION

### Corresponding Author

*E-mail: gamboa@irsamc.ups-tlse.fr.

### Notes

The authors declare no competing financial interest.

## ACKNOWLEDGMENTS

Computational support from the CALMIP/Toulouse facility is acknowledged as well as funding from ANR GASPARIM and CNRS/GDR2758.

## REFERENCES

(1) Bast, R.; Ekstrom, U.; Gao, B.; Helgaker, T.; Ruud, K.; Thorvaldsen, A. J. Phys. Chem. Chem. Phys. 2011, 13, 2627–2651.

(2) Mukamel, S. Nonlinear Optical Spectroscopy; Oxford University Press: Oxford, U.K., 1995.

(3) Barone, V. J. Chem. Phys. 2005, 122, 014108.

(4) Mills, I. M. Mol. Spectrosc. Mod. Res. 1972, 1, 115.

(5) Yagi, K.; Hirata, S.; Hirao, K. Phys. Chem. Chem. Phys. 2008, 10, 1781–1788.

(6) Carbonniere, P.; Pouchan, C. Theo. Chem. Acc. 2012, 131, 1–8.

(7) Hückel, E. H. Z. Phys. 1931, 70, 204–86.

(8) Montagnon, L.; Spiegelman, F. J. Chem. Phys. 2007, 127, 084111–13.

(9) Porezag, D.; Frauenheim, T.; Köhler, T.; Seifert, G.; Kaschner, R. Phys. Rev. B 1995, 51, 12947–12957.

(10) Seifert, G.; Porezag, D.; Frauenheim, T. Int. J. Quantum Chem. 1996, 58, 185–192.

(11) Elstner, M.; Porezag, D.; Jungnickel, G.; Elsner, J.; Haugk, M.; Frauenheim, T.; Suhai, S.; Seifert, G. Phys. Rev. B 1998, 58, 7260–7268.

(12) Oliveira, A.; Seifert, G.; Heine, T.; duarte, H. J. Braz. Chem. Soc. 2009, 20, 1193–1205.
(13) Frauenheim, T.; Seifert, G.; Elstner, M.; Hajnal, Z.; Jungnickel, G.; Porezag, D.; Suhai, S.; Scholz, R. *Phys. Stat. Solidi (b)* **2000**, 217, 41–62.

(14) Elstner, M. *Theo. Chem. Acc.* **2006**, 116, 316–325.

(15) Witek, H. A.; Irle, S.; Morokuma, K. *J. Chem. Phys.* **2004**, 121, 5163–5169.

(16) Witek, H. A.; Morokuma, K.; Stradomska, A. *J. Chem. Phys.* **2004**, 121, 5171–5178.

(17) Witek, H. A.; Morokuma, K.; Stradomska, A. *J. Theo. Comp. Chem.* **2005**, 04, 639–655.

(18) Elstner, M.; Hobza, P.; Frauenheim, T.; Suhai, S.; Kaxiras, E. *J. Chem. Phys.* **2001**, 114, 5149–5155.

(19) Zhechkov, L.; Heine, T.; Patchovskii, S.; Seifert, G.; Duarte, H. *J. Chem. Theor. Comput.* **2005**, 1, 841–847.

(20) Rapacioli, M.; Spiegelman, F.; Talbi, D.; Mineva, T.; Goursot, A.; Heine, T.; Seifert, G. *J. Chem. Phys.* **2009**, 130, 244304–10.

(21) Gaus, M.; Qiang, C.; Elstner, M. *J. Chem. Theo. Comp.* **2011**, 7, 931–948.

(22) Frauenheim, T.; Weich, F.; Köhler, T.; Uhlmann, S.; Porezag, D.; Seifert, G. *Phys. Rev. B* **1995**, 52, 11492–11501.

(23) Trani, F.; Barone, V. *J. Chem. Theo. Comp.* **2011**, 7, 713–719.

(24) Sherman, R. P.; Grinter, R. *J. Mol. Struct.: Theochem* **1986**, 135, 127–133.

(25) Masmoudi, M.; Massat, C.; Poteau, R. *Appl. Math. Comp. Sci.* **1996**, 6, 263–275.

(26) Wolff, S. *Int. J. Quantum Chem.* **2005**, 104, 645–659.

(27) Kristensen, K.; Jørgensen, P.; Thorvaldsen, A. J.; Trygve, H. *J. Chem. Phys.* **2008**, 129, 1–6.

(28) Heine, T.; Rapacioli, M.; Patchkovskii, S.; Frenzel, J.; Koster, A.; Calaminici, P.; Duarte, H. A.; Escalante, S.; Flores-Moreno, R.; Goursot, A.; Reveles, J.; Salahub, D.; Vela, A. *deMon-Nano Experiment*. http://physics.jacobs-university.de/theine/research/deMon/; accessed 2013.

(29) Yang, Y.; Yu, H.; Uork, D.; Cui, Q.; Elstner, M. *J. Phys. Chem. A* **2007**, 111, 10861–10873.

(30) Rapacioli, M.; Spiegelman, F.; Scemama, A.; Mirtschink, A. *J. Chem. Theor. Comput.* **2011**, 7, 44–55.