Tight-Binding Approximations to Time-Dependent Density Functional Theory— a fast approach for the calculation of electronically excited states

Tight-Binding Approximations to Time-Dependent Density Functional Theory
— a fast approach for the calculation of electronically excited states

Robert Rüger rueger@scm.com Scientific Computing & Modelling NV, De Boelelaan 1083, 1081 HV Amsterdam, The Netherlands Department of Theoretical Chemistry, Vrije Universiteit Amsterdam, De Boelelaan 1083, 1081 HV Amsterdam, The Netherlands Wilhelm-Ostwald-Institut für Physikalische und Theoretische Chemie, Linnéstr. 2, 04103 Leipzig, Germany    Erik van Lenthe Scientific Computing & Modelling NV, De Boelelaan 1083, 1081 HV Amsterdam, The Netherlands    Thomas Heine Wilhelm-Ostwald-Institut für Physikalische und Theoretische Chemie, Linnéstr. 2, 04103 Leipzig, Germany    Lucas Visscher Department of Theoretical Chemistry, Vrije Universiteit Amsterdam, De Boelelaan 1083, 1081 HV Amsterdam, The Netherlands
July 7, 2019
Abstract

We propose a new method of calculating electronically excited states that combines a density functional theory (DFT) based ground state calculation with a linear response treatment that employs approximations used in the time-dependent density functional based tight binding (TD-DFTB) approach. The new method termed TD-DFT+TB does not rely on the DFTB parametrization and is therefore applicable to systems involving all combinations of elements. We show that the new method yields UV/Vis absorption spectra that are in excellent agreement with computationally much more expensive time-dependent density functional theory (TD-DFT) calculations. Errors in vertical excitation energies are reduced by a factor of two compared to TD-DFTB.

pacs:
31.15.ee

I Introduction

Owing its success to the good compromise between accuracy and computational cost, density functional theory (DFT) based on the theorems Hohenberg and Kohn (1964) by Hohenberg and Kohn and employing the Kohn-Sham ansatz Kohn and Sham (1965) for the kinetic energy has become the most widely used method in both quantum chemistry and solid state theory over the last few decades. Rooted in the Kohn-Sham DFT framework, density functional based tight binding (DFTB) Porezag et al. (1995); Seifert, Porezag, and Frauenheim (1996) has been developed as a computationally very efficient approximation to DFT for systems too large to be treated with its parent method. DFTB’s efficiency stems from the use of an optimized minimum valence orbital basis that reduces the linear algebra operations, and a two center-approximation for the Kohn-Sham potential that allows precalculation and storage of integrals using the Slater-Koster technique Slater and Koster (1954). The self-consistent charge extension (SCC-DFTB, recently also called DFTB2) Elstner et al. (1998) accounts for density fluctuations and improves results for systems with polar bonds. A further extension known as DFTB3 Gaus, Cui, and Elstner (2011) has been developed to improve the description of hydrogen-bonded complexes and proton affinities. While DFTB is much more efficient than DFT it requires careful parametrization for all involved elements in order to yield accurate results. The limited availability of these parameters has historically slowed down the adoption of DFTB, but general purpose parameter sets covering large parts of the periodic table have recently become available Gaus, Goez, and Elstner (2013); Gaus et al. (2014); Lu et al. (2015); Kubillus et al. (2015); Wahiduzzaman et al. (2013); Oliveira, Philipsen, and Heine (2015).

As the Hohenberg-Kohn theorems only concern the ground state of a system, density functional theory is not applicable to the broad class of problems involving electronically excited states. The foundation for the excited state extension known as time-dependent density functional theory (TD-DFT) was later laid by Runge and Gross, who generalized Hohenberg and Kohn’s theorems to time-dependent external potentials Runge and Gross (1984). Based on their work, Casida calculated the linear response of the electron density to a perturbation in the external potential and from this derived an eigenvalue equation in the space of single orbital transitions from which the excited states of the electronic system can be obtained Casida (1995). In the field of quantum chemistry, Casida’s TD-DFT approach is today probably the most widely used method for the calculation of excited state properties. A recent review of TD-DFT can be found in reference 16.

An excited state calculation using TD-DFT is computationally quite demanding, much more so than the underlying DFT calculation of the ground state. For many systems it is therefore feasible to calculate the ground state, while a calculation of excited states is computationally out of reach. Various ways to reduce the computational complexity of TD-DFT have been put forward; based for example on partitioning into subsystems Casida and Wesołowski (2003); Neugebauer (2007), neglect of terms Hirata and Head-Gordon (1999); Grimme (2013); Risthaus, Hansen, and Grimme (2014), truncation of the single orbital transition space Grimme (2013); Rüger et al. (2015) and approximation of integrals Grimme (2013); Risthaus, Hansen, and Grimme (2014); Bannwarth and Grimme (2014); Niehaus et al. (2001a); Domínguez et al. (2013); Trani et al. (2011). Among the methods that approximate integrals is time-dependent density functional theory based tight binding (TD-DFTB) Niehaus et al. (2001a); Domínguez et al. (2013); Trani et al. (2011) developed by Niehaus et al., which builds on a DFTB ground state calculation and translates Casida’s linear response approach to the framework of DFTB. It has been found to yield very good results for transitions Niehaus et al. (2001a), making it especially suitable for the calculation of UV/Vis absorption spectra Rüger et al. (2015). Being computationally much cheaper than a full TD-DFT calculation and applicable to very large systems, TD-DFTB has been used in a variety of applications Joswig et al. (2003); Goswami et al. (2006); Frenzel, Joswig, and Seifert (2007); Li et al. (2007); Wang et al. (2007a, b); Li et al. (2008); Mitrić et al. (2009); Zhang et al. (2012); Fan et al. (2014) in which TD-DFT would not have been feasible. A review of TD-DFTB can be found in reference 37.

TD-DFTB inherits certain limitations from the DFTB ground state calculation it is based on: The electronic structure from DFTB, which is the basis of the excited state calculation, is of limited accuracy compared to a DFT calculation with a reasonable choice of exchange-correlation functional and orbital basis. Furthermore, whereas TD-DFT can be used for any system, historically the applicability of TD-DFTB was restricted to systems involving only elements for which DFTB parameters are available. With the development of the QUASINANO parameters, which are available throughout the periodic table, by Wahiduzzaman et al. Wahiduzzaman et al. (2013) this drawback was removed for the electronic part (and thus for TD-DFTB), even though TD-DFTB still requires a careful performance validation for the target system class.

These limitations are insofar particularly unfortunate as it is neither the calculation of the ground state’s electronic structure that is the computational bottleneck, nor requires the application of tight-binding approximations to the TD-DFT concept within Casida’s formulation any parameterization effort. In this article we introduce TD-DFT+TB, a new method for calculating electronically excited states that combines a DFT ground state with a linear response treatment that employs approximations similar to the ones used in TD-DFTB. We show that the cost of this calculation is approximately the same as of a ground state DFT calculation, and the accuracy of the excited state properties is much better than TD-DFTB.

The approach proposed in this article is inspired by and closely related to the sTDA Grimme (2013); Risthaus, Hansen, and Grimme (2014) and sTD-DFT Bannwarth and Grimme (2014) methods developed by Grimme et al., which also use a DFT ground state calculation and make TD-DFTB like approximations in Casida’s formalism. The main difference compared to our approach is, that these methods are based on DFT with a hybrid exchange-correlation functional Becke (1993). Hybrid functionals are usually employed to correct the underestimated charge-transfer excitation energies in TD-DFT with local functionals. However, we argue that for the calculation of optical absorption spectra, underestimated charge-transfer excitation energies only constitute a technical problem that can be solved conceptually easier and computationally more efficient by employing a physically motivated truncation of the single orbital transition space Rüger et al. (2015). Furthermore, it was recently shown Baerends, Gritsenko, and van Meer (2013); van Meer, Gritsenko, and Baerends (2014) by Baerends, Gritsenko, and van Meer that excitations loose their single orbital transition character with the admixture of Hartree-Fock exact exchange, which complicates the interpretation of the results. We therefore believe that the calculation of excited states of large systems should also be approached from a pure density functional standpoint. Both sTDA and sTD-DFT have been developed for hybrid functionals and contain free parameters that have been fitted to yield good results between 20% and 60% exact exchange. As such, they are not intended for and can not directly be used with local functionals. We believe that TD-DFT+TB as proposed in this article complements sTDA and sTD-DFT by making approximate TD-DFT methods also available for local exchange-correlation functionals.

The remainder of this article is organized as follows. In section II we recapitulate the most important equations from DFT and DFTB as well as their linear response extensions in order to set the stage for section III, in which we motivate and introduce TD-DFT+TB. We will also discuss its relation to other approximate TD-DFT methods, such as TD-DFTB, sTDA and sTD-DFT. In section IV we evaluate the accuracy and performance of the new method by calculating vertical excitation energies for a benchmark set of molecules as well as the UV/Vis absorption spectra of selected compounds. Section V summarizes our results and concludes the article.

Ii Review of methods

In order to establish the notation for the remainder of this article, this section contains a short summary of DFT and DFTB as well as their linear response extensions.

ii.1 Molecular orbitals from DFT(B)

Electronic structure calculations of molecular systems typically use atom centered basis functions , so that the molecular orbitals  can be written as

(1)

The basis functions are composed of primitives that may be Gaussian, Slater, numerical or any other functions that are centered at the atomic positions. For DFT the size of the basis is variable and can within the limits of computational affordability be chosen according to the desired accuracy, while DFTB on the other hand typically uses an optimized minimum valence orbital basis that is fixed during the DFTB parameter creation, and can not be changed at run-time.

The expansion coefficients  of the molecular orbitals are obtained by solving the secular equations

(2)

Here,

(3)

is the overlap between basis functions. In DFT, the Hamiltonian matrix elements  are calculated as

(4)

where

(5)

is the Kohn-Sham effective potential Kohn and Sham (1965), consisting of the external potential, an electrostatic term and the so-called exchange-correlation potential. Note that the effective potential  depends on the molecular orbitals themselves through their electronic density , so that equation (2) has to be solved self-consistently.

DFTB avoids the evaluation of integrals at run-time by replacing the actual density by a trial density . This trial density is a superposition of atomic contributions which are optimized within the parameterization process. Within the DFTB-inherent two-center approximation the effective potential is constructed by superposing two spherical atomic effective potentials Porezag et al. (1995); Seifert, Porezag, and Frauenheim (1996); Wahiduzzaman et al. (2013) or trial densities Elstner et al. (1998); Gaus, Cui, and Elstner (2011), which allows the Hamiltonian and overlap matrix elements to be precalculated and stored using the Slater-Koster technique Slater and Koster (1954). This shifts the computational bottleneck from the calculation of matrix elements to linear algebra operations which are dominated by the diagonalization. Together with the small matrix sizes due to the minimal valence basis set, DFTB is computationally extremely efficient.

ii.2 Excited states and absorption spectra from TD-DFT

Once the electronic structure of the ground state has been determined, excited states can be calculated using Casida’s linear response approach Casida (1995), which casts the problem of calculating excitation energies and excited states into an eigenvalue equation in the  dimensional space of single orbital transitions. The eigenvalue problem can be written as

(6)

where is the vertical excitation energy of the -th excited state. We adopt the usual convention of using the indexes  for occupied and  for virtual orbitals. The elements  correspond to the contribution of the transition from the occupied orbital  to the virtual orbital  and can be used to construct an approximate excited state wavefunction  from the Slater determinant  of the occupied Kohn-Sham orbitals Casida (1995).

(7)

Here we use  for the difference in orbital energy between the involved Kohn-Sham orbitals. The elements of the matrix  are given by

(8)

and looking back at equation (6), it is easy to see that it is the so-called coupling matrix  that shifts the excitation energies  away from the orbital energy differences . The coupling matrix depends on the multiplicity of the calculated excitations. For the sake of clarity we will restrict our discussion to singlet excitations for the moment. Triplet excitations pose no additional problems and their calculation will be discussed later. For singlet excitations in TD-DFT the elements of the coupling matrix are given by

(9)

where the kernel

(10)

incorporates both a Coulomb term and the second derivative of the DFT exchange-correlation functional .

For the prediction of photon absorption spectra it is necessary to calculate both the excitation energies  as well as the corresponding transition dipole moments . Excitation energies are immediately obtained as the eigenvalues of Casida’s equation (6), and using the eigenvector elements  the transition dipole moments  can be calculated as a linear combination of the transition dipole moments  of the single orbital transitions.

(11)

Here the transition dipole moments  of the single orbital transitions are calculated as

(12)

In order to make a connection to experimentally measured quantities, the theoretically calculated oscillator strengths  and excitation energies  can be related Mulliken (1939) to the molar absorptivity

(13)

Here is a normalized, typically peaked function that models the experimental line broadening. Both Gaussian and Lorentzian functions are common choices for .

It would be beyond the scope of this article to go into further details on the properties and problems of TD-DFT. A recent review of the strengths and weaknesses of TD-DFT in general can be found in reference 16. There is also an excellent book on TD-DFT, see reference 42. The method put forward in this article presents an approximation to TD-DFT and we therefore consider TD-DFT with a GGA exchange-correlation functional as the reference method.

ii.3 TD-DFTB as an approximation to TD-DFT

The calculation of the TD-DFT coupling matrix elements involves expensive two-center integrals and even though highly optimized implementations are available van Gisbergen, Snijders, and Baerends (1999), evaluating the integrals is still the computational bottleneck of the method. In order to make density functional based excited state calculations applicable to larger systems, Niehaus et al. have put forward TD-DFTB Niehaus et al. (2001a); Domínguez et al. (2013), which builds on top of a DFTB ground state calculation and uses DFTB-like approximations for the coupling matrix elements: The transition density in equation (9) is subjected to a multipole expansion truncated at first order (monopole approximation)

(14)

where is a spherically symmetric function centered on atom . This allows the singlet-singlet coupling matrix elements in TD-DFTB to be written as

(15)

The so-called atomic transition charges  are calculated from the molecular orbital coefficient and overlap matrix  and  using Mulliken population analysis Mulliken (1955).

(16)

While should in principle be calculated as a two-center integral over the product of the atom centered functions and the kernel from equation (10), it is in practice approximated as a function

(17)

of the internuclear distance and the chemical hardness  and  of atom  and  respectively, converging to the Coulomb interaction between two point charges for long distances Niehaus et al. (2001a); Elstner et al. (1998). The required atomic chemical hardness is not a free, tunable parameter, but rather an inherent property of the atoms themselves. There is, however, some freedom in the choice of the method used to obtain these values, e.g. from atomic DFT calculations by application of Janak’s theorem Mineva and Heine (2004, 2006); Elstner et al. (1998), or using a phenomenological model Ghosh and Islam (2009).

So far our discussion has been restricted to singlet-singlet excitations. For the calculation of singlet-triplet excitations the only change required is in the coupling matrix elements, which for singlet-triplet excitations in TD-DFTB are given by

(18)

Here the so called magnetic Hubbard parameters  are defined as

(19)

and can be calculated from atomic DFT calculations just like the chemical hardnesses.

In addition to the approximation of the coupling matrix, TD-DFTB also approximates the transition dipole moments  of the single orbital transitions. With the monopole approximation of the transition density from equation (14) the transition dipole moments of the single orbital transitions are easily written as Niehaus et al. (2001a)

(20)

One rather obvious limitation of the monopole approximation in equation (14) is that basis functions  and  residing on the same atom  do not contribute to the atomic transition charge . This leads to vanishing (or underestimated) transition charges for excitations involving localized molecular orbitals  and , such as and promotions. Due to the vanishing coupling matrix elements  these excitations are then predicted to be pure single orbital transitions  with an excitation energy  exactly. Furthermore, their transition dipole moment  is predicted to be zero. This failure has recently been corrected by Domínguez et al. through inclusion of one-center integrals of the exchange type Domínguez et al. (2013). However, this so-called on-site correction to TD-DFTB is fairly involved and we will restrict our discussion to TD-DFTB in its original formulation Niehaus et al. (2001a).

In summary, TD-DFTB is an approximation to TD-DFT that uses molecular orbitals obtained from a DFTB ground state calculation and approximates the coupling matrix and single orbital transition dipole moments in order to avoid integral evaluation at run-time. For a recent review of TD-DFTB we would like to refer the reader to reference 37.

Iii Td-Dft+tb

iii.1 Motivation and introduction

In subsection II.3 we have outlined how the TD-DFTB coupling matrix can be derived from its TD-DFT counterpart by making a monopole approximation for the transition density. While this is how TD-DFTB was originally introduced Niehaus et al. (2001a), it is interesting to note that the same equations can also be obtained as the linear response of the SCC-DFTB Hamiltonian Domínguez et al. (2013), just like TD-DFT was obtained from the linear response of DFT Casida (1995). In this sense all of the approximations that go into TD-DFTB have been done at the ground state level, and the subsequent excited state calculation is merely done consistently with the already present approximations.

This brings up an interesting question: Would more accurate results be obtained if the approximation was delayed until the linear response treatment? Or in other words, would it be better to do an approximate linear response of the DFT Hamiltonian than to look at the exact linear response of the SCC-DFTB Hamiltonian? In this article we want to propose to do TD-DFTB-like approximations in a linear response excited state calculation based on a DFT ground state. We will henceforth refer to this approach as TD-DFT+TB. The relationship between the different methods is illustrated in figure 1.

The basic idea of TD-DFT+TB is to use the molecular orbitals obtained from a DFT ground state calculation as input to an excited state calculation with the TD-DFTB coupling matrix from equation 15. Technically, this is very easy to do: Looking back at subsection II.3 it is evident that the only information needed about the ground state is the overlap matrix , the coefficient matrix  as well as the orbital energies  and occupations. Additionally, the information about which atom  the basis function  is centered on is also needed for the population analysis. However, all of this information could also be provided by a DFT instead of a DFTB ground state calculation.

One important thing to note is that the basis sets used in DFT are typically larger than the minimal basis set used in DFTB. In fact, the pre-optimized DFTB ground state densities are typically of higher quality compared to those obtained using a minimum basis set DFT approach. Therefore, it is important to employ a DFT basis that gives a sufficiently accurate ground state, even though this leads to more virtual orbitals and hence a larger coupling matrix in TD-DFT+TB compared to TD-DFTB.

A problem associated with the larger basis set in DFT is that the Mulliken population analysis Mulliken (1955) used in TD-DFTB for the calculation of the atomic transition charges  is known to become unstable for large basis sets, especially if diffuse basis functions are included. While Mulliken analysis is working sufficiently well for the minimal atomic orbital basis set used in TD-DFTB, we have found that for a basis of TZP quality Mulliken transition charges only poorly represent the transition density. This was also observed by Grimme, who instead proposed Grimme (2013) to use Löwdin population analysis Löwdin (1950) for which the atomic transition charges are calculated as

(21)

We have indeed found Löwdin transition charges to be much more reliable than the ones obtained from Mulliken population analysis, and therefore use Löwdin analysis as the default method of calculating transition charges in TD-DFT+TB. Benchmark results for both Mulliken and Löwdin transition charges can be found in section IV.1.

Figure 1: Relationships among the different computational methods.

In TD-DFTB the atomic transition charges  are used in both the approximation of the  coupling matrix elements as well as the approximation of the  single orbital transition dipole moments . While the former is what makes TD-DFTB so efficient, the latter has mostly technical reasons: The Slater-Koster files used in DFTB contain the matrix elements, but not the basis functions themselves, making it impossible to evaluate integrals over molecular orbitals at run time. This is a rather unpleasant deficiency introduced by the traditional DFTB Slater-Koster implementation. However, it is not a deficiency of the method itself, and with knowledge of the atomic orbitals obtained during the parameterization process and in combination with a suitable integral engine any expectation value can be calculated correctly using DFTB. This has been demonstrated for other properties, e.g. for NMR chemical shifts Heine et al. (1999), and will be applied to the transition dipole matrix elements here.

For TD-DFT+TB, we have two possibilities to compute the transition dipole moments:

  1. The simplest, and most approximate, way is to use the point charge approximation as done in TD-DFTB. This approximation would be most attractive if TD-DFT+TB would be used by employing two independent codes (one for DFT and one for TD-DFTB). However, in our present case there is negligible computational performance gain, so we do not follow this line.

  2. We calculate the transition dipole moments directly from the DFT molecular orbitals. Thus, we calculate the unapproximated single orbital transition dipole moments  from equation (12). In different words: We avoid here those TD-DFTB approximations that have been due to restrictions imposed by the Slater-Koster-type implementation and which are not resulting in significant performance gain.

One particularly attractive feature of TD-DFT+TB is that it does not rely at all on the DFTB parametrization. The only parameters used for the construction of the TD-DFTB coupling matrix are the chemical hardness  (for singlet excitations) and the magnetic Hubbard  parameter (for triplet excitations). These are just physical properties of the atoms that can be calculated and tabulated for the entire periodic table. We use the chemical hardness as tabulated by Ghosh and Islam Ghosh and Islam (2009) and have calculated the values for the magnetic Hubbard parameter  using the same details as specified earlier Wahiduzzaman et al. (2013). The numerical values of are given in table 5. All other parameters entering DFTB which are needed for describing the ground state, i.e. the form of the basis functions, the effective potential, and the repulsive potential needed for calculating the total energy and its gradients are not needed to build the TD-DFTB coupling matrix. TD-DFT+TB is therefore directly applicable to systems containing any combination of elements without the need of further parameterization.

TD-DFT TD-DFT+TB TD-DFTB
Molecular orbitals from DFT from DFTB
Coupling matrix  Eq. (9) Eq. (15)
Atomic transition
charges 
not used Eq. (21) Eq. (16)
Single orbital transition
dipole moments 
Eq. (12) Eq. (20)
Chemical hardness  &
Magnetic hubbard 
not used
precalculated by DFT
for spherical atoms111included in the DFTB parameters files in case of TD-DFTB
Table 1: Comparison of the methods.

In summary, TD-DFT+TB can be interpreted either as applying DFTB approximations to the Casida equations, or, equivalently, as TD-DFTB based on molecular orbitals from DFT. Technical choices are the calculation of charges and transition dipole moments. We propose to employ Löwdin instead of Mulliken atomic transition charges, and DFT transition dipole moments, but other options are definitely possible. A summarizing comparison of TD-DFT, TD-DFTB and TD-DFT+TB is given in table 1.

iii.2 Relation to other methods

TD-DFT+TB as introduced in the last subsection is quite closely related to the sTDA Grimme (2013); Risthaus, Hansen, and Grimme (2014) and sTD-DFT Bannwarth and Grimme (2014) methods developed by Grimme and coworkers: These methods also use molecular orbitals from a DFT calculation and use the same atomic monopole approximation for the transition density (which was originally introduced with TD-DFTB) in order to avoid the calculation of integrals. The major difference is that TD-DFT+TB is a pure density functional approach, while sTDA and sTD-DFT use hybrid exchange-correlation functionals Becke (1993) in both the calculation of the ground state and the excited states.

The primary reason why hybrid functionals with a fraction of exact Hartree-Fock exchange are often used in TD-DFT is that local functionals are known to drastically underestimate the excitation energies of charge-transfer excitations. It was shown Gritsenko and Baerends (2004); Baerends, Gritsenko, and van Meer (2013); van Meer, Gritsenko, and Baerends (2014) that this failure can be traced back to the different meaning of virtual orbital energies in Kohn-Sham DFT and Hartree-Fock: In DFT the virtual orbitals represent excited electrons interacting with other electrons, while in Hartree-Fock the virtual orbitals experience interaction with  electrons, so that they represent added rather than excited electrons. In other words, the Kohn-Sham HOMO-LUMO gap corresponds to the optical gap, while the Hartree-Fock HOMO-LUMO gap corresponds the fundamental gap, which is the difference between ionization energy and electron affinity. It is easy to see why this leads to underestimated charge-transfer excitation energies in DFT: If occupied and virtual orbital involved in a transition are localized on different fragments of the system, the transfered electron is essentially added to the acceptor fragment and its energy is determined by the acceptor’s fundamental gap, not by its optical gap. The fundamental gap is always larger than the optical gap, the difference being the interaction between the excited electron and the hole in its (now unoccupied) original orbital Baerends, Gritsenko, and van Meer (2013); Bredas (2014). In summary, local excitations are well described in Kohn-Sham DFT with local functionals, while charge-transfer excitations profit from admixture of exact exchange. The so-called range-separated hybrid functionals Iikura et al. (2001); Henderson, Janesko, and Scuseria (2008); Baer, Livshits, and Salzner (2010), where the amount of exact exchange increases with electron-electron distance, reflect this.

It is interesting to note though, that charge-transfer excitations typically have very small oscillator strengths. Looking at equation (12) it is easy to see that the transition dipole moment is zero if the involved orbitals  and  have no significant overlap, as is the case for charge-transfer excitations. So even though charge-transfer excitation energies are severely underestimated in DFT with local functionals, the obtained absorption spectra are usually not affected. There is, however, a technical problem associated with the underestimated charge transfer excitation energies for the specific application of calculating optical absorption spectra: Since, the matrix  that has to be diagonalized in Casida’s equation (6) is extremely large, it is typically diagonalized using iterative eigensolvers that only calculate the few lowest eigenvectors. If large numbers of spurious charge-transfer excitations are now predicted at much too low energies, many more excitations have to be calculated in order to cover the relevant energy range. This drastically slows down the calculation even though the spurious charge-transfer excitations do not noticeably affect the obtained absorption spectra. Grimme cites this issue as the main reason for the use of hybrid functionals in sTDA and sTD-DFT. However, as we have recently shown the problem can also be solved by intensity selection Rüger et al. (2015), that is by simply neglecting single orbital transitions with small transition dipole moments. This does not correct the energy of charge-transfer excitations but instead removes them from the spectrum altogether, leading to both a smaller number of excitations that have to be calculated as well as an overall smaller matrix  due to the reduced number of single orbital transitions. While hybrid functionals are likely unavoidable if one needs accurate charge-transfer excitation energies, we believe that for the specific application of calculating absorption spectra, intensity selection is a much simpler and computationally more efficient alternative to hybrid functionals.

Furthermore, we would like to point out that while the use of hybrid functionals cures the charge-transfer problem, it introduces other problems that are not present in pure density functional approaches: As pointed out by Baerends, Gritsenko, and van Meer virtual orbitals from Kohn-Sham DFT with local functionals represent excited electrons interacting with other electrons. The coupling matrix in equation (8) is usually small compared to the orbital energy differences  on the diagonal, making orbital energy differences  an excellent approximation to excitation energies  Zhang and Musgrave (2007). Furthermore the excitations are often dominated by just one single orbital transition , which drastically simplifies their interpretation Baerends, Gritsenko, and van Meer (2013); van Meer, Gritsenko, and Baerends (2014). Hartree-Fock virtual orbitals on the other hand represent added electrons, so that their orbital energy differences  have little relation to excitation energies  and are in fact much larger. It is actually not uncommon for the Hartree-Fock LUMO to be unbound with an orbital energy of . Furthermore, as the Hartree-Fock virtual orbitals interact with other electrons instead of , they are much more diffuse than in DFT. The Hartree-Fock virtuals are less suitable for the description of excited electrons and in general more of them are needed for the description of an excitation, meaning that excitations often lose the single orbital transition character they have in DFT, making their interpretation much more difficult Baerends, Gritsenko, and van Meer (2013); van Meer, Gritsenko, and Baerends (2014). These problems are less severe if the employed exchange-correlation functional only has a small fraction of exact exchange. It is, however, important to be aware of the fact that certain undesirable properties of time-dependent Hartree-Fock are reintroduced into TD-DFT if hybrid functionals are used.

In summary, we believe that there are good reasons to also approach excited state calculations for large systems from a pure density functional standpoint.

Iv Method evaluation

iv.1 Vertical excitation energies

In order to assess the accuracy of TD-DFT+TB we have calculated the lowest few excitation energies for the 28 molecules containing 1st and 2nd period elements in the benchmark set developed by Schreiber et al. Schreiber et al. (2008). For a direct comparison we have done the same calculations with TD-DFTB using the mio-1-1 set of parameters Elstner et al. (1998); Niehaus et al. (2001b); Yang et al. (2008); Gaus, Cui, and Elstner (2011). We use TD-DFT results as the reference against which TD-DFTB and TD-DFT+TB are compared. Both TD-DFT and TD-DFT+TB results were obtained using the PBE exchange-correlation functional Perdew, Burke, and Ernzerhof (1996) and a TZP basis set.

Note that some excitations in the benchmark set by by Schreiber et al. Schreiber et al. (2008) have significant double excitation character and are hence difficult to describe with conventional TD-DFT. See reference 60 and 61 for a detailed discussion and possible solutions to this problem. However, for the purpose of comparing TD-DFT+TB and TD-DFTB to TD-DFT this does not play a role, as all methods are equally affected by this issue.

The calculated vertical excitation energies and the root-mean-square deviation (RMSD) compared to TD-DFT are shown for the individual molecules in figure 2 and figure 3 for singlet-singlet and singlet-triplet excitations, respectively.

Figure 2: Vertical singlet-singlet excitation energies (left ordinate) for the molecules from Schreiber et al.’s test set Schreiber et al. (2008). The bars at the bottom represent the RMSD in vertical excitation energies compared to TD-DFT (to scale with the right ordinate).
Figure 3: Vertical singlet-triplet excitation energies (left ordinate) for the molecules from Schreiber et al.’s test set Schreiber et al. (2008). The bars at the bottom represent the RMSD in vertical excitation energies compared to TD-DFT (to scale with the right ordinate).

The calculated RMSD of all excitations in all molecules are shown in table 2. Following Casida et al.’s recommendation Casida et al. (1998) we only considered excitations that have an excitation energy  and no large contributions from transitions into unbound virtual orbitals with . In cases where the number of excitations satisfying these criteria differs between the methods, we compare the lowest common number of excitations. For singlet-singlet excitations in ethene and furan none of the calculated excitations satisfies both of Casida et al.’s criteria, so that these two molecules had to be excluded from the calculation of the overall RMSD for singlet-singlet excitations.

Compared to normal TD-DFTB, TD-DFT+TB is closer to TD-DFT for both singlet-singlet and singlet-triplet excitations. For singlet-singlet excitations switching from TD-DFTB to TD-DFT+TB reduces the RMSD by a factor of two from 0.301eV to 0.153eV. It is known that TD-DFTB is more accurate for singlet-singlet excitations than for singlet-triplet excitations Niehaus et al. (2001a) for which we calculated an RMSD of 0.489eV. We observe the same behavior for TD-DFT+TB, although with an RMSD of 0.215eV the difference in accuracy between singlet-singlet and singlet-triplet excitations is slightly smaller.

Multiplicity TD-DFTB TD-DFT+TB
singlet-singlet 0.30eV 0.15eV
singlet-triplet 0.49eV 0.22eV
Table 2: Root-mean-square deviation in vertical excitation energies of TD-DFTB and TD-DFT+TB for Schreiber et al.’s test set. TD-DFT is used as the reference method.
transition TD-DFT TD-DFT+TB TD-DFTB
() 13.47 0.32 11.53 1.94 11.53 0.64 11.53 0.00 12.73 0.00 12.73 0.00
14.94 0.34 9.55 4.94 12.99 0.98 9.55 3.44 13.90 0.88 10.19 3.71
Table 3: Dipole allowed transitions in N. All energies in eV. is the orbital energy difference for the dominant single orbital transition.

Note that for the calculation of the RMSD we have simply compared the excitation energies from the different methods according to their order in energy. We have not attempted to compensate for the fact that two excited states might switch in energy ordering when going from TD-DFT to one of the approximate methods. While this does not affect the comparison between TD-DFTB and TD-DFT+TB, the absolute errors in table 2 will be slightly underestimated and one should be careful when comparing them to the literature.

As mentioned in section III we have also run TD-DFT+TB calculations for Schreiber et al.’s test set using Mulliken instead of Löwdin population analysis for the calculation of the atomic transition charges. We found an RMSD of 0.449eV in vertical singlet-singlet excitation energies, which is three times larger than the 0.153eV obtained with Löwdin charges, indicating that Mulliken transition charges do not accurately model the transition density for the relatively large TZP basis set used. For singlet-triplet excitations we have furthermore found that unphysical transition charges sometimes lead to negative eigenvalues in equation (6) and hence imaginary excitation energies.

iv.2 Oscillator strengths and absorption spectra

In the last subsection we looked exclusively at vertical excitation energies. However, for the application of calculating UV/Vis absorption spectra both excitation energies and oscillator strengths have to be calculated.

iv.2.1 N

One difference between TD-DFTB and TD-DFT+TB is that the latter does not use the atomic transition charges for the calculation of the single orbital transition dipole moments. To illustrate the effect of this we have calculated the lowest excitations in N with a nuclear distance of 1.106Å. The results are shown in table 3. According to TD-DFT there are two dipole allowed transitions: A  state consisting mainly of a transition, and a  state dominated by a transition. Note that even though both of them have excitation energies we have found them to be well described and largely basis set independent due to the fact that they do not have contributions from transitions into unbound virtual orbitals. The  state is reasonably well described by both TD-DFTB and TD-DFT+TB, who both predict it to be dipole-allowed. Both methods underestimate the vertical excitation energy  of the  state, with the TD-DFTB energy being closer to the TD-DFT reference. However, this is mostly due to the larger orbital energy difference  in DFTB compared to DFT, since the coupling matrix induced shift is similar for both TD-DFTB and TD-DFT+TB. The transition into the  state is less well described with the approximate methods. TD-DFT predicts a coupling matrix induced shift  of almost 2eV while both approximate methods produce exactly a single orbital transition with . This is due to the atomic transition charges’ inability to model local transitions as mentioned in subsection II.3. Since TD-DFTB also uses the transition charges for the transition dipole moments, it incorrectly predicts the transition into the  state state to be dipole-forbidden. This is not the case in TD-DFT+TB, so that the method can at least be used to identify dipole-allowed and transitions, even though their excitation energies will be less accurate than those of  transitions.

However, for the large systems such approximate method are typically used for,  transitions usually have the largest oscillator strengths, so that TD-DFTB and TD-DFT+TB’s problems with localized transitions often do not noticeably affect the calculated absorption spectra.

iv.2.2 Fullerene C

As an example for the calculation of absorption spectra, we have calculated the UV/Vis spectrum of the C fullerene. This was one of the example systems in the original TD-DFTB article and also makes a good technical benchmark as almost a thousand excitations have to be calculated to cover the relevant energy range. For TD-DFT and TD-DFT+TB we used a TZP basis and the PBE exchange-correlation functional Perdew, Burke, and Ernzerhof (1996). TD-DFTB calculations were performed with both the 3ob-3-1 parameter set Gaus, Goez, and Elstner (2013); Gaus et al. (2014); Lu et al. (2015); Kubillus et al. (2015) and the QUASINANO2013 set by Wahiduzzaman et al. Wahiduzzaman et al. (2013). For calculations with the 3ob-3-1 parameter set the ground state calculation was performed at the DFTB3 Gaus, Cui, and Elstner (2011) level of theory. Conceptually this is slightly inconsistent, as the calculation of the excited states is based on the linear response of a Hamiltonian different from the one used for the calculation of the ground state. However, since the DFTB3 orbitals are generally of better quality than DFTB2 orbitals, this gives rather good results in practice.

The calculated spectra are shown in figure 4.

Figure 4: Absorption spectrum of the C fullerene. Experimental gas phase absorption spectrum from reference 63. Note that the authors quote a 100% uncertainty in the absolute absorption cross sections due to the vapor pressure–temperature relation. Theoretical spectra have been broadened with a Gaussian.

TD-DFT+TB reproduces the TD-DFT reference spectrum almost perfectly. TD-DFTB with the 3ob-3-1 parameter performs very well below 5.5eV but underestimates the intensity of an excitation seen at 5.9eV in the TD-DFT spectrum. All three spectra qualitatively reproduce the series of absorption bands of increasing intensity seen in the experimental spectrum Yasumatsu et al. (1996). However, the theoretical spectra are redshifted compared to experiment. Absolute intensities should not be compared to experiment, as the experimentally measured cross sections have an uncertainty of due to the sensitive vapor pressure–temperature relation of fullerenes Abrefah et al. (1992).

The TD-DFTB spectrum calculated with the QUASINANO2013.1 parameters shows a substantial blue-shift compared to the other methods and the experimental reference. However, the shape of the spectrum with its three bands of increasing oscillator strength is reasonably well described. The origin of the blue-shift can be traced back to differences in the Kohn-Sham orbital energies: DFT and DFTB with the 3ob-3-1 parameters show a HOMO-LUMO gap of about 1.6eV, while DFTB with the QUASINANO2013.1 parameters produces a gap of 2.3eV. Keeping in mind that the HOMO-LUMO gap in DFT represents the optical gap Zhang and Musgrave (2007); Baerends, Gritsenko, and van Meer (2013); van Meer, Gritsenko, and Baerends (2014), it is easy to understand why the QUASINANO2013.1 parameters predict overall larger excitation energies. The reason for the larger orbital energy differences with the QUASINANO2013.1 parameters is that they were optimized to reproduce band structures in solids Wahiduzzaman et al. (2013), for which relatively tight confinement potentials are required in the atomic calculations. However, the additional potential leads to increased orbital energy differences through quantum confinement and produces systematically overestimated excitation energies and blue-shifted absorption spectra. This illustrates how strongly TD-DFTB results can depend on details of the DFTB parametrization; a problem that does not exist in TD-DFT+TB.

Timings for the calculation of the C absorption spectra are shown in table 4.

TD-DFT TD-DFT+TB TD-DFTB
ground state 4min 38s 4min 33s < 1s
excited states 19h 37min 11min 35s 1min 26s
Table 4: Timings for the calculation of the 988 lowest singlet-singlet excitations in the C fullerene. The obtained spectra are shown in figure 4. All calculations were performed on an Intel Core i7-4770 processor.

The benchmark TD-DFT calculation take almost 20 hours on a recent workstation computer, only 5 minutes of which are spent calculating the ground state. With TD-DFT+TB the total wall-time decreases to about 16 minutes, which is a speedup by a factor of 73 compared to TD-DFT. With a total wall-time of less than 90 seconds TD-DFTB still much faster that TD-DFT+TB. The DFTB ground state calculation takes less than a second and is therefore completely negligible compared to the 5 minutes for the DFT ground state in TD-DFT+TB. Furthermore, due to the minimal basis set the space of single orbital transitions is much smaller in TD-DFTB, leading to a smaller matrix to be diagonalized: For DFT with a TZP basis set there are single orbital transitions, whereas DFTB only has transitions.

iv.2.3 Chlorophyll A

We have also calculated the UV/Vis absorption spectrum of chlorophyll A. Due to the magnesium ion at the center of the chlorin ring, DFTB parameters that allow calculations of chlorophyll have only recently become available Wahiduzzaman et al. (2013); Lu et al. (2015). The calculated spectra are shown in figure 5. The agreement between TD-DFT and TD-DFT+TB is again almost perfect throughout the entire energy range and both methods show the well-known and Soret absorption bands around 1.95eV and 2.8eV, respectively. The spectrum obtained with TD-DFTB and the 3ob-3-1 parameters is very close to TD-DFT below 3.2eV, but differs somewhat beyond that. All three methods reproduce the essential features of the experimental absorption spectrum Strain, Thomas, and Katz (1963); Du et al. (1998), although the energy gap between and Soret band is slightly underestimated. Note, however, that the experimental spectrum was recorded in solution, while our calculation corresponds to absorption in the gas phase. In the region above the theoretical spectra show more structure than the relatively flat experimental spectrum, which we attribute to the neglect of vibrational broadening in the theoretical spectra. With the QUASINANO2013 parameters we again observe a blue-shift of the entire spectrum, so that and Soret band are predicted at 2.5eV and 3.6eV respectively.

Figure 5: Absorption spectrum of chlorophyll A. Experimental spectrum measured in diethyl ether prepared by Scott Prahl based on reference 65 and 66. Theoretical spectra have been broadened with a Gaussian.

iv.2.4 Ir(ppy)

Our last example calculation is the UV/Vis absorption spectrum of fac-Ir(ppy) (an abbreviation for fac-Tris(2-phenylpyridine)iridium), a compound that is of interest in the context of highly efficient organic light emitting diodes Baldo et al. (1999). The calculated absorption spectra are shown in figure 6.

Figure 6: Absorption spectrum of Ir(ppy). Experimental spectrum measured in dichloromethane from reference 68. Note that the experimental reference does not give absolute absorptivities. The experimental spectrum was therefore scaled to reproduce the TD-DFT value at the peak just above . Theoretical spectra have been broadened with a Gaussian.

Below 4.4eV TD-DFT+TB agrees well with the TD-DFT reference spectrum. Beyond that energy range the oscillator strengths from TD-DFT+TB seem to be overestimated, so that the predicted absorption is overall too strong. Both methods reproduce the principal features of the experimentally measured spectrum Fine et al. (2012), though the absorption spectra are slightly redshifted compared to experiment. Note that the experimental spectrum was measured in solution, while our calculations correspond to gas phase absorption. Absolute experimental absorption coefficients were not given in reference 68 and can hence not be compared to theory. Due to the iridium atom at the center of the complex, DFTB calculations of Ir(ppy) can at the moment only be performed with the QUASINANO2013 set of parameters. For both the fullerene and the chlorophyll we had observed a blue-shift in the spectra calculated with these parameters, which is again the case for Ir(ppy).

V Conclusion

In summary we have presented a new method for calculating electronic excitations that combines molecular orbitals from a DFT ground state calculation with TD-DFTB like approximations for the coupling matrix from Casida’s linear response formalism. We have shown that the new method named TD-DFT+TB improves vertical excitations energies compared to TD-DFTB and yields electronic absorption spectra that almost perfectly agree with computationally much more costly TD-DFT calculations. In contrast to TD-DFTB, TD-DFT+TB does not rely on DFTB parametrization and is therefore applicable to molecular systems containing any combination of elements.

The new method is very easy to implement into existing DFT codes that already have support for TD-DFT, since it is essentially only a simplification of the coupling matrix. Alternatively it could also very easily be supported by a standalone DFTB implementation with TD-DFTB support: Instead of calculating the molecular orbitals using DFTB, one could read orbitals calculated by an external DFT code from disk and use them as input for the TD-DFTB calculation. While both approaches are viable, we believe that direct integration into a DFT code is the more user-friendly alternative. We have integrated TD-DFT+TB in this way into the 2016 release of the ADF modeling suite te Velde et al. (2001).

Our method is implemented and can be used for a wide range of systems where TD-DFT is computationally unfeasible. However, there is still room for improvements, and we are currently working on several enhancements and further validations. For heavier elements we are currently investigating whether orbital-dependent hardness parameters give superior performance compared to the presently used atomic ones. We are further assessing the performance of the approach by comparing smaller (DZP) and larger (TZ2P) basis sets.

Looking at the bigger picture, there are by now several related methods for the calculation of exited state properties of large systems, namely TD-DFTB, TD-DFT+TB, sTDA and sTD-DFT. It would be desirable to benchmark all these different methods against experimental data in order to be able to give clear recommendations to end users, regarding applicability and accuracy of the various methods. Based on the discussion in section III.2, we for example expect intensity selected TD-DFT+TB to be very suitable for the calculation and analysis of absorption spectra, while sTD-DFT should yield generally more accurate excitation energies. We believe that a consistently done benchmark study including all the different methods could serve to give end users the right tools for the right applications and would in general make such methods more approachable for non-expert users.

Acknowledgements.
The research leading to these results has received funding from the European Union’s Seventh Framework Programme (FP7-PEOPLE-2012-ITN) under project PROPAGATE, GA 316897.

References

Element [Ha]
H 1 -0.0717
He 2 -0.0865
Li 3 -0.0198
Be 4 -0.0230
B 5 -0.0196
C 6 -0.0226
N 7 -0.0254
O 8 -0.0278
F 9 -0.0298
Ne 10 -0.0317
Na 11 -0.0152
Mg 12 -0.0166
Al 13 -0.0140
Si 14 -0.0144
P 15 -0.0149
S 16 -0.0155
Cl 17 -0.0161
Ar 18 -0.0166
K 19 -0.0107
Ca 20 -0.0120
Sc 21 -0.0124
Ti 22 -0.0138
V 23 -0.0141
Cr 24 -0.0138
Mn 25 -0.0150
Fe 26 -0.0154
Co 27 -0.0158
Ni 28 -0.0168
Cu 29 -0.0171
Zn 30 -0.0169
Ga 31 -0.0134
Ge 32 -0.0136
As 33 -0.0136
Se 34 -0.0137
Br 35 -0.0138
Kr 36 -0.0138
Rb 37 -0.0096
Sr 38 -0.0107
Y 39 -0.0097
Zr 40 -0.0107
Nb 41 -0.0113
Mo 42 -0.0125
Tc 43 -0.0127
Ru 44 -0.0132
Rh 45 -0.0134
Element [Ha]
Pd 46 -0.0136
Ag 47 -0.0137
Cd 48 -0.0138
In 49 -0.0115
Sn 50 -0.0117
Sb 51 -0.0116
Te 52 -0.0115
I 53 -0.0114
Xe 54 -0.0114
Cs 55 -0.0083
Ba 56 -0.0094
La 57 -0.0089
Ce 58 -0.0090
Pr 59 -0.0111
Nd 60 -0.0116
Pm 61 -0.0120
Sm 62 -0.0124
Eu 63 -0.0127
Gd 64 -0.0091
Tb 65 -0.0132
Dy 66 -0.0134
Ho 67 -0.0137
Er 68 -0.0139
Tm 69 -0.0141
Yb 70 -0.0142
Lu 71 -0.0090
Hf 72 -0.0098
Ta 73 -0.0104
W 74 -0.0107
Re 75 -0.0109
Os 76 -0.0111
Ir 77 -0.0112
Pt 78 -0.0113
Au 79 -0.0108
Hg 80 -0.0114
Tl 81 -0.0107
Pb 82 -0.0110
Bi 83 -0.0109
Po 84 -0.0108
At 85 -0.0107
Rn 86 -0.0106
Fr 87 -0.0082
Ra 88 -0.0092
Ac 89 -0.0080
Th 90 -0.0084
Table 5: Values for the magnetic Hubbard parameters .
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
""
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
   
Add comment
Cancel
Loading ...
84152
This is a comment super asjknd jkasnjk adsnkj
Upvote
Downvote
""
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters
Submit
Cancel

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test
Test description