TightBinding Approximations to TimeDependent Density Functional Theory
— a fast approach for the calculation of electronically excited states
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 timedependent density functional based tight binding (TDDFTB) approach. The new method termed TDDFT+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 timedependent density functional theory (TDDFT) calculations. Errors in vertical excitation energies are reduced by a factor of two compared to TDDFTB.
pacs:
31.15.eeI 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 KohnSham 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 KohnSham 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 centerapproximation for the KohnSham potential that allows precalculation and storage of integrals using the SlaterKoster technique Slater and Koster (1954). The selfconsistent charge extension (SCCDFTB, 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 hydrogenbonded 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 HohenbergKohn 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 timedependent density functional theory (TDDFT) was later laid by Runge and Gross, who generalized Hohenberg and Kohn’s theorems to timedependent 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 TDDFT approach is today probably the most widely used method for the calculation of excited state properties. A recent review of TDDFT can be found in reference 16.
An excited state calculation using TDDFT 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 TDDFT have been put forward; based for example on partitioning into subsystems Casida and Wesołowski (2003); Neugebauer (2007), neglect of terms Hirata and HeadGordon (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 timedependent density functional theory based tight binding (TDDFTB) 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 TDDFT calculation and applicable to very large systems, TDDFTB 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 TDDFT would not have been feasible. A review of TDDFTB can be found in reference 37.
TDDFTB 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 exchangecorrelation functional and orbital basis. Furthermore, whereas TDDFT can be used for any system, historically the applicability of TDDFTB 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 TDDFTB), even though TDDFTB 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 tightbinding approximations to the TDDFT concept within Casida’s formulation any parameterization effort. In this article we introduce TDDFT+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 TDDFTB. 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 TDDFTB.
The approach proposed in this article is inspired by and closely related to the sTDA Grimme (2013); Risthaus, Hansen, and Grimme (2014) and sTDDFT Bannwarth and Grimme (2014) methods developed by Grimme et al., which also use a DFT ground state calculation and make TDDFTB like approximations in Casida’s formalism. The main difference compared to our approach is, that these methods are based on DFT with a hybrid exchangecorrelation functional Becke (1993). Hybrid functionals are usually employed to correct the underestimated chargetransfer excitation energies in TDDFT with local functionals. However, we argue that for the calculation of optical absorption spectra, underestimated chargetransfer 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 HartreeFock 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 sTDDFT 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 TDDFT+TB as proposed in this article complements sTDA and sTDDFT by making approximate TDDFT methods also available for local exchangecorrelation 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 TDDFT+TB. We will also discuss its relation to other approximate TDDFT methods, such as TDDFTB, sTDA and sTDDFT. 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 runtime.
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 KohnSham effective potential Kohn and Sham (1965), consisting of the external potential, an electrostatic term and the socalled exchangecorrelation potential. Note that the effective potential depends on the molecular orbitals themselves through their electronic density , so that equation (2) has to be solved selfconsistently.
DFTB avoids the evaluation of integrals at runtime 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 DFTBinherent twocenter 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 SlaterKoster 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 TDDFT
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 KohnSham orbitals Casida (1995).
(7) 
Here we use for the difference in orbital energy between the involved KohnSham 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 socalled 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 TDDFT 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 exchangecorrelation 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 TDDFT. A recent review of the strengths and weaknesses of TDDFT in general can be found in reference 16. There is also an excellent book on TDDFT, see reference 42. The method put forward in this article presents an approximation to TDDFT and we therefore consider TDDFT with a GGA exchangecorrelation functional as the reference method.
ii.3 TDDFTB as an approximation to TDDFT
The calculation of the TDDFT coupling matrix elements involves expensive twocenter 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 TDDFTB Niehaus et al. (2001a); Domínguez et al. (2013), which builds on top of a DFTB ground state calculation and uses DFTBlike 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 singletsinglet coupling matrix elements in TDDFTB to be written as
(15) 
The socalled 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 twocenter 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 singletsinglet excitations. For the calculation of singlettriplet excitations the only change required is in the coupling matrix elements, which for singlettriplet excitations in TDDFTB 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, TDDFTB 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 onecenter integrals of the exchange type Domínguez et al. (2013). However, this socalled onsite correction to TDDFTB is fairly involved and we will restrict our discussion to TDDFTB in its original formulation Niehaus et al. (2001a).
In summary, TDDFTB is an approximation to TDDFT 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 runtime. For a recent review of TDDFTB we would like to refer the reader to reference 37.
Iii TdDft+tb
iii.1 Motivation and introduction
In subsection II.3 we have outlined how the TDDFTB coupling matrix can be derived from its TDDFT counterpart by making a monopole approximation for the transition density. While this is how TDDFTB 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 SCCDFTB Hamiltonian Domínguez et al. (2013), just like TDDFT was obtained from the linear response of DFT Casida (1995). In this sense all of the approximations that go into TDDFTB 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 SCCDFTB Hamiltonian? In this article we want to propose to do TDDFTBlike approximations in a linear response excited state calculation based on a DFT ground state. We will henceforth refer to this approach as TDDFT+TB. The relationship between the different methods is illustrated in figure 1.
The basic idea of TDDFT+TB is to use the molecular orbitals obtained from a DFT ground state calculation as input to an excited state calculation with the TDDFTB 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 preoptimized 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 TDDFT+TB compared to TDDFTB.
A problem associated with the larger basis set in DFT is that the Mulliken population analysis Mulliken (1955) used in TDDFTB 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 TDDFTB, 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 TDDFT+TB. Benchmark results for both Mulliken and Löwdin transition charges can be found in section IV.1.
In TDDFTB 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 TDDFTB so efficient, the latter has mostly technical reasons: The SlaterKoster 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 SlaterKoster 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 TDDFT+TB, we have two possibilities to compute the transition dipole moments:

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

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 TDDFTB approximations that have been due to restrictions imposed by the SlaterKostertype implementation and which are not resulting in significant performance gain.
One particularly attractive feature of TDDFT+TB is that it does not rely at all on the DFTB parametrization. The only parameters used for the construction of the TDDFTB 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 TDDFTB coupling matrix. TDDFT+TB is therefore directly applicable to systems containing any combination of elements without the need of further parameterization.
TDDFT  TDDFT+TB  TDDFTB  

Molecular orbitals  from DFT  from DFTB  
Coupling matrix  Eq. (9)  Eq. (15)  

not used  Eq. (21)  Eq. (16)  

Eq. (12)  Eq. (20)  

not used 

In summary, TDDFT+TB can be interpreted either as applying DFTB approximations to the Casida equations, or, equivalently, as TDDFTB 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 TDDFT, TDDFTB and TDDFT+TB is given in table 1.
iii.2 Relation to other methods
TDDFT+TB as introduced in the last subsection is quite closely related to the sTDA Grimme (2013); Risthaus, Hansen, and Grimme (2014) and sTDDFT 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 TDDFTB) in order to avoid the calculation of integrals. The major difference is that TDDFT+TB is a pure density functional approach, while sTDA and sTDDFT use hybrid exchangecorrelation 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 HartreeFock exchange are often used in TDDFT is that local functionals are known to drastically underestimate the excitation energies of chargetransfer 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 KohnSham DFT and HartreeFock: In DFT the virtual orbitals represent excited electrons interacting with other electrons, while in HartreeFock the virtual orbitals experience interaction with electrons, so that they represent added rather than excited electrons. In other words, the KohnSham HOMOLUMO gap corresponds to the optical gap, while the HartreeFock HOMOLUMO 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 chargetransfer 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 KohnSham DFT with local functionals, while chargetransfer excitations profit from admixture of exact exchange. The socalled rangeseparated hybrid functionals Iikura et al. (2001); Henderson, Janesko, and Scuseria (2008); Baer, Livshits, and Salzner (2010), where the amount of exact exchange increases with electronelectron distance, reflect this.
It is interesting to note though, that chargetransfer 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 chargetransfer excitations. So even though chargetransfer 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 chargetransfer 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 chargetransfer 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 sTDDFT. 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 chargetransfer 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 chargetransfer 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 chargetransfer 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 KohnSham 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). HartreeFock 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 HartreeFock LUMO to be unbound with an orbital energy of . Furthermore, as the HartreeFock virtual orbitals interact with other electrons instead of , they are much more diffuse than in DFT. The HartreeFock 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 exchangecorrelation functional only has a small fraction of exact exchange. It is, however, important to be aware of the fact that certain undesirable properties of timedependent HartreeFock are reintroduced into TDDFT 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 TDDFT+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 TDDFTB using the mio11 set of parameters Elstner et al. (1998); Niehaus et al. (2001b); Yang et al. (2008); Gaus, Cui, and Elstner (2011). We use TDDFT results as the reference against which TDDFTB and TDDFT+TB are compared. Both TDDFT and TDDFT+TB results were obtained using the PBE exchangecorrelation 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 TDDFT. See reference 60 and 61 for a detailed discussion and possible solutions to this problem. However, for the purpose of comparing TDDFT+TB and TDDFTB to TDDFT this does not play a role, as all methods are equally affected by this issue.
The calculated vertical excitation energies and the rootmeansquare deviation (RMSD) compared to TDDFT are shown for the individual molecules in figure 2 and figure 3 for singletsinglet and singlettriplet excitations, respectively.
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 singletsinglet 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 singletsinglet excitations.
Compared to normal TDDFTB, TDDFT+TB is closer to TDDFT for both singletsinglet and singlettriplet excitations. For singletsinglet excitations switching from TDDFTB to TDDFT+TB reduces the RMSD by a factor of two from 0.301eV to 0.153eV. It is known that TDDFTB is more accurate for singletsinglet excitations than for singlettriplet excitations Niehaus et al. (2001a) for which we calculated an RMSD of 0.489eV. We observe the same behavior for TDDFT+TB, although with an RMSD of 0.215eV the difference in accuracy between singletsinglet and singlettriplet excitations is slightly smaller.
Multiplicity  TDDFTB  TDDFT+TB 

singletsinglet  0.30eV  0.15eV 
singlettriplet  0.49eV  0.22eV 
transition  TDDFT  TDDFT+TB  TDDFTB  

()  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 
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 TDDFT to one of the approximate methods. While this does not affect the comparison between TDDFTB and TDDFT+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 TDDFT+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 singletsinglet 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 singlettriplet 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 TDDFTB and TDDFT+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 TDDFT 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 TDDFTB and TDDFT+TB, who both predict it to be dipoleallowed. Both methods underestimate the vertical excitation energy of the state, with the TDDFTB energy being closer to the TDDFT 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 TDDFTB and TDDFT+TB. The transition into the state is less well described with the approximate methods. TDDFT 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 TDDFTB also uses the transition charges for the transition dipole moments, it incorrectly predicts the transition into the state state to be dipoleforbidden. This is not the case in TDDFT+TB, so that the method can at least be used to identify dipoleallowed 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 TDDFTB and TDDFT+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 TDDFTB article and also makes a good technical benchmark as almost a thousand excitations have to be calculated to cover the relevant energy range. For TDDFT and TDDFT+TB we used a TZP basis and the PBE exchangecorrelation functional Perdew, Burke, and Ernzerhof (1996). TDDFTB calculations were performed with both the 3ob31 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 3ob31 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.
TDDFT+TB reproduces the TDDFT reference spectrum almost perfectly. TDDFTB with the 3ob31 parameter performs very well below 5.5eV but underestimates the intensity of an excitation seen at 5.9eV in the TDDFT 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 TDDFTB spectrum calculated with the QUASINANO2013.1 parameters shows a substantial blueshift 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 blueshift can be traced back to differences in the KohnSham orbital energies: DFT and DFTB with the 3ob31 parameters show a HOMOLUMO gap of about 1.6eV, while DFTB with the QUASINANO2013.1 parameters produces a gap of 2.3eV. Keeping in mind that the HOMOLUMO 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 blueshifted absorption spectra. This illustrates how strongly TDDFTB results can depend on details of the DFTB parametrization; a problem that does not exist in TDDFT+TB.
Timings for the calculation of the C absorption spectra are shown in table 4.
TDDFT  TDDFT+TB  TDDFTB  

ground state  4min 38s  4min 33s  < 1s 
excited states  19h 37min  11min 35s  1min 26s 
The benchmark TDDFT calculation take almost 20 hours on a recent workstation computer, only 5 minutes of which are spent calculating the ground state. With TDDFT+TB the total walltime decreases to about 16 minutes, which is a speedup by a factor of 73 compared to TDDFT. With a total walltime of less than 90 seconds TDDFTB still much faster that TDDFT+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 TDDFT+TB. Furthermore, due to the minimal basis set the space of single orbital transitions is much smaller in TDDFTB, 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 TDDFT and TDDFT+TB is again almost perfect throughout the entire energy range and both methods show the wellknown and Soret absorption bands around 1.95eV and 2.8eV, respectively. The spectrum obtained with TDDFTB and the 3ob31 parameters is very close to TDDFT 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 blueshift of the entire spectrum, so that and Soret band are predicted at 2.5eV and 3.6eV respectively.
iv.2.4 Ir(ppy)
Our last example calculation is the UV/Vis absorption spectrum of facIr(ppy) (an abbreviation for facTris(2phenylpyridine)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.
Below 4.4eV TDDFT+TB agrees well with the TDDFT reference spectrum. Beyond that energy range the oscillator strengths from TDDFT+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 blueshift 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 TDDFTB like approximations for the coupling matrix from Casida’s linear response formalism. We have shown that the new method named TDDFT+TB improves vertical excitations energies compared to TDDFTB and yields electronic absorption spectra that almost perfectly agree with computationally much more costly TDDFT calculations. In contrast to TDDFTB, TDDFT+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 TDDFT, 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 TDDFTB 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 TDDFTB calculation. While both approaches are viable, we believe that direct integration into a DFT code is the more userfriendly alternative. We have integrated TDDFT+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 TDDFT 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 orbitaldependent 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 TDDFTB, TDDFT+TB, sTDA and sTDDFT. 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 TDDFT+TB to be very suitable for the calculation and analysis of absorption spectra, while sTDDFT 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 nonexpert users.
Acknowledgements.
The research leading to these results has received funding from the European Union’s Seventh Framework Programme (FP7PEOPLE2012ITN) under project PROPAGATE, GA 316897.References
 Hohenberg and Kohn (1964) P. Hohenberg and W. Kohn, Phys. Rev. 136, B864 (1964).
 Kohn and Sham (1965) W. Kohn and L. J. Sham, Phys. Rev. 140, A1133 (1965).
 Porezag et al. (1995) D. Porezag, T. Frauenheim, T. Köhler, G. Seifert, and R. Kaschner, Phys. Rev. B 51, 12947 (1995).
 Seifert, Porezag, and Frauenheim (1996) G. Seifert, D. Porezag, and T. Frauenheim, Int. J. Quantum Chem. 58, 185 (1996).
 Slater and Koster (1954) J. C. Slater and G. F. Koster, Phys. Rev. 94, 1498 (1954).
 Elstner et al. (1998) M. Elstner, D. Porezag, G. Jungnickel, J. Elsner, M. Haugk, T. Frauenheim, S. Suhai, and G. Seifert, Phys. Rev. B 58, 7260 (1998).
 Gaus, Cui, and Elstner (2011) M. Gaus, Q. Cui, and M. Elstner, J. Chem. Theory Comput. 7, 931 (2011).
 Gaus, Goez, and Elstner (2013) M. Gaus, A. Goez, and M. Elstner, J. Chem. Theory Comput. 9, 338 (2013).
 Gaus et al. (2014) M. Gaus, X. Lu, M. Elstner, and Q. Cui, J. Chem. Theory Comput. 10, 1518 (2014).
 Lu et al. (2015) X. Lu, M. Gaus, M. Elstner, and Q. Cui, J. Phys. Chem. B 119, 1062 (2015).
 Kubillus et al. (2015) M. Kubillus, T. Kubař, M. Gaus, J. Řezáč, and M. Elstner, J. Chem. Theory Comput. 11, 332 (2015).
 Wahiduzzaman et al. (2013) M. Wahiduzzaman, A. F. Oliveira, P. Philipsen, L. Zhechkov, E. van Lenthe, H. A. Witek, and T. Heine, J. Chem. Theory Comput. 9, 4006 (2013).
 Oliveira, Philipsen, and Heine (2015) A. F. Oliveira, P. Philipsen, and T. Heine, J. Chem. Theory Comput. 11, 5209 (2015).
 Runge and Gross (1984) E. Runge and E. K. U. Gross, Phys. Rev. Lett. 52, 997 (1984).
 Casida (1995) M. E. Casida, in Recent Advances in Density Functional Methods (1995) Chap. 5, pp. 155–192.
 Casida and HuixRotllant (2012) M. E. Casida and M. HuixRotllant, Annu. Rev. Phys. Chem. 63, 287 (2012).
 Casida and Wesołowski (2003) M. E. Casida and T. A. Wesołowski, Int. J. Quantum Chem. 96, 577 (2003).
 Neugebauer (2007) J. Neugebauer, J. Chem. Phys. 126, 134116 (2007).
 Hirata and HeadGordon (1999) S. Hirata and M. HeadGordon, Chem. Phys. Lett. 314, 291 (1999).
 Grimme (2013) S. Grimme, J. Chem. Phys. 138, 244104 (2013).
 Risthaus, Hansen, and Grimme (2014) T. Risthaus, A. Hansen, and S. Grimme, Phys. Chem. Chem. Phys. 16, 14408 (2014).
 Rüger et al. (2015) R. Rüger, E. van Lenthe, Y. Lu, J. Frenzel, T. Heine, and L. Visscher, J. Chem. Theory Comput. 11, 157 (2015).
 Bannwarth and Grimme (2014) C. Bannwarth and S. Grimme, Comput. Theor. Chem. 1040–1041, 45 (2014).
 Niehaus et al. (2001a) T. A. Niehaus, S. Suhai, F. Della Sala, P. Lugli, M. Elstner, G. Seifert, and T. Frauenheim, Phys. Rev. B 63, 085108 (2001a).
 Domínguez et al. (2013) A. Domínguez, B. Aradi, T. Frauenheim, V. Lutsker, and T. A. Niehaus, J. Chem. Theory Comput. 9, 4901 (2013).
 Trani et al. (2011) F. Trani, G. Scalmani, G. Zheng, I. Carnimeo, M. J. Frisch, and V. Barone, J. Chem. Theory Comput. 7, 3304 (2011).
 Joswig et al. (2003) J.O. Joswig, G. Seifert, T. A. Niehaus, and M. Springborg, J. Phys. Chem. B 107, 2897 (2003).
 Goswami et al. (2006) B. Goswami, S. Pal, P. Sarkar, G. Seifert, and M. Springborg, Phys. Rev. B 73, 205312 (2006).
 Frenzel, Joswig, and Seifert (2007) J. Frenzel, J.O. Joswig, and G. Seifert, J. Phys. Chem. C 111, 10761 (2007).
 Li et al. (2007) Q. S. Li, R. Q. Zhang, T. A. Niehaus, T. Frauenheim, and S. T. Lee, J. Chem. Theory Comput. 3, 1518 (2007).
 Wang et al. (2007a) X. Wang, Zhang, Niehaus, and T. Frauenheim, J. Phys. Chem. C 111, 2394 (2007a).
 Wang et al. (2007b) X. Wang, R. Q. Zhang, S. T. Lee, T. A. Niehaus, and T. Frauenheim, Appl. Phys. Lett. 90, 123116 (2007b).
 Li et al. (2008) Q. S. Li, R. Q. Zhang, S. T. Lee, T. A. Niehaus, and T. Frauenheim, J. Chem. Phys. 128, 244714 (2008).
 Mitrić et al. (2009) R. Mitrić, U. Werner, M. Wohlgemuth, G. Seifert, and V. BonačićKoutecký, J. Phys. Chem. A 113, 12700 (2009).
 Zhang et al. (2012) R.Q. Zhang, A. De Sarkar, T. A. Niehaus, and T. Frauenheim, Phys. Status Solidi B 249, 401 (2012).
 Fan et al. (2014) G.H. Fan, X. Li, J.Y. Liu, and G.Z. He, Comp. Theor. Chem. 1030, 17 (2014).
 Niehaus (2009) T. A. Niehaus, J. Mol. Struc.: THEOCHEM 914, 38 (2009).
 Becke (1993) A. D. Becke, J. Chem. Phys. 98, 1372 (1993).
 Baerends, Gritsenko, and van Meer (2013) E. J. Baerends, O. V. Gritsenko, and R. van Meer, Phys. Chem. Chem. Phys. 15, 16408 (2013).
 van Meer, Gritsenko, and Baerends (2014) R. van Meer, O. V. Gritsenko, and E. J. Baerends, J. Chem. Theory Comput. 10, 4432 (2014).
 Mulliken (1939) R. S. Mulliken, J. Chem. Phys. 7, 14 (1939).
 Marques et al. (2012) M. A. Marques, N. T. Maitra, F. M. Nogueira, E. Gross, and A. Rubio, eds., Fundamentals of TimeDependent Density Functional Theory (Springer Berlin Heidelberg, 2012).
 van Gisbergen, Snijders, and Baerends (1999) S. J. A. van Gisbergen, J. Snijders, and E. J. Baerends, Comput. Phys. Commun. 118, 119 (1999).
 Mulliken (1955) R. S. Mulliken, J. Chem. Phys. 23, 1833 (1955).
 Mineva and Heine (2004) T. Mineva and T. Heine, J. Phys. Chem. A 108, 11086 (2004).
 Mineva and Heine (2006) T. Mineva and T. Heine, Int. J. Quantum Chem. 106, 1396 (2006).
 Ghosh and Islam (2009) D. C. Ghosh and N. Islam, Int. J. Quantum Chem. 110, 1206 (2009).
 Löwdin (1950) P.O. Löwdin, J. Chem. Phys. 18, 365 (1950).
 Heine et al. (1999) T. Heine, G. Seifert, P. W. Fowler, and F. Zerbetto, J. Phys. Chem. A 103, 8738 (1999).
 Gritsenko and Baerends (2004) O. Gritsenko and E. J. Baerends, J. Chem. Phys. 121, 655 (2004).
 Bredas (2014) J.L. Bredas, Mater. Horiz. 1, 17 (2014).
 Iikura et al. (2001) H. Iikura, T. Tsuneda, T. Yanai, and K. Hirao, J. Chem. Phys. 115, 3540 (2001).
 Henderson, Janesko, and Scuseria (2008) T. M. Henderson, B. G. Janesko, and G. E. Scuseria, J. Phys. Chem. A 112, 12530 (2008).
 Baer, Livshits, and Salzner (2010) R. Baer, E. Livshits, and U. Salzner, Annual Review of Physical Chemistry 61, 85 (2010).
 Zhang and Musgrave (2007) G. Zhang and C. B. Musgrave, J. Phys. Chem. A 111, 1554 (2007).
 Schreiber et al. (2008) M. Schreiber, M. R. SilvaJunior, S. P. A. Sauer, and W. Thiel, J. Chem. Phys. 128, 134110 (2008).
 Niehaus et al. (2001b) T. A. Niehaus, M. Elstner, T. Frauenheim, and S. Suhai, J. Mol. Struc. (THEOCHEM) 541, 185 (2001b).
 Yang et al. (2008) Y. Yang, H. Yu, D. York, M. Elstner, and Q. Cui, J. Chem. Theory Comput. 4, 2067 (2008).
 Perdew, Burke, and Ernzerhof (1996) J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996).
 HuixRotllant et al. (2011) M. HuixRotllant, A. Ipatov, A. Rubio, and M. E. Casida, Chem. Phys. 391, 120 (2011).
 Casida and HuixRotllant (2015) M. E. Casida and M. HuixRotllant, in DensityFunctional Methods for Excited States (Springer Science + Business Media, 2015) pp. 1–60.
 Casida et al. (1998) M. E. Casida, C. Jamorski, K. C. Casida, and D. R. Salahub, J. Chem. Phys. 108, 4439 (1998).
 Yasumatsu et al. (1996) H. Yasumatsu, T. Kondow, H. Kitagawa, K. Tabayashi, and K. Shobatake, J. Chem. Phys. 104, 899 (1996).
 Abrefah et al. (1992) J. Abrefah, D. R. Olander, M. Balooch, and W. J. Siekhaus, Appl. Phys. Lett. 60, 1313 (1992).
 Strain, Thomas, and Katz (1963) H. H. Strain, M. R. Thomas, and J. J. Katz, Biochimica et Biophysica Acta 75, 306 (1963).
 Du et al. (1998) H. Du, R.C. A. Fuh, J. Li, L. A. Corkan, and J. S. Lindsey, Photochemistry and Photobiology 68, 141 (1998).
 Baldo et al. (1999) M. A. Baldo, S. Lamansky, P. E. Burrows, M. E. Thompson, and S. R. Forrest, Appl. Phys. Lett. 75, 4 (1999).
 Fine et al. (2012) J. Fine, K. Diri, A. Krylov, C. Nemirow, Z. Lu, and C. Wittig, Molecular Physics 110, 1849 (2012).
 te Velde et al. (2001) G. te Velde, F. M. Bickelhaupt, E. J. Baerends, C. Fonseca Guerra, S. J. A. van Gisbergen, J. G. Snijders, and T. Ziegler, J. Comput. Chem. 22, 931 (2001).
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 