Nonadiabatic Dynamics of Cycloparaphenylenes with TD-DFTB

Article Cite This: J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

pubs.acs.org/JCTC

Nonadiabatic Dynamics of Cycloparaphenylenes with TD-DFTB Surface Hopping Ljiljana Stojanović,*,† Saadullah G. Aziz,‡ Rifaat H. Hilal,§,‡ Felix Plasser,∥ Thomas A. Niehaus,⊥ and Mario Barbatti*,† †

Aix Marseille Univ, CNRS, ICR, Marseille, France Chemistry Department, Faculty of Science, King Abdulaziz University, Jeddah B.O. 208203, Saudi Arabia § Chemistry Department, Faculty of Science, Cairo University, Giza, Egypt ∥ Institute of Theoretical Chemistry, Faculty of Chemistry, University of Vienna, Währinger Straße 17, 1090 Vienna, Austria ⊥ Univ Lyon, Université Claude Bernard Lyon 1, CNRS, Institut Lumière Matière, F-69622 Lyon, France ‡

S Supporting Information *

ABSTRACT: We implemented a version of the decoherencecorrected fewest switches surface hopping based on linearresponse time-dependent density functional tight binding (TDDFTB), enhanced by transition density analysis. The method has been tested for the gas-phase relaxation dynamics of two cycloparaphenylene molecules, [8]CPP and [10]CPP, explaining some important features of their nonadiabatic dynamics, such as the origin of their long fluorescence lifetimes (related to the slow radiative emission from the S1 state) and the trend of increasing the fluorescence rate with the molecular size (related to an increase in the S1−S0 energy gaps and oscillator strengths in the larger molecule). The quality of the TD-DFTB electronic structure information was assessed through four quantities: excitation energies; charge-transfer (CT) numbers, which estimate the charge transfer character of states; participation ratio (PR), which describes delocalization of electronic density; and participation ratio of natural transition orbitals (PRNTO), which describes the multiconfigurational character of states. These quantities were computed during dynamics and recomputed for the same geometries with the higher-level long-range-corrected TD-LC-DFTB and a lower-level single-determinant approximation for the excited states, SD-(LC)-DFTB. Taking TD-LC-DFTB as the standard, TD-DFTB underestimates the excitation energies by ∼0.5 eV and overestimates CT and PR. SD-DFTB underestimates excitation energies and overestimates CT to the same extent that TD-DFTB does, but it predicts reasonable PR distributions. SD-LC-DFTB leads to an extreme overestimation of the excitation energies by ∼3 eV, overestimates the charge transfer character of the state, but predicts the PR values very close to those obtained with TD-LC-DFTB. The density functional tight binding (DFTB)16 method bridges the gap between first-principles electronic structure methods and empirical tight-binding approaches. DFTB can be seen as an approximation of density functional theory (DFT) at the generalized-gradient approximation (GGA) level. The basic algorithms implemented in the DFTB+ code,17 which we use here, are the non-self-consistent charge (non-SCC)18,19 and SCC-DFTB (DFTB2 and DFTB3)16 for the electronic ground state. The treatment of dispersion interactions (+D) is possible in the DFTB+ code via Lennard-Jones potentials,20 with the Slater−Kirkwood polarizable atomic model,21 and on the level of the DFT-D3 method.22 Electronic excitations can be calculated with linear-response time-dependent approach (LR TD-DFTB).23 A recent extension of the method now also allows for computations using long-range-corrected (LC) functionals,24,25 which provide a better description of state

1. INTRODUCTION Nonadiabatic effects influence the underlying relaxation mechanisms of a range of processes that occur after photoexcitation in systems of all sizes, spanning from single molecules to solid-state materials. Some of these nonadiabatic processes such as charge and energy transfer, electron−hole separation, singlet fission, etc.1−5 have attracted much attention, because of their potential to aid the development of photonic and photovoltaic materials. Therefore, many research groups have been aiming to develop and extend suitable theoretical algorithms for simulation of nonadiabatic processes in extended systems.6−8 Since the dimensionality of extended systems is a major limiting factor in their excited-state dynamics simulations, the methods that have been used for this purpose usually rely on mixed quantum-classical nonadiabatic dynamics approaches, such as the trajectory surface hopping in its different variants.6,7,9 These methods are often associated with fast semiempirical algorithms for electronic structure computations.7,10−15 © XXXX American Chemical Society

Received: September 27, 2017 Published: November 15, 2017 A

DOI: 10.1021/acs.jctc.7b01000 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation

units bonded in para positions to form a ring. They have been synthesized by several experimental groups37 and, since then, studied in detail, both experimentally38−41 and theoretically.42−44 They have received much attention, because of their potential use as precursors in the synthesis of armchair carbon nanotubes45 and as efficient emitters of tunable wavelengths.46 CPPs show unusual optical properties, markedly distinct from the linear paraphenylenes,47 such as sizeindependent absorption wavelengths, blue shifting of emission peaks with an increase in the ring size,39 and double-peaked emission spectra,39 among others. CPP excitations are characterized by long relaxation lifetimes in the nanosecond domain.41,48 Small CPPs are also interesting, because of their large quantum yields for singlet oxygen production in the reaction with triplet oxygen.41,48 As we shall discuss, nonadiabatic dynamics with LR TDDFTB explore states with somewhat large CT character. This is not surprising, given that conventional DFTB is parametrized for PBE, where exactly this same type of problem is expected to occur.49 Recent results for (LR) TD-LC-DFTB make clear that this level provides correct energy character, fixing the main deficiencies of conventional DFTB.25 Since long-range corrected functionals are indispensible for the correct description of charge-transfer states, but the analytical gradients for long-range corrected TD-DFTB are under development at the moment, we resort to a multilayer dynamics scheme, where the conventional simulations are performed with TD-DFTB and the TD-LC-DFTB computations of the excited states populated during dynamics are performed a posteriori. This approximation precludes us from drawing conclusions about the lifetimes of states, because of the mismatch between TDDFTB and TD-LC-DFTB gradients. Still, plenty of important insights about relaxation mechanisms could be deduced, primarily because the long-enough dynamics provides a map of the excited states hypersurfaces, particularly the S 1 hypersurface, and, thus, fluorescence and internal conversion processes can be studied, as well as the properties of occurring excitons. The Newton-X/DFTB+/TheoDORE interface that we present in this paper also allows analysis of the transition density matrices of the states populated during dynamics. We utilized this feature to perform the statistical analysis of excitation energies, charge transfer, density delocalization, and multiconfigurational character for a large ensemble of points generated during the dynamics, based on different approximations of the S1 state wave functions, (LR) TD-(LC)-DFTB and SD-(LC)-DFTB. This comparison allows us to derive a comprehensive picture about what types of deviations in the electronic densities are caused by an exploration of the conformational space by GGA functionals and single determinant approximations. This analysis confirmed that the full LR TD-DFTB representation of the wave function is a better approximation for the excited states than the singledeterminant one, often applied in the simulations of nanosized systems.10 This paper is organized as follows. In section 2, “Method Description”, we briefly review the fewest-switches surface hopping and TD-DFTB methods, providing details on the interfaces between DFTB+, Newton-X, and TheoDORE codes. Section 3, “Computational Details” contains the technical information about the methods used for the actual calculations. In section 4, “Results and Discussion”, we first discuss the reliability of the description of the excited states of the [8]CPP

localization and charge transfer excitations. The low computational cost of these algorithms renders levels up to LR TD(LC)-DFTB3 (+D3) suitable for simulations of excited states in extended systems. In this work, we implemented an interface between DFTB +17 and Newton-X,26,27 which allowed us to run decoherencecorrected fewest switches trajectory surface hopping (FSSH) with the TD-DFTB method23 available in the recently released version (v 1.3) of the DFTB+ code. There are reports of few other implementations of nonadiabatic dynamics with mixed quantum-classical methods based on DFTB. Jakowski and Morokuma implemented Ehrenfest dynamics based on Liouville-von Neumann equation for density matrices propagation and made several applications with the SCC-DFTB method.28 Ehrenfest dynamics based on real-time Kohn−Sham DFTB has been described in refs 14 and 29. Surface hopping dynamics, which are also based on real-time Kohn−Sham DFTB, is discussed in refs 10 and 15. In these surface hopping implementations, the excited states are represented by single Kohn−Sham determinants that are based on stationary orbitals (SD-DFTB), which may lead to an inaccurate description of adiabatic crossings.30 Surface hopping with a more robust description of the excited states should be obtained via linearresponse TD-DFTB. Nonadiabatic dynamics with linearresponse TD-DFTB, which has also been implemented by Mitrić et al.,11 based on an earlier version of DFTB, is the basis of our developments here. Our surface hopping based on linear-response TD-DFTB implementation relies on a general interface between the newest versions of the DFTB+ and Newton-X codes, which are both publicly available free of charge. It enables diverse new features to improve the description of ground-state and excitedstate energies, such as the DFTB third-order expansion (DFTB3 model)31 and onsite corrections for TD-DFTB:32 to accelerate the dynamics, such as fast nonadiabatic coupling calculations based on the recently proposed orbital derivative method;33 and to enhance the data analysis, via a direct module for advanced transition density analysis based on the program TheoDORE.34 Moreover, this interface is ready to profit from analytical energy gradients for TD-DFTB, based on long-rangecorrected functionals,24,25 which are currently under development. The FSSH/TD-DFTB dynamics (on the GGA level of functional) has a potential for the simulation of processes where the charge-transfer states are not involved in dynamics (e.g., relaxation processes of small molecules in solutions, where both solute and optically inactive solvent molecules are treated with TD-DFTB; see, e.g., ref 11), or in the cases where the occurring charge-transfer states are of short distance (for instance, in the process of hot electron injection in photovoltaic materials; see, e.g., ref 35). Even in cases where TD-DFTB is not expected to work well, mainly due to artifacts in the description of charge-transfer states, a multilayer approach may be devised,36 where dynamics is propagated with TD-DFTB and single points are recalculated at a higher level. Naturally, in such a case, time-dependent properties (population lifetimes, for instance) may not be adequately described, while statistical distributions of the excited state in equilibrium may still provide relevant information. As a test case for this new implementation, we choose to study the nonadiabatic dynamics of two cycloparaphenylenes (CPPs): [8]CPP and [10]CPP. Cycloparaphenylenes are hoopshaped conjugated molecules composed of repeating phenylene B


Article

Journal of Chemical Theory and Computation ⎡ −2Δt ⎤ PM → J = max⎢0, σ Re( c M ̅ * cJ̅ )⎥ 2 MJ ⎣ | cM̅ | ⎦

and [10]CPP on the TD-DFTB level, comparing the vertical excitations with TD-DFT results calculated with three different functionals. Next, we discuss the TD-DFTB excited-state dynamics of CPPs. Finally, we statistically analyze the quality of the TD-DFTB electronic structure information produced during the dynamics and evaluate the SD-(LC)-DFTB levels.

where Δt is the time step for the integration of the TDSE. Thus, once the initial conditions (positions and velocities) are created, to perform the FSSH simulations, one must solve the electronic TDSE at each nuclear position given by Newton’s equations to obtain the excitation energies and gradients, and to evaluate time-dependent NACs. These quantities are interpolated between two classical time steps to be used in the integration of the TDSE. Our FSSH/TD-DFTB implementation profits of the algorithmic infrastructure already implemented in Newton-X. In the standard way, time-dependent NACs are computed by using a finite differences method, as proposed by HammesSchiffer and Tully.52 In this approach, TD-NACs are given in terms of excited-state wave function overlaps between two subsequent nuclear time steps,

2. METHOD DESCRIPTION 2.1. Fewest-Switches Trajectory Surface Hopping (FSSH). FSSH is a mixed quantum-classical approach to quantum dynamics, which approximates the wave packet motion by a swarm of independent nuclear trajectories, each one evolving classically on a single electronic state.50 Nonadiabatic transitions between states are introduced by a stochastic algorithm that may instantaneously switch the surface on which the trajectory is propagated, based on evaluating of hopping probabilities.50 The hopping probabilities in FSSH are obtained by expanding the time-dependent electronic wave function as Ψ(r; R(t )) =

∑ cK(t )ψK(r; R(t )) K

σJK (t ) ≈ (1)

dcJ dt

= EJ cJ − iℏ ∑ σJK cK K

SJK (t ) = ⟨ψJ(t − Δt )|ψK (t )⟩

ψJ

sij(t ) = ⟨χi (t − Δt )|χj (t )⟩

(2)

(3)

The coefficients cI are further corrected for decoherence effects with the simplified decay of mixing method (SDM),51 according to cK̅ = cK e−Δt / τKM

∀K≠M

⎡ c cM̅ = M ⎢1 − |cM| ⎢⎣

2⎥

∑ K ≠M

⎤1/2

| cK̅ |

⎥⎦

EDFTB[ρ] = (4)

E 0[ρ0 ] + E1[ρ0 , δρ] + E2[ρ0 , (δρ)2 ] + E3[ρ0 , (δρ)3 ] + ···  

In these equations, M is the active state and the decoherence time τKM is given by the phenomenological equation −1 |EK − EM | ⎛ α⎞ 1 = ⎜1 + ⎟ τKM ℏ Kn ⎠ ⎝

(9)

Instead of the integration of the TDSE equation, having SJK computed with either method described above, the timedependent coefficients cK can also be obtained via the local diabatization method55 that was proposed by Granucci and Persico.13 In this case, the coefficients are obtained by a unitary transformation between two time steps, with the help of an adiabatic-to-diabatic transformation matrix obtained under the condition σJK = 0. 2.2. The Basics of the TD-DFTB Method. The basic idea of the density-functional-based tight binding (DFTB) method is to perform a Volterra expansion of the Kohn−Sham (KS) DFT total energy functional around a reference electronic density (ρ0),16

∂ψK ∂t

(8)

The wave function overlaps (SJK) can be evaluated based on Slater determinants (the determinant derivative method).53 The TD-NACs can be alternatively evaluated with the recently implemented orbital derivative method54 developed by Izmaylov’s group,33 which relies on time derivatives that are based on molecular orbitals. In this case, the calculation of the SJK overlaps is skipped, and the NACs are directly computed in terms of molecular orbital (or KS orbital) overlaps:

where EJ is the energy of the electronic state J and the bra-ket in the second term of the right side represents the time-derivative nonadiabatic coupling (NAC) between states J and K:

σJK =

1 [3SJK (t ) − 3SKJ(t ) − SJK (t − Δt ) + SKJ(t − Δt )] 4Δt (7)

where

where ψK(r;R(t)) is the adiabatic electronic wave function for state K, which is a function of the electron coordinates r and parametrically dependent on the classical nuclear coordinate R; cK(t) represents complex time-dependent expansion coefficients. The nuclear trajectories R(t) are propagated by integrating Newton’s equations on the active state M. Simultaneously to the propagation of Newton’s equations, a local-approximation of the time-dependent electronic Schrödinger equation (TDSE) is integrated to obtain the expansion coefficients cK(t). In the adiabatic representation, it is given by iℏ

(6)

DFTB1   DFTB2   DFTB3

(10)

(5)

where ρ0 is the sum of neutral atomic densities. Within the DFTB+ code, the second-order (SCC-DFTB or DFTB2) expansion16 and the third-order (DFTB3) expansion56 of the density functional have been implemented. The total energy within DFTB is derived by employing several approximations.16 A minimal basis set of localized

where EI is the potential energy of state I, Kn is the classical kinetic energy of the nuclei, and α is a parameter whose recommended value is 0.1 hartree. The FSSH probability50 between the active state M and another state K at each time step is computed as C


Article


can be calculated in terms of the wave function overlaps SJK (eq 8) or molecular (or KS) orbital overlaps sij (eq 9), depending on the approach chosen. The algorithm for computation of these quantities has been already implemented in the NewtonX code for the calculation of the TD-DFT NACs.53 Here, we only describe them briefly to point out the approximations used in the case of TD-DFTB. The Casida approximationwhich states that the excitedstate wave function can be represented using a CIS-like expansion, based on the singly excited Slater determinantsis applied.58 This approximation was previously utilized in the implementation of the surface hopping dynamics with the TDDFT method.59−61 The number of the singly excited Slater determinants that is obtained by applying a single excitation operator to the ground-state Slater determinant is Ndet = NoccNvirt, and the form of the excited-state wave function is

atomic orbitals, obtained by solving the DFT equations for atoms, is used in the expansion of KS orbitals. The matrix elements of the part of the Hamiltonian that is dependent only on the reference density and orbitals are evaluated based on DFT computations. Its diagonal elements are simply equal to the computed atomic orbital energies, and the nondiagonal elements are precomputed for pairs of elements for a range of internuclear distances and tabulated in so-called Slater−Koster (SK) files.16 Overlaps of atomic orbitals are similarly precomputed and stored.16 The terms of the energy expression, depending on the electronic density fluctuation δρ, are simplified using the Mulliken approximation by means of atomic point charges.16 Remaining integrals are parametrized using Hubbard-like parameters for each atom computed with DFT. Furthermore, the nuclear repulsion term and the terms that are dependent only on the reference density are approximated as a sum of short-range repulsive potentials that are dependent only on diatomic distance.16,56,57 The LR TD-DFTB method applies a similar procedure as LR TD-DFT in the computation of excitation energies,23 which are given as the eigenvalues ω of the equation ⎡ A B ⎤⎡ X ⎤ ⎡ 1 0 ⎤⎡ X ⎤ ⎥⎦⎢⎣ ⎥⎦ = ω⎢⎣ ⎢⎣ ⎥⎢ ⎥ 0 −1⎦⎣ Y ⎦ B* A* Y

|ψJ ⟩ =

(13)

i,a

where the CIS coefficients X + Y are the eigenvectors of eq 11 and |Φia⟩ are the singly excited Slater determinants. The overlap of the excited-state wave functions represented in the basis of Slater determinants at time steps t and t + Δt is ⟨ψ I(t )|ψ J(t + Δt )⟩

(11)

=

where the matrix elements are written in terms of orbital transitions from occupied orbitals (i,j) into virtual orbitals (a,b) as

∑ (X + Y )iaI*(X + Y )Jjb ⟨Φia(t )|Φjb(t + Δt )⟩ (14)

ia , jb

The overlap of two Slater determinants is computed as the determinant of the overlap matrix of KS orbitals χi obtained at two subsequent time steps.62 Furthermore, the overlaps of KS orbitals (eq 9) are represented in terms of overlaps of atomic orbitals ϕμ at two subsequent geometries from the dynamics:

Aia , jb = δijδab(εa − εi) + (ia||jb) Bia , jb = (ia||bj)

∑ (X + Y )iaJ |Φia⟩

(12)

In these equations, (ia∥jb) are two-electron integrals that involve both Coulomb and exchange-correlation kernels and ε represents the orbital energies from a previous DFTB groundstate calculation. In the framework of TD-DFTB, the twoelectron integrals are calculated by applying the generalized Mulliken approximation, according to which the transition densities between different Kohn−Sham orbitals are represented by point charges.23 The LR TD-DFTB method has recently been expanded beyond the generalized Mulliken approximation by introducing onsite exchange-like integrals.32 This correction improves the description of nπ* and σπ* states, which were not well-described in the original version of LR TD-DFTB.32 Moreover, long-range-corrected (LC) exchange-correlation functionals have also been implemented and benchmarked in the DFTB method (LC-DFTB)24 and in the LR TD-DFTB method for excited-state simulations (TD-LC-DFTB).24 Benchmark calculations on a set of small chromophores showed that the TD-LC-DFTB method can treat charge transfer excitations accurately. For this class of excitations, DFTB (like any DFT method) using the generalized-gradient approximation (GGA) underestimates excitation energies strongly.25 Although the lack of analytical excited-state gradients within TD-LC-DFTB precludes the use of this approach in surface hopping simulations at the moment, we will still use it here as benchmark results for single-point calculations. 2.3. Time-Dependent Wave Function Overlaps with TD-DFTB. As discussed above, the core quantity to compute the hopping probabilities are determined by using NACs, which

⟨χi (t )|χj′ (t + Δt )⟩ =

∑ cμ*i(t )cνj(t + Δt )⟨ϕμ(t )|ϕv′(t + Δt )⟩ μν

(15)

These atomic orbitals do not form an orthonormal set, because they are computed at subsequent geometries in the dynamics, which differ slightly from each other. As already mentioned, the atomic orbital overlaps are precomputed in DFTB at different internuclear distances and stored in Slater−Koster (SK) files. Thus, during an SCC-DFTB computation, these integrals are not calculated, but they are read from the SK files. During the dynamics, some of the atom pairs from two subsequent geometries might become displaced, with respect to each other, by very short distances, which are not covered within the precomputed atomic orbital overlaps, because they are not relevant for simulations of single molecules. In this case, we evaluate NACs by applying a linear extrapolation using the computed values from the last two time steps. This problem does not happen with the third-order DFTB and the 3ob SK set optimized for it, because atomic orbital overlaps in the standard 3ob set are computed for interatomic distances larger than 0.02 Bohr. This problem might occur in simulations with the second-order DFTB, where the mio Slater-Koster set should be used, in which the atomic orbital overlaps are defined for minimum interatomic distance of 0.4 Bohr. Therefore, we suggest using the third-order DFTB and the 3ob basis set to circumvent the extrapolation of the nonadiabatic couplings. Since DFTB is a valence-electron-only method, core occupied orbitals do not appear in the construction of Slater determinants. This fact accelerates the computation of wave D


Article

Journal of Chemical Theory and Computation function overlaps for two reasons: (1) the dimension of the Slater determinants is reduced; and (2) the number of Slater determinants is decreased from NoccNvirt to (Nocc − Ncore)Nvirt. Hence, the formal scaling of the determinant derivative method is reduced from Nocc5Nvirt5 to (Nocc − Ncore)5 Nvirt2 and the scaling of the orbital overlap method is reduced from NoccNvirt2 to (Nocc − Ncore)Nvirt2. We tested the nonadiabatic couplings between S3 and S2 states in one sequence of the TD-DFTB surface hopping dynamics of thiophene by comparing the computed NACs between these states with the values computed from the TDPBE wave functions imposing the geometries from the TDDFTB dynamics. In the case of the TD-DFT wave functions, all occupied orbitals are used in the construction of ground-state and excited Slater determinants. In both cases, the NACs were computed with the local diabatization method.55 The values of NACs obtained from TD-DFTB and TD-PBE wave functions are in excellent agreement (Figure 1), which indicates that the

J Ω AB =

1 2

∑ (Dab0J )2 a∈A b∈B

(16)

where a and b are molecular orbitals belonging to fragments A and B, and D0J ab is the element of the transition density matrix in the localized basis, for the transition from the ground state 0 into the excited state J.63 The charge transfer properties of the excited states are described with the charge transfer (CT) index, which is defined as a sum of the nondiagonal ΩJAB elements divided by their total sum (Ω), CTJ =

1 Ω

∑

J Ω AB

A B≠A

(17)

They have values between 0, in the case of localized excitations (Frenkel excitons), and 1, in the case of charge-transfer excitations. The participation ratio (PR) defines the number of fragments over which an exciton is delocalized. It is defined as63 PR J =

⎤ Ω2 ⎡⎢ 1 1 ⎥ + 2 2 2 ⎢⎣ ∑A (∑B ΩAB) ∑B (∑A ΩAB) ⎥⎦

(18)

where the first term quantifies the extent of the hole delocalization and the second term quantifies the extent of the electron delocalization over fragments. A compact description of the excited states, in terms of a relatively small number of transitions, can be achieved in the basis of natural transition orbitals (NTOs).68 The NTOs are obtained by a singular value decomposition of the transition density matrix represented in a molecular orbital basis, D0,J.

Figure 1. Comparison of nonadiabatic couplings (a.u.) between S3 and S2 states in the sequence of a trajectory of thiophene computed based on the TD-DFTB and TD-PBE wave functions.

D0, J = U diag( λ1 ,

λ 2 , ...)VT

(19)

where the matrices U and V collect the sets of hole and electron orbitals, respectively, and λ i are the weights of the corresponding transitions between them.63 The participation ratio of natural transition orbitals (PRNTO), which is defined as

neglect of core molecular orbitals does not significantly affect the excited-state overlaps. Moreover, TD-DFTB excited-state wave functions seem to approximate their TD-DFT counterparts accurately also in regions of strong nonadiabatic coupling. 2.4. Excited-State Transition Density Matrix Analysis. To analyze the nature of the excited states occurring during TD-DFTB surface hopping dynamics, we interfaced the TheoDORE program with the DFTB+ program. TheoDORE is a program suite intended for the electronic wave function analysis based on state and transition density matrices between electronic states.63−65 The implemented methods are particularly suitable for obtaining insight into excitonic structure in large conjugated systems66,67 and are, therefore, a natural choice for analyzing CPPs. The most important part of this analysis refers to the determination of charge-transfer character of the excited states and their spatial localization properties by computing charge transfer and fragment participation ratio numbers. In addition, the TheoDORE program also implements the computation of the natural transition orbitals of the excited states and their participation ratio of natural transition orbitals. Segmenting a molecule into suitable fragments, one can evaluate the extent of charge transfer between fragments upon excitation by computing the value

PRNTO =

2 (∑i λi) ∑i λi 2

(20)

quantifies the number of transitions between NTOs necessary to describe the excited state. This quantity was initially introduced as a collectivity measure in ref 69, and its physical meaning has recently been elucidated in the context of quantum entanglement.70 Its upper bound is the number of nonzero singular values (NSVD) of D0,J, and it is equal to this number in the case when all nonzero singular values have the same magnitude. Oppositely, the lowest value of PRNTO, reached in the case when the excited state can be represented by only one transition between NTOs, is 1.63

3. COMPUTATIONAL DETAILS We created 500 sets of initial conditions (ICs) (geometries and momenta) using a harmonic-oscillator Wigner distribution,71 based on normal modes in the ground states of [8]CPP and [10]CPP. The vertical excitation energies and oscillator strengths for the transitions to the six lowest-lying states were computed with the TD-DFTB method with the third-order expansion of density32 and 3ob-1−1 SK set72 (Cartesian E


Article


For the comparison with the TD-DFTB results, we also computed vertical excitation energies of both molecules with the long-range corrected TD-LC-DFTB method implemented in the DFTB+ code and with TD-DFT based on three different functionals (PBE,77 B3LYP,78,79 and CAM-B3LYP80), using the GAMESS-US code.81,82

coordinates of the ground state minima of [8]CPP and [10]CPP are given in section SI-I in the Supporting Information. The dispersion interactions were included via a Lennard-Jones potential21 with parameters adopted from the universal force field (UFF).20 The sampled geometries are used in the simulations of the photoabsorption spectra with the nuclear ensemble approach73 employing TD-DFTB excitation energies and oscillator strengths. The simulated spectra of both molecules are featured by broad bands peaked in both cases at ∼2.95 eV (see section SI-II in the Supporting Information). Nonadiabatic excited-state dynamics simulations were performed on TD-DFTB potential energy surfaces. The sets of ICs for dynamics propagation were selected from the initial set of 500 ICs starting from the bright S6 and S5 states applying a very narrow energy window (2.95 ± 0.05 eV) around the maxima of the absorption spectra. This procedure produced 72 ICs in the S5 state and 35 ICs in the S6 state for [8]CPP. We selected 35 ICs from the S5 state and 17 ICs from the S6 state for dynamics propagation. In the case of [10]CPP, 64 ICs in the S5 state and 32 initial conditions in the S6 state are filtered, from which we selected 32 ICs in the S5 and 16 ICs in the S6 state for dynamics propagation. We limited the number of trajectories to ∼50 per molecule, while keeping in mind that the 3 ps propagation of larger number of trajectories would require substantial computational time. The TD-DFTB dynamics (GGA functional) of [8]CPP lasted ∼10 CPU minutes per time step (on an Intel(R) Xeon 3.0 GHz processor), from which the electronic structure and gradients computations elapsed ∼90% of the time and ∼10% of the time was spent in nonadiabatic couplings computation. For the comparison, the same simulations with the TD-PBE/6-31G method require ∼110 CPU minutes for the electronic structure and gradient computations and ∼2 min for nonadiabatic coupling computations on the same processor. Nonadiabatic couplings were computed in both simulations with the overlap derivative method. Thus, the construction of excited-state functions in the valence-only approximation leads to ∼50% time savings for computations of nonadiabatic couplings of [8]CPP with the OD method. In the case of determinant derivative methods, the time savings would be even more pronounced, because of the fifth-order scaling with the number of occupied orbitals used in the Slater determinant constructions. Nonadiabatic events between excited states were taken into account by FSSH corrected for decoherence with the SDM method.51 Since TD-DFTB, as a single-reference method, suffers from the same problems as TD-DFT in the description of S1−S0 NAC topology,74 we introduced a condition according to which the internal conversion between S1 and S0 states is assumed, and trajectories are stopped, whenever the energy gap between these states drops below 0.15 eV. Newton’s equations of motion were integrated using the velocity Verlet algorithm75 with a time step of 0.5 fs. The TDSE was integrated using the fifth-order Butcher’s algorithm76 with a time step of 0.025 fs, applying interpolated electronic properties between classical steps. The nonadiabatic couplings between excited states were evaluated numerically using the orbital derivative method.33,54 The dynamics was propagated for 3 ps. The excited states occurring during dynamics were analyzed based on their transition density matrices. For this purpose, the molecules are fragmented into phenylene rings. CT, PR, and PRNTO values for the active excited state during the trajectories were computed using the TheoDORE program.

4. RESULTS AND DISCUSSION 4.1. Ground-State and Excited-State Geometries. The optimized geometries of the ground state and the first excited states of [8]CPP and [10]CPP are shown in Figure 2. The

Figure 2. Structures of [8]CPP and [10]CPP in the ground state and the first excited states optimized with (TD-)DFTB.

molecules are highly symmetric in their ground states, having the symmetries of the D4d ([8]CPP) and D5d ([10]CPP) point groups. Their characteristic geometrical features are influenced by three different effects: backbone strain, steric effects, and competition between π-conjugation and quinoidalization.40 Their ground states are featured by alternant dihedral angles between subsequent phenylene rings with values of ±32.1° in [8]CPP and ±35.5° in [10]CPP. These values are in excellent agreement with those obtained in TD-CAM-B3LYP optimization (32° and 36°, respectively).43 The optimized symmetric local minima of the S1 states are characterized by planarization of phenylene rings with decreased alternant dihedral angles, obtaining the values ±9.9° and ±20.2° in [8]CPP and [10]CPP, respectively. Another important distinction between S0 and S1 geometries is in the differences between double intraring and single inter-ring C−C bond lengths. These differences decrease in the S1 states, in comparison to the ground statesin [8]CPP, they are 0.07 in S0 and 0.032 Å in S1; in [10]CPP, they are 0.085 and 0.045 Å. These features indicate that quinoidalization, which is accompanied by ring planarization, is favored in the S1 states, compared to the ground states, where π-conjugation has a larger effect on the structures. Quinoidalization is more pronounced in the smaller molecule, because of the influence of larger ring strain, which is reflected in decreased dihedral angles and differences between the lengths of single and double C−C bonds. These findings are consistent with the experimental results obtained with Raman spectroscopy40 and geometrical observations from TD-CAM-B3LYP optimizations.43 F


Article


Table 1. Computed and Experimental Vertical Excitation Energies and Oscillator Strengths of Excited States of [8]CPP and [10]CPP TD-DFTB E (eV)

a

TD-PBE f

E (eV)

TD-B3LYP f

E (eV)

2.443 2.842 2.842 3.215 3.215 3.270

0.000 0.006 0.006 0.860 0.860 0.000

2.336 2.775 2.775 3.214 3.214 3.273

0.000 0.000 0.000 0.946 0.946 0.000

2.879 3.563 3.563 3.675 3.675 3.957

2.669 2.907 2.907 3.139 3.139 3.168

0.000 0.008 0.008 0.855 0.855 0.000

2.496 2.758 2.758 3.067 3.067 3.069

0.000 0.001 0.001 0.987 0.987 0.000

3.114 3.618 3.618 3.647 3.647 3.869

TD-CAM-B3LYP f

[8]CPP 0.0 0.007 0.007 1.322 1.322 0.000 [10]CPP 0.000 1.415 1.415 0.275 0.275 0.000

TD-LC-DFTB

Exp.

E (eV)

f

E (eV)

f

E (eV)

3.506 4.177 4.177 4.758 4.758 4.792

0.000 1.595 1.595 0.020 0.020 0.000

3.663 4.282 4.282 5.002 5.002 5.138

0.000 1.152 1.152 0.000 0.000 0.000

2.64a 3.67a

3.731 4.191 4.191 4.767 4.767 4.838

0.000 2.245 2.245 0.000 0.000 0.000

3.949 4.355 4.355 4.905 4.905 5.263

0.000 2.106 2.106 0.000 0.000 0.000

3.65a

Data taken from ref 39.

4.2. Vertical Excitation Energies. The vertical excitation energies and oscillator strengths of the first six excited states of [8]CPP and [10]CPP are compared with the corresponding TD-DFT values obtained with the PBE, B3LYP, and CAMB3LYP functionals in Table 1. Among the first six excited states at the ground-state geometries, both molecules have two pairs of doubly degenerate states: S2/S3 and S4/S5. One of the pairs, which corresponds to transitions from HOMO to LUMO + 1 and LUMO + 2, is characterized by large oscillator strengths. The other one, corresponding to transitions from HOMO − 1 and HOMO − 2 to LUMO, is characterized by small or vanishing oscillator strengths. TD-DFTB, TD-PBE, and TD-B3LYP predict the same state order in [8]CPP, according to which bright degenerate states are placed above dark ones. TD-CAMB3LYP and LC-TD-DFTB invert this order, predicting the bright states to have lower energies. In the case of [10]CPP, TD-DFTB and TD-PBE predict intense transitions to degenerate S4 and S5 states, whereas the S2 and S3 states have low oscillator strengths. This order is inverted with the range-separated functionals and also with the hybrid functional, which predicts that the energies of the close-lying pairs of degenerate states differ only by 0.03 eV. In comparison to the experimental absorption peak positions, the most reliable values of vertical excitation energies are obtained with the B3LYP functional. The TD-CAM-B3LYP excitation energies overestimate excitation energies by ∼0.5 eV in comparison to the absorption band peaks. It has already been noticed that range-separated TD-DFT methods have a tendency to blue-shift excitation energies of conjugated molecules, compared to the experimental absorption peak positions.43,83 Compared to the experimental position of the tail of absorption spectra of [8]CPP39 (which arises from the excitation to the S1 state), the TD-DFTB underestimates the excitation energy by ∼0.2 eV. The TD-DFTB energies corresponding to the highest oscillator strength transitions to S4 and S5 states are underestimated by 0.45 ([8]CPP) and 0.52 eV ([10]CPP), compared to the positions of the band peaks in the experimental absorption spectra. In the case of [8]CPP, TD-DFTB predicts the same order of the excited states as that obtained with TD-B3LYP. In the case of [10]CPP, it inverts the

order of the close-lying double-degenerate states. The TDDFTB energies are consistent with TD-PBE, and the TD-LCDFTB energies are consistent with the corresponding TDCAM-B3LYP values. According to the TD-DFTB results, the energy of S4/S5 excitation is independent of the system size, and the S0/S1 transition is blue-shifted when the ring size is increased. These two features agree with the experimental results.39,46 The initial conditions created using a Wigner distribution have broken symmetries, compared to the ground-state equilibrium geometry. For this reason, transitions to the S6 state, which are forbidden at the equilibrium geometry, become possible in the ensemble of geometries. The oscillator strength of the S6 state becomes comparable to the oscillator strengths of the degenerate S4 and S5 states. These three states are closelying, and they can exchange order at different geometries. Thus, we selected initial conditions starting from the S5 and S6 adiabatic states, including all six states in the dynamics simulations. The algorithm for the initial condition selection in the Newton-X code computes the probabilities for selecting an initial condition based on the oscillator strength for the selected state, which ensures that, at all selected geometries, the S5 and S6 states are bright states. 4.3. Excited-State Population Dynamics. After the photoexcitation, the molecules undergo ultrafast internal conversion from the initially populated S6 and S5 states via the S4 state to the lower excited states. This process is shown in Figure 3, which reports the fraction of trajectories in each electronic state in each time step. For both molecules, after this initial period taking ∼100 fs, the population remains mostly shared between S1 (∼40%) and S2 + S3 (∼40%). The occupations of S4 and S5 + S6 remain at ∼10% each until the end of the simulations. The ground state is not reached in any of the trajectories during the 3 ps simulations. Generally, the trajectories have a tendency to relax to the S1. However, because of the lack of dissipative modes to the environment, once the molecules relax to their S1 states, part of the population is transferred back to higher excited states, when they reach regions of strong nonadiabatic couplings between S1 and S2. Furthermore, the S3 state can be easily repopulated from the close-lying S2 state. This recurrent back-transfer G


Article


phosphorescence spectra, they confirmed that this reaction produces singlet oxygen, and from the triplet oxygen phosphorescence decay profile, they determined the quantum yield of the singlet oxygen generation. They assumed that the quantum yield for singlet oxygen generation is equivalent to the quantum yield for intersystem crossing (ΦISC = ΦO2). The rationale behind this hypothesis was that, within the range of applied triplet oxygen concentrations, the reaction of triplet quenching with oxygen is several orders of magnitude faster than competing phosphorescence process, meaning that the efficiency of triplet quenching is 100%.41 The quantum yield of internal conversion was then estimated as ΦIC = 1 − ΦISC − Φfl. From the quantum yields of these three processes, one can conclude that internal conversion (ΦIC = 0.48, Table 2) and intersystem crossing (ΦISC = 0.44) are the main deactivation mechanisms in [8]CPP, whereas fluorescence has negligible quantum yield (Φfl = 0.08). In [10]CPP, fluorescence becomes the dominant relaxation mechanism with a quantum yield of 0.46, whereas the quantum yield of intersystem crossing is significantly decreased (ΦISC = 0.17), compared to its value in [8]CPP. Internal conversion remains an important deactivation pathway in [10]CPP, with a quantum yield of 0.37. A recent wave packet dynamics simulation84 of the relaxation process in CPPs within a triplet manifold, which included T3 − T1 states, revealed the ultrafast internal conversion to the T1 state. Strong Jahn−Teller interactions occurring between doubly degenerate T3 and T2 states, which induce this ultrafast relaxation, promote indirect intersystem crossing from S1 to T1 state (S1 → T3/T2 → T1). Although the intersystem crossing and phosphorescence were proven to be important, we do not consider this reaction channel, because our dynamics did not include triplet states, and we limit the discussion to fluorescence and internal conversion. The fluorescence lifetimes are estimated from the measured fluorescence rate constants39 (as τfl = 1/kfl). The fluorescence lifetime of [8]CPP (204.3 ns) is ∼15 times longer than of [10]CPP (14.3 ns) (see Table 2). This shortening with the molecular size can be rationalized based on the dynamics results. The distributions of the excitation energies and oscillator strengths in the S1 states of both molecules are shown in Figure 4. The S1 excitation energy (Ev) distribution is centered at 1.5 eV in [8]CPP and it increases to 1.8 eV in [10]CPP, indicating the local minima on the S1 surfaces. The S1 → S0 oscillator strength (f) in [10]CPP is distributed toward higher values, compared to the that observed in [8]CPP. Their mean values also increase with the molecular size, from f [8]CPP = 0.04 to f [10]CPP = 0.10. These two trends contribute to decreasing the fluorescence lifetime (τfl) in [10]CPP, in

Figure 3. Fractional occupations of the excited states of [8]CPP (left) and [10]CPP (right) during dynamics.

explains the oscillatory behavior of the average population of these states during the dynamics (see Figure 3). Surface hopping dynamics for [n]CPPs have been reported by Adamska et al.43 for several ring sizes (n) between 9 and 16. Their simulations were based on configuration interaction singles (CIS) on top of a semiempirical Hamiltonian, propagated for 500 fs. Similar to our results, they also found that internal conversion to the S1 state is reached on the ultrafast scale. However, in their case, once the S1 is populated, the trajectory remains there, and the S1 state reaches 100% occupation within 200 fs. The main reason for this difference is that Adamska and co-workers propagated the trajectories using constant-temperature Langevin dynamics with solvation,7 while we employed a micro-Canonical ensemble in the gas phase. The nuclear equations of motion, in their case, include a velocity-dependent friction term, which emulates the kinetic energy transfer to the solvent. There is still a third set of dynamics simulations for [6]CPP, [8]CPP, and [10]CPP reported in ref44. These simulations, based on wave packet propagation for 200 fs, were performed involving four excited states computed on a 26-dimensional model Hamiltonian parametrized at the TD-B3LYP level. 4.4. Decay Mechanisms. After the S1 state is reached during dynamics, three kinetically competitive decay processes from S1 to the ground state may occur: fluorescence, internal conversion, and intersystem crossing to the T1 state,41 followed by phosphorescence to the S0 state. Available experimental and computed quantum yields, rate constants, and lifetimes are summarized in Table 2. Fujitsuka et al. estimated fluorescence quantum yields, reaction rates, and lifetimes from fluorescence intensity decay experiments39 and determined triplet lifetimes from phosphorescence decay experiments.41 Moreover, they estimated quantum yields and rate constants of the S1/T1 intersystem crossing and S1/S0 internal conversion from experiments of triplet quenching with triplet oxygen. Based on the

Table 2. Experimental Quantum Yields (Φ), Decay Rates (k), and Lifetimes (τ) of Fluorescence (fl), Internal Conversion (IC), Intersystem Crossing (ISC), and Phosphorescence (ph) of [8]CPP and [10]CPP Quantum Yields

Decay Rates −1

−1

Lifetimes −1

molecule

Φfl

ΦIC

ΦISC

kfl (× 10 s )

kIC (× 10 s )

kISC (× 10 s )

τfl (ns)

τph (μs)

[8]CPP [10]CPP

0.084a 0.46a

0.475b 0.368b

0.440b 0.171b

0.48a,b 7.00a,b

2.7c 5.5c

2.5b 2.6b

208.3d 14.3d

60b 58b

7

7

7

a Data taken from ref 39. (The absorption spectra of CPPs were recorded in THF.) bData taken from ref 41. ΦISC is assumed to be the quantum yield of singlet oxygen (see text). (The experiments are performed in THF.) cComputed as kIC = (kfl + kISC)/(ΦIC−1 − 1) applying the experimental rate constants for fluorescence and intersystem crossing and experimental quantum yield for internal conversion. dInternal fluorescence lifetimes computed as τfl = 1/kfl.

H


Article


Figure 4. Distributions of the S1 excitation energies (left) and oscillator strengths (right) of [8]CPP and [10]CPP during dynamics in single trajectories. Figure 5. Relative activation of different normal modes at the geometry near the S2/S1 conical intersection of [8]CPP.

comparison to the [8]CPP, because τfl is proportional to f−1Ev−2.73 Using the mean values for the excitation energies and oscillator strengths, we may estimate that the ratio between fluorescence lifetimes is 3.6. This value is underestimated, compared to the experimental value (∼15), but correctly indicates the qualitative trend. The disagreement between these two values could be partially attributed to environmental effects, since the experimental rate constant measurements were done in tetrahydrofuran, whereas the simulations were performed in the gas phase. The fluorescence lifetime shortening from [8]CPP to [10]CPP should directly lead to an increase of the fluorescence quantum yield, because it hinders competing processes. The experimental ratio between the fluorescence quantum yields is 5.5 (see Table 2). The increase in the oscillator strength of the S1 state in [10]CPP, in comparison to [8]CPP, has been explained by decreasing the energy gap between S1 and S2 states, which increases the vibronic coupling between them44 and by localization of excitations in the S1 states in larger CPPs, which is accompanied by asymmetrical geometrical distortions leading to an increase in the oscillator strength.43 4.5. Vibrational Dynamics. We have analyzed the vibrational relaxation using normal-mode analysis.85,86 The geometries at the S2/S1 hopping point for all trajectories of [8]CPP were collected, and each of them was then aligned with the [8]CPP S0 minimum. The alignment of structures is achieved by removing the translational degrees of freedom by coinciding their centers of masses. The structures are then superimposed by rotation within their center of masses systems. The optimal rotation matrix that maximizes the overlap between two geometries is computed applying Kabsch’s algorithm.87 After being aligned, the Cartesian difference between the S2/S1 hopping geometry and the reference geometry was computed and projected on the normal mode basis calculated for the ground state. A map of the active normal modes at the S2/S1 hopping point was built by averaging the computed activations of normal modes over all S2/S1 hopping points. The relative activations of normal modes at the hopping points, with respect to the ground-state geometry, are shown in Figure 5. They identify the main vibrational motions driving the excited-state relaxation. The most active are low-wavenumber torsional modes (υ < 120 cm−1) and radial breathing modes (RBM, 150 cm−1 < υ < 220 cm−1). Several normal modes in the region between 230 cm−1 and 400 cm−1 feature small activation

levels. They involve bending of Cipso−Cipso bonds, with respect to the ring planes. The radial flexural modes (RFM, υ ≈ 500 cm−1) also show small activation levels. They correspond to ring deformations induced by Cipso−Cipso vibrations outward from the ring and Cortho−Cortho vibrations inward toward the ring. Modes with the wavenumbers between 800 cm−1 and 900 cm−1 have medium activation levels. They are mostly ring deformation modes with out-of-plane C−H bending. The modes with wavenumbers of ∼1200 cm−1, also with medium activation level, are in-plane C−H bending modes. In the region between 1430 cm−1 and 1790 cm−1, a series of C−C stretching modes show small activation. They are characterized in more detail below. CPPs, as well as carbon nanotubes, have three different types of C−C stretching modes: G1, G2, and Gio vibrations.88 The G1 and G2 modes correspond to the circumferential stretching of the ortho−ortho and ipso−ipso C−C bonds. In the G1 modes, all phenylene rings experience identical in-phase motions with the same magnitudes, whereas in the G2 modes, the adjacent phenylene rings oscillate out of phase. The G2 vibrations have higher frequencies, compared to the G1 mode. The Gio modes correspond to ipso−ortho C−C bond vibrations and have lower frequencies, compared to G1 and G2 modes. The Gio vibrations also have two components, corresponding to inphase and out-of-phase (we denote them as Gio′ ) motions in adjacent rings. Gio (1430−1480 cm−1) and G2 (∼1790 cm−1) modes show small activation levels at the S2/S1 crossing. Rigid excited-state potential energy profiles along all S2/S1 active normal modes discussed above are computed with TDDFTB at geometries obtained by elongation along the normal modes. For most of the modes, the excited states remain approximately parallel to each other along the displacement coordinate (see section SI-III in the Supporting Information). Nevertheless, there are three exceptions: the energy profiles along two Gio modes (ν170 = 1432.5 cm−1 and ν174 = 1448.7 cm−1) and the G2 mode (ν201 = 1786.0 cm−1). These three C− C stretching modes induce a decrease in the S2−S1 gap (see Figure 6). For example, the geometrical distortion along mode ν174 for 0.5ν174, with respect to the equilibrium geometry (where ν174 represents the normalized normal mode), reduces the S2 − S1 gap from 0.3 eV to 0.15 eV. The small S2 − S1 gaps indicate the vicinity of a conical intersection between these two states, which could be obtained by relaxing the other I


Article


4.6. Charge Transfer and Exciton (De)localization in the S1 State. The excitation energies, charge transfer (CT) numbers, and participation ratio (PR) computed with TDDFTB along a single trajectory of [8]CPP and [10]CPP, considering only time steps when the molecules are in the S1 states, are shown in Figure 7. To verify the reliability of TDDFTB in the description of energies, charge transfer, and localization properties of the excited states during the dynamics, we recomputed the excited states, CT, and PR numbers with TD-LC-DFTB at the same geometries. These values are also plotted in Figure 7. As expected from our analysis of the vertical excitation energies, the distributions of the excitation energies show significant differences between S1 energies computed with the two methods. The TD-DFTB energies are red-shifted, in comparison to the TD-LC-DFTB energies, by ∼1 eV. Nevertheless, not all disagreement between the results is caused by TD-DFTB inaccuracy. We have already discussed that range-separated functionals had a tendency to blue-shift the excitation energies of CPPs, with respect to the experimental absorption peak for ∼0.5 eV, whereas the TDDFTB had a tendency to red-shift these energies for a similar magnitude (see the discussion of Table 1).

Figure 6. Excitation energies of S1 − S6 states of [8]CPP at the geometries obtained by distorting the molecules along three C−C stretching modes, ν170 (Gio, 1432.5 cm−1), ν174 (Gio, 1448.7 cm−1), and ν201 (G2, 1786.0 cm−1). The x-axis represents a fraction of a normalized normal mode for which the geometry is distorted, with respect to the equilibrium geometry.

coordinates along the stretching modes. The activation of ν170 and ν174 also induces S3/S4 and S5/S6 crossings. These series of state crossings induced by Gio and G2 imply that these modes play a central role for the relaxation of the excited states down to S1, and for the eventual repopulation of the higher excited states.

Figure 7. Distributions of the excitation energies (left), charge transfer numbers (middle), and participation ratio numbers (right) in the S1 state of [8]CPP (top of the figure) and [10]CPP (bottom of the figure), during dynamics of a single trajectory considering the TD-(LC)-DFTB and SD(LC)-DFTB methods. J


Article


We wanted to verify how the single-determinant (SD) approximation of the excited-state wave functions10 affects the description of their localization and CT properties. We represented the S1 states in the form of single determinants, selecting the one with the largest linear-response coefficients obtained either with TD-DFTB (rendering SD-DFTB) or TDLC-DFTB (rendering SD-LC-DFTB). The excitation energies in the SD approximation correspond to simple KS orbital energy gaps. CT and PR were computed with the renormalized SD wave functions. The three quantities for [8]CPP and [10]CPP are shown in Figure 7. In both cases, the energies of the S1 state at the SD-DFTB approximation are very close to the TD-DFTB values, whereas the SD-LC-DFTB excitation energies are overestimated by ∼3 eV in comparison to the LCTD-DFTB values. CT and PR are both insensitive to range separation, and SD-DFTB and SD-LC-DFTB provide very similar distributions. In the SD approximation, CT has a sharp distribution peaked at ∼0.85. (This result is similar to that obtained with TD-DFTB). However, the results for PR are somewhat surprising. For both molecules, PR distributions obtained with both SD approximations are in reasonable agreement with that from TD-LC-DFTB, our highest level. Some insights into the differences between compositions of the TD-DFTB and TD-LC-DFTB wave functions that induce distinct CT and exciton localization properties can be gained by analyzing their PRNTO distributions (Figure 8). The PRNTO

Moreover, TD-DFTB overestimates CT numbers of the S1 state of both molecules, in comparison with the CT numbers obtained from TD-LC-DFTB computations. According to the TD-LC-DFTB, the S1 states are obtained by mixing of Frenkel and CT excitons, with the CT distributions peaked at ∼0.65. On the other hand, TD-DFTB predicts almost pure CT character of the S1 states, with CT distributions peaked at ∼0.9 in both molecules. This indicates that TD-DFTB, similar to that observed for TD-DFT at the level of GGA functionals, overestimates the CT character of excited states. It has been confirmed, for instance, on a set of conjugated systems89 that TD-DFT with the GGA-level functionals predicts spurious CT character for the lowest excited state, whereas predominantly local character is obtained at the ADC(2) level and with TDDFT with hybrid and long-range-corrected functionals. The PR distributions calculated with TD-DFTB indicate that the S1 states of both molecules are delocalized over 8 fragments in [8]CPP and ∼9 fragments in [10]CPP (see Figure 7). The inclusion of long-range correction to the functional leads to some localization of excitons, in comparison to the TD-DFTB prediction, with the S1 state of [8]CPP mostly delocalized over 6−8 fragments and that of [10]CPP showing a broad distribution between 4 and 9 fragments and peaking at 6. (It has been already observed that the localization of excitons can be adequately described only with range-corrected hybrid DFT models.90,91) According to the TD-LC-DFTB, the broad PR distribution for [10]CPP implies that, in this molecule, both localized and delocalized excitons occur in the first excited state. The maximum of the PR distribution peaked at 6 fragments indicates a more pronounced exciton localization than in the case of the smaller molecule. This result is in agreement with the dynamics results of Adamska et al.,43 who noticed that a localization of excitons, absent in [9]CPP, occurs in larger CPPs (n ≥ 12), where they are localized on 5−7 rings. This exciton localization is attributed to its self-trapping induced by geometrical distortions. In addition, Nishihara et al.38 observed, based on the temperature dependence of the fluorescence lifetime, that the excitons of larger CPPs (n = 12) are localized at temperatures less than room temperature and delocalized at higher temperatures. Note that, even though, based on energy criterium, the B3LYP functional provides the most reliable description on the excited states of CPPs, it does not correctly describe the delocalization properties of the exciton in the S1 state, which is an important feature of dynamics. (See also discussion in ref 92.) We have optimized the geometries corresponding to the localized exciton in [10]CPP, starting from several different geometries taken from dynamics at which the system was in the S1 state and the exciton was localized. The optimizations were performed with TD-B3LYP and TD-CAMB3LYP, using the 631G basis set and with a small optimization radius. TDCAMB3LYP yielded nonsymmetric structure and localized S1 state (PR ≈ 7), whereas TD-B3LYP always converged to the completely symmetrical structure and completely delocalized S1 state (PR ≈ 10). This finding indicates that TD-B3LYP predicts that the localization of excitons in the S1 state is not the most favorable energetically, and probably it would be shown in dynamics performed based on TD-B3LYP energies and gradients. Most probably in the TD-B3LYP dynamics, the system would converge to the symmetrical structures and delocalized excitons in the S1 state.

Figure 8. Distributions of the PRNTO of the S1 states reached during dynamics in single trajectories of [8]CPP (left) and [10]CPP (right).

distributions show that the TD-DFTB wave functions, built via eq 13, are represented mainly by a single transition, whereas the TD-LC-DFTB wave functions are represented by more than one important transition between electron/hole NTO pairs. Thus, it is clear that the LC correction not only changes the orbitals (via the action of the HF exchange) but also changes the distribution of determinants contributing to the state description. Based on this constatation, we are in a position to understand the relationships between ΔE, CT, and PR with the different levels. TD-DFTB and SD-DFTB do not show a significant difference in the excitation energies, because the TD-DFTB wave functions are dominated by one transition. In addition, the two-electron integrals in eq 10 (which are included in TDDFTB but not SD-DFTB) are usually small, compared to the HOMO − LUMO gap for local/semilocal functionals. In contrast, when LC corrections are considered, the SD approximation is an inadequate representation of the TD-LCDFTB result, composed of many determinants. For functionals K


Article


(2) The experimentally observed trend of increasing the fluorescence rate constant in [10]CPP in comparison to [8]CPP is reproduced by the TD-DFTB dynamics. It originates from a simultaneous increase of the S0 − S1 energy gap and oscillator strength in [10]CPP, in comparison to its values in [8]CPP. (3) The vibrational dynamics driving the relaxation in the manifold of excited states is enabled by low-amplitude activation of C−C stretching modes, which induces a series of state crossings. These same modes may be responsible for the back transfer in the gas phase. By performing the analysis of the S1 states based on their transition density matrices, we concluded that the correct description of the charge-transfer and localization properties of the S1 states can be achieved only with long-range corrected (LC) TD-DFTB. According to TD-LC-DFTB, the S1 states are mixtures of the Frenkel and CT excitons, with the CT numbers peaking at ∼0.65. The TD-DFTB method has a tendency to overestimate their CT character, predicting the S1 states to have almost pure CT character. According to the TD-LC-DFTB transition densities, the S1 state is delocalized in the smaller molecule, whereas it is localized in the larger one. The feature of density localization in larger [n]CPPs (n ≥ 10), which has been observed earlier from the results of dynamics simulation43 and also experimentally, based on the temperature dependence of fluorescence lifetime38 (for n = 12), could not be reproduced from the TD-DFTB transition densities. We have also tested the performance of single-determinant (SD) approximations of the wave function with DFTB and LCDFTB. The results show that SD-DFTB provides a description of the excited states with quality similar to that of TD-DFTB. However, SD-LC-DFTB deteriorates the excitation energies. The density localization is well-described in both cases. Nevertheless, because this localization is based on a singledeterminant representation, it may fail in cases where the actual density is delocalized due to multideterminantal effects.

involving a fraction of Hartree−Fock exchange, the mentioned two-electron integrals are typically large. These two effects lead to the strong deviation in the excitation energy between the two levels. The CT distributions computed based on TD-DFTB and SD-DFTB both predict almost-pure CT character of the S1 states, implying that the highest-weight transitions between electron/hole NTO pairs are the CT transitions. The inclusion of the LC correction in TD-LC-DFTB reduces the CT value through the multideterminantal contribution to the density. Nevertheless, the LC correction in SD-DFTB does not change the CT values, in comparison to SD-DFTB, because, by construction, only one determinant is used. The situation of PR is more subtle than that of ΔE and CT. Although the TD-DFTB state is dominated by a single determinant, other determinants still contribute to the final result. In addition, these secondary contributions are responsible for donating the delocalized character to the state, reflected on the large PR. In SD-DFTB, however, only the main determinant is considered, resulting in more-localized densities and smaller PR. The inclusion of LC correction in TD-LC-DFTB leads to the same effectPR reductionbut for a different reason, because, now, the reduced PR is obtained as result of a multideterminantal composition. As far as we can see, the good agreement of PR computed with TD-LC-DFTB and SD-DFTB was just coincidence. If, for instance, the multideterminantal TD-LC-DFTB results were delocalized, the SD-DFTB would still predict a small PR. Once more, SD-LCDFTB and SD-DFTB show similar PR values, because only a single determinant is analyzed.

5. CONCLUSIONS In the present study, we implemented the decoherencecorrected fewest-switches surface hopping based on LR TDDFTB, enhanced by transition density analysis through interfaces between Newton-X, DFTB+, and TheoDORE programs. This implementation was applied to study the excited states and relaxation dynamics of two cycloparaphenylenes: [8]CPP and [10]CPP. The vertical excitation energies of [8]CPP and [10]CPP computed with the TD-DFTB and LC-TD-DFTB were compared with the TD-DFT energies obtained with the PBE, B3LYP, and CAM-B3LYP functionals. TD-DFTB predicts the same ordering of the states as TD-B3LYP for [8]CPP, which gives the closest excitation energies to the experimental values. For [10]CPP, the order of the close-lying pairs of degenerate states obtained with TD-DFTB does not correspond to that predicted with TD-B3LYP. From the comparison of energies, we can conclude that TD-DFTB systematically underestimates the excitation energies by ∼ 0.5 eV, which is a typical feature of DFT methods at the GGA functional level. However, it describes two important features of the absorption spectra well: the most intense peaks are independent of system size, and the S0 − S1 energy gap is blue-shifted with the increase of system size. Based on the results of the TD-DFTB surface hopping dynamics and normal-mode analysis, we can draw the following conclusions: (1) Both molecules relax to S1 within 100 fs. Up to 3 ps, there is no indication of further relaxation to S0. In the gas phase, repopulation of higher excited states, especially S2 and S3, does not allow for full population transfer to S1.

■

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jctc.7b01000. Cartesian coordinates, simulated spectrum, and potential energy curves (PDF)

■

AUTHOR INFORMATION

Corresponding Authors

* E-mail: [email protected]. *E-mail: [email protected]. Website: www.barbatti. org. ORCID

Ljiljana Stojanović: 0000-0003-2821-5110 Felix Plasser: 0000-0003-0751-148X Thomas A. Niehaus: 0000-0001-9576-7658 Mario Barbatti: 0000-0001-9336-6607 Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS This Project was funded by the Deanship of Scientific Research (DSR) King Abdulaziz University, Jeddah, under Grant No. 43130-35-RG. Therefore, the authors acknowledge DSR support L


Article


functional tight-binding method for simulations of complex materials properties. Phys. Rev. B: Condens. Matter Mater. Phys. 1998, 58, 7260− 7268. (17) Aradi, B.; Hourahine, B.; Frauenheim, T. DFTB+, a Sparse Matrix-Based Implementation of the DFTB Method. J. Phys. Chem. A 2007, 111, 5678−5684. (18) Porezag, D.; Frauenheim, T.; Kö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−12957. (19) Seifert, G.; Porezag, D.; Frauenheim, T. Calculations of molecules, clusters, and solids with a simplified LCAO-DFT-LDA scheme. Int. J. Quantum Chem. 1996, 58, 185−192. (20) Rappe, A. K.; Casewit, C. J.; Colwell, K. S.; Goddard, W. A.; Skiff, W. M. UFF, a full periodic table force field for molecular mechanics and molecular dynamics simulations. J. Am. Chem. Soc. 1992, 114, 10024−10035. (21) Elstner, M.; Hobza, P.; Frauenheim, T.; Suhai, S.; Kaxiras, E. Hydrogen bonding and stacking interactions of nucleic acid base pairs: A density-functional-theory based treatment. J. Chem. Phys. 2001, 114, 5149−5155. (22) Grimme, S.; Antony, J.; Ehrlich, S.; Krieg, H. A consistent and accurate ab initio parametrization of density functional dispersion correction (DFT-D) for the 94 elements H-Pu. J. Chem. Phys. 2010, 132, 154104. (23) 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. (24) Lutsker, V.; Aradi, B.; Niehaus, T. A. Implementation and benchmark of a long-range corrected functional in the density functional based tight-binding method. J. Chem. Phys. 2015, 143, 184107. (25) Kranz, J. J.; Elstner, M.; Aradi, B.; Frauenheim, T.; Lutsker, V.; Garcia, A. D.; Niehaus, T. A. Time-Dependent Extension of the LongRange Corrected Density Functional Based Tight-Binding Method. J. Chem. Theory Comput. 2017, 13, 1737−1747. (26) Barbatti, M.; Granucci, G.; Ruckenbauer, M.; Plasser, F.; Crespo-Otero, R.; Pittner, J.; Persico, M.; Lischka, H. NEWTON-X: A package for Newtonian Dynamics Close to the Crossing Seam (v. 2.0). Available via the Internet at www.newtonx.org, accessed 2017. (27) Barbatti, M.; Ruckenbauer, M.; Plasser, F.; Pittner, J.; Granucci, G.; Persico, M.; Lischka, H. Newton-X: A Surface-Hopping Program for Nonadiabatic Molecular Dynamics. WIREs: Comp. Mol. Sci. 2014, 4, 26−33. (28) Jakowski, J.; Morokuma, K. Liouville−von Neumann molecular dynamics. J. Chem. Phys. 2009, 130, 224106. (29) Lei, Y.; Yuan, S.; Dou, Y.; Wang, Y.; Wen, Z. Detailed Dynamics of the Nonradiative Deactivation of Adenine: A Semiclassical Dynamics Study. J. Phys. Chem. A 2008, 112, 8497−8504. (30) Maitra, N. T. On correlated electron-nuclear dynamics using time-dependent density functional theory. J. Chem. Phys. 2006, 125, 014110. (31) Yang, Y.; Yu, H.; York, D.; Cui, Q.; Elstner, M. Extension of the Self-Consistent-Charge Density-Functional Tight-Binding Method: Third-Order Expansion of the Density Functional Theory Total Energy and Introduction of a Modified Effective Coulomb Interaction. J. Phys. Chem. A 2007, 111, 10861−10873. (32) 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. (33) Ryabinkin, I. G.; Nagesh, J.; Izmaylov, A. F. Fast Numerical Evaluation of Time-Derivative Nonadiabatic Couplings for Mixed Quantum−Classical Methods. J. Phys. Chem. Lett. 2015, 6, 4200−4203. (34) Plasser, F. TheoDORE 1.4: A package for theoretical density, orbital relaxation, and exciton analysis, http://theodore-qc.sourceforge. net, 2017.

for Scientific Research. Also, the authors appreciate the kind cooperation provided by the Deanship of Scientific Research (DSR), King Abdulaziz University. The work of M.B. and L.S. was supported by Excellence Initiative of Aix-Marseille University (A*MIDEX) and the project Equip@Meso (ANR10-EQPX-29-01), both funded by the French Government “Investissements d’Avenir” program. M.B. also acknowledges funding from HPC resources from GENCI-CINES (Grant No. 2017-A0010810012). T.N. would like to thank the Laboratoire d’Excellence iMUST for financial support.

■

REFERENCES

(1) Zimmerman, P. M.; Zhang, Z.; Musgrave, C. B. Singlet fission in pentacene through multi-exciton quantum states. Nat. Chem. 2010, 2, 648−652. (2) Bredas, J. L.; Beljonne, D.; Coropceanu, V.; Cornil, J. Chargetransfer and energy-transfer processes in pi-conjugated oligomers and polymers: A molecular picture. Chem. Rev. 2004, 104, 4971−5003. (3) Thorsmølle, V. K.; Averitt, R. D.; Demsar, J.; Smith, D. L.; Tretiak, S.; Martin, R. L.; Chi, X.; Crone, B. K.; Ramirez, A. P.; Taylor, A. J. Morphology Effectively Controls Singlet-Triplet Exciton Relaxation and Charge Transport in Organic Semiconductors. Phys. Rev. Lett. 2009, 102, 017401. (4) Landry, B. R.; Subotnik, J. E. Quantifying the Lifetime of Triplet Energy Transfer Processes in Organic Chromophores: A Case Study of 4-(2-Naphthylmethyl)benzaldehyde. J. Chem. Theory Comput. 2014, 10, 4253−4263. (5) Kilina, S.; Kilin, D.; Tretiak, S. Light-Driven and Phonon-Assisted Dynamics in Organic and Semiconductor Nanostructures. Chem. Rev. 2015, 115, 5929−5978. (6) Akimov, A. V.; Prezhdo, O. V. The PYXAID Program for NonAdiabatic Molecular Dynamics in Condensed Matter Systems. J. Chem. Theory Comput. 2013, 9, 4959−4972. (7) Nelson, T.; Fernandez-Alberti, S.; Chernyak, V.; Roitberg, A. E.; Tretiak, S. Nonadiabatic Excited-State Molecular Dynamics Modeling of Photoinduced Dynamics in Conjugated Molecules. J. Phys. Chem. B 2011, 115, 5402−5414. (8) Jain, A.; Alguire, E.; Subotnik, J. E. An Efficient, Augmented Surface Hopping Algorithm That Includes Decoherence for Use in Large-Scale Simulations. J. Chem. Theory Comput. 2016, 12, 5256− 5268. (9) Wang, L.; Long, R.; Prezhdo, O. V. Time-Domain Ab Initio Modeling of Photoinduced Dynamics at Nanoscale Interfaces. Annu. Rev. Phys. Chem. 2015, 66, 549−579. (10) Pal, S.; Trivedi, D. J.; Akimov, A. V.; Aradi, B.; Frauenheim, T.; Prezhdo, O. V. Nonadiabatic Molecular Dynamics for Thousand Atom Systems: A Tight-Binding Approach toward PYXAID. J. Chem. Theory Comput. 2016, 12, 1436−1448. (11) Mitrić, R.; Werner, U.; Wohlgemuth, M.; Seifert, G.; BonačićKoutecký, V. Nonadiabatic Dynamics within Time-Dependent Density Functional Tight Binding Method. J. Phys. Chem. A 2009, 113, 12700− 12705. (12) Fabiano, E.; Keal, T. W.; Thiel, W. Implementation of surface hopping molecular dynamics using semiempirical methods. Chem. Phys. 2008, 349, 334−347. (13) Granucci, G.; Persico, M.; Toniolo, A. Direct semiclassical simulation of photochemical processes with semiempirical wave functions. J. Chem. Phys. 2001, 114, 10608−10615. (14) Gao, X.; Geng, H.; Peng, Q.; Ren, J.; Yi, Y.; Wang, D.; Shuai, Z. Nonadiabatic Molecular Dynamics Modeling of the Intrachain Charge Transport in Conjugated Diketopyrrolo-pyrrole Polymers. J. Phys. Chem. C 2014, 118, 6631−6640. (15) Gao, X.; Peng, Q.; Niu, Y.; Wang, D.; Shuai, Z. Theoretical insight into the aggregation induced emission phenomena of diphenyldibenzofulvene: A nonadiabatic molecular dynamics study. Phys. Chem. Chem. Phys. 2012, 14, 14207−14216. (16) Elstner, M.; Porezag, D.; Jungnickel, G.; Elsner, J.; Haugk, M.; Frauenheim, T.; Suhai, S.; Seifert, G. Self-consistent-charge densityM


Article

Journal of Chemical Theory and Computation (35) Hartland, G. V.; Besteiro, L. V.; Johns, P.; Govorov, A. O. What’s so Hot about Electrons in Metal Nanoparticles? ACS Energy Lett. 2017, 2, 1641−1653. (36) Grimme, S.; Bannwarth, C.; Shushkov, P. A Robust and Accurate Tight-Binding Quantum Chemical Method for Structures, Vibrational Frequencies, and Noncovalent Interactions of Large Molecular Systems Parametrized for All spd-Block Elements (Z = 1−86). J. Chem. Theory Comput. 2017, 13, 1989−2009. (37) Jasti, R.; Bhattacharjee, J.; Neaton, J. B.; Bertozzi, C. R. Synthesis, Characterization, and Theory of [9]-, [12]-, and [18]Cycloparaphenylene: Carbon Nanohoop Structures. J. Am. Chem. Soc. 2008, 130, 17646−17647. (38) Nishihara, T.; Segawa, Y.; Itami, K.; Kanemitsu, Y. Exciton recombination dynamics in nanoring cycloparaphenylenes. Chem. Sci. 2014, 5, 2293−2296. (39) Fujitsuka, M.; Cho, D. W.; Iwamoto, T.; Yamago, S.; Majima, T. Size-dependent fluorescence properties of [n]cycloparaphenylenes (n = 8−13), hoop-shaped [π]-conjugated molecules. Phys. Chem. Chem. Phys. 2012, 14, 14585−14588. (40) Fujitsuka, M.; Iwamoto, T.; Kayahara, E.; Yamago, S.; Majima, T. Enhancement of the Quinoidal Character for Smaller [n]Cycloparaphenylenes Probed by Raman Spectroscopy. ChemPhysChem 2013, 14, 1570−1572. (41) Fujitsuka, M.; Lu, C.; Iwamoto, T.; Kayahara, E.; Yamago, S.; Majima, T. Properties of Triplet-Excited [n]Cycloparaphenylenes (n = 8−12): Excitation Energies Lower than Those of Linear Oligomers and Polymers. J. Phys. Chem. A 2014, 118, 4527−4532. (42) Camacho, C.; Niehaus, T. A.; Itami, K.; Irle, S. Origin of the size-dependent fluorescence blueshift in [n]cycloparaphenylenes. Chem. Sci. 2013, 4, 187−195. (43) Adamska, L.; Nayyar, I.; Chen, H.; Swan, A. K.; Oldani, N.; Fernandez-Alberti, S.; Golder, M. R.; Jasti, R.; Doorn, S. K.; Tretiak, S. Self-Trapping of Excitons, Violation of Condon Approximation, and Efficient Fluorescence in Conjugated Cycloparaphenylenes. Nano Lett. 2014, 14, 6539−6546. (44) Reddy, V. S.; Camacho, C.; Xia, J.; Jasti, R.; Irle, S. Quantum Dynamics Simulations Reveal Vibronic Effects on the Optical Properties of [n]Cycloparaphenylenes. J. Chem. Theory Comput. 2014, 10, 4025−4036. (45) Omachi, H.; Nakayama, T.; Takahashi, E.; Segawa, Y.; Itami, K. Initiation of carbon nanotube growth by well-defined carbon nanorings. Nat. Chem. 2013, 5, 572−576. (46) Segawa, Y.; Fukazawa, A.; Matsuura, S.; Omachi, H.; Yamaguchi, S.; Irle, S.; Itami, K. Combined experimental and theoretical studies on the photophysical properties of cycloparaphenylenes. Org. Biomol. Chem. 2012, 10, 5979−5984. (47) Lukes, V.; Solc, R.; Barbatti, M.; Elstner, M.; Lischka, H.; Kauffmann, H.-F. Torsional potentials and full-dimensional simulation of electronic absorption and fluorescence spectra of para-phenylene oligomers using the semiempirical self-consistent charge densityfunctional tight binding approach. J. Chem. Phys. 2008, 129, 164905− 12. (48) Hines, D. A.; Darzi, E. R.; Jasti, R.; Kamat, P. V. Carbon Nanohoops: Excited Singlet and Triplet Behavior of [9]- and [12]Cycloparaphenylene. J. Phys. Chem. A 2014, 118, 1595−1600. (49) Dreuw, A.; Head-Gordon, M. Single-Reference ab Initio Methods for the Calculation of Excited States of Large Molecules. Chem. Rev. 2005, 105, 4009−4037. (50) Tully, J. C. Molecular-Dynamics with Electronic-Transitions. J. Chem. Phys. 1990, 93, 1061−1071. (51) Granucci, G.; Persico, M. Critical Appraisal of the Fewest Switches Algorithm for Surface Hopping. J. Chem. Phys. 2007, 126, 134114−11. (52) Hammes-Schiffer, S.; Tully, J. C. Proton transfer in solution: Molecular dynamics with quantum transitions. J. Chem. Phys. 1994, 101, 4657−4667. (53) Barbatti, M.; Pittner, J.; Pederzoli, M.; Werner, U.; Mitrić, R.; Bonačić-Koutecký, V.; Lischka, H. Non-Adiabatic Dynamics of

Pyrrole: Dependence of Deactivation Mechanisms on the Excitation Energy. Chem. Phys. 2010, 375, 26−34. (54) Stojanović, L.; Bai, S.; Nagesh, J.; Izmaylov, A.; Crespo-Otero, R.; Lischka, H.; Barbatti, M. New Insights into the State Trapping of UV-Excited Thymine. Molecules 2016, 21, 1603. (55) Plasser, F.; Granucci, G.; Pittner, J.; Barbatti, M.; Persico, M.; Lischka, H. Surface hopping dynamics using a locally diabatic formalism: Charge transfer in the ethylene dimer cation and excited state dynamics in the 2-pyridone dimer. J. Chem. Phys. 2012, 137, 22A514−13. (56) 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. (57) Darghouth, A.; Casida, M.; Taouali, W.; Alimi, K.; Ljungberg, M.; Koval, P.; Sánchez-Portal, D.; Foerster, D. Assessment of DensityFunctional Tight-Binding Ionization Potentials and Electron Affinities of Molecules of Interest for Organic Solar Cells Against FirstPrinciples GW Calculations. Computation 2015, 3, 616. (58) Casida, M. In Recent Advances in Density Functional Methods, Part I; Chong, D., Ed.; World Scientific: Singapore, 1995; p 155. (59) Plasser, F.; Crespo-Otero, R.; Pederzoli, M.; Pittner, J.; Lischka, H.; Barbatti, M. Surface Hopping Dynamics with Correlated SingleReference Methods: 9H-Adenine as a Case Study. J. Chem. Theory Comput. 2014, 10, 1395−1405. (60) Mitrić, R.; Werner, U.; Bonačić-Koutecký, V. Nonadiabatic dynamics and simulation of time resolved photoelectron spectra within time-dependent density functional theory: Ultrafast photoswitching in benzylideneaniline. J. Chem. Phys. 2008, 129, 164118. (61) Tapavicza, E.; Tavernelli, I.; Rothlisberger, U. Trajectory surface hopping within linear response time-dependent density-functional theory. Phys. Rev. Lett. 2007, 98, 023001−4. (62) Löwdin, P.-O. Quantum Theory of Many-Particle Systems. I. Physical Interpretations by Means of Density Matrices, Natural SpinOrbitals, and Convergence Problems in the Method of Configurational Interaction. Phys. Rev. 1955, 97, 1474. (63) Plasser, F.; Lischka, H. Analysis of Excitonic and Charge Transfer Interactions from Quantum Chemical Calculations. J. Chem. Theory Comput. 2012, 8, 2777−2789. (64) Plasser, F.; Bäppler, S. A.; Wormit, M.; Dreuw, A. New tools for the systematic analysis and visualization of electronic excitations. II. Applications. J. Chem. Phys. 2014, 141, 024107. (65) Plasser, F.; Wormit, M.; Dreuw, A. New tools for the systematic analysis and visualization of electronic excitations. I. Formalism. J. Chem. Phys. 2014, 141, 024106. (66) Mewes, S. A.; Mewes, J.-M.; Dreuw, A.; Plasser, F. Excitons in poly(para phenylene vinylene): a quantum-chemical perspective based on high-level ab initio calculations. Phys. Chem. Chem. Phys. 2016, 18, 2548−2563. (67) Panda, A. N.; Plasser, F.; Aquino, A. J. A.; Burghardt, I.; Lischka, H. Electronically Excited States in Poly(p-phenylenevinylene): Vertical Excitations and Torsional Potentials from High-Level Ab Initio Calculations. J. Phys. Chem. A 2013, 117, 2181−2189. (68) Martin, R. L. Natural transition orbitals. J. Chem. Phys. 2003, 118, 4775−4777. (69) Luzanov, A. V.; Sukhorukov, A. A.; Umanskii, V. É. Application of transition density matrix for analysis of excited states. Theor. Exp. Chem. 1976, 10, 354−361. (70) Plasser, F. Entanglement entropy of electronic excitations. J. Chem. Phys. 2016, 144, 194107. (71) Wigner, E. On the Quantum Correction For Thermodynamic Equilibrium. Phys. Rev. 1932, 40, 749−759. (72) Gaus, M.; Goez, A.; Elstner, M. Parametrization and Benchmark of DFTB3 for Organic Molecules. J. Chem. Theory Comput. 2013, 9, 338−354. (73) Crespo-Otero, R.; Barbatti, M. Spectrum Simulation and Decomposition with Nuclear Ensemble: Formal Derivation and Application to Benzene, Furan and 2-Phenylfuran. Theor. Chem. Acc. 2012, 131, 1237. N


Article

Journal of Chemical Theory and Computation (74) Levine, B. G.; Ko, C.; Quenneville, J.; Martínez, T. J. Conical intersections and double excitations in time-dependent density functional theory. Mol. Phys. 2006, 104, 1039−1051. (75) Swope, W. C.; Andersen, H. C.; Berens, P. H.; Wilson, K. R. A Computer-Simulation Method for the Calculation of EquilibriumConstants for the Formation of Physical Clusters of Molecules Application to Small Water Clusters. J. Chem. Phys. 1982, 76, 637− 649. (76) Butcher, J. C. A Modified Multistep Method for the Numerical Integration of Ordinary Differential Equations. J. Assoc. Comput. Mach. 1965, 12, 124−135. (77) Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized gradient approximation made simple. Phys. Rev. Lett. 1996, 77, 3865−3868. (78) Becke, A. D. Density-functional thermochemistry. III. The role of exact exchange. J. Chem. Phys. 1993, 98, 5648−5652. (79) Lee, C.; Yang, W.; Parr, R. G. Development of the Colle-Salvetti correlation-energy formula into a functional of the electron density. Phys. Rev. B: Condens. Matter Mater. Phys. 1988, 37, 785−789. (80) Yanai, T.; Tew, D. P.; Handy, N. C. A New Hybrid ExchangeCorrelation Functional Using the Coulomb-Attenuating Method (CAM-B3LYP). Chem. Phys. Lett. 2004, 393, 51−57. (81) Schmidt, M. W.; Baldridge, K. K.; Boatz, J. A.; Elbert, S. T.; Gordon, M. S.; Jensen, J. H.; Koseki, S.; Matsunaga, N.; Nguyen, K. A.; Su, S.; et al. General atomic and molecular electronic structure system. J. Comput. Chem. 1993, 14, 1347−1363. (82) Gordon, M. S.; Schmidt, M. W. Advances in electronic structure theory: GAMESS a decade later. In Theory and Applications of Computational Chemistry: The First Forty Years, Dykstra, C. E.; Frenking, G.; Kim, K. S.; Scuseria, G. E., Eds.; Elsevier: Amsterdam, 2005; Chapter 41, pp 1167−1189. (83) Nayyar, I. H.; Batista, E. R.; Tretiak, S.; Saxena, A.; Smith, D. L.; Martin, R. L. Role of Geometric Distortion and Polarization in Localizing Electronic Excitations in Conjugated Polymers. J. Chem. Theory Comput. 2013, 9, 1144−1154. (84) Reddy, V. S.; Irle, S. Indirect Intersystem Crossing (S1 → T3/ T2 → T1) Promoted by the Jahn−Teller Effect in Cycloparaphenylenes. J. Chem. Theory Comput. 2017, 13, 4944−4949. (85) Plasser, F.; Barbatti, M.; Aquino, A. J. A.; Lischka, H. ExcitedState Diproton Transfer in [2,2′-Bipyridyl]-3,3′-diol: the Mechanism Is Sequential, Not Concerted. J. Phys. Chem. A 2009, 113, 8490−8499. (86) Kurtz, L.; Hofmann, A.; de Vivie-Riedle, R. Ground state normal mode analysis: Linking excited state dynamics and experimental observables. J. Chem. Phys. 2001, 114, 6151−6159. (87) Kabsch, W. A solution for the best rotation to relate two sets of vectors. Acta Crystallogr., Sect. A: Cryst. Phys., Diffr., Theor. Gen. Crystallogr. 1976, 32, 922−923. (88) Pena-Alvarez, M.; Qiu, L.; Taravillo, M.; Baonza, V. G.; Delgado, M. C. R.; Yamago, S.; Jasti, R.; Navarrete, J. T. L.; Casado, J.; Kertesz, M. From linear to cyclic oligoparaphenylenes: electronic and molecular changes traced in the vibrational Raman spectra and reformulation of the bond length alternation pattern. Phys. Chem. Chem. Phys. 2016, 18, 11683−11692. (89) Mewes, S. A.; Plasser, F.; Dreuw, A. Universal Exciton Size in Organic Polymers is Determined by Nonlocal Orbital Exchange in Time-Dependent Density Functional Theory. J. Phys. Chem. Lett. 2017, 8, 1205−1210. (90) Körzdörfer, T.; Sears, J. S.; Sutton, C.; Brédas, J.-L. Long-range corrected hybrid functionals for π-conjugated systems: Dependence of the range-separation parameter on conjugation length. J. Chem. Phys. 2011, 135, 204107. (91) Nayyar, I. H.; Batista, E. R.; Tretiak, S.; Saxena, A.; Smith, D. L.; Martin, R. L. Localization of Electronic Excitations in Conjugated Polymers Studied by DFT. J. Phys. Chem. Lett. 2011, 2, 566−571. (92) Sundholm, D.; Taubert, S.; Pichierri, F. Phys. Chem. Chem. Phys. 2010, 12, 2751−2757.

O


Nonadiabatic Dynamics of Cycloparaphenylenes with TD-DFTB

Recommend Documents