Temperature Dependence of Radiative and Nonradiative Rates from

General Approach to Coupled Reactive Smoluchowski Equations: Integration and Application of Discrete Variable Representation and Generalized Coordinat...
18 downloads 0 Views 2MB Size
Article pubs.acs.org/JCTC

Temperature Dependence of Radiative and Nonradiative Rates from Time-Dependent Correlation Function Methods Shiladitya Banerjee,*,†,‡ Alberto Baiardi,† Julien Bloino,†,‡ and Vincenzo Barone† †

Scuola Normale Superiore, piazza dei Cavalieri 7, I-56126 Pisa, Italy Consiglio Nazionale delle Ricerche, Istituto di Chimica dei Composti OrganoMetallici (ICCOM-CNR), UOS di Pisa, Area della Ricerca CNR, Via G. Moruzzi 1, I-56124 Pisa, Italy



ABSTRACT: The temperature dependence of the rate constants in radiative and nonradiative decays from excited electronic states has been studied using a time-dependent correlation function approach in the framework of the adiabatic representation and the harmonic oscillator approximation. The present work analyzes the vibrational aspect of the processes, which gives rise to the temperature dependence, with the inclusion of mode-mixing, as well as of frequency change effects. The temperature dependence of the rate constants shows a contrasting nature, depending on whether the process has been addressed within the Franck−Condon approximation or beyond it. The calculation of the Duschinsky matrix and the shift vector between the normal modes of the two states can be done in Cartesian and/or internal coordinates, depending on the flexibility of the investigated molecule. A new computational code has been developed to calculate the rates of intersystem crossing, internal conversion, and fluorescence for selected molecules as functions of temperature.

1. INTRODUCTION Following an electronic excitation, there are multiple pathways for the decay of excited states, which can be broadly classified into radiative and nonradiative processes. The radiative decay pathway from an excited singlet state to the ground electronic state, i.e., fluorescence, usually has a lifetime ranging from 10−9 to 10−6 seconds.1 On the other hand, the molecule can relax from one excited singlet state to another by the nonradiative process of internal conversion (IC), which is very fast, occurring typically within 10−14 to 10−10 seconds.1 A singlet excited state can also decay nonradiatively within 10−10 to 10−8 seconds1 by the process of intersystem crossing to a near-lying triplet state, which can subsequently decay to the ground electronic state by the very slow radiative process of phosphorescence, whose lifetime can range from 10−8 to 103 seconds.1 While fluorescence is probably the most interesting and widely studied process, for molecules exhibiting a multitude of near-lying excited states, the nonradiative processes might become significant decay channels prior to or competing with fluorescence. Additionally, it has been observed that for systems where the optical excitation to the first excited state is symmetry-forbidden, following an excitation to the first bright state, the molecule relaxes to the first excited state by internal conversion and then fluoresces to the ground state, a process © 2015 American Chemical Society

made possible by symmetry relaxation after the nonradiative decay.2 The decay rate of excited states has always been an interesting subject of study from early times.3−5 Fluorescence rates and quantum yields have been measured and calculated in numerous cases for diverse classes of molecules.6−8 In recent years, nonradiative processes have also gathered interest among experimentalists9 as well as among theoreticians.10,11 The rate constant of internal conversion can be evaluated through timeindependent12 and time-dependent formalisms.10 The timedependent formalisms can be traced back to 1955, when Kubo proposed a formalism for calculating the temperature dependence of rate constants of radiative and nonradiative transitions.13 In particular, Lin and co-workers have developed a time-dependent formulation based on the generating function approach.10,11,14 Their work has demonstrated the reliability and feasibility of the time-dependent approach to estimate the rates of such processes. In ref 15, Lin explains a formalism for the calculation of absorption and emission lineshapes from a time-dependent approach based on generating functions. For the calculation of the rate of intersystem crossing, a timedependent formalism has been developed and implemented in Received: October 27, 2015 Published: December 18, 2015 774

DOI: 10.1021/acs.jctc.5b01017 J. Chem. Theory Comput. 2016, 12, 774−786

Article

Journal of Chemical Theory and Computation

is always a possibility that following an excitation to the first available bright state, the molecule can decay nonradiatively to the first excited state (S1) via internal conversion, following which fluorescence to the ground state (S0) can take place, the transition now being allowed due to symmetry-relaxation during the process of internal conversion.2 This led to the speculation of using one of the diamondoid molecules as an exemplary case for calculating the rate constant of internal conversion. The molecule studied here was [123]-tetramantane,35 according to the nomenclature employing the Balaban− Schleyer convention.36 Flexible molecules, often encountered in experimental or theoretical studies, can undergo significant structural changes during the transition from one electronic state to another. Treating such systems in terms of Cartesian coordinates often leads to an abruptly decaying correlation function in the timedependent approach,34 which in turn incorrectly leads to a smooth, broad spectrum without any vibrational fine structure. In order to solve this problem, an efficient internal coordinatebased approach has been developed recently by some of us.37 In that work, the definition of the Duschinsky transformation has been modified to support internal coordinates, thus surpassing the limitations of Cartesian coordinates and obtaining improved and more accurate spectra (see ref 37 for more details on the approach). This tool has been used here to calculate the rate of triplet to singlet nonradiative decay for the flexible molecule biphenyl using the internal coordinate-based approach. As another illustration of the calculation of the triplet to singlet nonradiative rate constant using the time-dependent approach, the bis(terpyridine)-ruthenium(II) complex, [Ru(tpy)2]2+, was studied. The vibrationally resolved T2 → S0 phosphorescence spectrum of this molecule in acetonitrile solvent (S0 being the ground electronic state and T2 the second excited triplet state), computed recently38 at the B3PW91/ LANL2DZ39−41 level of theory, reproduced the line-shape of the experimental phosphorescence spectrum very accurately, thus confirming the applicability of the protocol for taking into account the vibronic effects correctly. The same method and level of theory have been used in the present work to estimate the rate of the nonradiative T2 → S0 transition in this complex. Previous theoretical investigations42 based on the statespecific solvation scheme43−45 within the polarizable continuum model (PCM)46/time-dependent density functional theory (TD-DFT) framework have shown that the solvent reorganization energy provides a good estimation of the polar contributions to the inhomogeneous broadening in solution. This, in turn, allows one to obtain accurate absorption lineshapes of molecules in solution. Here, we have used the value of the solvent reorganization energy calculated at the ground state equilibrium geometry to estimate the value of the electronic coupling between the S0 and T2 states. This paper has been organized in the following way: In Section 2, we present the underlying theory for our implementation and the analytical expressions for the calculation of the rate constants. The computational details are presented in Section 3, while Section 4 contains the validations for the code and results on the temperature dependence of the rate constants of radiative and nonradiative decay processes in certain systems of interest. Finally, Section 5 summarizes and concludes this work.

the quantum-chemical package VIBES16 by Marian and coworkers. In both above-mentioned time-dependent formalisms, the effects of mode mixing between the vibrational normal modes of the ground and excited states (the so-called Duschinsky effect17) have been included (however, not in ref 15) together with temperature effects. In the present work, a combination of these approaches has been implemented, with the aim of providing a general tool for the calculation of the appropriate time-dependent autocorrelation function involved in different types of radiative and nonradiative transitions at various temperatures. The formalism can be broadly applied at two levels: the Franck−Condon (FC) regime18−20 and beyond the Franck−Condon regime. The formalism for the calculation of the rate of intersystem crossing falls under the first category, while that for the rate of internal conversion belongs to the second. For fluorescence, the dependence of the transition electric dipole moment on the nuclear coordinates is neglected within the FC approximation, while it is considered in the Herzberg−Teller (HT) regime.21 Heller and co-workers have also developed time-dependent formalisms for the calculation of vibronic spectra within the FC and HT regimes22 using a semi-classical wave packet propagation approach, which has often found its way into quantum chemical packages such as ORCA.23 Another type of nonradiative process involves the transfer of charge or excitation energy from one molecule (donor) to another (acceptor). These processes, namely, electron transfer,24 proton transfer,25 and excitation energy transfer are treated in a different theoretical framework, using classical Förster’s theory.26 The electronic coupling depends on the interaction of the charges or dipoles and their separation.24 The vibronic contribution is usually determined from the spectral overlap between the donor-emission and acceptor-absorption spectra. These processes are not addressed in this article and are being studied now for a future work. The time-dependent approach for the calculation of rate constants of internal conversion has been investigated for an isomer of tetramantane, a molecule which belongs to the class of novel, saturated, cyclic hydrocarbons called diamondoids.27,28 The hardness, inertness, and high thermal conductivity of diamondoids have made them subjects of widespread industrial use.29,30 Additionally, diamondoids have very attractive optical properties, e.g., they absorb light around and above 6 eV and show excellent fluorescence in the solid and gas phases.27,31 This makes them useful as monochromatic electron-emitters. Due to the diversity in their shape, size, and composition, a large number of artificial diamondoids with tunable optical properties have been synthesized32,33 over the past decade. For example, with an adequate number, position, and type of substituent(s), the fluorescence can be shifted to the visible region, therefore, making them useful as biological marker-dyes.27 For these reasons, in the recent past, diamondoids have been subjects of large-scale experimental and theoretical investigations. A systematic study of the vibrationally resolved absorption, emission, and resonance Raman spectra of a large number of naturally occurring and artificial diamondoids has been performed recently.34,35 Comparisons to experiments have been made wherever possible, and good agreements have been observed. Due to their high symmetry, diamondoids often have closely spaced excited states.2,35 For the same reason, optical excitation from the ground to the first excited electronic state is forbidden in certain diamondoids.2,35 As a result, there 775

DOI: 10.1021/acs.jctc.5b01017 J. Chem. Theory Comput. 2016, 12, 774−786

Article

Journal of Chemical Theory and Computation

2. THEORY 2.1. Calculation of Rate Constants and Their Temperature Dependence. The general expressions for the radiative and nonradiative processes can be broadly classified into two categories: zeroth order (at the FC level of approximation) and first order (where the first-order derivatives of certain relevant physical parameters have to be considered). In the adiabatic regime, the rate constant for intersystem crossing can be expressed by Fermi’s golden rule as47 kISC =

is calculated self-consistently with the excited state solute electron density. The rate constants for radiative processes like fluorescence and phosphorescence can also be expressed in a way similar to eq 1, but an additional dependence on the emitted light frequency needs to be introduced in the rate expression. For example, the rate of fluorescence can be expressed within the FC approximation as48 kF =

2π 2 V f (ΔE , T ) ℏ

∑ Z −1e−βE ⟨Φ̅ v ̅ |Φv ⟩2 δ(ΔE + Ev v̅

v ,v

− Ev̅)

I(ω) = IFC(ω) =

kISC

4ω3 2 V f (ΔE , T ) 3ℏc 3

(6)

with ω being the frequency of the emitted photon. In this case, the coupling element V is given by the magnitude of the transition electric dipole moment μ (vectors and matrices have been expressed in bold throughout the paper). Going beyond the Franck−Condon approximation, the electronic coupling V can be expressed as a Taylor series around its value at the reference geometry, usually truncated after the first order, generally referred to as the HT term.

(2)

where |Φ̅ v⟩̅ and |Φv ⟩ denote the initial and final vibrational wave functions, respectively. Here and throughout the paper, final state quantities are represented by double overbars and initial state ones with a single overbar. Also, β = (1/kBT) (kB being the Boltzmann constant), and Z is the thermal vibrational partition function of the initial state. E v and Ev̅ are the vibrational energies in the respective electronic states. The harmonic approximation has been used throughout. Within the Condon approximation, the coupling element V is usually evaluated at a particular reference geometry and assumed to be constant. Expressing the Dirac delta function as the Fourier transform of a time (t)-dependent correlation function, the rate constant can be expressed as16 R2 = 2 Z −1 ℏ

(5)

where the fluorescence intensity I(ω) is given by

where the general coupling term V is in this case the spin−orbit coupling between the singlet and triplet electronic states, and ΔE represents the adiabatic energy difference between the two states. The temperature (T) dependence of the rate constant originates from the Franck−Condon weighted density of states, f(ΔE, T) and is therefore dependent only on the vibrational coordinates. f (ΔE , T ) =

∫ I(ω)dω

(1)

⎛ ∂μ ⎞ ⎟⎟ (Q k − Q 0) + ... k = 1 ⎝ ∂Q k ⎠ N

μ = μ0 +

∑ ⎜⎜

(7)

Q0

where the transition dipole moment has been expanded as a Taylor series with respect to the final state normal coordinates around its equilibrium value (μ0) at a reference geometry Q 0 . For convenience, from now on, the first-order derivative ⎛ ∂μ ⎞ ⎜ ⎟ with respect to the kth normal mode will be written as ⎝ ∂Q k ⎠Q 0



∫−∞ dt e

iΔEt

χISC (t )

μk. Within this level of approximation, I(ω) cannot be expressed by eq 6 but as

(3)

where R is the spin−orbit coupling constant, and χISC(t) is a time-dependent correlation function. Additionally, for nonradiative transitions in solutions, the spin−orbit coupling has been approximated by the inhomogeneous solvent-broadening, which in turn is estimated by the solvent reorganization energy in going from one electronic state to the other. For a process involving the ground and one excited state of a molecule, the solvent reorganization energy is defined as42,43 λssp = ΔEneq − ΔEeq

I(ω) = IFC(ω) + IFC/HT(ω) + IHT(ω)

(8)

where IFC/HT(ω) =

4ω3 3ℏc 3

∑ Z −1e−βE δ(ΔE + Ev v̅

− Ev̅)

v ,v

∑ μ0 ·μk ⟨Φ̅ v ̅ |Φv ⟩⟨Φ̅ v ̅ |Q k|Φv ⟩ k

(4)

(9)

and

where ΔEneq and ΔEeq represent the vertical transition energies at the equilibrium geometry of the ground state within the nonequilibrium and equilibrium solvation models, respectively, computed using a state-specific solvation scheme.42,43 In the latter case, all the degrees of freedom of the solvent are allowed to come into equilibrium with the electron density in the excited state, whereas in the former case only the fast degrees of freedom of the solvent are allowed to do so.43 It should be noted that in the state-specific scheme of solvation, in contrast to the usual linear-response approach, the solvent reaction field

IHT(ω) =

4ω3 3ℏc 3

∑ Z −1e−βE δ(ΔE + Ev v̅

− Ev̅)

v ,v

∑ μk ·μm ⟨Φ̅ v ̅ |Q k|Φv ⟩⟨Φ̅ v ̅ |Q m|Φv ⟩ k ,m

(10)

In the time domain, eq 8 can be expressed as the sum of the Fourier transform of three time-dependent correlation functions 776

DOI: 10.1021/acs.jctc.5b01017 J. Chem. Theory Comput. 2016, 12, 774−786

Article

Journal of Chemical Theory and Computation I(ω) = +

+

2ω3 2 −1 μ Z 3π ℏc 3 0

The electronic coupling, i.e., the nonadiabatic coupling vector can be expressed as48,49



∫−∞ dt eiΔEtχ(t )

2ω3 3π ℏc 3

∑ μ0 ·μk Z −1 ∫

2ω3 3π ℏc 3

∑ μk ·μm Z −1 ∫



−∞

k



−∞

k ,m

⟨Ψelec|(∂V /∂Q k)|Ψ̅elec⟩ ⟨Ψelec|Pk̂ |Ψ̅elec⟩ = (E ̅ − E )

dt eiΔEt χFC/HT; k (t )

Here, (∂V /∂Q k) is the gradient of the Coulomb potential of the initial state with respect to the kth normal coordinate of the final state. E̅ and E are the total energies of the initial and the final electronic states, respectively. Applying the Condon approximation, the gradients and the energies are evaluated at the lower state equilibrium geometry. 2.2. Analytical Expressions for Correlation Functions: Franck−Condon Level. As discussed in more details in refs 50 and 51, the time-dependent correlation function of eq 19 is computed by solving Gaussian type integrals, obtained using Mehler’s formula for the density matrix of an oscillator,50 tracing over the final state electronic coordinates. The resulting autocorrelation function can be expressed as48

dt eiΔEt χHT; k , m (t ) (11)

In eq 12, χ(t) denotes the correlation function at the Franck− Condon level (similar to χISC(t) in eq 3). χFC/HT;k(t) and χHT;k,m(t) refer to the mixed Franck−Condon/Herzberg−Teller and pure Herzberg−Teller contributions, respectively. The rate expression of IC also goes beyond the Franck− Condon regime. Here, both the electronic and the vibrational wave functions of the two electronic states are coupled by the nuclear momentum operator. Considering the process of internal conversion between an initial state |Ψ̅ elecΦ̅ v⟩̅ and a final state |ΨelecΦv ⟩, the interaction between the two states is governed by the nonadiabatic interaction facilitated by the nuclear-momentum operator48 2

 = −ℏ

∑ ⟨ΨelecΦv | k

∂Ψ̅elec ∂Φ̅ v ̅ ∂Q k ∂Q k

χ (t ) =

∂Ψ̅elec

k

∂Q k

⟩⟨Φv |



∂Φ̅ v ̅ ∂Q k

(21)

(12)

where K is the shift vector between the normal modes of the two states, and J is the Duschinsky matrix, defined by the Duschinsky transformation, which describes the mixing of the normal modes in the initial and final states17,52



Q̅ = JQ + K

(13)

∑ kIC;k ,m

A = a + JT aJ̅

(14)

k ,m

B = b + JT bJ̅

where the rate constant with contributions from the kth and mth vibrational normal modes are given by48 kIC; k , m

2π = R k , m ∑ Z −1e−βE v ̅Pk , mδ(ΔE + E v − E v ̅ ) ℏ v ,v

(15)

(16)

is the nonadiabatic coupling vector and Pk , m

= ⟨Φv |Pk̂ |Φ̅ v ̅ ⟩⟨Φv |Pm̂ |Φ̅ v ̅ ⟩

akk ̅ =

ω̅ k sin(ℏτ ̅ω̅ k )

b kk ̅ =

ω̅ k tan(ℏτ ̅ω̅ k)

akk =

ωk sin(ℏτ ωk)

b kk =

ωk tan(ℏτ ωk)

∂ ∂Q k

(17)

τ̅ = − τ =

1 R k , mZ −1 ℏ2

t ℏ

(25)

ω̅ k and ωk are the vibrational frequencies of the k normal mode in the two states. The correlation functions have been damped with Lorentzians of suitable widths to include the effects of homogeneous broadening. 2.3. Analytical Expressions for Correlation Functions: Beyond the Franck−Condon Level. To calculate the correlation function for the rate of internal conversion, χk,m(t), the following matrices and vector are additionally introduced48

(18)



∫−∞ dt eiΔEtχk ,m (t )

i t − ℏ kBT

th

the nuclear momentum operator for the kth normal mode. Within the time domain, the rate constant can be expressed as48 kIC; k , m =

(24)

where

is the coupling between the vibrational wave functions of the two states, with Pk̂ = −iℏ

(23)

The elements of the diagonal matrices a,̅ a , b̅ and b are defined in the following way:

Here, R k , m = ⟨Ψelec|Pk̂ |Ψ̅elec⟩⟨Ψelec|Pm̂ |Ψ̅elec⟩

(22)

The matrices A and B, which have been introduced to obtain compact expressions for the correlation functions, are defined as

The overall rate is now obtained by adding the contributions of the different pairs of normal modes

kIC =

det( aa̅ ) det(B) det(B − AB−1A) × exp[(i /ℏ)(KT(b̅ − a̅ )J(B − A)−1(b − a)JT K)]

Applying the Condon approximation, the expression simplifies to  = −ℏ2 ∑ ⟨Ψelec|

(20)

(19)

where χk,m(t) is the corresponding correlation function. 777

DOI: 10.1021/acs.jctc.5b01017 J. Chem. Theory Comput. 2016, 12, 774−786

Article

Journal of Chemical Theory and Computation χFCHT; k , m (t ) = χ (t )[iℏTr(GkmS−1) + (S−1F)T

E = b̅ − a̅ F = ( KTEJ KTEJ ) S=

Gkm(S−1F) − HTk S−1F]

⎛ B −A ⎞ ⎜ ⎟ ⎝− A B ⎠

This is similar to the IC correlation function given by eq 27, but the forms of the matrix Gkm and vector Hk will differ;48 the first two terms in the product give the pure Herzberg−Teller contribution, while the last term gives the contribution of the mixed Franck−Condon and Herzberg−Teller term. In this case,

(26)

The correlation function is expressed as a product of the FC correlation function and terms coming from the coupling of the vibrational wave functions by the nuclear momentum vector χk , m (t ) = χ (t )[iℏTr(GkmS−1) + (S−1F)T Gkm(S−1F) −

HTkmS−1F]

(27)

Gkm

The matrix Gkm and vector Hkm used in eq 27 are defined in terms of other matrices and vectors, in the following way: Gkm

(32)

(33)

kBT

becomes undefined. The diagonal matrix b̅ is substituted by the diagonal matrix γ ̅ whose elements are given by52

(29)

In eqs 28 and (29),

γkk̅ =

... 0

G22,km =

... ⎞ ... ⎟ ⎟ ... ⎟⎟ ... ⎠

For calculations at 0 K, the previous formalism needs to be modified based on first-order Taylor expansion of trigonometric functions52 because the term 1 appearing in eq 25

⎛ H1,km ⎞ ⎟⎟ Hkm = ⎜⎜ ⎝ H 2,km ⎠

G21,km =

... 01, N + m ... ... ... 1k , N + m ... ...

Hk = ( 01 ... 1k ... 02N )T

(28)

and

G12,km =

⎛ 01,1 01,2 ⎜ ... ... =⎜ ⎜⎜ 0k ,1 0k ,2 ⎝ ... ...

and

⎛ G11,km G12,km ⎞ ⎟⎟ = ⎜⎜ ⎝G21,km G22,km ⎠

G11,km =

(31)

⎛ ⎞ ⎜ ⎟ ⎜ ⎟ T ⎜−b kk (J aJ̅ )m ⎟ ⎜ ⎟ 0 ⎜ ⎟ ⎝ ⎠ ... ... ⎛ ⎞ ⎜ ⎟ 0 ⎜ ⎟ T ⎜ b kk (J bJ̅ )m ⎟ ⎜ ⎟ 0 ⎜ ⎟ ⎝ ⎠ ... ... ⎛ ⎞ ⎜ ⎟ 0 ⎜ ⎟ T ⎜ akk(J aJ̅ )m ⎟ ⎜ ⎟ 0 ⎜ ⎟ ⎝ ⎠ ... ... ⎛ ⎞ ⎜ ⎟ 0 ⎜ ⎟ T ⎜−akk(J bJ̅ )m ⎟ ⎜ ⎟ 0 ⎜ ⎟ ⎝ ⎠ ...

⎛ ∂s ⎞ Bij = ⎜⎜ i ⎟⎟ ⎝ ∂xj ⎠eq

(35)

Here, si refers to the ith internal coordinate, and xj is the jth Cartesian coordinate. The matrix (and other) notations in this part of the text are independent of those used previously. B is involved in the transformation of the energy Hessian from Cartesian (Hx) to internal (Hs) coordinates. For adiabatic models, the transformation is the following:

Hs = BT HxB

H1,km = ( ... b kk (K EJ)m 0 ... ) T

(34)

In addition, a ̅ is set to zero everywhere, except in the determinant prefactor of eq 21, where it is set to γ.̅ 2.4. Extension to Internal Coordinates. For the computations in internal coordinates, the approach described in ref 37 has been used, which is based on the extension of the Duschinsky transformation proposed by Reimers.53 This approach relies on the Wilson theory, where the internal coordinates s are expanded as a Taylor series up to the rst-order in terms of the Cartesian coordinates. The matrix of the firstorder derivatives is usually referred to as Wilson matrix B, where

T

T

iω̅ k ℏ

T

H 2,km = ( ... −akk(K EJ)m 0 ... )

(36)

Once Hs is built, the normal modes in internal coordinates (Ls) can be computed as eigenvectors of GHs, where the Wilson matrix G is defined as follows:

(30)

Here, terms like (JT b̅ J)m refer to the mth row of the matrix (JT b̅ J) or other corresponding matrices. (JT b̅ J)m is then multiplied by b kk and inserted into the kth row of G12,km in eq 30. The other rows of G12,km are null. The total correlation function at the HT level, (χFCHT;k,m(t)) is given as48

G = BM−1BT

(37)

with M being the diagonal matrix of the atomic masses. Unlike in Cartesian coordinates, when internal coordinates are used, Ls is not an orthogonal matrix. Nevertheless, for adiabatic models, the Duschinsky transformation can be generalized as follows: 778

DOI: 10.1021/acs.jctc.5b01017 J. Chem. Theory Comput. 2016, 12, 774−786

Article

Journal of Chemical Theory and Computation

Js = L̅ −s 1Ls K s = L̅ −s 1(s − s ̅ )

total propagation time of 100 ps was used with a time step of 0.1 fs for the calculation of autocorrelation functions. To ensure their decay to zero, they have been damped by applying Lorentzians with full width at half-maximum (fwhm) between 50 and 70 cm−1 depending on the system investigated. All calculations were done with the adiabatic Hessian (AH) model at the FC and FCHT levels. Following the notation adopted by some of us,52 the methods employed to do the calculations will be labeled AH|FC and AH|FCHT in the following. In all cases, the rate constants have been calculated between T = 0 and T = 300 K with a step of 50 K in order to study their variation with temperature.

(38)

The correlation functions can then be calculated by the formalism described before using Js and Ks instead of J and K, respectively.

3. COMPUTATIONAL DETAILS The theory presented above has been implemented in a standalone program (written in Fortran) which can interact directly with the GAUSSIAN package54 to extract the vibrational frequencies, Duschinsky matrix, and shift vector. Extension to support other electronic structure calculation programs is straightforward. This tool can also be used to combine data from multiple sources, especially to compensate missing data from a particular package, such as nonadiabatic couplings. The FFTW (Fastest Fourier Transform in the West) package55 has been used to carry out Fourier transforms. Geometry optimizations and force constant calculations have been carried out at the DFT level (and TD-DFT for excited states) using GAUSSIAN in all cases. Most calculations were performed using the B3LYP functional56 in conjunction with the 6-31G(d),57 TZVP,58 and SNSD59 basis sets. The last is a double-ζ quality basis set, supplemented by diffuse functions. For bis(terpyridine)-ruthenium(II) ([Ru(tpy) 2 ] 2+ ), the B3PW91 hybrid functional39,40 was used with the LANL2DZ basis set, and the corresponding pseudopotential for describing the inner electrons of Ru and augmented with polarization functions on C, N, O, and Ru. Solvent effects were simulated with the polarizable continuum model (PCM). Its integral equation formalism (IEF-PCM60) with the default parameters of GAUSSIAN was used to account for the bulk solvent effects of acetonitrile (dielectric constant of 35.69) and n-hexane (dielectric constant of 1.88). Geometry optimizations were performed with tight convergence criteria (maximum displacements and forces smaller than 6 × 10−5 Å and 1.5 × 10−5 Hartree/Bohr respectively) with an ultrafine integration grid (99 radial shells and 590 angular points per shell). Force constants were computed at the equilibrium geometry of each state. For excited electronic states, they are computed by numerical differentiation of analytical gradients using steps of 0.01 Å along each Cartesian coordinate. First derivatives of the electric transition dipole moments needed for HT calculations are readily available during the process. Spin−orbit and nonadiabatic couplings were calculated using the CASSCF (complete active space self-consistent field61) level of theory. For spin−orbit couplings, the implementation in GAUSSIAN was used, while nonadiabatic couplings were obtained from MOLPRO.62 It is worth reminding that highaccuracy computation of the electronic coupling is not the focus of this work, so for certain complex systems, the electronic couplings have been estimated using approximate approaches. For the vibronic calculations, as mentioned before, the Duschinsky matrix and shift vector were directly extracted from GAUSSIAN. For the internal coordinate description, delocalized internal coordinates (DICs)63 have been used as a nonredundant set of internal coordinates. As already discussed in ref 37, DICs give more reliable results with respect to the widely used Z-matrix coordinates in a fully automated procedure, so they have been chosen here as well. Then, a

4. RESULTS 4.1. Time-Dependent Formalism at Franck−Condon Level of Approximation: Rate of Intersystem Crossing in Phenalenone. In order to validate our implementation, the rate of intersystem crossing between the S1 and T1 states of phenalenone was computed and compared to the results obtained by Tatchen and co-workers,16 who used a timedependent correlation function approach based on the generating functions, without taking temperature effects into account. To remain as close as possible to their protocol, the same electronic structure calculations were performed. Figure 1 shows the ground state equilibrium geometry of phenalenone.

Figure 1. Ground state equilibrium geometry of phenalenone at the B3LYP/TZVP level of theory. Brown and white balls represent C and H atoms, respectively. Red ball represents the O atom.

The value of the spin−orbit coupling (SOC) was taken from ref 16, where the authors used the SPOCK-CI kit64 interfaced with TURBOMOLE65 for the calculation. For a better comparison of our results, the value of the adiabatic energy difference between the S1 and the T1 states calculated at the DFT/MRCI (multireference configuration interaction) level of theory66 was also taken from the same reference. The correlation functions were damped using Lorentzians with fwhm of 10 cm−1 as done in ref 16. Comparison between the rate constants obtained from our calculations (also at 0 K) and those of ref 16 have been summarized in Table 1. The comparison shows a very good agreement, validating our implementation for the calculation of the rates of ISC at the FC level of approximation. 4.2. Time-Dependent Formalism beyond Franck− Condon Level: Internal Conversion and Herzberg−Teller Effect on the Fluorescence of Azulene. The rate of S1 → S0 internal conversion in azulene was calculated at the B3LYP/6779

DOI: 10.1021/acs.jctc.5b01017 J. Chem. Theory Comput. 2016, 12, 774−786

Article

Journal of Chemical Theory and Computation Table 1. Comparison between Our Calculated Values of Rate Constants of ISC from the S1 State to the Two Triplet Spin Microstates T1x and T1y and Those of Ref 16

a

states involved

(SOC)2 (cm−2)

kISC (1/s)a

kISC (1/s)b

S1 → T1x S1 → T1y

1823.0 234.0

1.8 × 1010 5.7 × 109

2.0 × 1010 5.8 × 109

Ref 16. bThis work. Figure 3. Superimposed equilibrium geometries of S0 (red) and T1 (blue) states of biphenyl. In the S0 state, the two phenyl rings are tilted with respect to each other by 40°.

31G(d) level of theory. This was chosen to match as closely as possible the calculations of Lin and co-workers48 who used the B3LYP/SV(P) level of theory. The SV(P) basis set,67 while directly available in TURBOMOLE but not in GAUSSIAN, is qualitatively similar to the 6-31G(d) basis set.65 Figure 2 shows the ground state equilibrium geometry of azulene.

must be used in order to describe correctly vibronic effects.37,69 In order to highlight the difference of the Duschinsky transformation using the two different coordinate systems for the T1 ← S0 transition, in Figure 4, a graphical representation of

Figure 2. Ground state equilibrium geometry of azulene at the B3LYP/6-31G(d) level of theory. Brown and white balls represent C and H atoms, respectively.

The total rate constant of IC in azulene at 300 K, as the sum of the contribution of the individual modes calculated according to eqs 14 and 19, is equal to 1.41 × 1011 s−1. This is a slight improvement (considering the experimental value68 of 5.0 × 1011 s−1) over the value of 0.23 × 1011 s−1, obtained by Lin and co-workers. The S1 → S0 radiative decay rate of azulene was calculated at both Franck−Condon and Herzberg−Teller levels of approximation. At the Franck−Condon level, the calculated value of the rate constant was 1.07 × 106 s−1 at 300 K, slightly higher than the experimental value48 of 1.0 × 106 s−1. The rate constant, evaluated at an interval of 50 K in the range of 0 to 300 K, was found to be nearly independent of temperature, showing a slight decrease from 1.12 × 106 s−1 at 0 K to 1.07 × 106 s−1 at 300 K. The transition dipole moment derivatives were negligible for most of the modes, the largest term being approximately 0.008 atomic units. In spite of this, inclusion of Herzberg−Teller contributions increased the rate constant obtained at the Franck−Condon level; at 300 K, kF for S1 → S0 fluorescence was found to be 1.7 × 106 s−1. 4.3. Rate of T1 → S0 Intersystem Crossing in Biphenyl: An Internal Coordinate Approach. As already mentioned above, biphenyl has been used to test the reliability of our implementation in internal coordinates. In fact, the angle between the two phenyl rings is 40° at the equilibrium geometry of the S0 state (Figure 3) and near null for the T1 state, and therefore, a large-amplitude distortion is involved in the electronic transition. As a consequence, internal coordinates

Figure 4. Comparison of the shift vector (left) and Duschinsky matrix (right) between the normal modes of the T1 and S0 states of biphenyl using Cartesian (top) and internal (bottom) coordinates. Blue bars correspond to positive displacements and red bars to negative ones.

Kx, Jx, Ks, and Js, calculated from the B3LYP/SNSD-optimized S0 and T1 states, is reported. Even if the overall structure of the shift vector is similar, in the Cartesian representation, the shift of the high-energy modes is significantly larger. The Duschinsky matrix is similar in both cases, so the more pronounced effect of switching from the Cartesian to internal coordinates is on the shift vector. This effect is caused by the poor description of torsion as linear combination of Cartesian coordinates, which leads to an incorrect definition of the Duschinsky transformation. In fact, as shown in Figure 5, the phosphorescence spectrum computed at the FC level (using GAUSSIAN) for this transition is reproduced correctly only using internal coordinates, whereas significant discrepancies are present when Cartesian coordinates are employed. 780

DOI: 10.1021/acs.jctc.5b01017 J. Chem. Theory Comput. 2016, 12, 774−786

Article

Journal of Chemical Theory and Computation

Figure 6. Vibrational normal modes which have been neglected for the spectral simulations according to the reduced dimensionality scheme followed here for biphenyl. The modes have been numbered in order of increasing vibrational frequencies.

Figure 5. Comparison of the T1 → S0 phosphorescence spectrum of biphenyl computed at the TI AH|FC level (stick) in Cartesian (top) and DICs (bottom) with the experimental one (dashed curve).

Let us remark that the simulation of the phosphorescence spectrum has been carried out using the reduced dimensionality scheme presented in ref 70 to neglect the contributions due to the torsional normal mode. In the Cartesian case, also an inplane ring deformation mode (a graphical representation of this mode is given in Figure 6), which is strongly mixed to the torsion through the Duschinsky matrix, has been neglected. It is noteworthy that the main difference between the two spectra is due to the presence of intense vibronic transition also in the low-energy range of the spectrum (19000−21000 cm−1). The presence of those additional bands can be understood better from the plot of the Sharp and Rosenstock C matrix,71 which is reported in Figure 7. C is involved in the simulation of vibronic spectra within the Sharp and Rosenstock TI formulation and, more specifically, controls the contribution from mixing of the final state modes in the FC integral. As shown in Figure 7, in the Cartesian representation, the extent of coupling is much higher, which results in some combination bands contributing to the spectrum. However, when calculated using the internal coordinates, the extent of coupling is much lower, which improves the spectral quality and the intensities. The spin−orbit coupling calculated at the CASSCF(6,6)/6311+G* level of theory using GAUSSIAN was 0.2 cm−1. The convergence of the active space, which consisted of three highest occupied π and three lowest unoccupied π* orbitals, was checked. Henceforth, the rate constant for the T1 → S0

Figure 7. Sharp and Rosenstock C matrix from Cartesian (top) and internal (bottom) coordinates based on AH|FC calculations. The same normalization factor of 0.33, corresponding to the highest magnitude of the two coupling matrices, has been used to facilitate visual comparisons. The representation is obtained by calculating the square of the elements of C and normalizing them with respect to the maximum value of C2ij. A shade of gray is associated with each element (i, j) in the figure based on the equivalence (0, white; 1, black).

781

DOI: 10.1021/acs.jctc.5b01017 J. Chem. Theory Comput. 2016, 12, 774−786

Article

Journal of Chemical Theory and Computation nonradiative decay at 300 K was found to be 4.64 × 104 s−1. The calculated value can be compared to the experimentally determined T1 decay rate of 2 × 104 s−1 at 300 K in methylcyclohexane−isopentane, as reported in ref 72 or (0.74− 0.80) × 104 s−1 in n-hexane at 300 K, as in ref 73. Additionally, we predict that the rate constant decreases with temperature, as shown in Figure 8.

Table 2. Comparison of Vertical Excitation Energies at S0 Minimum Geometry Using State-Specific PCM-Solvation Model for S0 → T2 Transition transition

ΔEeq (eV)

ΔEneq (eV)

S0 →T2

1.8959

2.1688

As mentioned above, explicit calculation of the spin−orbit coupling was not performed in this case. Instead, the coupling between the S0 and the T2 states was approximated as the solvent reorganization energy for the S0 → T2 transition at the S0 equilibrium geometry (computation of the T2 → S0 emission solvent reorganization energy has to be done with a ground state DFT calculation instead of the state specific, TD-DFT scheme, which results in an inaccurate estimation of the solvent reorganization energy42 and is, therefore, not used here). The rate of the T2 → S0 transition was then calculated by substituting R with the solvent reorganization energy (0.29 eV) in eq 3. Table 3 presents the variation of the rate constant Table 3. Comparison of Variation of kISC with Temperature in Gas Phase and in Solution

Figure 8. Temperature dependence of the rate constant for the T1 → S0 intersystem crossing of biphenyl in gas phase. The rate constant decreases slowly but steadily with temperature.

For a closer comparison to the experimental rate constant given in ref 73, the ground state was reoptimized in n-hexane solution. The spin−orbit coupling calculated at the new S0 equilibrium geometry changed to 0.1 cm−1. The same Duschinsky matrix and shift vector as used for the gas phase calculations (solvent effects on the normal modes and vibrational frequencies are neglected, considering the very low polarity of the solvent) were used to calculate the correlation function. kISC was calculated to be 1.16 × 104 s−1 at 300 K. This is in better agreement with the experimentally determined value of the rate constant. 4.4. Rate of Nonradiative T2 → S0 Transition in [Ru(tpy)2]2+. Figure 9 shows the ground state equilibrium geometry of [Ru(tpy)2]2+ in acetonitrile solution. Table 2 shows the vertical excitation energies calculated using the equilibrium and nonequilibrium models of statespecific solvation.

phase

temperature (K)

kISC (× 1014 s−1)

gas

0 50 100 150 200 250 300 0 50 100 150 200 250 300

0.91 0.81 0.65 0.52 0.42 0.34 0.28 1.50 1.50 1.45 1.33 1.18 1.03 0.89

solution (acetonitrile)

of the T2 → S0 transition as a function of temperature. Since the electronic factor remains constant, the variation of the rate constant reflects solely the vibronic effects. Figure 10 shows that the rate constant decreases with an increase in temper-

Figure 10. Temperature dependence of the rate constant for the T2 → S0 nonradiative transition of [Ru(tpy)2]2+ in gas phase (black) and acetonitrile solution (red) calculated at the B3PW91/LANL2DZ level of theory.

Figure 9. Ground state equilibrium geometry of [Ru(tpy)2]2+ in acetonitrile solution. Brown and white balls represent C and H atoms, respectively, while blue balls represent N. The large yellow ball represents Ru. 782

DOI: 10.1021/acs.jctc.5b01017 J. Chem. Theory Comput. 2016, 12, 774−786

Article

Journal of Chemical Theory and Computation

case of [123]-tetramantane is equal to 144. At 300 K, the value of the rate constant was found to be 1.44 × 1010 s−1. Figure 12 shows the variation of kIC with temperature. As shown in the plot, the rate constant increases gradually with

ature, as already observed for biphenyl in Figure 8. Results based on electronic structure calculations in the gas phase also showed a similar behavior (Figure 10), the values of the rate constants being higher than in solution. A similar, nonArrhenius behavior of the rate constant has also been observed in the calculated rates of electron transfer between bacteriopheophytin and ubiquinone within the FC regime by Borrelli and Peluso.47 As a general observation, kISC is predicted to be very high for the T2 → S0 transition, but it should be kept in mind that an explicit calculation of the spin−orbit coupling was not performed here. 4.5. Rate of Internal Conversion in [123]-Tetramantane and Its Variation with Temperature. The B3LYP/ TZVP-optimized ground state geometry of [123]-tetramantane is shown in Figure 11.

Figure 12. Temperature dependence of the rate constant for the S4 → S1 internal conversion of [123]-tetramantane in the gas phase.

temperature in the given range, the relative increase from 0 to 300 K being about 55%. Once again, the increase can be ascribed to vibronic effects since the nonadiabatic coupling vector does not vary with temperature. The variation of kIC with temperature has a contrasting nature compared to that of kISC. This might stem from the fact that for the former the rate expression goes beyond the Franck−Condon level of approximation. As shown in eq 27, it consists of the product of the FC correlation function and a difference of the contributions from the two HT-like terms. This might be the cause for the change in the behavior of the overall function with temperature. Comparing to the experimentally determined2 excited state decay rate (∼2.4 × 109 s−1) of [123]-tetramantane, it can be predicted that internal conversion might be an important decay channel for the excited state relaxation of diamondoids, especially in the presence of a large number of close excited electronic states. We observe a similar but smaller temperature-dependent increase in the rate constant beyond the FC regime in our calculated values of the rate constant of fluorescence of azulene when HT contributions are taken into account, i.e., the rate constant for the S1 → S0 fluorescence increased from 1.63 × 106 s−1 at 0 K to 1.70 × 106 s−1 at 300 K.

Figure 11. Optimized ground state geometry of [123]-tetramantane at the B3LYP/TZVP level of theory. Brown and white balls represent C and H atoms, respectively.

The analysis of the electronic excitation energies predicted the fourth excited singlet state, S4, to be the optical transition with the highest oscillator strength, about seven times more intense than the S0 → S1 transition. The transition to S1 is dominated by a HOMO−LUMO excitation, while that to S4 involves a dominant HOMO → (LUMO+1) transition. Additionally, the S1 and S4 states are relatively close in energy (vertical gap of 0.45 eV at the S0 equilibrium geometry) compared to the S1−S0 pair. This hints at the possibility that following an electronic transition to the S4 state the molecule decays to the S1 state by internal conversion, followed by fluorescence emission to the S0 state. Therefore, we studied the rate constant of internal conversion from the S4 to the S1 state in this work. The nonadiabatic coupling vector was calculated using CASSCF(6,5)/6-311G* at the equilibrium geometry of the lower state. The active space included the three highest occupied and the two lowest unoccupied molecular orbitals. kIC has been calculated using eqs 14 and 19, where only the sum over the diagonal elements kIC;k,k has been taken to estimate the total rate constant. An earlier theoretical study on azulene48 has shown that about 96% of the total rate constant comes from the diagonal (kIC;k,k) terms. Preliminary tests over the strongly coupled modes of this system confirmed this trend, so this simplification was used here as well. This is done in order to save an N-fold increase in computation time, N being the number of vibrational normal modes of the system, which in

5. CONCLUSION A state of the art time-dependent formalism has been implemented, validated, and applied to calculate the rates of radiative and nonradiative processes in medium-to-large-size molecules. The selected case studies range from novel molecules with increasing, large-scale industrial applications like diamondoids to transition metal complexes used in solar cells and for biological applications. The variation of the rate constants with temperature has been studied for each example. Our tool has been designed to support calculations in internal coordinates (in addition to those in Cartesian coordinates), which makes it suitable to treat flexible systems like biomolecules, where nonradiative transitions often play a significant role in ruling photochemical and photophysical processes.74 783

DOI: 10.1021/acs.jctc.5b01017 J. Chem. Theory Comput. 2016, 12, 774−786

Article

Journal of Chemical Theory and Computation

singlet electronic states of ethylene. J. Chem. Phys. 1998, 108, 2044− 2055. (11) Lin, S. H. Rate of Interconversion of Electronic and Vibrational Energy. J. Chem. Phys. 1966, 44, 3759−3767. (12) Rashev, S. Calculation of internal conversion rate constants of single vibronic levels in S1 benzene. J. Chem. Phys. 1994, 101, 6632− 6639. (13) Kubo, R.; Toyozawa, Y. Application of the Method of Generating Function to Radiative and Non-Radiative Transitions of a Trapped Electron in a Crystal. Prog. Theor. Phys. 1955, 13, 160−182. (14) Lin, S. H. Energy Gap Law and Franck-Condon Factor in Radiationless Transitions. J. Chem. Phys. 1970, 53, 3766−3767. (15) Lin, S. H. Spectral Bandshape of Absorption and Emission of Molecules in Dense Media. Theoret. chim. Acta (Berl.) 1968, 10, 301− 310. (16) Etinski, M.; Tatchen, J.; Marian, C. M. Time-dependent approaches for the calculation of intersystem crossing rates. J. Chem. Phys. 2011, 134, 154105. (17) Duschinsky, F. On the Interpretation of Electronic Spectra of Polyatomic Molecules. Acta Physicochim. URSS 1937, 7, 551−566. (18) Franck, J. Elementary processes of photochemical reactions. Trans. Faraday Soc. 1926, 21, 536−542. (19) Condon, E. U. A theory of intensity distribution in band systems. Phys. Rev. 1926, 28, 1182−1201. (20) Condon, E. U. Nuclear motion associated with electronic transitions in diatomic molecules. Phys. Rev. 1928, 32, 858−872. (21) Herzberg, G.; Teller, E. Schwingungsstruktur der Elektronenübergänge bei mehratomigen Molekülen. Z. Phys. Chem. (Leipz.) 1933, 21, 410−446. (22) Reimers, J. R.; Wilson, K. R.; Heller, E. J. Complex time dependent wave packet technique for thermal equilibrium systems: Electronic spectra. J. Chem. Phys. 1983, 79, 4749−4757. (23) Neese, F. ORCA, an ab initio, density functional and semiempirical program package; University of Bonn: Bonn, Germany, 2007. (24) Jortner, J. Temperature dependent activation energy for electron transfer between biological molecules. J. Chem. Phys. 1976, 64, 4860− 4867. (25) Warshel, A. Dynamics of reactions in polar solvents. Semiclassical trajectory studies of electron-transfer and proton-transfer reactions. J. Phys. Chem. 1982, 86, 2218−2224. (26) Fö r ster, Th. 10th Spiers Memorial Lecture. Transfer Mechanisms of Electronic Excitation. Discuss. Faraday Soc. 1959, 27, 7−17. (27) Landt, L.; Klünder, K.; Dahl, J. E.; Carlson, R. M. K.; Möller, T.; Bostedt, C. Optical Response of Diamond Nanocrystals as a Function of Particle Size, Shape, and Symmetry. Phys. Rev. Lett. 2009, 103, 047402. (28) Landt, L.; Staiger, M.; Wolter, D.; Klünder, K.; Zimmermann, P.; Willey, T. M.; van Buuren, T.; Brehmer, D.; Schreiner, P. R.; Tkachenko, B. A.; Fokin, A. A.; Möller, T.; Bostedt, C. The influence of a single thiol group on the electronic and optical properties of the smallest diamondoid adamantane. J. Chem. Phys. 2010, 132, 024710. (29) Schwertfeger, H.; Fokin, A. A.; Schreiner, P. R. Diamonds are a Chemist’s Best Friend: Diamondoid Chemistry Beyond Adamantane. Angew. Chem., Int. Ed. 2008, 47, 1022−1036. (30) Yang, W. L.; Fabbri, J. D.; Willey, T. M.; Lee, J. R. I.; Dahl, J. E.; Carlson, R. M. K.; Schreiner, P. R.; Fokin, A. A.; Tkachenko, B. A.; Fokina, N. A.; Meevasana, W.; Mannella, N.; Tanaka, K.; Zhou, X. J.; van Buuren, T.; Kelly, M. A.; Hussain, Z.; Melosh, N. A.; Shen, Z. X. Monochromatic Electron Photoemission from Diamondoid Monolayers. Science 2007, 316, 1460−1462. (31) Steglich, M.; Huisken, F.; Dahl, J. E.; Carlson, R. M. K.; Henning, T. Electronic spectroscopy of FUV-irradiated diamondoids: a combined experimental and theoretical study. Astrophys. J. 2011, 729, 91. (32) Fokin, A. A.; Butova, E. D.; Chernish, L. V.; Fokina, N. A.; Dahl, J. E. P.; Carlson, R. M. K.; Schreiner, P. R. Simple Preparation of

As a general observation, it is seen that, within the Franck− Condon regime, the rate constants show a non-Arrhenius decrease with temperature, whereas beyond the Franck− Condon regime, they increase with temperature. Since the electronic coupling does not vary with temperature, the temperature dependence of the rate constants can be ascribed to vibronic effects, which show different behavior within and beyond the FC approximation. Due to the singularity of the nonadiabatic coupling vectors near crossing points of surfaces,75 an alternative diabatic approach can be used to calculate the rate of internal conversion in molecules.76 The diabatic coupling and its variation with the nuclear coordinates will be used in the near future in conjunction with the time-dependent approach at the Franck−Condon and Herzberg−Teller levels of approximation to estimate the rates of internal conversion as an extension of the present work using a similar recipe as for the radiative rate.



AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Notes

The authors declare no competing financial interest.



ACKNOWLEDGMENTS The research leading to these results has received funding from the European Research Council under the European Union’s Seventh Framework Programme (FP/2007-2013)/ERC Grant Agreement n. [320951]. The work was also supported by the Italian MIUR (FIRB 2012: “Progettazione di materiali nanoeterogenei per la conversione di energia solare”, prot.: RBFR122HFZ). We thank Prof. Peter Saalfrank (Potsdam, Germany) and Dr. Jan Boyke Schönborn (Potsdam) for many fruitful discussions. The high performance computation of the DREAMS center (http://dreamshpc.sns.it) and fruitful discussions with Dr. Camille Latouche are also acknowledged.



REFERENCES

(1) Jaffe, H. H.; Miller, A. L. The fates of electronic excitation energy. J. Chem. Educ. 1966, 43, 469−473. (2) Richter, R.; Wolter, D.; Zimmermann, T.; Landt, L.; Knecht, A.; Heidrich, C.; Merli, A.; Dopfer, O.; Reiss, P.; Ehresmann, A.; Petersen, J.; Dahl, J. E.; Carlson, R. M. K.; Bostedt, C.; Möller, T.; Mitrić, R.; Rander, T. Size and shape dependent photoluminescence and excited state decay rates of diamondoids. Phys. Chem. Chem. Phys. 2014, 16, 3070−3076. (3) Borissow, P. Journ. Russ. Phys-Chem. 1906, 37, 249. (4) Vavilov, S. I. Z. Phys. 1924, 22, 266. (5) Crosby, G. A.; Demas, J. A. Measurement of photoluminescence quantum yields. J. Phys. Chem. 1971, 75, 991−1024. (6) Perrin, F. Polarisation of fluorescence and mean life of excited molecules. J. Phys. Radium 1926, 7, 390−401. (7) Ross, J. A.; Jameson, D. M. Time-resolved methods in biophysics. 8. Frequency domain fluorometry: applications to intrinsic protein fluorescence. Photochem. Photobiol. Sci. 2008, 7, 1301−1312. (8) Hammer, M.; Schweitzer, D.; Richter, S.; Königsdörffer, E. Sodium fluorescein as a retinal pH indicator? Physiol. Meas. 2005, 26, N9−N12. (9) Ricci, M.; Bradforth, S. E.; Jimenez, R.; Fleming, G. R. Internal conversion and energy transfer dynamics of spheroidene in solution and in the LH-1 and LH-2 light-harvesting complexes. Chem. Phys. Lett. 1996, 259, 381−390. (10) Hayashi, M.; Mebel, A. M.; Liang, K. K.; Lin, S. H. Ab initio calculations of radiationless transitions between excited and ground 784

DOI: 10.1021/acs.jctc.5b01017 J. Chem. Theory Comput. 2016, 12, 774−786

Article

Journal of Chemical Theory and Computation Diamondoid 1,3-Dienes via Oxetane Ring Opening. Org. Lett. 2007, 9, 2541−2544. (33) Landt, L.; Bostedt, C.; Wolter, D.; Möller, T.; Dahl, J. E. P.; Carlson, R. M. K.; Tkachenko, B. A.; Fokin, A. A.; Schreiner, P. R.; Kulesza, A.; Mitrić, R.; Bonačić-Koutecký, V. Experimental and theoretical study of the absorption properties of thiolated diamondoids. J. Chem. Phys. 2010, 132, 144305. (34) Banerjee, S.; Stüker, T.; Saalfrank, P. Vibrationally resolved optical spectra of modified diamondoids obtained from timedependent correlation function methods. Phys. Chem. Chem. Phys. 2015, 17, 19656−19669. (35) Banerjee, S.; Saalfrank, P. Vibrationally resolved absorption, emission and resonance Raman spectra of diamondoids: a study based on time-dependent correlation functions. Phys. Chem. Chem. Phys. 2014, 16, 144−158. (36) Balaban, A. T.; Von Rage Schleyer, P. Systematic classification and nomenclature of diamond hydrocarbons-I: Graph-theoretical enumeration of polymantanes. Tetrahedron 1978, 34, 3599−3609. (37) Baiardi, A.; Bloino, J.; Barone, V. Accurate Simulation of Resonance-Raman Spectra of Flexible Molecules: An Internal Coordinates Approach. J. Chem. Theory Comput. 2015, 11, 3267− 3280. (38) Latouche, C.; Baiardi, A.; Barone, V. Virtual Eyes Designed for Quantitative Spectroscopy of Inorganic Complexes: Vibronic Signatures in the Phosphorescence Spectra of Terpyridine Derivatives. J. Phys. Chem. B 2015, 119, 7253−7257. (39) Perdew, J. P. Density-functional approximation for the correlation energy of the inhomogeneous electron gas. Phys. Rev. B: Condens. Matter Mater. Phys. 1986, 33, 8822−8824. (40) Perdew, J. P.; Burke, K.; Wang, Y. Generalized gradient approximation for the exchange-correlation hole of a many-electron system. Phys. Rev. B: Condens. Matter Mater. Phys. 1996, 54, 16533− 16539. (41) Becke, A. D. Density-functional thermochemistry. III. The role of exact exchange. J. Chem. Phys. 1993, 98, 5648−5652. (42) Ferrer, F. J. A.; Improta, R.; Santoro, F.; Barone, V. Computing the inhomogeneous broadening of electronic transitions in solution: a first-principle quantum mechanical approach. Phys. Chem. Chem. Phys. 2011, 13, 17007−17012. (43) Improta, R.; Barone, V.; Scalmani, G.; Frisch, M. J. A statespecific polarizable continuum model time dependent density functional theory method for excited state calculations in solution. J. Chem. Phys. 2006, 125, 054103. (44) Barone, V.; Cossi, M.; Tomasi, J. A new definition of cavities for the computation of solvation free energies by the polarizable continuum model. J. Chem. Phys. 1997, 107, 3210−3221. (45) Cossi, M.; Scalmani, G.; Rega, N.; Barone, V. New developments in the polarizable continuum model for quantum mechanical and classical calculations on molecules in solution. J. Chem. Phys. 2002, 117, 43−54. (46) Tomasi, J.; Mennucci, B.; Cammi, R. Quantum Mechanical Continuum Solvation Models. Chem. Rev. 2005, 105, 2999−3094. (47) Borrelli, R.; Peluso, A. The temperature dependence of radiationless transition rates from ab initio computations. Phys. Chem. Chem. Phys. 2011, 13, 4420−4426. (48) Niu, Y.; Peng, Q.; Deng, C.; Gao, X.; Shuai, Z. Theory of Excited State Decays and Optical Spectra: Application to Polyatomic Molecules. J. Phys. Chem. A 2010, 114, 7817−7831. (49) Peng, Q.; Niu, Y.; Deng, C.; Shuai, Z. Vibration correlation function formalism of radiative and non-radiative rates for complex molecules. Chem. Phys. 2010, 370, 215−222. (50) O’Rourke, R. C. Absorption of Light by Trapped electrons. Phys. Rev. 1953, 91, 265−270. (51) Niu, Y. L.; Peng, Q.; Shuai, Z. G. Promoting mode-free formalism for excited state radiationless decay process with Duschinsky rotation effect. Sci. China, Ser. B: Chem. 2008, 51, 1153−1158. (52) Baiardi, A.; Bloino, J.; Barone, V. General Time Dependent Approach to Vibronic Spectroscopy Including Franck-Condon,

Herzberg-Teller, and Duschinsky Effects. J. Chem. Theory Comput. 2013, 9, 4097−4115. (53) Reimers, J. R. A practical method for the use of curvilinear coordinates in calculations of normal-mode-projected displacements and Duschinsky rotation matrices for large molecules. J. Chem. Phys. 2001, 115, 9103−9109. (54) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Mennucci, B.; Petersson, G. A.; Nakatsuji, H.; Caricato, M.; Li, X.; Hratchian, H. P.; Bloino, J.; Janesko, B. G.; Izmaylov, A. F.; Lipparini, F.; Zheng, G.; Sonnenberg, J. L.; Liang, W.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Vreven, T.; Throssell, K.; Montgomery, J. A.Jr; Peralta, J. E.; Ogliaro, F.; Bearpark, M.; Heyd, J. J.; Brothers, E.; Kudin, K. N.; Staroverov, V. N.; Keith, T.; Kobayashi, R.; Normand, J.; Raghavachari, K.; Rendell, A.; Burant, J. C.; Iyengar, S. S.; Tomasi, J.; Cossi, M.; Rega, N.; Millam, J. M.; Klene, M.; Knox, J. E.; Cross, J. B.; Bakken, V.; Adamo, C.; Jaramillo, J.; Gomperts, R.; Stratmann, R. E.; Yazyev, O.; Austin, A. J.; Cammi, R.; Pomelli, C.; Ochterski, J. W.; Martin, R. L.; Morokuma, K.; Zakrzewski, V. G.; Voth, G. A.; Salvador, P.; Dannenberg, J. J.; Dapprich, S.; Parandekar, P. V.; Mayhall, N. J.; Daniels, A. D.; Farkas, O.; Foresman, J. B.; Ortiz, J. V.; Cioslowski, J.; Fox, D. J. Gaussian Development Version, Revision I.03; Gaussian, Inc., Wallingford, CT, 2014. (55) Frigo, M.; Johnson, S. G. Proc. IEEE 2005, 93, 216−231. (56) 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. (57) Hariharan, P. C.; Pople, J. A. Influence of polarization function on molecular-orbital hydrogenation energies. Theor. Chem. Acc. 1973, 28, 213−222. (58) Schäfer, A.; Huber, C.; Ahlrichs, R. Fully optimized contracted Gaussian basis sets of triple zeta valence quality for atoms Li to Kr. J. Chem. Phys. 1994, 100, 5829−5835. (59) Barone, V.; Biczysko, M.; Bloino, J. Fully anharmonic IR and Raman spectra of medium-size molecular systems: accuracy and interpretation. Phys. Chem. Chem. Phys. 2014, 16, 1759−1787. (60) Cances, E.; Mennucci, B.; Tomasi, J. J. A new integral equation formalism for the polarizable continuum model: Theoretical background and applications to isotropic and anisotropic dielectrics. J. Chem. Phys. 1997, 107, 3032−3041. (61) Jensen, F. Introduction to Computational Chemistry; Wiley and Sons Ltd., 2007. (62) Werner, H. J.; Knowles, P. J.; Knizia, G.; Manby, F. R.; Schütz, M.; Celani, P.; Korona, T.; Lindh, R.; Mitrushenkov, A.; Rauhut, G.; Shamasundar, K. R.; Adler, T. B.; Amos, R. D.; Bernhardsson, A.; Berning, A.; Cooper, D. L.; Deegan, M. J. O.; Dobbyn, A. J.; Eckert, F.; Goll, E.; Hampel, C.; Hesselmann, A.; Hetzer, G.; Hrenar, T.; Jansen, G.; Köppl, C.; Liu, Y.; Lloyd, A. W.; Mata, R. A.; May, A. J.; McNicholas, S. J.; Meyer, W.; Mura, M. E.; Nicklass, A.; O’Neill, D. P.; Palmieri, P.; Peng, D.; Pflüger, K.; Pitzer, R.; Reiher, M.; Shiozaki, T.; Stoll, H.; Stone, A. J.; Tarroni, R.; Thorsteinsson, T.; Wang, M.; MOLPRO, version 2012.1, a package of ab initio programs, 2012. (63) Baker, J.; Kessi, A.; Delley, B. The generation and use of delocalized internal coordinates in geometry optimization. J. Chem. Phys. 1996, 105, 192−212. (64) Kleinschmidt, M.; Tatchen, J.; Marian, C. M. SPOCK.CI: A multireference spin-orbit configuration interaction method for large molecules. J. Chem. Phys. 2006, 124, 124101. (65) Furche, F.; Ahlrichs, R.; Hättig, C.; Klopper, W.; Sierka, M.; Weigend, F. Turbomole. WIREs Comput. Mol. Sci. 2014, 4, 91−100. (66) Grimme, S.; Waletzke, M. A combination of Kohn-Sham density functional theory and multi-reference configuration interaction methods. J. Chem. Phys. 1999, 111, 5645−5655. (67) Schäfer, A.; Horn, H.; Ahlrichs, R. Fully optimized contracted Gaussian basis sets for Atoms Li to Kr. J. Chem. Phys. 1992, 97, 2571− 2577. 785

DOI: 10.1021/acs.jctc.5b01017 J. Chem. Theory Comput. 2016, 12, 774−786

Article

Journal of Chemical Theory and Computation (68) Gustav, K.; Storch, M. Non-radiative deactivation of molecules II. Theoretical study of the internal conversion rates in azulene. Int. J. Quantum Chem. 1990, 38, 1−10. (69) Cerezo, J.; Zúniga, J.; Requena, A.; Á vila Ferrer, F. J.; Santoro, F. Harmonic Models in Cartesian and Internal Coordinates to Simulate the Absorption Spectra of Carotenoids at Finite Temperatures. J. Chem. Theory Comput. 2013, 9, 4947−4958. (70) Pedone, A.; Bloino, J.; Barone, V. Role of Host-Guest Interactions in Tuning the Optical Properties of Coumarin Derivatives Incorporated in MCM-41: A TD-DFT Investigation. J. Phys. Chem. C 2012, 116, 17807−17818. (71) Santoro, F.; Improta, R.; Lami, A.; Bloino, J.; Barone, V. Effective method to compute Franck-Condon integrals for optical spectra of large molecules in solution. J. Chem. Phys. 2007, 126, 84509. (72) Zimmerman, H. E.; Kamm, K. S.; Werthemann, D. P. Mechanisms of electron demotion. Direct measurement of internal conversion and intersystem crossing rates. Mechanistic organic photochemistry. J. Am. Chem. Soc. 1975, 97, 3718−3725. (73) Heinzelmann, W.; Labhart, H. Triplet-triplet spectra and triplet quantum yields of some aromatic hydrocarbons in liquid solution. Chem. Phys. Lett. 1969, 4, 20−24. (74) Sobolewski, A. L.; Domcke, W.; Dedonder-Lardeux, C.; Jouvet, C. Excited-state hydrogen detachment and hydrogen transfer driven by repulsive 1πσ* states: A new paradigm for nonradiative decay in aromatic biomolecules. Phys. Chem. Chem. Phys. 2002, 4, 1093−1100. (75) Nakamura, H.; Truhlar, D. G. Direct diabatization of electronic states by the fourfold way. II. Dynamical correlation and rearrangement processes. J. Chem. Phys. 2002, 117, 5576−5593. (76) Islampour, R.; Miralinaghi, M. Theoretical Study of Internal Conversion Decay Rates Associated with the Three Lowest Singlet Electronic States in Pyrazine. J. Phys. Chem. A 2009, 113, 2340−2349.

786

DOI: 10.1021/acs.jctc.5b01017 J. Chem. Theory Comput. 2016, 12, 774−786