The internal conversion (IC) process from S_{1} to S_{0} and the intersystem crossing (ISC) transition from T_{1} to S_{0} are two essential processes in functional molecular material design. Despite their importance, it is currently impossible to evaluate the rate of these processes for a large set of molecules and, therefore, perform high-throughput virtual screening in large-scale data to gain more physical insight. In this work, we explore possible approaches to accelerate the calculations of IC and ISC rates based on a systematic reduction of the number of modes included in the computation and the study of the importance of the different parameters and the influence of their accuracy on the final result. The results reproduce the experimental trends with systematic errors that are ultimately due to the approximations of the theory. We noted that plausible results for ISC in planar molecules are only obtained by including the effect of Hertzberg–Teller coupling. Our method establishes the feasibility and expected accuracy of the computation of nonradiative rates in the virtual screening of molecular materials.

## INTRODUCTION

The nonradiative internal conversion (IC) process from the first singlet excited state S_{1} to the ground state S_{0} and the intersystem crossing (ISC) process from the first triplet excited state T_{1} to S_{0} are two of the most widely studied processes in organic photochemistry with a large set of experimental data available, which can be used to validate the theory and understand its limitation.^{1–6} A molecule in the S_{1} state can generally be de-excited via fluorescence (S_{1} → S_{0} + *h*ν), undergo ISC processes (e.g., S_{1} → T_{1}) or the IC process (S_{1} → S_{0}). The molecule in the T_{1} state can also relax to the ground state via the ISC process (T_{1} → S_{0}). The rates of these processes play crucial roles in understanding and designing functional molecular materials with desirable optical properties. For example, thermally activated delayed fluorescence (TADF) materials with delayed fluorescence can be designed by decreasing the IC rate.^{7} Strong fluorescence quantum yield of organic light-emitting diode (OLED) requires a small singlet-triplet gap and a weak IC rate.^{8} Feng *et al.*^{9} raised a strategy to achieve solution-phase organic room-temperature phosphorescence by suppressing the ISC between the T_{1} state and S_{0} states through spatial separation. Nonradiative relaxation channels are detrimental for all devices converting solar radiation to chemical fuels (photocatalysis^{10}) or electrical power (photovoltaic^{11}).

A promising way to design and develop functional organic materials is high-throughput virtual screening (HTVS),^{12–15} where properties of interest are computed for large sets (currently 10^{4}–10^{5}) of candidate molecules. For example, thousands of novel TADF OLED molecules across the visible spectrum were successfully identified through a large-scale computer-driven search.^{16,17} Hundreds of singlet fission materials outside the previous design principles were discovered by screening a database of experimentally known molecules.^{18} New insights into design rules were gained from the statistical study and machine learning of organic acceptors for photovoltaic devices.^{19–21} An important observation that can be made on the basis of recent literature is that, while nonradiative transitions play essential roles in most of the relevant applications, they are generally excluded from the screening process. There are clearly many theoretical works^{22–35} in this area spanning many decades, but they tend to focus on fewer examples or families of the related compound and have not yet been employed to determine the photophysical rate of large datasets. Different considerations become important in the applications to HTVS, such as the relative costs of the various components of the calculation and the sensitivity and robustness of the results to the detail of the electronic structure calculations (density functional, basis set, number of modes included). A study of these aspects is a prerequisite for incorporating the computation of nonradiative rates in HTVS.

An approach that is potentially amenable to be used in conjunction with HTVS is that based on time-dependent perturbation theory (Fermi’s golden rules), where the transitions are seen as population transfer between initial and final vibronic states, assuming that the potential energy surfaces of the relevant electronic states can be considered harmonic. This framework is routinely used to interpret phenomena like TADF^{36} or quantitatively predict the light-emitting efficiency.^{37} There are two general approaches for performing the calculation of the Franck–Condon terms. One is the time-independent sum-over-states method, which considers all possible transitions between the initial and final vibronic state.^{38–41} This approach is straightforward and easy to implement. The other is the time-dependent method, where the delta functions in Fermi’s golden rules can be treated as the Fourier transform of autocorrelation functions.^{30,42–44} Unlike the original time-independent approach, the correlation function approach does not suffer from the problem of the huge number of transitions in large-sized systems. However, general and robust codes may need to avoid numerical instability of time integral. In any case, this component of the methodology is not the rate-determining step of the calculation. It is also typically used for a few systems at the time, where it is possible to control all the model parameters. To expand similar methods to larger datasets, one needs to automatize the construction of the vibronic model and the computation of nonadiabatic coupling for IC or spin–orbit coupling for the ISC process, respectively, possibly studying the effect of reducing the accuracy of the quantum chemical calculations on the resulting rate. As usual, when the dataset of interest increases, it becomes possible to introduce an empirical correction to the systematic errors of the methodology.^{18} In the case of IC and ISC, a well-known limitation is the anharmonicity, which can enhance the rates significantly, especially for systems with a large energy gap between the initial and final electronic states (e.g., above 20 000 cm^{−1}).^{33,45} Finally, there will be aspects that are likely outside the scope of HTVS, such as relaxation mechanism involving large amplitude motions or conical intersection.^{46–50} Some of these channels, like isomerization, are suppressed in larger rigid molecules in the condensed phase but are still possible. In the context of materials discovery, they should be studied with conventional approaches after the range of candidates to be included is reduced by larger-scale HTVS.

In this work, the IC and ISC rates are calculated using Fermi’s golden rule, coupled with electronic structure calculations of realistic molecules and the harmonic assumption. Several approximations are introduced to systematically reduce the number of vibrational modes (and vibronic states) to make it possible to evaluate nonradiative rates rapidly for a potentially large number of molecules. This work focuses on validating the methodology and correcting systematic errors by comparison to a dataset of a few tens of experimentally determined nonradiative rates.

## METHODOLOGY AND COMPUTATIONAL DETAILS

### Nonradiative rate

The theoretical background is outlined briefly to provide the notation required to describe the approximations and the results.^{30,32,34} The transition rate from one initial state to several final states can be calculated using Fermi’s golden rule,^{51}

where *ℏ* is the Planck constant, $Vfi=fV\u0302i$ is the interaction between a final state $f$ and the initial state $i$, and $\delta E$ is the Dirac delta function.

For the IC process, the nuclear kinetic energy operator $T\u0302N$ is regarded as the interaction $V\u0302$ in Eq. (1). Within the harmonic approximation, it can be expressed as

where *ω*_{j} and *Q*_{j} are the frequency and the dimensionless coordinate of the normal mode *j*, respectively.

By using the vibronic state representation under the Born–Oppenheimer approximation, the interaction between a final state and the initial state can be rewritten as

where $\varphi i/fq,Qj$ is the electronic part of the initial/final vibronic state with ** q** being the coordinate of electrons and $\u220fj=13N\u22126\Lambda wji/fi/fQj$ is the vibrational part given by a product of 3

*N*− 6 (or 3

*N*− 5 for linear systems) vibrational states $\Lambda wji/fi/fQj$. For vibrational states, the harmonic approximation is applied based on the potential energy surface of the electronic state $\varphi i/f$, and $wji/f$ is the occupation number of the normal mode

*j*. For convenience, a vibronic state $\varphi aq,Qj\u220fj=13N\u22126\Lambda wjaaQj$ is labeled as $awja$.

By ignoring high-order terms of the operator $T\u0302N$, the interaction in Eq. (3) can be approximated as

In this work, we focused on the IC process from S_{1} to S_{0}. The index *a* of vibronic states $awja$ can be set as S_{1} and S_{0} for the initial and final states, respectively. It is also reasonable to assume that the vibrational wavefunctions of the initial state are all ground states at low temperature (i.e., $wji=0$for all *j*), while the occupation number set of vibration modes could be arbitrary for final states and is indicated as $wj$ (in the Table S1, we show that the activation energy is very large, and this is a reasonable approximation). The initial state and a final state are labeled as $S10$ and $S0wj$ individually, and Eq. (4) can then be expressed as

where $\u210f\omega j\varphi S0\u2202\u2202Qj\varphi S1$ is the nonadiabatic coupling *V*_{NA,j} between S_{0} and S_{1} for the normal mode *j*. Squares of vibrational integrals $\Lambda wjS0Qj\u2202\u2202Qj\Lambda 0S1Qj$ and $\Lambda wjS0Qj|\Lambda 0S1Qj$ can be computed analytically under the harmonic approximation,

Here, *y*_{j} is the Huang–Rhys (HR) factor of the normal mode *j* for the transition from S_{1} to S_{0}. It could be computed by projecting the total structure difference between S_{1} and S_{0} equilibrium geometries onto the displacement vector of the normal mode *j*.^{52,53}

For the delta function in Eq. (1), the energy gap between two states *E*_{f} − *E*_{i} can be calculated as the sum of the electronic energy gap $ES0\u2212ES1$ and the nuclear energy difference $\u2211j\u210fwj\omega j$ under the harmonic approximation. The overall expression of the IC rate can, thus, be written as

It is commonly assumed that the modes inducing the transition (with nonzero nonadiabatic couplings) and accepting the electronic energy (with nonzero HR factors) are distinct.^{54} If *V*_{NA,j} is nonzero, we therefore assume *y*_{j} = 0 for a mode *j*, and the only value of *w*_{j} contributing to the summation in Eq. (8) is *w*_{j} = 1 for which $\Lambda wjS0\u2202\u2202Qj\Lambda 0S12$ does not vanish (1/2). By discarding the zero values of the nuclear integrals in Eqs. (6) and (7) for some particular set of occupation numbers [e.g., *w*_{j} ≠ 1 with *V*_{NA,j} ≠ 0 for Eq. (6), *w*_{j} ≠ 0 with *y*_{j} ≠ 0 for Eq. (7)], the rate in Eq. (8) can be simplified as

This equation can be further simplified by replacing *ℏω*_{j} with an average *ℏω*_{NA} weighted by the nonzero nonadiabatic couplings,

with the effective nonadiabatic coupling $VNAeff$ being

and

In Eq. (10), the last term is an energy-shifted Franck–Condon weighted density (FCWD), which can be defined as

Then, Eq. (10) can be expressed compactly as a product of the effective nonadiabatic coupling and an independent Franck–Condon term, i.e.,

The ISC rate from T_{1} to S_{0} can be calculated similarly by replacing the interaction operator with spin–orbit coupling (SOC) operator $H\u0302SOC$. The final formula is

### Experimental data

To determine the validity of the methodology, we considered several experimental datasets for IC and ISC rates. Experimental energy gaps are used as the electronic energy gap terms $ES1\u2212ES0$ and $ET1\u2212ES0$ for IC and ISC rates in Eqs. (14) and (15), respectively, as we only focus on the accuracy of coupling terms and the Franck–Condon factors in this work. For the calculation of hypothetical molecules, one should clearly also compute this energy term either using modern range separated functionals^{55–57} or a calibration procedure against experimental data.^{18} Experimental data of 12 molecules (shown in Fig. 1) for the IC process are obtained from Refs. 58–60 that focused on the discussion of the energy gap law.^{25} The experimental IC rates were calculated using the experimental fluorescence lifetime and quantum yield for naphthalene, anthracene, tetracene, and coronene.^{58} The experimental IC rates for azulene, azulene-d8, guaiazulene, trimethyl-azulene, phenylbenzoxalene, dibenzoxalene, bromodibenzoxalene, and chlorodibenzoxalene were calculated by using the experimental fluorescence lifetime.^{59,60} The dataset for the ISC processes contains 17 molecules (shown in Fig. 1), including naphthalene, biphenylene, biphenyl, fluorene, anthracene, phenanthrene, pyrene, fluoranthene, tetracene, chrysene, 3,4-benzphenanthrene, 1,2-benzanthracene, *p*-terphenyl, *m*-terphenyl, 1,2,5,6-dibenzanthracene, 1,3,5-triphenyl-benzene, and hexahelicene.^{22} Their ISC rates were computed using experimental T_{1} lifetime. It is worth noticing that the existing collections of nonradiative rates from literature (e.g., in Ref. 61) proved to be not immediately useful for this work because of several inconsistencies (mismatch with cited sources, quantum yields summing to values greater than 1) that impose the need for manual verification of each entry. The current dataset can be considered a validation dataset (it spans the relevant range of excitation energy) and the results presented here can be regarded as a prerequisite for the application of a similar methodology to a larger dataset.

### Computational details

To obtain FCWD terms in Eqs. (14) and (15), HF factors are calculated by projecting the geometry difference for the corresponding processes onto the displacement vectors of the normal modes of the final electronic state (i.e., S_{0} for both cases) under curvilinear coordinates.^{53} In practice, restricted DFT, unrestricted DFT, and TDDFT calculations are carried out using the Gaussian 16 software (Revision A.03)^{62} to obtain the optimized S_{0}, T_{1}, and S_{1} molecular structures, respectively. Normal modes information is also calculated by DFT in Gaussian 16 using the S_{0} geometry. The results are presented for calculations in a vacuum with the supplementary material (Table S2) reporting a modest effect of the solvent treated within the polarizable continuum model^{63,64} on these results. After stable numerical tests (Fig. S2a), the delta function in FCWD is simulated by a normalized Gaussian broadening function with a standard deviation of 200 cm^{−1}.

For the IC process, nonadiabatic couplings are computed by projecting the nonadiabatic coupling matrix elements onto the normal modes of S_{0}. This procedure is preferable to various diabatization schemes^{65} for possible applications in HTVS, where the quality of the diabatization cannot be verified automatically. The NAC matrix elements between S_{1} and S_{0} are calculated in the Gaussian 16 package using the TDDFT approach proposed by Send and Furche,^{66} which is the most common implementation of NAC matrix element computations within the DFT framework and has been used to study internal conversions between low-lying electronic excited states in organic systems.^{32,67–71} For the ISC process, the SOC between S_{0} and T_{1} is obtained using TDDFT without the Tamm–Dancoff approximation in the ORCA package.^{72,73} The S_{0} geometry is used for both coupling terms.

Since the B3LYP functional^{74–76} provides good energies and vibronic spectral shapes,^{77,78} it is adopted as the default functional for the calculations above. The 6-31G(d)^{79} basis set is used as default for all DFT/TDDFT calculations except the SOC calculation, for which a larger basis set def2-TZVP basis set^{80} is adopted. To balance the accuracy and efficiency of our approaches for HTVS, different basis sets (3-21G,^{81} 6-31G(d) and def2-TZVP) and the alternative M06-2X functional^{82} are also used to compute parameters for IC/ISC rate calculations.

### Reduction of modes set for the FCWD

Given that the size of vibronic modes set, $wk$ in Eq. (13) increases exponentially with increasing number of vibration modes, thus, evaluating the FCWD explicitly could be too expensive if the complete set of vibrational modes is used. At the same time, the evaluation of the nonradiative rates should be independent of the size of the molecule for the intended application in HTVS. We define a reduced effective set of normal modes and HR factors in the following way. We define the reorganization energy *λ*_{j}, the contribution to the total reorganization energy due to mode *j*, in terms of the mode frequency and HR factor,

We then define a spectral density $J\omega $ introducing a finite broadening *σ* for each mode,

Figure 2(a) illustrates the relationship between the discrete values of $\omega j,\lambda j$ and the continuous function $J\omega $ for a particular molecule using *σ* = 160 cm^{−1}. [See Figs. 4(a) and 4(b) and S1 for other IC and ISC molecules.]

This function can be fitted quite reasonably for molecular systems as the sum of a limited number of Gaussian functions (e.g., four functions) of varying widths,

Figure 1(b) shows the fitting results using one to four Gaussian functions. The effective mode set with four modes reproduces an acceptable fitting curve compared to the original spectral density. Besides the residual sum of squares of the fitting (Table S3), the total effective reorganization energy *λ*^{eff} can also be computed using the parameter set $\omega k\u2032,\lambda k\u2032$ to evaluate the fitting quality (Table S4 for both IC and ISC molecules) compared to the original total reorganization energy *λ*.

Important modes that promote nonradiative processes are those with strong HR factors (significant FC overlaps) and high energy (easy to overcome the energy gap). It is, therefore, vital to include in the reduced modes the weakly coupled modes in the region of the C–H stretching (∼3000 cm^{−1}). To capture the information at high frequency, one effective mode is initialized at a high-frequency region (>2000 cm^{−1}), while the other three modes are set arbitrarily (larger than 0) as an initial guess for the fitting process.

After the fitting procedure, parameters set $\omega k\u2032,\lambda k\u2032$ would be used to obtain frequencies and HR factors ($yk\u2032=\lambda k\u2032/\u210f\omega k\u2032$) of effective modes to estimate the FCWD in Eq. (13) (see Table S4). By comparing the FCWD using different effective sets to the benchmark result of the six-mode complete set, the results for naphthalene (IC process) in Fig. S2b suggest the validity of the reduction approach (also see Fig. S3 and Table S5 for tetracene and azulene).

## RESULTS AND DISCUSSION

### Internal conversion

With the B3LYP/6-31G(d) calculation level, the IC rates of 12 molecules are computed using our approach, as displayed in Table I and Fig. 3. The results show qualitative agreement with experimental data (e.g., the correct order of rates) and the energy gap law, except for coronene, molecule **12a** in Table I. The considerable underestimation of the rate for molecules with larger energy gaps is generally attributed to the anharmonicity effect.^{33,45}

. | $ES1\u2212ES0$ (cm^{−1})
. | Expt. $logkIC/s\u22121$ . | Cal. $logkIC/s\u22121$ . | λ (cm^{−1})
. | $VNAeff$ (cm^{−1})
. | ω_{NA} (cm^{−1})
. |
---|---|---|---|---|---|---|

1a | 14 660 | 11.85 | 9.52 | 3397 | 185 | 1188 |

2a | 14 620 | 11.85 | 10.07 | 4040 | 202 | 1545 |

3a | 13 900 | 12.36 | 10.13 | 3496 | 210 | 1521 |

4a | 15 900 | 11.25 | 8.91 | 3452 | 178 | 1192 |

5a | 15 000 | 12.31 | 10.43 | 5328 | 123 | 1559 |

6a | 16 300 | 11.89 | 9.66 | 4784 | 108 | 1543 |

7a | 15 900 | 12.02 | 9.66 | 4717 | 105 | 1351 |

8a | 16 000 | 12 | 9.72 | 4710 | 109 | 1521 |

9a | 32 200 | 5.3 | −4.21 | 2413 | 27 | 1807 |

10a | 26 700 | 5.95 | −2.34 | 1803 | 28 | 1821 |

11a | 21 200 | 7.48 | 0.37 | 1364 | 32 | 1716 |

12a | 23 800 | 5.78 | −6.88 | 552 | 151 | 1348 |

. | $ES1\u2212ES0$ (cm^{−1})
. | Expt. $logkIC/s\u22121$ . | Cal. $logkIC/s\u22121$ . | λ (cm^{−1})
. | $VNAeff$ (cm^{−1})
. | ω_{NA} (cm^{−1})
. |
---|---|---|---|---|---|---|

1a | 14 660 | 11.85 | 9.52 | 3397 | 185 | 1188 |

2a | 14 620 | 11.85 | 10.07 | 4040 | 202 | 1545 |

3a | 13 900 | 12.36 | 10.13 | 3496 | 210 | 1521 |

4a | 15 900 | 11.25 | 8.91 | 3452 | 178 | 1192 |

5a | 15 000 | 12.31 | 10.43 | 5328 | 123 | 1559 |

6a | 16 300 | 11.89 | 9.66 | 4784 | 108 | 1543 |

7a | 15 900 | 12.02 | 9.66 | 4717 | 105 | 1351 |

8a | 16 000 | 12 | 9.72 | 4710 | 109 | 1521 |

9a | 32 200 | 5.3 | −4.21 | 2413 | 27 | 1807 |

10a | 26 700 | 5.95 | −2.34 | 1803 | 28 | 1821 |

11a | 21 200 | 7.48 | 0.37 | 1364 | 32 | 1716 |

12a | 23 800 | 5.78 | −6.88 | 552 | 151 | 1348 |

According to Eq. (14), the IC rate is decided by the total reorganization energy *λ* and the nonadiabatic coupling $VNAeff$, which are also listed in Table I (and Fig. S4). The most significant deviation of coronene is probably due to the computed reorganization energy, which is underestimated using DFT while, experimentally,^{83} it is closer to the other molecules. One can observe from Table I and Fig. S4 an approximate decrease of the nonadiabatic coupling with increase of the energy gap (which could be anticipated from the Hellmann–Feynman formulation of the nonadiabatic coupling^{84}) of one order of magnitude for the dataset considered. The computation of $VNAeff$ for each molecule is therefore essential although clearly less crucial than the evaluation of the HR factors. As shown in Fig. S6, the level of theory does not affect the results of the predicted rate too significantly.

A possible method to further accelerate the calculation is derived from the observation that the spectral densities of organic molecules have similar shapes across a diverse set of molecules (this is illustrated in Fig. 4 and Fig. S1). It is, therefore, possible to define a single reference effective set of frequencies and HR factors $\omega keff,ykeff$ of a typical molecule with total reorganization energy *λ*^{ref}. Then, for any arbitrary molecule, one can define the set of effective frequencies and HR factors as $\omega k=\omega keff,yk=ykeff\lambda /\lambda ref$, which only requires the calculation of the total reorganization energy *λ* of this molecule (saving the calculation of the frequency). This approximation consists in assuming that a common spectral density can capture the phenomenology of organic conjugated molecules. Figure 4(c) shows that the rates computed with this approximation are similar to the original results in Table I, with the spectral density of guaiazulene chosen randomly as the reference spectral density (the results are similar if one performs the calculation with an average spectral density among those computed). These suggest that, in the context of high-throughput virtual screening, one can avoid the computation of the individual spectral density and, therefore, the calculation of the normal modes, which computationally costs as much as the optimization of the excited state (see Fig. 5).

Despite the qualitative agreement, there is a substantial deviation from the experimental rate that increases with increase in the energy gap, which is normally attributed to the anharmonic effect not captured by the model. Empirical corrections are possible, but they would require a large experimental dataset to become more robust. For example, a constant multiplicative correction would align very well with the experimental and computed data for molecules with an energy gap below 18 000 cm^{−1}. Still, a different correction is required for molecules with a higher energy gap, suggesting that a three-parameter empirical correction is probably needed.

Finally, as mentioned in the computational detail part, different calculation levels are performed to balance the accuracy and efficiency of our approach. Figure 5 shows the computation cost for each component by using four effective modes. The cost of optimization is the highest, and the percentage of optimization is over 50%. In contrast, the average cost for calculating the FC term for the IC and ISC rates using four modes is only several seconds, and the percentage of FC calculation cost is much less than 1%. Moreover, the NAC and SOC calculation cost is generally low, less than 16%. The cost of the frequency calculation increased from around 10% to 20% as the calculation level increased. At the same time, the rate results in Fig. S6 suggest that the calculated IC rate is less sensitive to the basis set. Therefore, the 3-21G basis set could be acceptable for virtual screening approaches 6-31G(d) would be a better choice if an accuracy comparison is needed in a narrow energy gap region (e.g., molecules around 15 000 cm^{−1} in our case) since the total reorganization energy can be computed well using 6-31G(d) (Fig. S5) instead of 3-21G. Similar results are seen with the functional M06-2X with a slightly larger deviation from the experimental trend.

### Intersystem crossing

The ISC rates of 17 molecules are calculated using Eq. (15), and a summary of quantum calculations results and experimental data is provided in Table II. However, the results of the ISC rate (“cal. $logkISC/s\u22121$” in Table II) are entirely unsatisfactory and not in agreement with either the experimental rates or the energy gap law (also see Fig. S7a). As shown in Table II, SOCs significantly vary from 0.1 to 10^{−5} cm^{−1} (also see Fig. S7b), which is the origin of the poor results.

. | $ET1\u2212ES0$ (cm^{−1})
. | Expt. $logkISC/s\u22121$ . | Cal. $logkISC/s\u22121$ . | Cal. $logkISC(HT)/s\u22121$ . | λ (cm^{−1})
. | V_{SOC} (cm^{−1})^{a}
. | $VSOCeff$ (cm^{−1})^{a}
. | ω_{SOC} (cm^{−1})
. |
---|---|---|---|---|---|---|---|---|

1b | 21 300 | −0.41 | −7.40 | 1.54 | 3655 | 1.9 × 10^{−5} | 0.39 | 763 |

2b | 19 000 | −0.17 | 1.96 | ⋯ | 4952 | 6.4 × 10^{−2} | ⋯ | ⋯ |

3b | 22 900 | −0.70 | 0.77 | ⋯ | 4606 | 2.2 × 10^{−1} | ⋯ | ⋯ |

4b | 23 800 | −0.77 | −3.58 | 0.43 | 3920 | 5.9 × 10^{−3} | 0.36 | 813 |

5b | 14 900 | 1.34 | −3.43 | 3.18 | 2633 | 2.1 × 10^{−4} | 0.34 | 677 |

6b | 21 600 | −0.60 | −1.79 | 1.99 | 3993 | 7.4 × 10^{−3} | 0.35 | 796 |

7b | 16 900 | 0.30 | −3.74 | 2.31 | 2715 | 4.1 × 10^{−4} | 0.27 | 749 |

8b | 18 500 | 0.08 | −5.64 | 1.94 | 3144 | 7.8 × 10^{−5} | 0.31 | 741 |

9b | 10 300 | 2.85 | −1.32 | 4.66 | 2087 | 3.3 × 10^{−4} | 0.32 | 645 |

10b | 20 000 | −0.44 | −4.60 | 1.24 | 3196 | 5.4 × 10^{−4} | 0.30 | 734 |

11b | 20 000 | −0.43 | −0.75 | ⋯ | 2773 | 2.5 × 10^{−1} | ⋯ | ⋯ |

12b | 16 500 | 0.52 | −1.90 | 2.35 | 2625 | 3.7 × 10^{−3} | 0.33 | 690 |

13b | 20 600 | −0.44 | 1.71 | ⋯ | 4255 | 3.5 × 10^{−1} | ⋯ | ⋯ |

14b | 22 700 | −0.68 | 0.98 | ⋯ | 4433 | 3.2 × 10^{−1} | ⋯ | ⋯ |

15b | 18 300 | −0.17 | −4.76 | 1.38 | 2605 | 4.2 × 10^{−4} | 0.31 | 714 |

16b | 22 600 | −0.77 | 0.63 | ⋯ | 4408 | 2.2 × 10^{−1} | ⋯ | ⋯ |

17b | 19 000 | −0.37 | −1.33 | 1.54 | 3542 | 1.9 × 10^{−2} | ⋯ | ⋯ |

. | $ET1\u2212ES0$ (cm^{−1})
. | Expt. $logkISC/s\u22121$ . | Cal. $logkISC/s\u22121$ . | Cal. $logkISC(HT)/s\u22121$ . | λ (cm^{−1})
. | V_{SOC} (cm^{−1})^{a}
. | $VSOCeff$ (cm^{−1})^{a}
. | ω_{SOC} (cm^{−1})
. |
---|---|---|---|---|---|---|---|---|

1b | 21 300 | −0.41 | −7.40 | 1.54 | 3655 | 1.9 × 10^{−5} | 0.39 | 763 |

2b | 19 000 | −0.17 | 1.96 | ⋯ | 4952 | 6.4 × 10^{−2} | ⋯ | ⋯ |

3b | 22 900 | −0.70 | 0.77 | ⋯ | 4606 | 2.2 × 10^{−1} | ⋯ | ⋯ |

4b | 23 800 | −0.77 | −3.58 | 0.43 | 3920 | 5.9 × 10^{−3} | 0.36 | 813 |

5b | 14 900 | 1.34 | −3.43 | 3.18 | 2633 | 2.1 × 10^{−4} | 0.34 | 677 |

6b | 21 600 | −0.60 | −1.79 | 1.99 | 3993 | 7.4 × 10^{−3} | 0.35 | 796 |

7b | 16 900 | 0.30 | −3.74 | 2.31 | 2715 | 4.1 × 10^{−4} | 0.27 | 749 |

8b | 18 500 | 0.08 | −5.64 | 1.94 | 3144 | 7.8 × 10^{−5} | 0.31 | 741 |

9b | 10 300 | 2.85 | −1.32 | 4.66 | 2087 | 3.3 × 10^{−4} | 0.32 | 645 |

10b | 20 000 | −0.44 | −4.60 | 1.24 | 3196 | 5.4 × 10^{−4} | 0.30 | 734 |

11b | 20 000 | −0.43 | −0.75 | ⋯ | 2773 | 2.5 × 10^{−1} | ⋯ | ⋯ |

12b | 16 500 | 0.52 | −1.90 | 2.35 | 2625 | 3.7 × 10^{−3} | 0.33 | 690 |

13b | 20 600 | −0.44 | 1.71 | ⋯ | 4255 | 3.5 × 10^{−1} | ⋯ | ⋯ |

14b | 22 700 | −0.68 | 0.98 | ⋯ | 4433 | 3.2 × 10^{−1} | ⋯ | ⋯ |

15b | 18 300 | −0.17 | −4.76 | 1.38 | 2605 | 4.2 × 10^{−4} | 0.31 | 714 |

16b | 22 600 | −0.77 | 0.63 | ⋯ | 4408 | 2.2 × 10^{−1} | ⋯ | ⋯ |

17b | 19 000 | −0.37 | −1.33 | 1.54 | 3542 | 1.9 × 10^{−2} | ⋯ | ⋯ |

^{a}

*V*_{SOC} and $VSOCeff$ are computed at the B3LYP/def2-TZVP and B3LYP/3-21G calculation level, respectively.

It is easy to see that all molecules with nearly vanishing SOC are planar hydrocarbons molecules (e.g., naphthalene, pyrene) due to the predominant *π* − *π** excitation character of the T_{1} states for these molecules. In other words, near-zero SOC is characteristic of planar organic molecules. Instead, high-order spin–orbit interaction coupled with nuclear motion (i.e., Herzberg–Teller coupling and nonadiabatic coupling) should be included for ISC rate calculations of these molecules.^{85–87}

Here, we propose an approximated method for SOC correction to improve our approach for predicting the ISC rate for planar molecules with negligible SOC, starting by introducing the geometry (*Q*) dependence of the SOC,

Considering out-of-plane vibration modes in the original vibrational set consisting of all normal modes that introduce *σ* − *π** excitation to T_{1} states, the SOC *V*_{SOC}(*Q*) can be expressed as

Since out-of-plane modes have zero HR factor, the final equation of the ISC rate$kISC(HT)$, corrected by Herzberg–Teller (HT) coupling, has a similar approximation to the one for IC rate calculation [Eq. (10)],

Here, the effective SOC $VSOCeff$ is written as

and $\u210f\omega \u0304SOC$ is the SOC-weighted energy shift,

After the correction for planar molecules with negligible SOC, the results show good improvement for rate predictions and follow the energy gap law well, as shown in Fig. 6. More importantly, the rates in Fig. 6 are corrected by including the effective SOC $VSOCeff$ calculated using a much smaller basis set (3-21G), which suggests that small basis sets are acceptable to capture $VSOCeff$ (also see Table S6). That is, this correction for planar molecules can still be acceptable for HTVS. Furthermore, all planar molecules in our calculation have similar $VSOCeff$ (∼0.32 cm^{−1}) and $\u210f\omega \u0304SOC$ (∼750 cm^{−1}), as shown in Table II and Fig. S8. This may suggest a further approximation in the ISC evaluation for planar organic molecules, consisting of using a fixed average $VSOCeff$ and $\u210f\omega \u0304SOC$ in Eq. (21).

In the case of molecules for which we do not perform the correction (nonplanar molecules or molecules with considerable SOC), the ISC rates via Eq. (15) show consistency with the trend of corrected rates, except for two nonplanar acenes (**11b** and **17b**). The underestimation of these two molecules would result from the small dihedral angles among rings (e.g., below 10°), making high-order vibronic coupling terms still contribute to the ISC process. We also note that, for four nonplanar molecules (biphenyl, *m*/*p*-terphenyl, and triphenylbenzene), low-frequency modes associated with the planarization of the molecule contribute significantly to geometry change during the molecular relaxation of the ISC process (Fig. S9). Within our framework, the vibronic manifold of this molecule will contain a low-frequency effective mode (below 200 cm^{−1}) with a large HR factor. This slows down the FCWD calculation without contributing much to the rate while suffering from large anharmonic effects. We have verified (Fig. S10) that this mode can be removed with marginal effect on the final results and, for practical application, one can set a lower boundary for the frequency of the normal modes to be included in the fitting process.

## CONCLUSION

In this work, a rapid calculation approach for IC and ISC rates is proposed to make the evaluation of such rates possible in conjunction with high-throughput virtual screening and contribute to organic materials discovery. The method is based on the standard time-dependent perturbation theory (Fermi’s golden rule), where the rate equation is expressed as a product of a squared coupling (nonadiabatic and spin–orbit for IC and ISC, respectively) and a Franck–Condon term. We noticed, however, that for ISC in planar molecules, where the spin–orbit coupling is very small, it is essential to incorporate a Herzberg–Teller correction. The work presented here can be considered a necessary stepping stone for establishing the feasibility and expected accuracy of nonradiative rate in the case of molecules that are hypothetical or for which there are no experimental data available.

We observed that the Franck–Condon term could be approximated rather effectively with a limited number of effective modes, which can be determined through a fitting of the spectral density. We have also observed that the spectral densities displayed by organic molecules tend to be similar, and it is even possible to ignore the difference between them and assume a common reference spectral density, which would avoid the need to evaluate normal modes for each molecule. The terms entering the evaluation of the rates (except the energy gap, not considered in this work) are not very sensitive to the basis set size and the functional used, suggesting that they can be evaluated within a few CPU hours for molecules of typical sizes and that they can be incorporated into virtual screening protocols.

While the experimental trends are well reproduced, there are important quantitative discrepancies when the computed and experimental rates are compared. These are ultimately due to the limitation of the harmonic approximation. In the context of virtual screening, it is not reasonable to introduce individual corrections for each molecule, but it may be possible to introduce corrections to the main systematic error based on a large dataset of experimental observations. As this requires a large set of homogeneous and validated data, one can anticipate that the incorporation of nonradiative rate prediction in materials discovery is likely to become a joint effort involving high-throughput experimentation and theory.

## SUPPLEMENTARY MATERIAL

The supplementary material file contains further details on the results of parameters for IC/ISC rate calculations, validation of our approach, and testing of different calculation levels.

## ACKNOWLEDGMENTS

We acknowledge the financial support from the China Scholarship Council (CSC) and the European Union (European Innovation Council, Project No. 101057564).

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**Lei Shi**: Data curation (equal); Formal analysis (equal); Investigation (equal); Resources (equal); Software (equal); Validation (equal); Visualization (lead); Writing – original draft (lead); Writing – review & editing (equal). **Xiaoyu Xie**: Methodology (equal); Software (supporting); Supervision (supporting); Writing – review & editing (equal). **Alessandro Troisi**: Conceptualization (lead); Funding acquisition (equal); Methodology (equal); Project administration (lead); Resources (lead); Supervision (lead); Writing – review & editing (equal).

## DATA AVAILABILITY

The data that support the findings of this study are available within the article and its supplementary material.