Predicting Phosphorescence Rates of Light Organic Molecules Using Time-Dependent Density Functional Theory and the Path Integral Approach to Dynamics.

In this work, we present a general method for predicting phosphorescence rates and spectra for molecules using time-dependent density functional theory (TD-DFT) and a path integral approach for the dynamics that relies on the harmonic oscillator approximation for the nuclear movement. We first discuss the theory involved in including spin-orbit coupling (SOC) among singlet and triplet excited states and then how to compute the corrected transition dipole moments and phosphorescence rates. We investigate the dependence of these rates on some TD-DFT parameters, such as the nature of the functional, the number of roots, and the Tamm-Dancoff approximation. After that, we evaluate the effect of different SOC integral schemes and show that our best method is applicable to a large number of systems with different excited state characters.

Introduction

Phosphorescent materials have been a major focus of research during the past few years, with applications in organic light-emitting diodes (OLEDs),[1−6] light-emitting electrochemical cells,[7,8] photovoltaic cells,[9,10] chemical sensors,[11−14] and bioimaging.[15−19] Because most molecular materials presenting this kind of long-lived emission contain heavy metals to increase spin–orbit coupling (SOC), the interest in purely organic molecules presenting triplet emission has also grown recently,[20,21] for economical and ecological reasons.

As they exhibit weak SOC, light organic molecules usually have rather slow phosphorescence rates, between 10–3 and 103 s–1, and nonradiative pathways prevail at room temperature.[22,23] Usually, these compounds become phosphorescent only at about 77 K, in organic glasses or as pure materials, when the intersystem crossing back to the ground state slows down, and the emission quantum yield then gets higher. Yet, some molecules do present high phosphorescence rates, even at room temperature, and these are currently under study for various applications.[24−26]

The theoretical prediction of these rates is rather challenging since it must account for relativity to allow for singlet and triplet mixing. The rate for a radiative transition between an initial state Ψi and a final state Ψf can be calculated from Fermi’s Golden Rule as[27]where ω is the frequency of the photon, n is the refractive index of the medium, ℏ is the planck constant divided by 2π, c is the speed of light, μ̂ is the dipole moment operator, and E is the energy of a given state. The rate thus depends on the transition dipole matrix element, and because nonrelativistic operators do not couple spin eigenstates, ⟨ΨT|μ̂|ΨS⟩ = 0 after spin integration and there can only be transitions between states of same multiplicity.

However, the rate is, in practice, nonzero due to the relativistic mixing of singlets and triplets. In nonrelativistic quantum mechanics, the SOC operator (ĤSOC) coupling states with different multiplicity can be added to the nonrelativistic Hamiltonian in an ad hoc manner. The SOC effects are usually accounted for in the literature by using first-order perturbation theory, with the perturbed triplet |ΨT1SOC⟩ given by[28]so that phosphorescence rates can be obtained from eq under the Franck–Condon assumption that the transition dipole moment is independent of nuclear coordinates. As there are three spin sublevels for the triplets, the final observed rate kPobs is calculated considering the Boltzmann population of each state[29]where ΔE12 is the energy difference between states 1 and 2 with the three sublevels labeled from 1 to 3 and k1, k2, and k3 correspond to the individual rate from each level.

This simple approach works in some cases for molecules with heavy atoms, where the SOC is large enough to induce a strong mixing and a large transition dipole moment,[4,29,30] but it fails when these matrix elements are small and vibronic coupling must be accounted for.[28] Also, the use of eq neglects the mixture of excited triplets with the ground singlet and the effects of triplet–triplet coupling.

In this work, we present a combination of methods to overcome these issues which is both fast and applicable to large systems. We suggest using Quasi-Degenerate Perturbation Theory (QDPT)[31] to calculate the mixing between singlets and triplets obtained from time-dependent density functional theory (TD-DFT),[32,33] with or without the Tamm–Dancoff approximation (TDA),[34] to get a more complete picture of the SOC mixing. The Golden Rule rate equation is then solved in the time domain by using the path integral solution for the multidimensional harmonic oscillator. In contrast to full wave packet propagation schemes like multiconfiguration time-dependent Hartree (MCTDH),[35,36] the time evolution of a correlation function can be solved analytically if the harmonic approximation holds. In particular, we propose to use our recently published implementation of the path integral approach[37] (often called time-dependent or semiclassical approach[38−42]) to predict the rate from eq , since it accounts for temperature effects, vibronic coupling, and the Duschinsky rotation between modes in an efficient way.

We first analyze the effect of some parameters from TD-DFT and the SOC integrals on a small subset of molecules and then show that this method can be applied to a larger number of molecules with different excited state characters, including solvent effects.

Theory

SOC Matrix Elements and Coupled States

In QDPT, the mixed ground and excited states can be obtained by diagonalizing the approximate relativistic Hamiltonian,in the basis of the unmixed states.[31] Here, ĤBO is the Born–Oppenheimer Hamiltonian and ĤSOC can be written in a second quantized form as[43]

The second quantized spin operators read as follows[44]where a and b are the creation operators for an alpha electron and a beta electron, respectively, and p and q denote arbitrary molecular orbitals. In eq , the SOC integrals in the spherical basis are given by

These SOC matrix elements can be calculated using the Wigner–Eckart theorem for spherical tensor operators[43,44] of rank 1,The Clebsch–Gordan coefficients (CGC) in the brackets only couple states for which |S′ – S| is at most 1 and M′ + m = M. The reduced matrix element can be calculated aswhere |l| is at most one, depending on the states coupled, and C is a simple constant depending on S and determined as the inverse CGC for S = M and S′ = M′ = S – l. Thus, it is only necessary to compute the SOC integrals z and the transition spin density matrices on the right-hand side of the above equation.

In TD-DFT, assuming a closed shell reference, the spin adapted singlet and triplet state vectors can be calculated as single excitationswith X̂ being a general notation for the normalized versions of eqs , 7, and 8 for triplet excitations as well as for the singlet excitation operator obtained fromafter normalization so that ⟨0|X̂†X̂|0⟩ = 1. In effect, the operators in eqs , 7, and 8 are multiplied by , while the one in eq is divided by . Here and in the following, the labels i, j, ... refer to occupied molecular orbitals and the labels a, b, ..., to virtual molecular orbitals, as opposed to the general labels p, q, ..., which may refer to either of these. The symbol x in eq is a generic notation for the amplitude corresponding to the excitation operator X̂. For singlets, we will use the notation x = s while for triplets the convention x = t will be preferred. Technically, the triplet equations are solved only once in TD-DFT/TDA since the three sublevels yield identical results so long as the Hamiltonian does not affect the spin functions. This result is then used to compute the reduced matrix elements from eq . Since only singlet and triplet states will be considered in the following, a new notation is introduced in which |0⟩ ≡ |Ψ00,0⟩ denotes the ground state and |S⟩ ≡ |Ψ0,0⟩ and |T⟩ ≡ |Ψ1,⟩ denote the singlets and triplets, with s and t being the corresponding excitation amplitudes. Since the SOC operator is Hermitian, it is enough to compute the matrix elements for which S′ ≤ S. Then, the surviving transition densities can be written asThe delta symbols in the above equations merely indicate which contributions survive for a specific choice of the generic labels p and q. The only necessary coefficients C in eq are then C+1 = 1 and . In full TD-DFT, the SOC matrix elements are computed from the renormalized |X + Y⟩ vectors (such that ⟨X + Y|X + Y⟩ = 1), which are not necessarily normalized if hybrid functionals are used.[45] Note that it is also possible to handle singly excited open shell systems in a spin adapted manner, as discussed elsewhere.[31]

After the matrix for ĤREL is built and diagonalized, we obtain the mixed singlet–triplet SOC states:with the complex eigenvector coefficients C for singlets and D for triplets. If solvent effects are to be included, one can use the linear response conductor-like polarizable continuum model (LR-CPCM)[46,47] to approximately correct the energies of the TD-DFT states before constructing the ĤREL matrix.

Finally, the transition dipole matrix elements can be computed. As the ground state only couples through SOC with the triplets, its SOC-corrected wave function is given byand the transition dipoles from ground to excited states can be obtained fromwhere μ denotes the dipole intergrals, while the SOC one body reduced density matrix Γ is given byThis in turn contains contributions from unmixed quantities, including the ground state densitythe ground state to excited singlet transition densityand the triplet–triplet transition densitykeeping in mind that only transition densities between triplet states of the same M value will survive after spin integration and that these yield identical results independent of M.

Path Integral Approach for the Rates

In order to compute the rates of phosphorescence, we solve eq for the first three mixed (triplet dominated) excited states, according to our previously described method.[37] This is based on the idea first explored by Lin[41,48] in the 1970s that one can use the Fourier Transform (FT) representation of the delta functionto switch the problem from the energy to the time domain. The formulation was explored further by Tannor and Heller in the 80s, see refs (36) and (42) and references therein, and was recently revisited by Tatchen and Pollak[38] and began to be applied in different contexts.[39,40] Using this formulation, it is possible to calculate kP from the FT of a correlation function χ(t) that is computed from the path integral of the multidimensional harmonic oscillator:with α being a collection of parameters. As we have shown, this formulation can be implemented very efficiently and yields as a result both the rate and the emission spectrum.[37]

The correlation function itself depends on the electronic transition dipole μe, which in the present case is calculated from eq and its derivatives,as well as on the functionswhich contain information on the vibrational coupling at the Franck–Condon (FC) and Herzberg–Teller (HT) levels. Here Ĥ is the vibrational Hamiltonian, Qk is a normal coordinate, and the bar denotes the coordinates of the final state. The derivatives in eq are calculated numerically, while the functions ρ are evaluated as discussed in our earlier work. These functions depend on the Dushinsky rotation and displacement matrices that connect initial and final normal coordinates, and they also account for temperature dependence which is introduced by summing over all initial states weighed by a Boltzmann factor.

Using this method to solve eq together with the SOC-corrected transition dipole moment matrix elements, we can now predict phosphorescence rates (and intrinsic lifetimes) including the solvent, Duschinsky, and the Herzberg–Teller effects.

Computational Details

All calculations were carried out using the development verison of the ORCA software package.[49] Structures were optimized using the B3LYP,[50,51] BP86,[52] or WB97X-D3[53,54] functionals and the def2-TZVP(-F) basis set.[55] For ground states, the optimization followed the restricted Kohn–Sham (KS) DFT process, while for excited state singlets, TD-DFT optimizations were carried out using the restricted KS determinant as reference. For triplet states, unrestricted KS ground state calculations were performed rather than computing the triplet excited states from TD-DFT with a restricted KS reference. In order to accelerate the computation of two electron integrals, the resolution of identity approximation was used for the Coulomb part (RIJ) and the chain of spheres algorithm for the exchange part (COSX), with the corresponding auxiliary basis and grid settings.[56,57] The DFT grid was set to GRID5, and the COSX grid was GRIDX5, with all the other parameters chosen as default unless otherwise stated. When LR-CPCM was used, the solvent of choice was ethanol to be consistent with the reference data.[58] Unless otherwise mentioned, the spin–orbit integrals were calculated using the RI-SOMF(1X)[59] approximation.

For the excited states, TD-DFT with or without TDA was employed.[60] The convergence criteria for both the SCF and geometry optimizations was set to TIGHT, with a convergence threshold of 10–8 Hartree in the former and 10–6 Hartree in the latter case. No optimized structure presented negative frequencies. It was assumed that the geometries and Hessians of the triplet spin sublevels were the same and, to calculate the rates, eq was used although no significant zero-field splitting was encountered. During the calculation of SOC integrals using mean-field approaches, the KS ground state density was used when necessary.

The individual rates were calculated using the ORCA_ESD module recently developed, using default settings and Duschinsky mixing, including our automatic selection of parameters for the FT step and normal modes described in terms of internal coordinates,[37] but with the temperature set to 77 K according to the reference data.[58] For calculating the numerical derivatives, it had to be ensured that the phases of the quantities involved correspond to each other in each single point calculation on the perturbed molecular geometries. In order to fix the phases of the MOs and the excited state eigenvectors, the largest element of each vector was set to be always positive. Similarly, the ĤSOC eigenvectors were multiplied by a phase such that the element with the largest absolute value was always real and positive. After each displacement, the vectors were adjusted in such a way that the same element remained positive and, in the SOC case, real and thus comparable to the reference.

To further accelerate the numerical calculation of the transition dipole moment derivatives and avoid using central differences in all displacements, we now calculate the difference between the norm of the matrix element after one step and the reference norm. If it is smaller than 5δq|μ⃗ref|, where δq is the step size (0.01 by default), the derivative is taken to be less relevant and it is calculated from a simple difference. Our testing shows that the results are virtually the same in all cases and one can avoid most of the necessary single point calculations.

Results and Discussion

For initial testing, a set of well studied molecules were chosen, which have different excited state characters, from purely π–π* to n–π* states. This set consists of anthracene (ANT), naphtalene (NAP), phenanthrene (PHE), carbazole (CAR), benzofuran (BFUR), quinoxaline (QNX), biacetyl (BIA), benzophenone (BZP), and anthrone (ATQ), all shown in Figure .

Figure 1

Set of molecules chosen for testing the proposed method.

These molecules present experimental rates varying from 0.01 to 340 s–1 [58] and are thus suitable for evaluating the range of applicability of our method. However, since the range is large and so are the errors in the predictions, we found that the ratio between the experimental and the theoretical rate was the most suitable quantity for comparison. Here, we define this ratio as R = max(ktheo,kexp)/min(ktheo,kexp) and will use the average of this ratio within a test set (R̅) as a measure of accuracy. This measure is equally sensitive to underestimating and overestimating the rate constant; e.g., the estimates 0.01 and 100 for the value 1 are both 2 orders of magnitude off, and hence, R = 100 for both. Note that the unsigned relative errors would be 0.99 and 99, respectively.

Dependence on TD-DFT Parameters

All geometries for the S0 and T1 states of the test set were optimized using BP86, B3LYP, and WB97X-D3 functionals, and its rates were calculated using TD-DFT excited states with LR-CPCM energy corrections for the solvent. At first, the number of TD-DFT roots were varied from 10 to 25, in order to evaluate the effect of both the nature of the functionals and the number of roots necessary to converge the results. It was found that in all cases they converged after coupling the 20 lowest energy excited singlets and the 20 lowest energy triplets, so we decided to include the 25 lowest energy roots for each multiplicity as the default to ensure convergence. Thus, the final SOC matrix involves 101 states, including the ground state, 25 singlet states, and 75 triplet states (25 for each sublevel).

When comparing the effect of different functionals, we found that B3LYP was the one with the smallest error, with R̅ = 1.98. For BP86, the error was determined to be R̅ = 5.69 while R̅ = 6.97 was found for the range-corrected WB97X-D3. In Figure , a comparison between the experimental and the calculated rates is presented, indicating the amount of the Herzberg–Teller effect for each molecule.

Figure 2

Comparison between different functionals, coupling the first 25 roots obtained from TD-DFT. The values under the labels are the relative Herzberg–Teller contributions for B3LYP using the RI-SOMF approach to obtain the SOC integrals.

These rates are typically obtained from the experimental lifetimes of the triplets (τT) and the quantum yields for phosphorescence and intersystem crossing (ΦP and ΦISC) as[22,58]There is a significant amount of error in the determination of all three, making it difficult to find even a well-defined experimental value. It is clear from Figure that, on average, B3LYP gives the best predictions, although all three functionals are well within the acceptable error range, with BP86 overestimating and WB97X-D3 underestimating the results. These rates are mostly due to the Herzberg–Teller (vibronic coupling) effect and would not be obtained if only the FC approximation was used. Temperature effects are also accounted for without any approximation here.

The effect of TDA and the solvent were also considered using B3LYP as the functional of choice; see Table . As can be seen, the solvent effect was not important for these molecules, since there were almost no charge-transfer states involved and the CPCM correction was small. On the other hand, the effect of TDA was more dramatic. This was due in particular to large deviations for CAR and QNX with RCAR = 13.7 and RQNX = 17.5. Looking closely at these results, it is possible to determine that the differences are caused by large derivatives of transition dipole moments over certain modes that are not encountered when using full TD-DFT. It is known that TDA yields triplet energies closer to the experiment,[61] while full TD-DFT is better for properties such as transition dipole moments[62,63] and it is the latter which makes a difference here. If one considers only the other molecules in this set, where there was no such error, the TDA approximation is actually the best, with an R̅ = 1.65. Even so, full TD-DFT seems to be more robust in general and will be used from here on.

Table 1

Average Errors Obtained from Full TD-DFT, TDA, and Gas Phase Calculations for the Test Set Shown in Figure a

	R̅	MSRE	MURE
TD-DFT	1.98	68.7	91.2
vacuum	2.15	71.9	103
TDA	4.76	362	374

Mean signed relative errors (MSRE) and mean unsigned relative errors (MURE) are given as percentages.

Effect of Different Schemes for SOC Integrals

Since there are different approaches to compute the SOC integrals, we also investigated the effect of these on the overall rates. The SOC operator ĤSOC can be included in nonrelativistic quantum mechanics using the Breit–Pauli (BP) approximation[59]where the SOC operator is composed of one electron terms (ĤSOC(1)) and two electron terms (ĤSOC(2)). In constructing a second quantized representation for these operators, it should be kept in mind that we are aiming for a one body description of the form shown in eq . Thus, the matrix elements z(− consist of a sum, z(− = h(− + g(−, of a true one electron term h(− and an effective one body term g(− arising from ĤSOC(2) under the mean field approximation. The explicit forms of these matrix elements arewhere the complex integrals are labeled as (1*1|g(1,2)|2*2) and h(− and g(− can be obtained following the recipe in eq from the components of the Cartesian vector operatorsHere, α is the fine structure constant, 1 and 2 denote electronic spatial coordinates, r1 = |r1| and r12 = |r12| denote the distance between electron 1 and nucleus A and between electron 1 and 2, respectively. The angular momentum operators l are given as l1 = r1 × p1 and l12 = r12 × p1, where p1 is the momentum operator belonging to electron 1. As we will treat SOC as a one electron property, the question remains in what precise form the two electron effects can be treated.

The simplest approximation for the BP SOC Hamiltonian is to use the one electron part and account for the other electrons through a shielding on the effective nuclear charge (i.e., replace Z with some effective charges Zeff). Here, we use the charges from Koseki et al.[64] and will refer to this approach as ZEFF. It is also possible to include the two electron terms approximately, using mean field approaches as in Hartree–Fock theory. We use the efficient mean field approach called RI-SOMF[59] that makes use of the RI scheme to accelerate Coulomb integrals and another mean field approach named AMFI-A,[43] by Schimmelpfennig, that uses precomputed atomic densities and only accounts for one-center integrals. A fourth option is to compute the BP SOC integrals similarly to the SOMF but including local DFT exchange and correlation, here named VEFF-2X,[59] which we will simply refer to as VEFF in this study. The results for the test set are presented in Figure .

Figure 3

Comparison between different approaches to compute the BP SOC integrals and their effect on the predicted phosphorescence rates using B3LYP.

As it can be seen directly from Figure , the RI-SOMF and VEFF approaches as well as the simplest ZEFF approach yield results which agree well with each other and with experiment. The comparatively mild variation in the rates predicted by these three methods can be attributed to differences in the quality of the SOC matrix elements and their derivatives. For PHE and NAP, VEFF and RISOMF differ quite substantially, both as far as the transition dipoles and their derivatives are concerned. However, both methods yield results in the same order of magnitude. The relatively large deviation of the AMFI-A values from the other results is connected to the one-center approach, which should be avoided for the present purposes. In fact, if the RI-SOMF calculations are carried out using the one-center approximation, the results are very close to the AMFI-A values.

A well-balanced method should produce both accurate SOC matrix elements and reliable SOC derivatives, since the quality of the latter has a significant impact on the calculated phosphorescence rates. It seems that ZEFF, although only containing one electron terms, produces quite accurate rates, being even slightly more accurate than RI-SOMF with an R̅ZEFF = 1.56. Even so, the mean field RI-SOMF approach should be more reliable if heavier atoms are also present and a better treatment of the two electron terms becomes necessary. In this sense, it is more general and it does not depend on the choice of the exchange-correlation functional or other empirical parameters either. We thus recommend RI-SOMF as the default approach, while retaining VEFF as a possibility and ZEFF as a cheaper alternative for weak coupling.

Results for a Larger Test Set

Finally, having defined a suitable set of default choices for our method, B3LYP functional, full TD-TDF with LR-CPCM corrections, and RI-SOMF SOC integrals, we investigate the accuracy of this method on a larger set. In addition to the ones portrayed in Figure , we add the well-known molecules presented in Figure , which are pyrene (PYR), biphenyl (BIP), benzonitrile (BNT), adenine (ADE), xanthone (XAN), and anthraquinone (ATQ).[58]

Figure 4

Molecules added to those in Figure to form the final evaluation set.

From these 15 molecules, 5 have only π–π* excited states with very weak SOC, 5 have well-defined n–π* states that give rise to larger SOC matrix elements, and 5 have some sort of mixed states. It can be seen in Figure that these differences give rise to phosphorescence rates that are up to 4 orders apart, e.g., between anthracene (ATQ) and anthraquinone (ATQ), as expected from the usual theoretical background of molecular photophysics.[22,23] Again, as in the previous case, most of these rates are due to vibronic coupling between the triplet and the ground state and would not be correct under the simpler FC assumption.

Figure 5

Predicted phosphorescence rates for all molecules presented in Figures and 4 using our best methods. Due to lack of complete data, for pyrene, we assumed that ΦISC = 1 – ΦF and, for adenine, ΦP = ΦISC when using eq (here, ΦF is the fluorescence quantum yield).

As explained in our previous work regarding fluorescence,[37] the emission rates can be also related to the experimental spectrum, and in Figure , we show the comparison between the experimental and theoretical phosphorescence spectra of biacetyl. As can be seen, there is a good correspondence in both the intensity and the position of the bands. In this case, about 50% of the intensity is due to the HT effect and the associated theoretical rate is in excellent agreement with the experiment.

Figure 6

Predicted phosphorescence spectrum and rate for biacetyl in ethanol.[65] The theoretical curve was blue-shifted by 2880 cm–1 to match the experimental 0–0 transition energy.

Conclusions

In this work, we have provided a derivation of the spin–orbit coupling between TD-DFT ground and excited states and the application of the corrected transition dipole moment matrix elements for the computation of phosphorescence rates and spectra. In order to solve the relevant rate equations, we use our previously developed approach for fluorescence, based on the path integral solution of the multidimensional harmonic oscillator model, and investigated the effects of different parameters on the predicted rates.

From our results, we conclude that the functional B3LYP yielded the best predictions, when using the full TD-DFT equations, at least when compared with BP86 and WB97X. The results obtained using the TDA approximation were also very good but had a large deviation for two test cases and are thus less robust. For simple organic molecules such as those explored here, spin–orbit coupling is weak and the Franck–Condon approximation is not sufficient to give qualitatively correct predictions. To achieve this, the Herzberg–Teller effects must be accounted for; otherwise, the results deviate enormously from the experiment. In terms of the spin–orbit integrals, we conclude that the effective nuclear charges as developed by Koseki et al.[64] and the mean-field integrals with or without local DFT contributions[59] compare well to the experiment, at least for light organic molecules. We recommend the RI-SOMF method in particular and advise against using the one-center approximation for calculating phosphorescence rates.

We finally show that, under the use of these optimal parameters, phosphorescence rates can be predicted for a series of diverse molecules that span almost 4 orders of magnitude in terms of the phosphorescence rates. Using these results, emission spectra can also be predicted in good agreement with the experiment.