Self-Consistent Optimization of Excited States within Density

Nov 20, 2015 - Department of Chemistry, Advanced Materials Science and Engineering Center, and Institute for Energy Studies, Western Washington Univer...
1 downloads 13 Views 958KB Size
Article pubs.acs.org/JCTC

Self-Consistent Optimization of Excited States within DensityFunctional Tight-Binding Tim Kowalczyk,*,†,‡ Khoa Le,‡ and Stephan Irle† †

Institute of Transformative Bio-Molecules (WPI-ITbM) and Department of Chemistry, Graduate School of Science, Nagoya University, Nagoya 464-8602, Japan ‡ Department of Chemistry, Advanced Materials Science and Engineering Center, and Institute for Energy Studies, Western Washington University, Bellingham, Washington 98225, United States S Supporting Information *

ABSTRACT: We present an implementation of energies and gradients for the ΔDFTB method, an analogue of Δ-self-consistent-field density functional theory (ΔSCF) within density-functional tight-binding, for the lowest singlet excited state of closed-shell molecules. Benchmarks of ΔDFTB excitation energies, optimized geometries, Stokes shifts, and vibrational frequencies reveal that ΔDFTB provides a qualitatively correct description of changes in molecular geometries and vibrational frequencies due to excited-state relaxation. The accuracy of ΔDFTB Stokes shifts is comparable to that of ΔSCF-DFT, and ΔDFTB performs similarly to ΔSCF with the PBE functional for vertical excitation energies of larger chromophores where the need for efficient excited-state methods is most urgent. We provide some justification for the use of an excited-state reference density in the DFTB expansion of the electronic energy and demonstrate that ΔDFTB preserves many of the properties of its parent ΔSCF approach. This implementation fills an important gap in the extended framework of DFTB, where access to excited states has been limited to the time-dependent linear-response approach, and affords access to rapid exploration of a valuable class of excited-state potential energy surfaces.



INTRODUCTION A variety of critical energy conversion phenomena are mediated at the molecular level by electronically excited states. Photoabsorption and charge separation in organic photovoltaics,1,2 photoprotection of biological tissue through internal conversion of pigments such as melanin,3 and the chemical storage of solar energy by photosynthetic4 or artificial means5−7 all proceed through vertical electronic excitation to a low-lying excited state, followed by evolution on the excited-state potential energy surface (PES). For modeling the real-time dynamics of excited-state processes, nonadiabatic molecular dynamics studies, often based on fewest switches surface hopping,8−11 have proliferated and represent significant progress. The success of these models hinges on the ability to efficiently explore both ground- and excited-state PES, including identification and characterization of local minima, transition states, and conical intersections.12 The practitioner’s options for characterization of the groundstate PES span a wide range of methods, reflecting the central significance of the ground-state electronic structure problem. From wave function-based correlation methods and approximate density functional theory (DFT) through more approximate but more computationally efficient semiempirical methods and molecular mechanics force fields, each set of approximations reflects a different balance of accuracy against computational cost. Semiempirical methods such as densityfunctional tight-binding (DFTB), while typically less accurate © XXXX American Chemical Society

for energies than ab initio and DFT approaches, nevertheless often deliver sufficiently accurate geometries13−15 to make them valuable tools for rapid sampling of the ground-state PES. Options for characterizing excited-state PES not only are fewer but also are typically both less accurate and more computationally expensive than their ground-state counterparts. For simulating processes that depend on the interplay of ground and excited states, theoretical models with the following traits are especially desirable: (1) Even-handed treatment of ground and excited states from a unified set of approximations.16−18 (2) Efficient evaluation of ground- and excited-state analytical gradients.19−21 (3) Preservation of accuracy away from equilibrium regions of the PES. The final trait is especially important for photochemical processes involving bond-breaking and -forming events such as excited-state proton transfer and photodegradation. Within the Kohn−Sham (KS) DFT framework, the linearresponse time-dependent DFT (TD-DFT) description of excited states22 is effective for medium-to-large molecules near their equilibrium structures. Analytical gradients are available for some commonly employed exchange-correlation Received: August 2, 2015

A

DOI: 10.1021/acs.jctc.5b00734 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation

Below, we briefly review essential features of the selfconsistent-charge density-functional tight-binding (SCCDFTB) method for ground states and discuss its extension to electronically excited states. Particular attention is devoted to the choice of reference density in the ΔDFTB formalism, for which we provide an empirical justification of the choice made in this work. Singlet vertical excitation energies from ΔDFTB are then compared to those from related methods and to reference values for a variety of organic chromophores. Taking advantage of the availability of analytical gradients, we also calculate optimized geometries, emission energies, Stokes shifts, and excited-state vibrational frequencies for representative test sets of organic molecules. Finally, we revisit the previously reported ΔDFTB calculation to resolve a question posed by the earlier study before summarizing our main conclusions and plans for future work.

functionals, but current realizations of TD-DFT require modifications such as range-separation23 or a spin-flip reference24,25 in order to correctly describe certain chargetransfer excitations26 and conical intersections,27 respectively. Some of these pathologies can alternatively be corrected by a multireference approach such as Grimme’s DFT/multireference configuration interaction (DFT/MRCI)28 or constrained DFT configuration interaction.27 The direct self-consistent-field (SCF) optimization of excited states in DFT (ΔSCF-DFT or ΔSCF) readily affords analytical gradients and treats ground and excited states on similar footing.29 It has proven to be successful in describing excitedstate PES topologies, showing a robustness against qualitatively incorrect crossings between different excited-state PES relative to TD-DFT. 30 The restricted open-shell Kohn−Sham (ROKS)31 and restricted ensemble Kohn−Sham (REKS)32 methods are similar in approach but diverge in their treatment of spin symmetry. An elegant approximation to the timeindependent excited-state DFT theory of Ayers, Levy, and Nagy33 was recently reported by Evangelista and co-workers.34 Their orthogonality-constrained DFT method (OCDFT) further elucidates approximations built into ΔSCF while introducing a more rigorous approach to spin-adaptation via orthogonality constraints. The constricted variational DFT approach of Ziegler and co-workers35 also provides a theoretical bridge between TD-DFT and ΔSCF.36 Semiempirical approaches to excited states include truncated configuration interaction on top of semiempirical molecular orbital methods, with the INDO/X model 37 and the orthogonalization-corrected OMx family of methods38 standing out as especially promising for valence excitation energies of organic molecules. Excited states within DFTB, in keeping with the method’s origins in DFT, can be calculated from a linearresponse time-dependent approach (TD-DFTB).39,40 Significant investment in the development of TD-DFTB has led to efficient implementations and extensions41,42 with accuracy approaching that of TD-DFT in limited circumstances.43 However, TD-DFTB inherits some of the pathologies of TDDFT mentioned above,43,44 limiting its utility in simulations requiring a robust description of charge-transfer excited states or conical intersections. One may ask whether a method analogous to ΔSCF-DFT, implemented in a DFTB framework, might offer certain advantages over TD-DFTB as has been demonstrated in the case of excited-state DFT. We are aware of only one previous, isolated example of such a ΔSCF-DFTB calculation, hereafter abbreviated ΔDFTB, in the literature. Elstner and co-workers compared excited-state PES scans along an isomerization coordinate of a protonated Schiff base model described by TD-DFTB and by a ΔDFTB approach.45 The excitation energy was underestimated relative to TD-DFTB and showed a different topology at the geometry of the S1/S0 conical intersection. These calculations did not account for geometry relaxation in the excited state. Given the growing need for fast and reliable excited-state molecular dynamics simulations, an assessment of the properties and performance of the ΔSCF approach within DFTB should be both valuable and timely. In this article, we examine the adaptation of the ΔSCF-DFT approach to the DFTB framework, analyzing and benchmarking its description of lowlying excited-state PES relative to alternative theoretical approaches and to experimental reference data.



THEORY ΔSCF Excited States in DFT. The self-consistent calculation of electronically excited states within KS DFT by the ΔSCF method predates the linear-response TD-DFT approach22,46−48 and has been a subject of renewed interest in recent years. Paired with modern exchange-correlation functionals, it has proven to be useful for the prediction of core,49 valence,29 and Rydberg excitation energies,50 including for molecules adsorbed on surfaces,51 with an accuracy often rivaling that of adiabatic TD-DFT when calculated with the same exchange-correlation functional.29,30 Here, we review the key ideas underpinning ΔSCF excited states to show how they can be adopted within the DFTB framework. The central premise of ΔSCF is the interpretation of the density constructed from KS spin-orbitals {ϕσi (r)}, which have been self-consistently optimized under the constraint of a particular non-Aufbau occupation pattern, occ ≠ {1, ..., N}, as an excited-state density, ρe (r) =

∑ ∑ σ

|ϕi σ (r)|2 (1)

i ∈ occ(σ )

Here, i indexes the spatial orbitals with spin label σ, and the occupied set of KS orbitals is selected according to the character of the excited state of interest. It is common for the lowest-lying singlet excited state S1 in organic materials to be dominated by a one-electron transition from the highest occupied molecular orbital (HOMO) to the lowest unoccupied molecular orbital (LUMO), corresponding to the occupation pattern occ(α) = {1, ..., N − 1, N + 1}

(2)

occ(β) = {1, ..., N }

(3)

We therefore emphasize excited states of HOMO → LUMO character in what follows; generalizations to other excited states29,36 are possible and worthy of exploration in future work. Excited-state densities obtained from a single set of ΔSCF KS orbitals {ϕσi }m are stationary densities within adiabatic TDDFT;29 however, they are also heavily spin-contaminated due to the arbitrary assignment of the transition to α (or β) orbitals. The subscript m on the set of orbitals is introduced to indicate that these orbitals constitute such a mixed-spin state. The effects of spin contamination can be alleviated at the level of the Hamiltonian (as in ROKS)31,52 or at the level of the energy through the Ziegler sum rule.34,46 The latter approach uses the B

DOI: 10.1021/acs.jctc.5b00734 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation energy difference between the spin-contaminated energy and the lowest triplet state energy (with its own set of optimized KS orbitals {ϕ̃ σi }t) to approximate the energy of the spin-pure S1 state,

EW [δm] =

(4)

E[ρ] ≈ E[ρ0 ] ⎡

2



δ E[ρ] ⎥ ∫ d3r′ ∫ d3r″⎢⎣ δρ ″δρ′ ⎦

δρ″δρ′ + O(N3) ρ″ 0 ρ′ 0

(5)

where ρ′ = ρ(r′) and ρ″ = ρ(r″). In the spin-polarized formalism, truncating the expansion at second order and introducing a set of confined Slater-type atomic orbitals lead to a tight-binding-like energy expression with additional secondorder terms Eγ and EW E SCC‐DFTB[ρ] = E 0[ρ0 ] + E rep[ρ0 ] + E γ [δρ] + EW [δm] (6)

The first term E is the sum over occupied orbitals of the corresponding eigenvalues of the Hamiltonian matrix constructed from ρ0 in a two-center approximation.59 The repulsive energy term Erep approximates the remaining Coulomb and exchange-correlation contributions to the total energy (sometimes referred to collectively as DFT double-counting terms). This term is approximated as a sum over pairs of short-ranged, atom-centered potentials that have been fit to reference DFT calculations.60 The second-order term Eγ accounts for charge transfer between any two atomic centers A and B,53 0

E γ [δρ] =

1 2

∑ δqA δqBγAB AB

A

ll ′

(8)

The atomic spin constants employed in this work are obtained from ref 63 and are also provided in the Supporting Information. SCC-DFTB employs a minimal set of Slater-type basis functions. All required integrals over these functions are precomputed and tabulated, resulting in substantial computational savings over KS DFT. The successes of ground-state SCC-DFTB and its third-order extension DFTB364,65 offer ample empirical evidence that the SAD reference density is a suitable zeroth-order density for the expansion in eq 5. The treatment of singlet excited states in ΔDFTB proceeds as follows. A HOMO → LUMO promotion in the α block of spin-orbitals affects both E0 and Erep by replacing the SAD reference density ρ0 with an excited-state reference density ρ0e constructed from the non-Aufbau set of occupied orbitals (eq 1). The second-order terms are calculated exactly as for the spin-unrestricted ground state, with the reference charges q0a now deriving from the excited-state reference density ρ0e . As in ground-state SCC-DFTB, iterative construction and diagonalization of the Hamiltonian relaxes the orbitals and density to charge self-consistency. A separate DFTB calculation of the triplet state followed by application of the Ziegler sum rule corrects the previously obtained mixed-spin ΔDFTB energy to obtain the energy of the singlet excited state, in direct analogy to the ΔSCF procedure in DFT. Reference Density for ΔDFTB. The centrality of the Taylor expansion to the SCC-DFTB formalism leads to a critical question about underlying assumptions in any selfconsistent DFTB treatment of excited states: Do the SAD reference density ρ0 and its excited-state counterpart ρ0e provide a suitable zeroth-order approximation for the determination of self-consistent excited states within DFTB? Given that ρ0 is constructed from neutral atomic densities in their ground electronic states, we might anticipate that ρ0 is more similar to the ground-state density ρg than to the density ρe of any bound excited state. Evidence from excited-state DFT suggests that ρ0e should nevertheless be a satisfactory zeroth-order approximation to ρe. In ΔSCF-DFT, non-Aufbau occupation of the ground-state optimized orbitals, yielding an initial guess density ρ*, is often sufficient to achieve SCF convergence to the target excited state. Since ρ0g and ρg in DFTB are similar enough to justify eq 5, we anticipate that ρ0e will also be similar to the density ρ* obtained from direct non-Aufbau occupation of the groundstate optimized orbitals in DFTB. Thus, lessons learned from ΔSCF-DFT suggest reason for optimism. A more definitive answer regarding the suitability of the Taylor expansion for ΔDFTB is among the goals of this investigation.

When calculated with standard hybrid exchange-correlation functionals, spin-purified ΔSCF vertical excitation energies of closed-shell organic molecules are roughly as accurate, on average, as those obtained from TD-DFT, although excitation energies from the two approaches may differ substantially from one another for a given molecule.29 Adaptation of the ΔSCF Approach to DFTB. As a semiempirical approximation derived from KS DFT, the ground-state SCC-DFTB model retains all features of KS DFT necessary for the formulation of a ΔDFTB approach to excited states. The formalism and working equations of SCCDFTB are well-documented and have been extensively reviewed;53−57 here, we endeavor to provide only the minimum of detail required to establish the approximations that define ΔDFTB excited states. At the heart of SCC-DFTB and its extensions is the assumption that a reference density ρ0, constructed as a superposition of atomic densities (SAD), provides a useful zeroth-order approximation to the true ground-state density ρ. The DFT energy can then be expanded in a Taylor series in fluctuations of the density away from ρ0,53,58

+

∑ ∑ δpAl δpAl′ W All′ WllA′

σ

E(S1) = 2E[{ϕi σ }m] − E[{ϕĩ }t ]

1 2



COMPUTATIONAL DETAILS For ease of comparison, vertical excitation energies from ΔDFTB, ΔSCF, and TD-DFT are calculated at the reference geometries provided in the corresponding vertical excitation energy benchmark data sets. For benchmark calculations that include geometry optimization, ΔDFTB excitation energies are reported at DFTB ground- and excited-state optimized geometries, whereas TD-DFT and ΔSCF excitation energies are calculated at the corresponding DFT ground- and excitedstate optimized geometries. Excited-state DFT calculations use the semilocal exchange-correlation functional PBE66 or the

(7)

where the net atomic charges δqA = qA − q0A are taken in practice as Mulliken charges calculated from the density and γAB is an effective electron−electron interaction potential. The other second-order term in eq 6 accounts for spin polarization of the SCC density, δm = δρα − δρβ, through the Mulliken spin 61,62 l,β populations δplA = δql,α A − δqA of each shell, C

DOI: 10.1021/acs.jctc.5b00734 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation hybrid exchange-correlation functionals B3LYP67 or PBE0,68,69 as indicated below. The 6-31G* basis set is used for all DFT calculations unless otherwise specified. ΔDFTB gradients with respect to nuclear coordinates X are calculated analytically via differentiation of eq 4, σ

∇X E(S1) = 2∇X E[{ϕi σ }m] − ∇X E[{ϕĩ }t ]

(9)

The ΔDFTB Hessian is evaluated via numerical differentiation of eq 9 for the computation of excited-state vibrational frequencies. ΔDFTB energies and gradients were implemented in a development version of the DFTB+ program package.70 All DFTB and ΔDFTB calculations reported below use the mio parameter set.53 For Stokes shifts, we also provide DFTB3 energies calculated with the 3OB parameter set64,71 in the Supporting Information, but the results obtained with the mio parameters are emphasized here to ensure an even-handed comparison of ΔDFTB and TD-DFTB energies. Ground-state DFTB calculations were performed in version 1.2 of DFTB+. All reference ground- and excited-state DFT calculations were performed with version 4.3 of the Q-Chem program package.72

Figure 1. (a) Illustration identifying the locations of hypothetical SAD (blue solid), ground-state DFTB (green dashed), and excited-state ΔDFTB (red dashed) densities in the density metric space of D’Amico and co-workers.73 (b, c) Calculated distances between reference (ρ0g), excited-state reference (ρ0e ), DFTB (ρg), and ΔDFTB (ρe) densities for formaldehyde at the DFTB ground-state-optimized geometry.



RESULTS AND DISCUSSION Suitability of the Reference Density. The SCC-DFTB energy expression (eq 5) is defined in terms of an expansion about a reference density for which a second-order correction is determined self-consistently. A variational approach to excited states like ΔDFTB likewise requires a reference density that is suitably similar to the SCC density for the expansion to be valid and useful. To provide some evidence that this condition is satisfied in ΔDFTB, it would be preferable to have a way of measuring the similarity of the reference density to the SCC density for both ground and excited states. Such a tool would allow us to determine, for example, whether the similarity between the ground-state reference density ρ0g and its SCC density ρg is stronger or weaker than the similarity between the excited-state reference density ρ0e and its SCC density ρe. A suitable measure of the similarity between different densities was recently introduced by Capelle and co-workers,73 who define a distance Dρ between any two electron densities ρ1 and ρ2 as Dρ(ρ1 , ρ2 ) =

∫ d3r|ρ2 (r) − ρ1(r)|

In the case of formaldehyde, Figure 1b shows that the SAD reference density ρ0g differs from the ground-state DFTB density ρg by 1.42 electrons according to the Capelle metric. The excited-state density ρe is a similar distance from ρ0g (1.36 electrons) and from ρg (1.43 electrons), reflecting its different orbital occupation pattern. The distance from the excited-state reference density to the excited-state density, Dρ(ρ0e , ρe) = 1.50, is only slightly larger than Dρ(ρ0g, ρg) and is similar to Dρ(ρ0g, ρ0e ). For our purposes, the precise numerical difference between these distances is inconsequential; the important conclusion is that if ρ0g is a suitable reference density for the calculation of ρg in DFTB, then ρ0e is also a suitable reference density for the calculation of ρe. This analysis was applied to several other small organic dyes (see Figure 3 and Supporting Information) with the same conclusion that the excited-state reference density is suitable for ΔDFTB. Excitation Energies. To characterize the predictions of ΔDFTB for vertical excitation energies, we calculated DFTB and ΔDFTB energies across the benchmark test set of 28 small organic molecules for which high-level reference calculations were reported by Thiel and co-workers.74 Only the lowest singlet excitation S1 ← S0 in each molecule is considered. As anticipated from the use of non-Aufbau occupations, the SCC equations of ΔDFTB can be difficult to converge in situations where frontier orbitals become degenerate. Consequently, for this benchmark data set, we exclude 4 of the 28 molecules in the full test set due to convergence challenges. Two of the excluded molecules (benzene and s-triazine) have both a degenerate HOMO and LUMO, whereas the other excluded molecules (thymine and uracil) exhibit orbital swapping during SCC iterations in the excited state. Such cases could be resolved with specialized convergence strategies like the maximum overlap method;75 here, we focus on the quality of converged ΔDFTB energies.

(10)

This definition equips the space of all densities with a valid metric, e.g., satisfying a triangle inequality, which permits a consistent measurement of distance between any two densities. N-electron densities are located on a hypersphere of radius N, so that two densities with the same N must have a distance within the range 0 ≤ Dρ ≤ 2N. Hypothetical locations of ρ0g, ρg, and ρe on a hypersphere of radius N are illustrated in Figure 1a. We measured Dρ between the four densities (ρ0g, ρg, ρ0e , ρe) for the case of a formaldehyde molecule in its ground-stateoptimized geometry. Real-space densities for the evaluation of eq 10 were obtained in cube file format from DFTB+ using the Waveplot program for several different grid sizes to ensure that the numerical integration of Dρ is converged. The results are portrayed in Figure 1b,c. Note that the space of all possible densities is of infinite dimension; consequently, there is no requirement that the triangles in Figure 1b,c coexist in the same plane. D

DOI: 10.1021/acs.jctc.5b00734 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation The performance of ΔDFTB relative to theoretical best estimates for the remaining 24 molecules is compared in Table 1 with that of TD-DFTB and TD-DFT with the pure PBE and Table 1. Deviation of Vertical Excitation Energies from Theoretical Best Estimates (in eV) for the Modified Thiel Test Set method

mean error

mean absolute error

RMS error

TD-DFTBa ΔDFTB (mixed) ΔDFTB (pure) TD-PBEb TD-PBE0b

−0.14 −0.75 −0.45 −0.41 −0.01

0.33 0.79 0.64 0.45 0.19

0.41 1.04 0.82 0.51 0.23

a

TD-DFTB excitation energies obtained from ref 43. bLinear response TD-DFT with the PBE and PBE0 functionals and the 6-31G* basis set.

hybrid PBE0 functionals. Both TD-DFTB and ΔDFTB tend to underestimate the excitation energies. ΔDFTB without spin purification, denoted “mixed” in Table 1, makes the most severe underestimation with a mean error (ME) of −0.75 eV; the spin purification correction reduces the size of the error, but it remains larger than that observed for TD-DFTB. While neither DFTB approach matches the performance of TD-DFT, the similarity of the mean absolute error (MAE) to the magnitude of the ME for both ΔDFTB methods suggests that the ΔDFTB error is mostly a systematic underestimation of the vertical excitation energy rather than random error. In that case, the shape of the ΔDFTB PES, and hence of excited-state geometries, may still be robust. This possibility will be explored in greater detail in the discussion of geometry optimization below. Figure 2 illustrates the deviations of ΔDFTB, TD-DFTB, and TD-DFT from theoretical best estimates for two subclasses of molecules in the Thiel set, the carbonyl-containing species and the heterocycles possessing one to two heteroatoms. For the π* ← n transitions in carbonyl compounds, TD-DFTB slightly outperforms spin-pure ΔDFTB with the exception of the simplest case, formaldehyde. The carbonyl compounds include a few of the exceptional cases for which spin-mixed ΔDFTB overestimates the vertical excitation energy. In these cases, spin purification actually results in a poorer estimate of ΔE: for molecules with E(S1) > E(T1), spin purification via the Ziegler sum rule increases the predicted ΔE. Both TD-DFTB and ΔDFTB are especially inaccurate for benzoquinone, for which TD-PBE0/6-31G* also yields a substantial error. Given their mediocre performance over the carbonyl compounds, ΔDFTB and TD-DFTB likely benefit from some error cancellation to achieve the satisfyingly small errors for the amides in Figure 2a. Among the heterocycles, ΔDFTB severely underestimates the vertical excitation energies of furan and pyrrole, which are π* ← π transitions, but it performs well for the π* ← n transitions in the diazines. In contrast, the TD-DFTB error is small for furan and pyrrole but substantial for pyridine. The error histograms of Figure 2 show that TD-DFTB and ΔDFTB errors are only weakly correlated. Thus, one cannot assume that any particular failure of TD-DFTB precludes the success of ΔDFTB, and vice versa. To determine whether the quality of the excited-state reference density is affecting the accuracy of the excitation energies, we carried out the density distance analysis described in the previous section on examples for which ΔDFTB

Figure 2. Histograms showing the deviation of TD-DFTB (green), ΔDFTB in its spin-mixed (blue) and spin-pure (purple) forms, and TD-PBE0 (gold) from theoretical best estimates. (a) Molecules possessing the carbonyl group. (b) Heterocycles with up to two heteroatoms.

performs particularly poorly (benzoquinone) and particularly well (pyridazine). The density distances are shown in Figure 3. Although Dρ(ρ0e , ρe) is larger for benzoquinone than for pyridazine, the corresponding ground-state distance Dρ(ρ0g, ρg) is also larger for benzoquinone. We find empirically that the distance between the reference and SCC-optimized densities grows with increasing total electron number (see additional examples in the Supporting Information), so the ratio of Dρ(ρ0e , ρe) to Dρ(ρ0g, ρg) is a more meaningful gauge of the suitability of the excited-state reference density than Dρ(ρ0e , ρe) alone. Thus, the fact that Dρ(ρ0e , ρe) ≈ Dρ(ρ0g, ρg) for both systems indicates that the quality of the excited-state reference density is not at fault for the poor performance of ΔDFTB for benzoquinone. The previously documented performance of ΔSCF-DFT for larger organic chromophores29 suggests that ΔDFTB may also

Figure 3. Distances between reference and SCC densities for the ground and excited states of (a) benzoquinone (b) pyridazine according to the metric defined by eq 10. E

DOI: 10.1021/acs.jctc.5b00734 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation

ΔPBE which can be achieved with the use of a hybrid functional like PBE0, it is likely that ongoing developments in ground-state DFTB44,76 will directly yield improved performance for ΔDFTB as well. Finally, a comparison of the statistics in Tables 1 and 2 supports the notion that ΔDFTB, like its parent method ΔSCF, shows improved performance for the larger organic chromophores where its computational efficiency is the most valuable. Geometry Optimization and Vibrational Frequencies. Simple and efficient access to analytical gradients is a strong point of ΔDFTB because it enables fast geometry optimization in the S1 excited state. Here, we characterize the accuracy of ΔDFTB-optimized geometries by comparison to reference calculations from higher-level theories. Transitions involving the π systems of the carbonyl CO group and of the CC skeleton are of central importance to the photophysics of organic molecules. Accordingly, we selected acetone and methylenecyclopropene (MCP), shown in Figure 4, as simple representatives for these broad classes of excitations. For each molecule, the ground- and excited-state geometries were optimized with SCC-DFTB and with spinpure ΔDFTB, respectively. Figure 4 presents DFTB/ΔDFTB optimized geometries of the ground and S1 excited states, comparing key geometrical features against previously reported geometries from multireference perturbation theory (CASPT2) and variational Monte Carlo (VMC) calculations.77 For the S1 state of MCP, the geometry is constrained to remain planar, in keeping with the procedure followed in the reference calculations. Excited-state bond lengths and angles from ΔDFTB differ from the CASPT2 and VMC references more substantially than their DFTB counterparts in the ground state. However, the differences are only a fraction of the total geometry change induced by the S1 ← S0 transition. For example, the CO bond length in the S1 state of acetone is about 3 pm longer according to ΔDFTB than according to either reference method. This error is larger than the corresponding DFTB error in the ground state but amounts to only about 20% of the total change in bond length upon excitation. ΔDFTB correctly predicts the direction of every heavy-atom bond length and angle change from S0 to S1 in both species. Thus, to the extent that DFTB itself is an accurate descriptor of the ground state, ΔDFTB appears to provide a qualitatively accurate description of the S1-optimized geometry. The geometry changes of MCP in Figure 4c,d further support this conclusion: although the errors in ΔDFTB bond lengths for this molecule are larger than those of the DFTB ground state, ΔDFTB correctly predicts that the cyclopropenyl CC bond is longer than the ethenyl CC bond in the excited state.

perform more reliably for such systems. To test this hypothesis and to compare the effects of spin purification in ΔSCF and in ΔDFTB, we calculated ΔDFTB vertical excitation energies for the larger-dye test set compiled in ref 29. Two of the 16 dyes (zinc phthalocyanine and a chlorinated cyanine dye) are excluded from the test set because the elements Zn and Cl are not part of the standard mio parameter set; one other (rhodamine 6G) is excluded because an SCC-converged excitation energy could not be obtained. The excitation energies for the other 13 dyes are compared against predictions from ΔSCF at the PBE/6-311+G* and PBE0/6-311+G* levels in Table 2. Table 2. Deviation of Vertical Excitation Energies from Theoretical Best Estimates (in eV) for the Large Organic Dye Test Set method

mean error

mean absolute error

RMS error

ΔDFTB (mixed) ΔDFTB (pure) ΔPBE (mixed)a ΔPBE (pure)a ΔPBE0 (mixed)a ΔPBE0 (pure)a

−0.68 −0.59 −0.79 −0.61 −0.47 −0.07

0.68 0.59 0.79 0.64 0.51 0.24

0.73 0.64 0.87 0.71 0.56 0.30

Statistics for the ΔSCF (PBE/6-311+G* and PBE0/6-311+G*) excitation energies of ref 29 are re-evaluated over the 13-dye subset for which ΔDFTB statistics are reported.

a

As for the molecules in the Thiel test set, ΔDFTB systematically underestimates vertical excitation energies as indicated by the fact that the ME and MAE have the same magnitude. ΔDFTB appears to inherit this underestimation from the PBE functional against which it is parametrized; the ME and MAE for ΔPBE reveal that this method also systematically underestimates ΔE. This underestimation can be offset by spin purification and exact exchange as illustrated for ΔPBE0, but it is only partially offset in ΔDFTB, which does not account for exact exchange in the formalism employed here. Still, spin purification of ΔPBE and of ΔDFTB excitation energies reduces the MAE and root-mean-square error (RMSE) over the test set. A comparison of the RMSE of spin-purified ΔDFTB (0.64 eV) over this test set to the RMSE of spin-purified ΔPBE (0.71 eV) reveals that ΔDFTB achieves a similar accuracy to its parent method ΔPBE over this test set. This performance is noteworthy for several reasons. First, it is remarkable that ΔDFTB predicts vertical excitation energies for these larger organic dyes with negligible loss of accuracy relative to a DFTbased approach with significantly higher computational overhead. Furthermore, in light of the significant improvement over

Figure 4. Comparison of optimized bond lengths (Å) and angles (deg) for the S0 and S1 states of acetone and methylenecyclopropene predicted by CASPT2, VMC, and ΔDFTB. (a) Acetone, S0. (b) Acetone, S1 (Cs symmetry). (c) Methylenecyclopropene, S0. (d) Methylenecyclopropene, S1. CASPT2 and VMC reference data are from ref 77. F

DOI: 10.1021/acs.jctc.5b00734 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation

changes (Table S3 in the Supporting Information) indicates that both factors contribute. The smaller errors for ΔDFTB Stokes shifts relative to vertical excitation energies corroborates the notion that the ΔDFTB excited-state PES is globally shifted to lower energies but otherwise quite parallel to the true S1 PES for modest geometrical distortions. To test the prediction of molecular properties calculated from ΔDFTB geometries, we acquired excited-state Hessians for S1 states of acetone, acetamide, and aniline by numerical differentiation of the ΔDFTB gradient. Table 4 compares the frequencies of selected vibrational modes predicted by ΔDFTB/DFTB against reference TD-DFT/DFT results calculated at the B3LYP/6-31G* level.78 ΔDFTB correctly predicts the direction of each frequency shift, and even the relative magnitudes of the shifts are in reasonable agreement with TD-DFT. Given the high computational cost of evaluating the TD-DFT Hessian for large molecules and recent successes in applying ΔSCF to excited-state vibrations,79,80 these smallmolecule benchmarks are an encouraging indicator that ΔDFTB will be useful for computationally efficient simulation of excited-state IR spectra. Photoisomerization. The availability of excited-state gradients in ΔDFTB allows us to revisit to a question posed in the first report of a ΔDFTB energy calculation.45 This investigation also provides an opportunity to study a nontrivial ΔDFTB PES in greater detail. Elstner and co-workers examined excited state cis−trans isomerization in models of the protonated Schiff base (PSB) of retinal45 using TD-DFTB, TD-DFT, and CASPT2. Figure 5a depicts the geometry of the PSB of 6-s-cis retinal at the initial 11-cis conformation (bottom, dihedral angle φ = 0°) and at the twisted conformation (top, φ = 90°) corresponding to a S1/S0 conical intersection. The current state of the art in DFT models of PSB isomerization is analyzed in ref 81. In the study of Elstner and co-workers, qualitative differences were observed in the isomerization pathways on the S1 PES predicted by the linear-response approaches (TD-DFTB, TDDFT) and the CASPT2 reference. These differences were traced back to the well-documented difficulty of TD-DFT to describe excitations with charge-transfer character26 and to overpolarization of the excited-state density.82 ΔSCF methods were cited as relatively immune to these problems. However, a scan of the central dihedral angle φ governing cis−trans isomerization, with all other internal degrees of freedom held fixed at their ground-state optimized values, suggested an unphysical energy gap instead of a conical intersection between S0 and S1 at φ = 90°. It was suggested by reference to a planewave ΔSCF-DFT study83 that this energy gap would persist in a fully relaxed scan, but the issue could not be definitively

Excited-state geometry optimization also affords access to emission energies, which can be used to estimate the gas-phase Stokes shift. Table 3 reports deviations of the predicted Stokes Table 3. Stokes Shifts (in eV) Predicted by Spin-Purified ΔB3LYP and ΔDFTB for Selected Organic Dyesa

a

Deviations from experimental reference values31 are given in parentheses.

shift from experiment for selected large organic dyes used in a previous benchmark of ROKS excited states.31 Of the nine dyes in the original set, three were excluded due to the presence of atoms outside of the mio parameter set; a fourth (coronene) possesses degenerate frontier orbitals and was excluded on account of the concomitant SCC convergence issues discussed previously. As for vertical excitation energies, ΔDFTB captures the essential trends in Stokes shifts and, with the exception of indigo, correctly ranks the magnitude of the Stokes shift in these five dyes. The Stokes shifts are overestimated in most cases, but not significantly more so than for ΔB3LYP. The overestimation is largest for dyes with the larger Stokes shifts and could arise from overestimation of either geometry relaxation in the excited state or the energetic costs of distorting the geometry in ground-state DFTB to the excitedstate minimum. Inspection of these contributing energy

Table 4. Selected Equilibrium Bond Lengths Req and Vibrational Frequencies ν in the S0 and S1 States Predicted by DFT/TDDFT and DFTB/ΔDFTBa acetone CO S0 (DFTB) S0 (DFT) S1 (ΔDFTB) S1 (TD-DFT) shift (DFTB) shift (DFT) a

acetamide CO

acetamide C−N

aniline C−N

Req (Å)

ν (cm−1)

Req (Å)

ν (cm−1)

Req (Å)

ν (cm−1)

Req (Å)

ν (cm−1)

1.208 1.217 1.378 1.313 +0.170 +0.096

1775 1776 1218 1311 −557 −465

1.223 1.221 1.425 1.294 +0.202 +0.073

1733 1813 1281 1492 −452 −321

1.383 1.369 1.415 1.479 +0.032 +0.110

1338 1360 1421 1394 +83 +34

1.385 1.400 1.353 1.379 −0.032 −0.021

1410 1316 1406 1284 −4 −32

DFT/TD-DFT results were reported in ref 78. G

DOI: 10.1021/acs.jctc.5b00734 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation

Figure 5. (a) ΔDFTB-optimized geometries of the protonated Schiff base of 6s-cis retinal along the 11-cis−trans isomerization pathway. Bottom and top structures correspond to dihedral angles of 0° and 90°, respectively. (b−d) Scans of the S0 and S1 PES along the 11-cis−trans dihedral angle described by DFTB (green +) and ΔDFTB (purple ×), respectively: (b) Fixed scan, with other degrees of freedom held at the optimized S0 geometry; (c) relaxed scan along the S0 PES; (d) relaxed scan along the S1 PES.

suitable zeroth-order density for the analogous expansion in ground-state DFTB. ΔDFTB with the mio parameter set typically underestimates excitation energies in the examples surveyed here. The spinpurification correction shifts energies higher, but this is not usually sufficient to overcome the underestimation. Computation of excitation energies for several larger dyes revealed that ΔDFTB excitation energies are comparable in accuracy to ΔPBE but are not currently competitive with ΔSCF using hybrid or more sophisticated exchange-correlation functionals. Excited-state geometries from ΔDFTB were found to be qualitatively reliable, with correct predictions of bond length and angle changes in the model compounds acetone and methylenecyclopropene. Changes in vibrational frequencies upon relaxation to the S1 minimum-energy structure are also in qualitative agreement with TD-DFT. In keeping with the apparent quality of ΔDFTB geometries, ΔDFTB is competitive with ΔB3LYP for the prediction of Stokes shifts in the cases examined here. Revisiting the application of excited-state DFTB to isomerization of protonated Schiff bases, we demonstrated that relaxed PES scans along the isomerization dihedral angle in the S1 excited state preserve the previously observed spurious energy gap between S0 and S1 where experiments and multireference ab initio calculations identify a conical intersection. Many of the pathologies of ΔDFTB that were identified in this investigation could be mitigated with modest advances in methodology. For example, SCC convergence difficulties are generally of the same type that arise in ΔSCF, where the maximum overlap method (MOM)75 has proven to be valuable in achieving SCF convergence. An implementation of MOM in DFTB consistent with ΔDFTB is under development to enable calculations on systems with frontier orbital degeneracies. The S1/S0 conical intersection in the protonated Schiff base model is

resolved without access to excited-state geometry optimization in the S1 state along the dihedral angle. Relaxed scans of the S0 and S1 PES along the dihedral angle φ are presented in Figure 5c,d; a fixed scan analogous to that of ref 45 is also provided in Figure 5b because our calculations were performed on the PSB of the full 6-s-cis retinal in lieu of a model PSB. All of the scans, and in particular the S1-relaxed scan in Figure 5d, show a significant energy gap at φ = 90°, consistent with the prediction of ref 45. The gap increases by roughly 0.6 eV, from 0.5 to 1.1 eV, upon relaxation of all other degrees of freedom on the S1 PES. This finding, while in disagreement with experiment, is consistent with previous ΔSCF-DFT studies83 and provides further evidence that qualitative features of PES in ΔDFTB will generally reflect their ΔSCF-DFT counterparts.



CONCLUSIONS We have introduced a theoretical framework and working equations for a ΔSCF approach to excited states in DFTB, which we call ΔDFTB. Iterative solution the SCC-DFTB equations with a predetermined non-Aufbau occupation of the orbitals converts an excited-state reference density to an optimized excited-state DFTB density. The Ziegler sum rule is employed to offset the error in the energy due to spin contamination. Calculation of the ΔDFTB energy and analytical gradient have been implemented within the DFTB+ program package. We used a recently introduced measure of distances between electronic densities to demonstrate that the excited-state reference density generated by non-Aufbau orbital occupation in ΔDFTB is a suitable reference density for the Taylor expansion of the energy. The argument relies on the assumption, demonstrated empirically in the success of ground-state DFTB, that the SAD reference density is a H

DOI: 10.1021/acs.jctc.5b00734 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation

dissipation within natural and artificial bacteriochlorophyll proteins. Nat. Commun. 2014, 5, 5287. (6) Sekar, N.; Ramasamy, R. P. Recent advances in photosynthetic energy conversion. J. Photochem. Photobiol., C 2015, 22, 19−33. (7) Ganesh, I. Solar fuels vis-a-vis electricity generation from sunlight: The current state-of-the-art. Renewable Sustainable Energy Rev. 2015, 44, 904−932. (8) Tully, J. C. Molecular dynamics with electronic transitions. J. Chem. Phys. 1990, 93, 1061. (9) Nelson, T.; Fernandez-Alberti, S.; Roitberg, A. E.; Tretiak, S. Nonadiabatic Excited-State Molecular Dynamics: Modeling Photophysics in Organic Conjugated Materials. Acc. Chem. Res. 2014, 47, 1155−1164. (10) Cui, G.; Thiel, W. Generalized trajectory surface-hopping method for internal conversion and intersystem crossing. J. Chem. Phys. 2014, 141, 124101. (11) Petit, A. S.; Subotnik, J. E. How to calculate linear absorption spectra with lifetime broadening using fewest switches surface hopping trajectories: A simple generalization of ground-state Kubo theory. J. Chem. Phys. 2014, 141, 014107. (12) Subotnik, J. E.; Alguire, E. C.; Ou, Q.; Landry, B. R.; Fatehi, S. The Requisite Electronic Structure Theory To Describe Photoexcited Nonadiabatic Dynamics: Nonadiabatic Derivative Couplings and Diabatic Electronic Couplings. Acc. Chem. Res. 2015, 48, 1340−1350. (13) Zheng, G.; Irle, S.; Morokuma, K. Performance of the DFTB method in comparison to DFT and semiempirical methods for geometries and energies of C20-C86 fullerene isomers. Chem. Phys. Lett. 2005, 412, 210−216. (14) Yilmazer, N. D.; Korth, M. Comparison of Molecular Mechanics, Semi-Empirical Quantum Mechanical, and Density Functional Theory Methods for Scoring Protein-Ligand Interactions. J. Phys. Chem. B 2013, 117, 8075−8084. (15) Cui, Q.; Elstner, M. Density functional tight binding: values of semi-empirical methods in an ab initio era. Phys. Chem. Chem. Phys. 2014, 16, 14368−14377. (16) Yost, S. R.; Kowalczyk, T.; Van Voorhis, T. A multireference perturbation method using non-orthogonal Hartree-Fock determinants for ground and excited states. J. Chem. Phys. 2013, 139, 174104. (17) Sundstrom, E. J.; Head-Gordon, M. Non-orthogonal configuration interaction for the calculation of multielectron excited states. J. Chem. Phys. 2014, 140, 114103. (18) Liu, X.; Subotnik, J. E. The Variationally Orbital-Adapted Configuration Interaction Singles (VOA-CIS) Approach to Electronically Excited States. J. Chem. Theory Comput. 2014, 10, 1004−1020. (19) Chiba, M.; Tsuneda, T.; Hirao, K. Excited state geometry optimizations by analytical energy gradient of long-range corrected time-dependent density functional theory. J. Chem. Phys. 2006, 124, 144106. (20) Zhang, D.; Peng, D.; Zhang, P.; Yang, W. Analytic gradients, geometry optimization and excited state potential energy surfaces from the particle-particle random phase approximation. Phys. Chem. Chem. Phys. 2015, 17, 1025−1038. (21) Kaduk, B.; Tsuchimochi, T.; Van Voorhis, T. Analytic energy gradients for constrained DFT-configuration interaction. J. Chem. Phys. 2014, 140, 18A503. (22) Casida, M. E. Time-dependent density-functional theory for molecules and molecular solids. J. Mol. Struct.: THEOCHEM 2009, 914, 3−18. (23) Tawada, Y.; Tsuneda, T.; Yanagisawa, S.; Yanai, T.; Hirao, K. A long-range-corrected time-dependent density functional theory. J. Chem. Phys. 2004, 120, 8425−8433. (24) Shao, Y.; Head-Gordon, M.; Krylov, A. I. The spin-flip approach within time-dependent density functional theory: Theory and applications to diradicals. J. Chem. Phys. 2003, 118, 4807. (25) Harabuchi, Y.; Maeda, S.; Taketsugu, T.; Minezawa, N.; Morokuma, K. Automated Search for Minimum Energy Conical Intersection Geometries between the Lowest Two Singlet States S0/ S1-MECIs by the Spin-Flip TDDFT Method. J. Chem. Theory Comput. 2013, 9, 4116−4123.

likely recoverable through a multireference configuration interaction built from constrained DFTB and ΔDFTB electronic states.27,84 More rigorous treatments of the spin adaptation, analogous to ROKS31 or OCDFT,34 are also likely to improve the convergence behavior of the ΔDFTB algorithm. It is conceivable that the development of DFTB parameter sets specialized for particular classes of excited states, by analogy to the development of specialized exchange-correlation functionals for excited states,33,85 could afford limited improvements in accuracy without a loss of generality. However, since the formulation and application of excited-state exchange-correlation functionals are unresolved matters, we advise the continued adoption of existing DFTB parameter sets for use in ΔDFTB. Given its promising description of excited-state geometries and opportunities for improved robustness, ΔDFTB stands to become an efficient semiempirical quantum chemical method for molecular dynamics near equilibrium on low-lying singlet excited states. Applications underway will be described in future work.



ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jctc.5b00734. Vertical excitation energies for the small and large organic molecule test sets; vertical excitation and emission energies and optimized geometries for molecules in the Stokes shift test set; density distances for selected molecules; atomic spin constants; and optimized geometries for ground and excited states of acetone, acetamide, and aniline (PDF)



AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Funding

T.K. gratefully acknowledges a JSPS postdoctoral fellowship and start-up support from Western Washington University. S.I. acknowledges support from a CREST grant by JST. Notes

The authors declare no competing financial interest.



DEDICATION This article is dedicated to the memory of Tom Ziegler, whose early and recent work in excited-state DFT provided inspiration for this study.



REFERENCES

(1) Janssen, R. A.; Nelson, J. Factors Limiting Device Efficiency in Organic Photovoltaics. Adv. Mater. 2013, 25, 1847−1858. (2) Brédas, J.-L.; Norton, J. E.; Cornil, J.; Coropceanu, V. Molecular Understanding of Organic Solar Cells: The Challenges. Acc. Chem. Res. 2009, 42, 1691−1699. (3) Piletic, I. R.; Matthews, T. E.; Warren, W. S. Probing Near Infrared Photo-Relaxation Pathways in Eumelanins and Pheomelanins. J. Phys. Chem. A 2010, 114, 11483−11491. (4) Fassioli, F.; Dinshaw, R.; Arpin, P. C.; Scholes, G. D. Photosynthetic light harvesting: excitons and coherence. J. R. Soc., Interface 2014, 11, 20130901. (5) Wahadoszamen, M.; Margalit, I.; Ara, A. M.; van Grondelle, R.; Noy, D. The role of charge-transfer states in energy transfer and I

DOI: 10.1021/acs.jctc.5b00734 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation (26) Dreuw, A.; Head-Gordon, M. Failure of Time-Dependent Density Functional Theory for Long-Range Charge-Transfer Excited States: The Zincbacteriochlorin-Bacteriochlorin and Bacteriochlorophyll-Spheroidene Complexes. J. Am. Chem. Soc. 2004, 126, 4007− 4016. (27) Kaduk, B.; Van Voorhis, T. Communication: Conical intersections using constrained density functional theoryconfiguration interaction. J. Chem. Phys. 2010, 133, 061102. (28) Grimme, S.; Waletzke, M. A combination of KohnSham density functional theory and multi-reference configuration interaction methods. J. Chem. Phys. 1999, 111, 5645−5655. (29) Kowalczyk, T.; Yost, S. R.; Van Voorhis, T. Assessment of the ΔSCF density functional theory approach for electronic excitations in organic dyes. J. Chem. Phys. 2011, 134, 054128. (30) Maurer, R. J.; Reuter, K. Assessing computationally efficient isomerization dynamics: SCF density-functional theory study of azobenzene molecular switching. J. Chem. Phys. 2011, 135, 224303. (31) Kowalczyk, T.; Tsuchimochi, T.; Chen, P.-T.; Top, L.; Van Voorhis, T. Excitation energies and Stokes shifts from a restricted open-shell Kohn-Sham approach. J. Chem. Phys. 2013, 138, 164101. (32) Filatov, M. Spin-restricted ensemble-referenced Kohn-Sham method: basic principles and application to strongly correlated ground and excited states of molecules. WIREs Comput. Mol. Sci. 2015, 5, 146−167. (33) Ayers, P. W.; Levy, M.; Nagy, A. Time-independent densityfunctional theory for excited states of Coulomb systems. Phys. Rev. A: At., Mol., Opt. Phys. 2012, 85, 042518. (34) Evangelista, F. A.; Shushkov, P.; Tully, J. C. Orthogonality Constrained Density Functional Theory for Electronic Excited States. J. Phys. Chem. A 2013, 117, 7378−7892. (35) Ziegler, T.; Krykunov, M.; Cullen, J. The implementation of a self-consistent constricted variational density functional theory for the description of excited states. J. Chem. Phys. 2012, 136, 124107. (36) Park, Y. C.; Krykunov, M.; Ziegler, T. On the relation between adiabatic time dependent density functional theory (TDDFT) and the ΔSCF-DFT method. Introducing a numerically stable ΔSCF-DFT scheme for local functionals based on constricted variational DFT. Mol. Phys. 2015, 113, 1636−1647. (37) Voityuk, A. INDO/X: A new semiempirical method for excited states of organic and biological molecules. J. Chem. Theory Comput. 2014, 10, 4950−4958. (38) Silva-Junior, M. R.; Thiel, W. Benchmark of Electronically Excited States for Semiempirical Methods: MNDO, AM1, PM3, OM1, OM2, OM3, INDO/S, and INDO/S2. J. Chem. Theory Comput. 2010, 6, 1546−1564. (39) Niehaus, T. A.; Suhai, S.; Della Sala, F.; Lugli, P.; Elstner, M.; Seifert, G.; Frauenheim, T. Tight-binding approach to time-dependent density-functional response theory. Phys. Rev. B: Condens. Matter Mater. Phys. 2001, 63, 085108. (40) Trani, F.; Scalmani, G.; Zheng, G.; Carnimeo, I.; Frisch, M. J.; Barone, V. Time-Dependent Density Functional Tight Binding: New Formulation and Benchmark of Excited States. J. Chem. Theory Comput. 2011, 7, 3304−3313. (41) Heringer, D.; Niehaus, T. A.; Wanko, M.; Frauenheim, T. Analytical excited state forces for the time-dependent densityfunctional tight-binding method. J. Comput. Chem. 2007, 28, 2589− 2601. (42) Barone, V.; Carnimeo, I.; Scalmani, G. Computational Spectroscopy of Large Systems in Solution: The DFTB/PCM and TD-DFTB/PCM Approach. J. Chem. Theory Comput. 2013, 9, 2052− 2071. (43) Domínguez, A.; Aradi, B.; Frauenheim, T.; Lutsker, V.; Niehaus, T. A. Extensions of the Time-Dependent Density Functional Based Tight-Binding Approach. J. Chem. Theory Comput. 2013, 9, 4901− 4914. (44) Lundberg, M.; Nishimoto, Y.; Irle, S. Delocalization errors in a hubbard-like model: Consequences for density-functional tightbinding calculations of molecular systems. Int. J. Quantum Chem. 2012, 112, 1701−1711.

(45) Wanko, M.; Garavelli, M.; Bernardi, F.; Niehaus, T. A.; Frauenheim, T.; Elstner, M. A global investigation of excited state surfaces within time-dependent density-functional response theory. J. Chem. Phys. 2004, 120, 1674−1692. (46) Ziegler, T.; Rauk, A.; Baerends, E. J. Calculation of multiplet energies by the Hartree-Fock-Slater method. Theor. Chim. Acta 1977, 43, 261−271. (47) von Barth, U. Local-density theory of multiplet structure. Phys. Rev. A: At., Mol., Opt. Phys. 1979, 20, 1693−1703. (48) Runge, E.; Gross, E. K. U. Density-Functional Theory for TimeDependent Systems. Phys. Rev. Lett. 1984, 52, 997−1000. (49) Besley, N. A.; Gilbert, A. T. B.; Gill, P. M. W. Self-consistentfield calculations of core excited states. J. Chem. Phys. 2009, 130, 124308. (50) Cheng, C.-L.; Wu, Q.; Van Voorhis, T. Rydberg energies using excited state density functional theory. J. Chem. Phys. 2008, 129, 124112. (51) Maurer, R. J.; Reuter, K. Excited-state potential-energy surfaces of metal-adsorbed organic molecules from linear expansion -selfconsistent field density-functional theory (SCF-DFT). J. Chem. Phys. 2013, 139, 014708. (52) Filatov, M.; Shaik, S. A spin-restricted ensemble-referenced Kohn-Sham method and its application to diradicaloid situations. Chem. Phys. Lett. 1999, 304, 429−437. (53) Elstner, M.; Porezag, D.; Jungnickel, G.; Elsner, J.; Haugk, M.; Frauenheim, T.; Suhai, S.; Seifert, G. Self-consistent-charge densityfunctional tight-binding method for simulations of complex materials properties. Phys. Rev. B: Condens. Matter Mater. Phys. 1998, 58, 7260. (54) Seifert, G. Tight-binding density functional theory: An Approximate Kohn-Sham DFT scheme. J. Phys. Chem. A 2007, 111, 5609−5613. (55) Elstner, M. SCC-DFTB: What is the proper degree of selfconsistency? J. Phys. Chem. A 2007, 111, 5614−5621. (56) Oliveira, A. F.; Seifert, G.; Heine, T.; Duarte, H. A. Densityfunctional based tight-binding: an approximate DFT method. J. Braz. Chem. Soc. 2009, 20, 1193−1205. (57) Seifert, G.; Joswig, O. Density-functional tight binding - an approximate density-functional theory method. WIREs Comput. Mol. Sci. 2012, 2, 456−465. (58) Foulkes, W. M. C.; Haydock, R. Tight-binding models and density-functional theory. Phys. Rev. B: Condens. Matter Mater. Phys. 1989, 39, 12520. (59) Porezag, D.; Frauenheim, T.; Köhler, T.; Seifert, G.; Kaschner, R. Construction of tight-binding-like potentials on the basis of densityfunctional theory: Application to carbon. Phys. Rev. B: Condens. Matter Mater. Phys. 1995, 51, 12947. (60) Gaus, M.; Chou, C.-P.; Witek, H.; Elstner, M. Automatized parametrization of SCC-DFTB Repulsive Potentials: Application to hydrocarbons. J. Phys. Chem. A 2009, 113, 11866−11881. (61) Köhler, C.; Seifert, G.; Gerstmann, U.; Elstner, M.; Overhof, H.; Frauenheim, T. Approximate density-functional calculations of spin densities in large molecular systems and complex solids. Phys. Chem. Chem. Phys. 2001, 3, 5109−5114. (62) Zheng, G.; Witek, H. A.; Bobadova-Parvanova, P.; Irle, S.; Musaev, D. G.; Prabhakar, R.; Morokuma, K.; Lundberg, M.; Elstner, M.; Köhler, C.; et al. Parameter Calibration of Transition-Metal Elements for the Spin-Polarized Self-Consistent-Charge DensityFunctional Tight-Binding (DFTB) Method: Sc, Ti, Fe, Co, and Ni. J. Chem. Theory Comput. 2007, 3, 1349−1367. (63) Köhler, C. Bercksichtigung von Spinpolarisationseffekten in einem dichtefunktional- basierten Ansatz. Ph.D. Thesis, Universität Paderborn, 2004. (64) Gaus, M.; Cui, Q.; Elstner, M. DFTB3: Extension of the SelfConsistent-Charge Density-Functional Tight-Binding Method (SCCDFTB). J. Chem. Theory Comput. 2011, 7, 931−948. (65) Gaus, M.; Goez, A.; Elstner, M. Parametrization and benchmark of DFTB3 for organic molecules. J. Chem. Theory Comput. 2013, 9, 338−354. J

DOI: 10.1021/acs.jctc.5b00734 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation (66) Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized gradient approximation made simple. Phys. Rev. Lett. 1996, 77, 3865−3868. (67) Becke, A. D. Density-functional Thermochemistry. III. The Role of Exact Exchange. J. Chem. Phys. 1993, 98, 5648−5652. (68) Perdew, J. P.; Ernzerhof, M.; Burke, K. Rationale for mixing exact exchange with density functional approximations. J. Chem. Phys. 1996, 105, 9982−9985. (69) Adamo, C.; Barone, V. Toward reliable density functional methods without adjustable parameters: The PBE0 model. J. Chem. Phys. 1999, 110, 6158−6170. (70) Aradi, B.; Hourahine, B.; Frauenheim, T. DFTB+, a Sparse Matrix-Based Implementation of the DFTB method. J. Phys. Chem. A 2007, 111, 5678−5684. (71) Kubillus, M.; Kubar̆, T.; Gaus, M.; R̆ ezác,̆ J.; Elstner, M. Parameterization of the DFTB3 Method for Br, Ca, Cl, F, I, K, and Na in Organic and Biological Systems. J. Chem. Theory Comput. 2015, 11, 332−342. (72) Shao, Y.; Gan, Z.; Epifanovsky, E.; Gilbert, A. T. B.; Wormit, M.; Kussmann, J.; Lange, A. W.; Behn, A.; Deng, J.; Feng, X.; et al. Advances in molecular quantum chemistry contained in the Q-Chem 4 program package. Mol. Phys. 2015, 113, 184−215. (73) D’Amico, I.; Coe, J. P.; França, V. V.; Capelle, K. Quantum Mechanics in Metric Space: Wave Functions and Their Densities. Phys. Rev. Lett. 2011, 106, 050401. (74) Schreiber, M.; Silva-Junior, M. R.; Sauer, S. P. A.; Thiel, W. Benchmarks for electronically excited states: CASPT2, CC2, CCSD, and CC3. J. Chem. Phys. 2008, 128, 134110. (75) Gilbert, A. T. B.; Besley, N. A.; Gill, P. M. W. Self-consistent field calculations of excited states using the maximum overlap method (MOM). J. Phys. Chem. A 2008, 112, 13164−13171. (76) Niehaus, T. A.; Della Sala, F. Range separated functionals in the density functional based tight-binding method: Formalism. Phys. Status Solidi B 2012, 249, 237−244. (77) Guareschi, R.; Filippi, C. Ground- and Excited-State Geometry Optimization of Small Organic Molecules with Quantum Monte Carlo. J. Chem. Theory Comput. 2013, 9, 5513−5525. (78) Zeng, Q.; Liu, J.; Liang, W. Molecular properties of excited electronic state: Formalism, implementation, and applications of analytical second energy derivatives within the framework of the timedependent density functional theory/molecular mechanics. J. Chem. Phys. 2014, 140, 18A506. (79) Hanson-Heine, M. W. D.; George, M. W.; Besley, N. A. Calculating excited state properties using Kohn-Sham density functional theory. J. Chem. Phys. 2013, 138, 064101. (80) Hanson-Heine, M. W. D.; Wriglesworth, A.; Uroos, M.; Calladine, J. A.; Murphy, T. S.; Hamilton, M.; Clark, I. P.; Towrie, M.; Dowden, J.; Besley, N. A.; et al. Calculating singlet excited states: Comparison with fast time-resolved infrared spectroscopy of coumarins. J. Chem. Phys. 2015, 142, 154119. (81) Valsson, O.; Filippi, C.; Casida, M. E. Regarding the use and misuse of retinal protonated Schiff base photochemistry as a test case for time-dependent density-functional theory. J. Chem. Phys. 2015, 142, 144104. (82) Gritsenko, O. V.; van Gisbergen, S. J. A.; Schipper, P. R. T.; Baerends, E. J. Origin of the field-counteracting term of the KohnSham exchange-correlation potential of molecular chains in an electric field. Phys. Rev. A: At., Mol., Opt. Phys. 2000, 62, 012507. (83) Molteni, C.; Frank, I.; Parrinello, M. An Excited State Density Functional Theory Study of the Rhodopsin Chromophore. J. Am. Chem. Soc. 1999, 121, 12177−12183. (84) Scholz, R.; Luschtinetz, R.; Seifert, G.; Jägeler-Hoheisel, T.; Körner, C.; Leo, K.; Rapacioli, M. Quantifying charge transfer energies at donoracceptor interfaces in small-molecule solar cells with constrained DFTB and spectroscopic methods. J. Phys.: Condens. Matter 2013, 25, 473201. (85) Samal, P.; Harbola, M. K. Exploring foundations of timeindependent density functional theory for excited states. J. Phys. B: At., Mol. Opt. Phys. 2006, 39, 4065.

K

DOI: 10.1021/acs.jctc.5b00734 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX