Vibronic Effects on Rates of Excitation Energy Transfer and Their

Mar 29, 2016 - ... the transition densities (determined using the regular TD-DFT formulation) of the .... Table 1 shows the computed values of the tot...
0 downloads 0 Views 1MB Size
Article pubs.acs.org/JCTC

Vibronic Effects on Rates of Excitation Energy Transfer and Their Temperature Dependence 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: Vibronic effects on rate constants governing excitation energy transfer between different electronic states have been studied within the adiabatic regime and the harmonic oscillator approximation, possibly including bulk solvent effects with the polarizable continuum model. A recent implementation in the GAUSSIAN package has been extended for this purpose, with the possibility of taking into account frequency shifts and mode mixing effects between ground and excited electronic states. The potentialities of the new computational tool have been analyzed for a number of case studies, including excimers of naphthalene and zinc porphyrin complexes. In all cases, the change of rate constants with temperature arising from the inclusion of vibronic effects is in good qualitative agreement with the available experimental data.

1. INTRODUCTION Photophysical and photochemical processes play a central role in molecular systems. Two of the most widely studied processes of this kind, which are also utilized in large scale industrial applications, are electronic absorption and emission (especially fluorescence). An interesting radiation-less transition connected to the above-mentioned optical processes, often taking place in biological systems, is known as the excitation energy transfer (EET).1 EET involves the transfer of energy from one electronic state of a molecule, the donor, to an electronic state of another one, the acceptor. The donor and acceptor can be different molecules, in which case the process is known as intermolecular EET; but often, they are two moieties of the same molecule, joined by a bridging unit, resulting in an intramolecular EET. Intermolecular EET can also take place between excited states of the monomeric units of a dimer; such dimers being often stabilized only in the excited states by what is known as an excimeric coupling between the states.1 Naphthalene and anthracene dimers represent paradigmatic examples of excimeric coupling and EET.1 Castanheira and Martinho2 have studied the solvatochromic shift of excimers of naphthalene and pyrene in different solvents at room temperature, showing that such excimers are stable in solutions. Intermolecular EET is very common in light harvesting complexes, like photosystems in bacteria and plants,3 where solar energy absorbed by the chromophores is transferred as excitation energy to the reaction centers by antenna proteins.3 Porphyrins are a class of molecules resembling the photosynthetic units of chlorophyll in addition to possessing interesting properties, like stable excited states,4 efficient excitation energy transfer,4,5 etc. In addition, they are good absorbers of visible radiation.4 These properties can be tuned by introduction of metal atoms like zinc in the core,6 similar to © 2016 American Chemical Society

the presence of magnesium in chlorophyll, making them useful as artificial light harvesting compounds.7 Therefore, the EET rates of zinc-porphyrin and related complexes have been subjects of several experimental and theoretical investigations. Intramolecular EET between triplet excited states of zincporphyrin (ZnP) or free base porphyrin (H2P) units connected by bridging phenylacetylene units have been studied by Eng and co-workers.8 The authors have performed a systematic study of the variation of the electronic coupling with the angle between the plane of the porphyrin rings and the bridges. The coupling was calculated as half the difference of excitation energy between the lowest singlet to triplet transition in the donor and acceptor units, which was then varied as a function of the cosine of the dihedral angle. Although the experimental rates were measured for ZnP-xB-H2P systems, i.e., EET between the first excited triplet (T1) states of different moieties, the theoretical calculations of electronic coupling were performed on symmetrical units joined by bridges, i.e. for ZnP-xB-ZnP and/or H2P-xB-H2P systems. The number of phenylacetylene units x was varied from 2 to 5. The alkyl and aryl substituents on the vinyl units of the porphyrin moieties were not considered in the theoretical simulations, because they were assumed to have negligible effects on the electronic coupling. Various models for the calculation of EET rates have been proposed, the most widely used being based on Förster’s theory. In this framework, the electronic coupling is proportional to the product of the transition dipole moments of the donor and acceptor and inversely proportional to the cube of the distance between them. The rate constant is calculated from Received: February 12, 2016 Published: March 29, 2016 2357

DOI: 10.1021/acs.jctc.6b00157 J. Chem. Theory Comput. 2016, 12, 2357−2365

Article

Journal of Chemical Theory and Computation

additional Coulomb screening effect, which is often approximated as the square of the refractive index n, i.e. 1 κμD μA R= 4π ϵ0 n2r 3 (3)

the electronic coupling and the so-called spectral overlap between the area-normalized donor-emission and acceptorabsorption spectra (not accounting for static disorders and coherence9). A more accurate yet feasible computational model of the electronic coupling beyond the point-dipole approximation has been recently formulated by Curutchet, Mennucci, and co-workers in terms of Coulomb, exchange-correlation, and overlap interactions between the donor and acceptor units.1 In most of those works, the spectral overlaps have been adapted from the experimental absorption and emission spectra or in some cases by vertical approaches not involving excited state optimization and/or force constant calculation. In this work, we have calculated the donor-emission and acceptorabsorption spectra by explicit evaluation of vibronic overlap integrals including frequency changes and mode mixing (the Duschinsky effect10) in the excited state, within the adiabatic regime. More general approaches which also take charge transfer effects into consideration, have been proposed, for instance, in refs 11, 12, and 13. Such effects might influence the rates in systems like H2P-xB-H2P where the donor and acceptor are linked by bridging atoms, but as we shall see later, the measured rate constants can be predicted with a fair degree of accuracy neglecting such effects, especially for the bridge lengths considered in the present work. The spectral overlap and the electronic coupling factor have been then used to calculate the rate of EET at different temperatures. The harmonic approximation has been used throughout, and the dependence of the electronic coupling on the nuclear coordinates has not been considered. The paper is organized as follows. The theory is outlined in section 2, while section 3 contains the computational details. The results are discussed in section 4, and section 5 summarizes the work.

Recently, an improvement over the above-mentioned scheme has been proposed by Mennucci and co-workers, to properly take into account the electron density and the effects of electron exchange and correlation. Casida’s equations16 generalized for two chromophores with a common excitation frequency are used, a condition which is of course fulfilled when the donor and acceptor are the same molecule, for instance, symmetrical units joined by a bridge or for dimers. The interaction between chromophores can be split into Coulomb (RCoul), exchange-correlation (RXC), and overlap (ROvlp) contributions1 R = RCoul + RXC − ROvlp Coul

2π 2 RS ℏ

RCoul = RXC =

1 κμD μA 4π ϵ0 r 3

∫r d rd r′ρA*(r) |r −1 r′| ρD (r′) ∫r d rd r′ρA*(r)gXC(r, r′)ρD (r′)

ROvlp = − ω0

∫r d rρA*(r)ρD (r)

(5)

Here, ρA(r) and ρD(r′) are the transition densities (determined using the regular TD-DFT formulation) of the noninteracting acceptor and donor respectively, and gXC(r, r′) is the exchangecorrelation kernel. ω0 is the resonant frequency of the noninteracting chromophoric units. For EET between excited states of nonidentical units, the average of the resonant frequencies is considered. Solvent effects are accounted for, by additional coupling terms. The explicit solvent contribution is a three-body Coulomb interaction between the two chromophores and the solvent. The reaction field of the solvent can also modify the other three terms (RCoul, RXC, and ROvlp), through what is known as its implicit contribution. For systems where the donor and acceptor moieties are quite close to each other, more than two electronic states might be involved in the process of energy transfer. Additionally, the overlap of the densities of the donor and acceptor fragments could become non-negligible. Electronic couplings calculated using a localized operator partitioning method17,18 address both issues. If the electronic coupling is not too strong, as the cases considered in the present work, the formalism reduces to the Förster limit.17 2.2. The Spectral Overlap. As mentioned before, the spectral overlap S is calculated as

(1)

where R is the electronic coupling, and S is the overlap integral between the normalized absorption and emission bandshapes of the acceptor and the donor, respectively. We will first focus on the electronic coupling and return to the overlap integral later. In its most simple form, R depends on the strength of the transition electric dipole moments of the acceptor (μA) and donor (μD), their relative orientation described by an orientation factor (κ) and the cube of the distance (r) between them15 R=

(4)

Ovlp

where R , R , and R are determined from the transition frequencies and transition densities of the noninteracting systems

2. THEORY 2.1. Rate Expression and Electronic Coupling. In the classical picture, the interaction between two oscillating dipoles gives rise to a progressive transfer of excitation energy from the initially oscillating donor-dipole to the initially static acceptordipole. This phenomenon is known as excitation energy transfer (EET). The quantitative formulation of the rate expression of EET was developed by Förster.14 According to Förster’s theory, the rate constant for EET is derived from the Fermi’s golden rule and has the form15 kEET =

XC

S=

∫ SA(ω)SD(ω)dω

(6)

where SA(ω) is the area-normalized absorption bandshape of the acceptor, while SD(ω) is the area-normalized emission bandshape of the donor, which can be connected to the molar absorptivity coefficient and the fluorescence intensity, respectively. In general, the spectra are computed from single electronic transitions, which are then broadened by means of symmetric

(2)

Here, ϵ0 is the permittivity of vacuum. The presence of solvent effects in the different components of R introduces an 2358

DOI: 10.1021/acs.jctc.6b00157 J. Chem. Theory Comput. 2016, 12, 2357−2365

Article

Journal of Chemical Theory and Computation

where g is the final state gradient, computed at the ground state geometry, and F is related to the Cartesian, excited state Hessian (FX ) calculated at the ground state geometry in the following way

Gaussian or Lorentzian distribution functions. However, as already shown by some of us,19 even low resolution spectra can exhibit an asymmetric shape due to the vibronic contributions. Thus, in order to get a reliable description of the electron transfer rate, it is necessary to correctly account for the spectral bandshapes and, hence, of the vibronic contributions. Extensive work on this subject has made possible the simulation of such spectra for medium- and even large-sized systems. We will present here some key points on their calculation. Vibronic absorption and emission spectra can be simulated by time independent (TI) and time dependent (TD) models. In the TI approaches, the spectrum is computed as a sum of individual vibronic transitions. On the other hand, in the TD approaches, the spectrum is obtained as the Fourier transform of the appropriate autocorrelation function. The main advantage of the TD approaches over the TI ones is the straightforward inclusion of temperature effects without additional computational cost. Another important aspect is that band assignment, which is a significant advantage of the TI models over the TD ones, is of no interest in the computation of kEET. In the following, the variation of kEET with temperature will be discussed; therefore, the TD-based approach proposed in our group will be used. Additional details can be found in refs 20 and 21. This TD approach relies on the harmonic approximation for both the initial and final potential energy surfaces (PES). In fact, within this approximation, the data needed to simulate the spectra are only the harmonic frequencies of the two states (Ω̅ and Ω ) together with their normal mode transformation matrices (L̅ and L ). Initial and final state parameters are denoted by single and double overbars, respectively. The normal modes are used to derive a linear relation between the normal coordinates of the states (Q̅ and Q ), usually known as the Duschinsky transformation Q̅ = JQ + K

F = L̅ −1M−1/2FX M−1/2(L̅ T)−1

However, this might lead to imaginary vibrational frequencies, as the Hessian is not computed at a stationary point of the PES and the harmonic approximation is not sufficient to treat those modes. In such cases, one can use the so-called linear coupling method (LCM) also known as the vertical gradient (VG) approach,21 where the Duschinsky matrix is assumed to be the identity matrix and the shift vector is approximated from the final-state gradient (g ), together with the normal modes and frequencies of the initial state −2

K = −Ω̅ L̅ TM−1/2g

(7)

J = L̅ TL (8)

where M is the diagonal matrix of the atomic masses, while X̅ and X denote the equilibrium structures of the Cartesian coordinates in the initial and final states, respectively. However, due to the size and flexibility of systems under investigation, the PES of the initial and final states may be shifted by a large extent, and, in such cases, the adiabatic model may face strong difficulties, and vertical approaches described below become more appropriate. In the vertical Hessian (VH) model, the Duschinsky matrix is obtained by diagonalizing the final state Hessian calculated at the initial state equilibrium geometry. The shift vector can be computed as −1

K = −F L̅ TM−1/2g

(11)

Due to the approximation of the final state frequencies and the lack of information about the shift between the minima of the two states, the vertical model often predicts incorrect band position or fine-structure.19 Furthermore, in the VG model, since the final state PES is assumed to be equal to the initial state PES, for the case of emission, one has to know the excited state minimum geometry and the corresponding normal modes. Thus, the computational gain which usually comes with a VG absorption spectrum is not achieved. To retain the same advantage, the excited state can be approximated by the shifted ground state as for absorption, the shifts being opposite in sign than in absorption. Vertical models usually have a tendency to compensate the inaccuracy of the adiabatic ones in describing large amplitude deformations.22 However, since all the systems studied in this work are semirigid, the adiabatic model will be used as a standard for all the vibronic simulations. In the presence of multiple electronic states vertical models are more consistent than their adiabatic counterparts. Finally, an important aspect, briefly discussed before, is connected to the simulation of the band broadening. Several strategies can be devised to include those effects. The most common one is to employ Lorentzian and Gaussian distribution functions with empirically set half widths at halfmaximum (HWHM). An interesting aspect of the TD approach on this matter is the possibility to combine easily both functions to obtain complex effects. Alternatively, more refined methods have been proposed in the literature to evaluate the best suited bandwidth. For polar solvents, the effect of inhomogeneous broadening can be calculated from the solvent relaxation energy,23 which is the difference between the excitation energies obtained from state specific nonequilibrium and equilibrium solvation models. In the first case, the slow degrees of freedom of the solvent are assumed to have not equilibrated with the solute density in the excited state, whereas full equilibration is taken into account for the latter case. The state specific solvent relaxation energy is calculated as23

Different methods can be used to derive expressions for the Duschinsky matrix, J and the shift vector, K, depending on the reference geometry used for the expansion of the final state PES. Within the so-called adiabatic models, the final state PES is expanded about its equilibrium geometry, and J and K are given by

K = L̅ TM−1/2(X − X̅ )

(10)

λ = λneq − λeq

(12)

where λneq and λeq are the state specific excitation energies within the nonequilibrium and equilibrium solvation regimes, respectively. The broadening is then simulated using a Gaussian whose half width at half-maximum can be calculated as23

(9) 2359

DOI: 10.1021/acs.jctc.6b00157 J. Chem. Theory Comput. 2016, 12, 2357−2365

Article

Journal of Chemical Theory and Computation

W = 4 ln(2)kBTλ

(except certain test cases mentioned in the relevant result sections).

(13)

thus accounting for solvent broadening at different temperatures (T). Otherwise, the HWHM may be set by a direct comparison of the simulated spectra with the experimental ones. In this work, we have illustrated one exemplary case for each choice, where we have used the solvent relaxation energy and direct comparison with experimental spectra to compute the broadened absorption and emission spectra. A recent work has shown that convolution of individually broadened line shapes for the solute (using empirical broadening) and the solvent (from the distribution of vertical transition energies taken at different instants of time during a molecular dynamics simulation of the solvent considering the solute to be fixed at its ground state minimum geometry)24 provides another way to predict accurate broadening parameters. However, since this would lead to increased computational cost, particularly for calculating kEET, where we need to compute two spectra for each system, we do not consider that approach in the present work.

4. RESULTS 4.1. Rate of Intermolecular EET: Naphthalene Dimer in the Gas Phase. In analogy with the calculations carried out by Mennucci and co-workers,1 the ground electronic state (S0) of naphthalene was optimized at the B3LYP/6-31+G(d) level of theory in the gas phase. The coupling between the first excited singlet (S1) state of the two monomers was calculated at the PBE0/6-311+G(d,p) level of theory, as done in ref 1. The distance between the monomers was taken as 4 Å. The separated, planar monomers are assumed to be in the eclipsed conformation, i.e., the distance of separation between them gives the distance between their centers (see Figure 1).

3. COMPUTATIONAL DETAILS Geometry optimizations and force constant calculations at the equilibrium geometries of the ground and excited states have been performed at the DFT and TD-DFT levels of theory using a development version of GAUSSIAN.25 The choice of functionals and basis sets was made in such a way that the calculations remained consistent with previous studies, especially where direct comparisons are sought. The B3LYP,26 CAM-B3LYP,27 and PBE028 functionals have been used for different systems, mostly in conjunction with the polarized, double-ζ 6-31G(d)29 basis set, possibly augmented with diffuse functions. Solvent effects (dichloromethane, with dielectric constant of 8.93 and refractive index of 1.42 at 293 K and cyclohexane, with dielectric constant of 2.02 and refractive index of 1.43) have been taken into account, using the integral equation formalism of the polarizable continuum model (IEF-PCM)30 as implemented in GAUSSIAN. The default GAUSSIAN parameters have been used for the convergence criteria and fineness of the grid for calculating integrals. Electronic couplings have been calculated at the TD-DFT level of theory. Since most of the test or practical cases involve intramolecular energy transfer between two fragments of the same molecule, the fragments are saturated at the connection point by adding phantom hydrogen atoms, and the couplings between the different excited states are calculated. The donor-emission and acceptor-absorption spectra are calculated from the same fragments as mentioned above. The GAUSSIAN vibronic code has been extended for the computation of spectral overlaps and therefore, kEET. In parallel, an independent, recent FORTRAN implementation of the TDapproach was used for calculating the spectra, followed by their overlap and hence, kEET. The code can be interfaced with any quantum chemical package which can calculate the Duschinsky matrix and shift vector between the normal modes of the two states. Temperature effects have been accounted for in the implemented code. A time-step of 0.1 fs has been used to compute the correlation functions for a total propagation time of 100 ps, and the fastest Fourier transformation in the West (FFTW31) package has been used for the Fourier transformations. The correlation functions have been damped with Lorentzians, whose HWHM lie between 80 and 100 cm−1

Figure 1. Eclipsed conformation of the naphthalene dimer. Brown balls represent carbons while white balls represent hydrogens.

Table 1 shows the computed values of the total coupling, including the comparisons with the values provided in Figure 4 of ref 1. Table 1. Comparison of Our Calculated Value of Total Coupling between the S1 States of the Monomers of Naphthalene Dimer, Separated by 4 Å, with the Calculated Value of Ref 1a distance (Å)

states involved

R (eV)

R (eV) (ref 1)

4

S1−S1

0.0338

0.03 −0.01

kEET (× 108 s−1)

gas solution (DCM) gas solution (DCM)

0 0 300 300

9.41 4.68 (0.50) 6.47 3.68 (0.57)

Table 4. Comparison of kEET (Ratio in Parentheses) and Coupling Values in Dichloromethane Solution at 0 K, Using Different Models for the Electronic Couplinga vibronic effects

coupling

R (eV)

kEET at 0 K, s−1

AH|FC

full dipole

0.029 0.025

4.68 × 108 3.22 × 108 (0.69)

a The full coupling refers to that calculated using eq 4, while the dipole coupling has been calculated using eq 2.

Table 2. Comparison of Calculated Values of Total and Explicit Solvent Mediated Coupling between the S1 States of the Monomers in Naphthalene Dimer in Dichloromethane Solution, at 4 Å Separation, with the Calculated Values of Ref 1a 4

temp (K)

have also been provided. A general observation is that kEET is larger in the gas phase than in DCM solution, and the variation of kEET with temperature is more pronounced in the gas phase. It is noteworthy that in a recent work,34 we observed a similar temperature dependence for the rate of intersystem crossing, and, in that case also, the variation was found to be more pronounced in the gas phase than in solution. It is seen that the ratio between the kEET computed in solution and in vacuum is lower than 1.0 as in Figure 6 of ref 1. However, in that work, the ratio did not depend on temperature (always 0.9 for 4 Å separation) because the spectral overlap of naphthalene was assumed to be the same in the gas and solution phases, hence the ratio became equal to the ratio of the electronic couplings in the two phases. Here, we see that when the spectral overlap is taken into account, the ratio of kEET changes considerably and varies with temperature. Considering only the ratio of the electronic coupling factor, kEET,DCM/kEET,gas amounts to 0.88, almost equal to the value of 0.9 as reported in ref 1. Dichloromethane being a moderately polar solvent, naphthalene dimer in dichloromethane was used as a test system to calculate the inhomogeneous broadening from the solvent relaxation energy (eqs 12 and 13). For naphthalene, the solvent relaxation energy is approximately 43.5 cm−1, which amounts to a half width at half-maximum of 317 cm−1 at 300 K. The corresponding kEET is 5.27 × 108 s−1, i.e., we observed an increase compared to the value of 3.68 × 108 s−1, though the order of magnitude remained unchanged. 4.2.1. The Dipole Approximation for Naphthalene Dimer in DCM. Table 4 shows the values of kEET using different

Figure 2. Normalized emission and absorption spectra of naphthalene at 300 K as functions of energy.

distance (Å)

phase

models for calculating the electronic coupling R. The orientation factor has been assumed to be 1.0 for calculating the dipole−dipole coupling, i.e., the dipole moments of the naphthalene monomers are parallel to each other, separated by 4 Å. The dipole coupling is comparable to the coupling calculated by eq 4; and the same applies to the corresponding rate constants. Although from the above observations, it seems that the dipole coupling can be confidently used for obtaining the rate constants, it has one important limitation. For example,

a

The reference values have been taken again from a plot provided in that reference, so only approximate values could be reported. 2361

DOI: 10.1021/acs.jctc.6b00157 J. Chem. Theory Comput. 2016, 12, 2357−2365

Article

Journal of Chemical Theory and Computation eq 4 allows also to calculate the coupling between triplet electronic states, whereas the dipole approximation fails in this case, because the transition moment for the S0 → Tn states vanishes, due to symmetry restrictions. 4.2.2. Determining the Broadening Factor: Naphthalene Dimer in Cyclohexane. Cyclohexane being a nonpolar solvent, the solvent relaxation energy for S0 → S1 excitation of naphthalene in cyclohexane solution is negligible. Therefore, a direct comparison to the available experimental absorption and/or emission spectra can be made to determine the broadening factor. Equilibrium geometries, normal modes, and vibrational frequencies of the ground and excited states, obtained at the PBE0/6-31+G(d) level of theory, in cyclohexane solution, have been used to calculate kEET. Comparison with experimental35 absorption and emission spectra showed that a Lorentzian with a HWHM of 260 cm−1 is suitable for reproducing the overall line shape and fine-structure. At 298 K, the calculated value of kEET in cyclohexane solution is 5.48 × 108 s−1. The electronic coupling in this case, obtained from eq 4, is approximately 0.029 eV, as in DCM solution (see Figure 4).

Figure 5. CAM-B3LYP/6-31G(d)-optimized ground state geometry of the complex studied here. The brown and white balls represent the carbon and hydrogen atoms, respectively, while the large yellow ball represents zinc. The blue balls represent N atoms.

Table 5. Vertical Transition Energies (TD-CAM-B3LYP/631G(d)) from the Ground State to the Third Singlet Excited State, for the Complex molecule

transition

excitation energy (eV)

zinc-porphyrin free-base porphyrin

S0 → S3 S0 → S3

3.54 3.45

cm−1 was taken from ref 36, and the rate-constant of EET was calculated as 4.35 × 108 s−1 at 300 K. Although experimental or theoretical results are not available for this system, to the best of our knowledge, this value can be compared to the measured kEET value of 2.80 × 108 s−1 for S1− S1 EET in a similar complex where the zinc-porphyrin and freebase porphyrin units (with additional substituents on the phenyl rings) are separated by two phenylacetylene bridges connected to either side of a phenyl group.37 4.3.2. Triplet−Triplet Coupling. We calculated kEET for the T1-T1 energy transfer in ZnP-2B-H2P. As shown in ref 8, to maximize the coupling between the zinc-porphyrin fragments, the dihedral angle between the planes of the porphyrin units and the bridging phenylacetylene units was taken to be around 60° (65° in our case at the optimized geometry). All electronic structure calculations were done at the B3LYP/6-31G(d) level of theory, which was also used to calculate the electronic coupling. The optimized ground state geometry is shown in Figure 6. Regarding the vibronic simulations, the AH model has been used to compute the absorption and emission spectra.

Figure 4. Digitalized experimental absorption spectrum of naphthalene compared with the calculated one. The measurement and simulation have been carried out at 298 K in cyclohexane solution. The theoretical spectrum has been blue-shifted by 17 nm.

4.3. Rate of Intramolecular EET in Porphyrin Complexes. 4.3.1. Singlet−Singlet Coupling. To study the rate of intramolecular excitation energy transfer between two fragments of a complex joined by a bridging unit, a zinc-bound porphyrin complex acting as the donor and free-base porphyrin acting as the acceptor have been considered. The electronic coupling between the S3 states of the donor and acceptor moieties, which are involved in the EET, was calculated recently by Mennucci and co-workers.36 We have selected the complex labeled as Mc in ref 36, i.e., the entire bridge has been considered in the calculation of the electronic coupling. The structure of the complex, optimized at the CAM-B3LYP/631G(d) level of theory, is shown in Figure 5. The S0 → S3 vertical excitation energies at the respective ground state equilibrium geometries (CAM-B3LYP/6-31G(d) as in ref 36) are shown in Table 5. Our calculated values show excellent agreement with the ones given in Table 2 of ref 36. The donor-emission and acceptor-absorption spectra were calculated, followed by the computation of spectral overlap. The S3−S3 coupling of 156.4

Figure 6. B3LYP/6-31G(d)-optimized ground state geometry of ZnP2B-H2P. The brown and white balls represent the carbon and hydrogen atoms, respectively, while the large yellow ball represents zinc. The blue balls represent N atoms.

The interchromophoric separation was approximately 12.5 Å. Our calculated electronic coupling for ZnP-2B-H2P was 36.2 cm−1, which can be compared to the value of 25 cm−1 for ZnP2B-Zn, calculated in ref 8. The spectral overlap between the ZnP emission and H2P absorption spectra (involving the T1 state of each unit) was calculated as 3.61 × 107 s−1 at 150 K (to reproduce the experimental conditions used in ref 8), in fair agreement with the experimental value of 2.0 × 107 s−1. At 0 K, 2362

DOI: 10.1021/acs.jctc.6b00157 J. Chem. Theory Comput. 2016, 12, 2357−2365

Article

Journal of Chemical Theory and Computation

Figure 7. Comparison of the overlap between the normalized donor emission and acceptor absorption spectra involving T1 and S0 states of ZnP-2BZnP (left) and ZnP-2B-H2P (right), at 150 K. As seen, the overlap is greater for the latter than the former, due to the red-shifted absorption spectrum of H2P.

kEET was computed to be 3.33 × 107 s−1. Although this is in contrast to the expected temperature dependence of kEET, it has to be noted that the moieties are different, so the effect of temperature on the vibronic contributions might increase the magnitude of the spectral overlap instead of decreasing it. We have next calculated kEET for the symmetrical ZnP-2BZnP system at 150 K. The electronic coupling remained essentially unaltered, but kEET decreased to 2.70 × 107 s−1. This can be rationalized in terms of the decreased spectral overlap between the acceptor-absorption and donor-emission spectra for the symmetrical ZnP-2B-ZnP complex. For the previous case of ZnP-2B-H2P, the absorption spectrum of H2P was redshifted compared to that of ZnP, resulting in a larger overlap with the emission spectrum of ZnP, as shown in Figure 7. The EET rate for ZnP-2B-H2P has been calculated also using the VG model. The absorption and emission spectra of H2P and ZnP are reported in Figure 8. It is noteworthy that the VG

above, the AH model can represent the most suitable approach for computing kEET at least for semirigid systems. 4.4. Rate Constants with Electronic Couplings Computed at Different Levels of Theory. Caricato and co-workers have shown in a recent work15 that the coupling between the donor and acceptor electronic states can be calculated at the CCSD (coupled cluster singles doubles38) level of theory using their implementation in GAUSSIAN. Solvent effects can be also taken into account when needed. The authors compared the coupling between different excited states involving π → π* transitions in indole dimer, for TDDFT(B3LYP) and CCSD excitations. Their findings indicate that for indole, the CCSD level of theory is more accurate in predicting the correct order of excited states and the associated transition electric dipole moments. A number of possible conformations of the dimer and various distances between the monomers have been studied in the simulations, and the dependence of the electronic coupling on those factors has been investigated. Here, we have calculated the rate of excitation energy transfer between the monomers in the indole dimer, corresponding to a separation of 6 Å. The vibronic effects were treated at the DFT/TD-DFT level. The ground state optimization and force constant calculation were performed at the B3LYP/aug-ccpVDZ39 level of theory to remain consistent with the calculations of ref 15. From the transition dipole moments and oscillator strengths corresponding to vertical TD-DFT calculations, it was found that the first excited state corresponds to the La state described in ref 15. The coupling between the S1 states of the monomers at 6 Å separation in the so-called stacked conformation was found to be 86.45 cm−1 at the TDB3LYP level. This agrees well with the calculated value of ref 15 (Figure 4). Finally, kEET at 0 K was found to be 1.82 × 108 s−1. A direct comparison with experiment is not possible here, but tryptophan units (whose main chromophoric parts are indole moieties), where EET is quite common,40 are known to have excited state lifetimes in the ns range.40 As shown in ref 15, at the CCSD level of theory, the coupling increases further, and averaging over different conformations increases the electronic coupling to about 130 cm−1, which would increase the value of kEET.

Figure 8. Normalized absorption spectrum (black) of H2P and emission spectrum (red) of ZnP, calculated using the VG approximation at 0 K.

spectra are significantly narrower than their AH counterparts, and, as a consequence, the overlap is decreased, and kEET at 0 K falls by an order of magnitude, to 2.19 × 106 s−1. This is due to the intrinsic approximation of the VG models, since the equilibrium position of the final state is extrapolated and mode mixing effects are neglected. The more refined VH model could not be employed here due to the presence of an imaginary frequency in the final state PES. Therefore, as mentioned

5. CONCLUSION Excitation energy transfer is an important nonradiative pathway for decay of excited states to other near-lying ones, along with 2363

DOI: 10.1021/acs.jctc.6b00157 J. Chem. Theory Comput. 2016, 12, 2357−2365

Journal of Chemical Theory and Computation



the other widely studied radiationless processes taking place near conical intersections, such as internal conversion facilitated by nonadiabatic coupling41,42 and intersystem crossing facilitated by spin−orbit coupling.43 Rates of intra- and intermolecular excitation energy transfer between two electronic states and their temperature dependence have been studied within Fermi’s golden rule framework, i.e., considering the first order rate constant to be proportional to the product of the square of the coupling between the electronic states and the overlap between the area-normalized donor-emission and acceptor-absorption spectra. The couplings between different electronic states (both singlet and triplet) of two chromophoric units, one acting as the donor and the other as the acceptor, have been calculated using a recent implementation in the GAUSSIAN code, with contributions from Coulomb, exchangecorrelation, and overlap interactions between the chromophoric units. Processes taking place in solution have been simulated with the IEF-PCM model, where the solvent’s implicit and explicit contributions influence the coupling. Dipole−dipole coupling can also be used as a measure of the electronic coupling between singlet states. The vibronic bandshapes have been computed at the harmonic level, within the TD approach, using both adiabatic and vertical models for the final state PES. Temperature effects have been taken into account and have been found to influence the magnitudes of the rate constants, the variation being more pronounced in the gas phase than in solution. Computation of the spectral overlap (and hence kEET) has also been added to the vibronic code already available in the GAUSSIAN package within both adiabatic and vertical models. The above formalism permits to build a fairly accurate tool for predicting the rates of excitation energy transfer for diverse classes of compounds, ranging from excimeric coupling in naphthalene dimers to bridged porphyrin and/or zincporphyrin complexes, which represent models for artificial photosynthesis. For small molecules, the CCSD level of theory, which is a very accurate yet expensive computational method, can be used for obtaining improved electronic couplings. The vibrational contributions should be included to complement the increased accuracy of electronic structure computations with the aim of predicting reliable rate constants. However, as shown from the previous results, for most practical cases, electronic couplings computed within the TD-DFT regime are sufficient to predict fairly accurate values of kEET when vibronic effects are taken into account.



Article

REFERENCES

(1) Iozzi, M. F.; Mennucci, B.; Tommasi, J.; Cammi, R. J. Chem. Phys. 2004, 120, 7029−7040. (2) Castanheira, E. M. S.; Martinho, J. M. G. J. Photochem. Photobiol., A 1994, 80, 151−156. (3) Mennucci, B.; Curutchet, C. Phys. Chem. Chem. Phys. 2011, 13, 11538−11550. (4) Peng, H. Q.; Niu, L. Y.; Wu, Y. Z.; Tung, C. H.; Yang, Q. Z. Chem. Rev. 2015, 115, 7502−7542. (5) Terazono, Y.; Kodis, G.; Bhushan, K.; Zaks, J.; Madden, C.; Moore, A. L.; Moore, T. A.; Fleming, G. R.; Gust, D. J. Am. Chem. Soc. 2011, 133, 2916−2922. (6) Tomizaki, K. Y.; Loewe, R. S.; Kirmaier, C.; Schwartz, J. K.; Retsek, J. L.; Bocian, D. F.; Holten, D.; Lindsen, J. S. J. Org. Chem. 2002, 67, 6519−6534. (7) Balaban, T. S.; Goddard, R.; Linke-Schaetzel, M.; Lehn, J. M. J. Am. Chem. Soc. 2003, 125, 4233−4239. (8) Eng, M. P.; Ljungdahl, T.; Martensson, J.; Albinsson, B. J. Phys. Chem. B 2006, 110, 6483−6491. (9) Scholes, G. D.; Jordanides, X. J.; Fleming, G. R. J. Phys. Chem. B 2001, 105, 1640−1651. (10) Duschinsky, F. Acta Physicochim. URSS 1937, 7, 551−566. (11) Pereverzev, A.; Bittner, E. R.; Burghardt, I. J. Chem. Phys. 2009, 131, 034104. (12) Izmaylov, A. J.; Mendive-Tapia, D.; Bearpark, M. J.; Robb, M. A.; Tully, J. C.; Frisch, M. J. J. Chem. Phys. 2011, 135, 234106. (13) Endicott, J. S.; Joubert-Doriol, L.; Izmaylov, A. F. J. Chem. Phys. 2014, 141, 034104. (14) Förster, Th. Ann. Phys. (Paris) 1948, 2, 55−75. (15) Caricato, M.; Curutchet, C.; Mennucci, B.; Scalmani, G. J. Chem. Theory Comput. 2015, 11, 5219−5228. (16) Casida, M. E.; Jamorski, C.; Casida, K. C.; Salahub, D. R. J. Chem. Phys. 1998, 108, 4439−4449. (17) Nagesh, J.; Izmaylov, A. F.; Brumer, P. J. Chem. Phys. 2015, 142, 084114. (18) Khan, Y.; Brumer, P. J. Chem. Phys. 2012, 137, 194112. (19) Pedone, A.; Bloino, J.; Barone, V. J. Phys. Chem. C 2012, 116, 17807−17818. (20) Niu, Y.; Peng, Q.; Deng, C.; Gao, X.; Shuai, Z. J. Phys. Chem. A 2010, 114, 7817−7831. (21) Baiardi, A.; Bloino, J.; Barone, V. J. Chem. Theory Comput. 2013, 9, 4097−4115. (22) Cerezo, J.; Zùñiga, J.; Requena, A.; Ferrer, F. J. A.; Santoro, F. J. Chem. Theory Comput. 2013, 9, 4947−4958. (23) Ferrer, F. J. A.; Improta, R.; Santoro, F.; Barone, V. Phys. Chem. Chem. Phys. 2011, 13, 17007−17012. (24) Cerezo, J.; Ferrer, F. J. A.; Prampolini, G.; Santoro, F. J. Chem. Theory Comput. 2015, 11, 5810−5825. (25) 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. (26) Lee, C.; Yang, W.; Parr, R. G. Phys. Rev. B: Condens. Matter Mater. Phys. 1988, 37, 785−789.

AUTHOR INFORMATION

Corresponding Authors

*E-mail: [email protected]. *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). The high performance computation of the DREAMS center (http://dreamshpc.sns.it) is acknowledged. 2364

DOI: 10.1021/acs.jctc.6b00157 J. Chem. Theory Comput. 2016, 12, 2357−2365

Article

Journal of Chemical Theory and Computation (27) Yanai, T.; Tew, D. P.; Handy, N. C. Chem. Phys. Lett. 2004, 393, 51−57. (28) Adamo, C.; Barone, V. J. Chem. Phys. 1999, 105, 9982−9985. (29) Hariharan, P. C.; Pople, J. A. Theor. Chem. Acc. 1973, 28, 213− 222. (30) Cances, E.; Mennucci, B.; Tomasi, J. J. J. Chem. Phys. 1997, 107, 3032−3041. (31) Frigo, M.; Johnson, S. G. Proc. IEEE 2005, 93 (2), 216−231. (32) Schael, F.; Rubin, M. B.; Speiser, S. J. Photochem. Photobiol., A 1998, 115, 99−108. (33) Borrelli, R.; Peluso, A. Phys. Chem. Chem. Phys. 2011, 13, 4420− 4426. (34) Banerjee, S.; Baiardi, A.; Bloino, J.; Barone, V. J. Chem. Theory Comput. 2016, 12, 774−786. (35) Maeda, H.; Maeda, T.; Mizuno, K. Molecules 2012, 17, 5108− 5125. (36) Caprasecca, S.; Curutchet, C.; Mennucci, B. J. Chem. Theory Comput. 2012, 8, 4462−4473. (37) Kilsa, K.; Kajanus, J.; Martensson, J.; Albinsson, B. J. Phys. Chem. B 1999, 103, 7329−7339. (38) Bartlett, R.; Musial, M. Rev. Mod. Phys. 2007, 79, 291−352. (39) Dunning, T. H., Jr J. Chem. Phys. 1989, 90, 1007−1023. (40) Liu, T.; Callis, P. R.; Hesp, B. H.; de Groot, M.; Buma, W. J.; Broos, J. J. Am. Chem. Soc. 2005, 127, 4104−4113. (41) Ryabinkin, I. G.; Nagesh, J.; Izmaylov, A. F. J. Phys. Chem. Lett. 2015, 6, 4200−4203. (42) Gherib, R.; Ryabinkin, I. G.; Izmaylov, A. F. J. Chem. Theory Comput. 2015, 11, 1375−1382. (43) Kleinschmidt, M.; Tatchen, J.; Marian, C. M. J. Chem. Phys. 2006, 124, 124101.

2365

DOI: 10.1021/acs.jctc.6b00157 J. Chem. Theory Comput. 2016, 12, 2357−2365