The internal conversion (IC) process from S1 to S0 and the intersystem crossing (ISC) transition from T1 to S0 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.

The nonradiative internal conversion (IC) process from the first singlet excited state S1 to the ground state S0 and the intersystem crossing (ISC) process from the first triplet excited state T1 to S0 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 S1 state can generally be de-excited via fluorescence (S1 → S0 + hν), undergo ISC processes (e.g., S1 → T1) or the IC process (S1 → S0). The molecule in the T1 state can also relax to the ground state via the ISC process (T1 → S0). 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 T1 state and S0 states through spatial separation. Nonradiative relaxation channels are detrimental for all devices converting solar radiation to chemical fuels (photocatalysis10) or electrical power (photovoltaic11).

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 104–105) 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 works22–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 TADF36 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.

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 

(1)

where is the Planck constant, Vfi=fV̂i is the interaction between a final state f and the initial state i, and δE is the Dirac delta function.

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

(2)

where ωj and Qj 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

(3)

where ϕi/fq,Qj is the electronic part of the initial/final vibronic state with q being the coordinate of electrons and j=13N6Λwji/fi/fQj is the vibrational part given by a product of 3N − 6 (or 3N − 5 for linear systems) vibrational states Λwji/fi/fQj. For vibrational states, the harmonic approximation is applied based on the potential energy surface of the electronic state ϕi/f, and wji/f is the occupation number of the normal mode j. For convenience, a vibronic state ϕaq,Qjj=13N6ΛwjaaQj is labeled as awja.

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

(4)

In this work, we focused on the IC process from S1 to S0. The index a of vibronic states awja can be set as S1 and S0 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=0for 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

(5)

where ωjϕS0QjϕS1 is the nonadiabatic coupling VNA,j between S0 and S1 for the normal mode j. Squares of vibrational integrals ΛwjS0QjQjΛ0S1Qj and ΛwjS0Qj|Λ0S1Qj can be computed analytically under the harmonic approximation,

(6)
(7)

Here, yj is the Huang–Rhys (HR) factor of the normal mode j for the transition from S1 to S0. It could be computed by projecting the total structure difference between S1 and S0 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 EfEi can be calculated as the sum of the electronic energy gap ES0ES1 and the nuclear energy difference jwjωj under the harmonic approximation. The overall expression of the IC rate can, thus, be written as

(8)

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 VNA,j is nonzero, we therefore assume yj = 0 for a mode j, and the only value of wj contributing to the summation in Eq. (8) is wj = 1 for which ΛwjS0QjΛ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., wj ≠ 1 with VNA,j ≠ 0 for Eq. (6), wj ≠ 0 with yj ≠ 0 for Eq. (7)], the rate in Eq. (8) can be simplified as

(9)

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

(10)

with the effective nonadiabatic coupling VNAeff being

(11)

and

(12)

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

(13)

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

(14)

The ISC rate from T1 to S0 can be calculated similarly by replacing the interaction operator with spin–orbit coupling (SOC) operator ĤSOC. The final formula is

(15)

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 ES1ES0 and ET1ES0 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 functionals55–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. 5860 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 T1 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.

FIG. 1.

Top: 12 molecules for the IC rate calculation in our work and their corresponding labels. Bottom: 17 molecules for the ISC rate calculation in our work and their corresponding labels.

FIG. 1.

Top: 12 molecules for the IC rate calculation in our work and their corresponding labels. Bottom: 17 molecules for the ISC rate calculation in our work and their corresponding labels.

Close modal

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., S0 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 S0, T1, and S1 molecular structures, respectively. Normal modes information is also calculated by DFT in Gaussian 16 using the S0 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 model63,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 S0. This procedure is preferable to various diabatization schemes65 for possible applications in HTVS, where the quality of the diabatization cannot be verified automatically. The NAC matrix elements between S1 and S0 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 S0 and T1 is obtained using TDDFT without the Tamm–Dancoff approximation in the ORCA package.72,73 The S0 geometry is used for both coupling terms.

Since the B3LYP functional74–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 set80 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 functional82 are also used to compute parameters for IC/ISC rate calculations.

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,

(16)

We then define a spectral density Jω introducing a finite broadening σ for each mode,

(17)

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

FIG. 2.

(a) The broadened spectral density plot (black line) and reorganization energy plot (orange points) for bromodibenzoxalene. (b) Fitting results with different effective mode numbers (1, 2, 3, and 4 modes) compared to the original spectra density (black line). The spectral density results around 3000 cm−1 are multiplied by 30 for clarity.

FIG. 2.

(a) The broadened spectral density plot (black line) and reorganization energy plot (orange points) for bromodibenzoxalene. (b) Fitting results with different effective mode numbers (1, 2, 3, and 4 modes) compared to the original spectra density (black line). The spectral density results around 3000 cm−1 are multiplied by 30 for clarity.

Close modal

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,

(18)

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 ωk,λk 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 ωk,λk would be used to obtain frequencies and HR factors (yk=λk/ωk) 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).

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

TABLE I.

Calculated IC rates and related parameters at the B3LYP/6-31G(d) calculation level.

ES1ES0 (cm−1)Expt. logkIC/s1Cal. logkIC/s1λ (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 
ES1ES0 (cm−1)Expt. logkIC/s1Cal. logkIC/s1λ (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 
FIG. 3.

Comparison of experimental (black squares) and calculated (red circles) IC rates for the 12 molecules. B3LYP/6-31G(d) is used for all quantum calculations.

FIG. 3.

Comparison of experimental (black squares) and calculated (red circles) IC rates for the 12 molecules. B3LYP/6-31G(d) is used for all quantum calculations.

Close modal

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 coupling84) 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 ω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 ωk=ωkeff,yk=ykeffλ/λ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).

FIG. 4.

(a) and (b) Normalized spectral densities of the reorganization energy of 12 molecules for the IC process (normalized by total reorganization energy λ). (c) The plot of IC rate vs energy gap for 12 molecules computed with molecule-specific spectral densities (red circles) and with a single rescaled spectral density (blue triangles). (d) The reference spectral density (guaiazulene molecule).

FIG. 4.

(a) and (b) Normalized spectral densities of the reorganization energy of 12 molecules for the IC process (normalized by total reorganization energy λ). (c) The plot of IC rate vs energy gap for 12 molecules computed with molecule-specific spectral densities (red circles) and with a single rescaled spectral density (blue triangles). (d) The reference spectral density (guaiazulene molecule).

Close modal
FIG. 5.

The average computation cost and IC calculation using (a) B3LYP/3-21G level, (c) B3LYP/6-31G(d) level, and (e) B3LYP/def2-TZVP level. The average computation cost of the ISC calculation in (b) B3LYP/3-21G level, (d) B3LYP/6-31 G(d) level, and (f) B3LYP/def2-TZVP level. The inner pie charts are the percentage of computation time of each step.

FIG. 5.

The average computation cost and IC calculation using (a) B3LYP/3-21G level, (c) B3LYP/6-31G(d) level, and (e) B3LYP/def2-TZVP level. The average computation cost of the ISC calculation in (b) B3LYP/3-21G level, (d) B3LYP/6-31 G(d) level, and (f) B3LYP/def2-TZVP level. The inner pie charts are the percentage of computation time of each step.

Close modal

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.

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/s1” 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.

TABLE II.

Calculated ISC rates and related parameters at the B3LYP/6-31G(d) calculation level.

ET1ES0 (cm−1)Expt. logkISC/s1Cal. logkISC/s1Cal. logkISC(HT)/s1λ (cm−1)VSOC (cm−1)aVSOCeff (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 ⋯ ⋯ 
ET1ES0 (cm−1)Expt. logkISC/s1Cal. logkISC/s1Cal. logkISC(HT)/s1λ (cm−1)VSOC (cm−1)aVSOCeff (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

VSOC 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 T1 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,

(19)

Considering out-of-plane vibration modes in the original vibrational set consisting of all normal modes that introduce σπ* excitation to T1 states, the SOC VSOC(Q) can be expressed as

(20)

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

(21)

Here, the effective SOC VSOCeff is written as

(22)

and ω̄SOC is the SOC-weighted energy shift,

(23)

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 ω̄SOC (∼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 ω̄SOC in Eq. (21).

FIG. 6.

Results of calculated ISC rate compared to experimental data (black circles). Blue squares (open and solid) correspond to the computed results via Eq. (15), and red stars show the corrected ISC rates via Eq. (21). Here, solid squares represent molecules for which the HT correction is not introduced: biphenyl, m/p-terphenyl, triphenylbenzene, nonplanar acene (e.g., hexahelicene), and biphenylene.

FIG. 6.

Results of calculated ISC rate compared to experimental data (black circles). Blue squares (open and solid) correspond to the computed results via Eq. (15), and red stars show the corrected ISC rates via Eq. (21). Here, solid squares represent molecules for which the HT correction is not introduced: biphenyl, m/p-terphenyl, triphenylbenzene, nonplanar acene (e.g., hexahelicene), and biphenylene.

Close modal

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.

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.

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.

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

The authors have no conflicts to disclose.

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).

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

1.
N. J.
Turro
,
Modern Molecular Photochemistry
(
University Science Books
,
1991
).
2.
C. M.
Marian
,
Wiley Interdiscip. Rev.: Comput. Mol. Sci.
2
,
187
(
2012
).
3.
S.
Dimitrov
,
B.
Schroeder
,
C.
Nielsen
,
H.
Bronstein
,
Z.
Fei
,
I.
Mcculloch
,
M.
Heeney
, and
J.
Durrant
,
Polymers
8
,
14
(
2016
).
4.
T. J.
Penfold
,
E.
Gindensperger
,
C.
Daniel
, and
C. M.
Marian
,
Chem. Rev.
118
,
6975
(
2018
).
5.
H. C.
Friedman
,
E. D.
Cosco
,
T. L.
Atallah
,
S.
Jia
,
E. M.
Sletten
, and
J. R.
Caram
,
Chem
7
,
3359
(
2021
).
7.
Q.
Peng
,
D.
Fan
,
R.
Duan
,
Y.
Yi
,
Y.
Niu
,
D.
Wang
, and
Z.
Shuai
,
J. Phys. Chem. C
121
,
13448
(
2017
).
8.
B.
Milián-Medina
and
J.
Gierschner
,
Org. Electron.
13
,
985
(
2012
).
9.
C.
Feng
,
S.
Li
,
L.
Fu
,
X.
Xiao
,
Z.
Xu
,
Q.
Liao
,
Y.
Wu
,
J.
Yao
, and
H.
Fu
,
J. Phys. Chem. Lett.
11
,
8246
(
2020
).
10.
X.
Li
,
P. M.
Maffettone
,
Y.
Che
,
T.
Liu
,
L.
Chen
, and
A. I.
Cooper
,
Chem. Sci.
12
,
10742
(
2021
).
11.
N.
Kaur
,
M.
Singh
,
D.
Pathak
,
T.
Wagner
, and
J. M.
Nunzi
,
Synth. Met.
190
,
20
(
2014
).
12.
P.
Friederich
,
A.
Fediai
,
S.
Kaiser
,
M.
Konrad
,
N.
Jung
, and
W.
Wenzel
,
Adv. Mater.
31
,
1970188
(
2019
).
13.
C. F.
Perkinson
,
D. P.
Tabor
,
M.
Einzinger
,
D.
Sheberla
,
H.
Utzat
,
T.-A.
Lin
,
D. N.
Congreve
,
M. G.
Bawendi
,
A.
Aspuru-Guzik
, and
M. A.
Baldo
,
J. Chem. Phys.
151
,
121102
(
2019
).
14.
Y.
Cui
,
P.
Zhu
,
X.
Liao
, and
Y.
Chen
,
J. Mater. Chem. C
8
,
15920
(
2020
).
15.
R.
Pollice
,
G.
dos Passos Gomes
,
M.
Aldeghi
,
R. J.
Hickman
,
M.
Krenn
,
C.
Lavigne
,
M.
Lindner-D’Addario
,
A.
Nigam
,
C. T.
Ser
,
Z.
Yao
, and
A.
Aspuru-Guzik
,
Acc. Chem. Res.
54
,
849
(
2021
).
16.
R.
Gómez-Bombarelli
,
J.
Aguilera-Iparraguirre
,
T. D.
Hirzel
,
D.
Duvenaud
,
D.
Maclaurin
,
M. A.
Blood-Forsythe
,
H. S.
Chae
,
M.
Einzinger
,
D.-G.
Ha
,
T.
Wu
,
G.
Markopoulos
,
S.
Jeon
,
H.
Kang
,
H.
Miyazaki
,
M.
Numata
,
S.
Kim
,
W.
Huang
,
S. I.
Hong
,
M.
Baldo
,
R. P.
Adams
, and
A.
Aspuru-Guzik
,
Nat. Mater.
15
,
1120
(
2016
).
17.
Ö. H.
Omar
,
T.
Nematiaram
,
A.
Troisi
, and
D.
Padula
,
Sci. Data
9
,
54
(
2022
).
18.
D.
Padula
,
Ö. H.
Omar
,
T.
Nematiaram
, and
A.
Troisi
,
Energy Environ. Sci.
12
,
2412
(
2019
).
19.
S. A.
Lopez
,
B.
Sanchez-Lengeling
,
J.
de Goes Soares
, and
A.
Aspuru-Guzik
,
Joule
1
,
857
(
2017
).
20.
W.
Sun
,
Y.
Zheng
,
K.
Yang
,
Q.
Zhang
,
A. A.
Shah
,
Z.
Wu
,
Y.
Sun
,
L.
Feng
,
D.
Chen
, and
Z.
Xiao
,
Sci. Adv.
5
,
eaay4275
(
2019
).
21.
W.
Li
,
H.
Ma
,
S.
Li
, and
J.
Ma
,
Chem. Sci.
12
,
14987
(
2021
).
22.
W.
Siebrand
,
J. Chem. Phys.
47
,
2411
(
1967
).
23.
W.
Siebrand
,
J. Chem. Phys.
46
,
440
(
1967
).
24.
W.
Siebrand
and
D. F.
Williams
,
J. Chem. Phys.
49
,
1860
(
1968
).
25.
R.
Englman
and
J.
Jortner
,
Mol. Phys.
18
,
145
(
1970
).
26.
A.
Nitzan
,
J.
Jortner
, and
P. M.
Rentzepis
,
Chem. Phys. Lett.
8
,
445
(
1971
).
27.
K. F.
Freed
,
Acc. Chem. Res.
11
,
74
(
1978
).
28.
R.
Islampour
and
M.
Miralinaghi
,
J. Phys. Chem. A
111
,
9454
(
2007
).
29.
S.
Rashev
and
D. C.
Moule
,
J. Chem. Phys.
130
,
134307
(
2009
).
30.
Y.
Niu
,
Q.
Peng
,
C.
Deng
,
X.
Gao
, and
Z.
Shuai
,
J. Phys. Chem. A
114
,
7817
(
2010
).
31.
I.
Gimondi
,
C.
Cavallotti
,
G.
Vanuzzo
,
N.
Balucani
, and
P.
Casavecchia
,
J. Phys. Chem. A
120
,
4619
(
2016
).
32.
R. R.
Valiev
,
V. N.
Cherepanov
,
G. V.
Baryshnikov
, and
D.
Sundholm
,
Phys. Chem. Chem. Phys.
20
,
6121
(
2018
).
33.
R. R.
Valiev
,
R. T.
Nasibullin
,
V. N.
Cherepanov
,
G. V.
Baryshnikov
,
D.
Sundholm
,
H.
Ågren
,
B. F.
Minaev
, and
T.
Kurtén
,
Phys. Chem. Chem. Phys.
22
,
22314
(
2020
).
34.
R. R.
Valiev
,
R. T.
Nasibullin
,
V. N.
Cherepanov
,
A.
Kurtsevich
,
D.
Sundholm
, and
T.
Kurtén
,
Phys. Chem. Chem. Phys.
23
,
6344
(
2021
).
35.
P.
Cui
and
Y.
Xue
,
Appl. Nanosci.
11
,
2837
(
2021
).
36.
K.
Shizu
,
M.
Uejima
,
H.
Nomura
,
T.
Sato
,
K.
Tanaka
,
H.
Kaji
, and
C.
Adachi
,
Phys. Rev. Appl.
3
,
014001
(
2015
).
37.
Z.
Shuai
,
D.
Wang
,
Q.
Peng
, and
H.
Geng
,
Acc. Chem. Res.
47
,
3301
(
2014
).
38.
T. E.
Sharp
and
H. M.
Rosenstock
,
J. Chem. Phys.
41
,
3453
(
1964
).
39.
R.
Botter
,
V. H.
Dibeler
,
J. A.
Walker
, and
H. M.
Rosenstock
,
J. Chem. Phys.
44
,
1271
(
1966
).
40.
V. I.
Baranov
,
L. A.
Gribov
, and
B. K.
Novosadov
,
J. Mol. Struct.
70
,
1
(
1981
).
41.
V. I.
Baranov
and
L. A.
Gribov
,
J. Mol. Struct.
70
,
31
(
1981
).
42.
S. H.
Lin
,
J. Chem. Phys.
44
,
3759
(
1966
).
43.
A.
Baiardi
,
J.
Bloino
, and
V.
Barone
,
J. Chem. Theory Comput.
9
,
4097
(
2013
).
44.
Z.
Shuai
,
Chin. J. Chem.
38
,
1223
(
2020
).
45.
A.
Humeniuk
,
M.
Bužančić
,
J.
Hoche
,
J.
Cerezo
,
R.
Mitrić
,
F.
Santoro
, and
V.
Bonačić-Koutecký
,
J. Chem. Phys.
152
,
054107
(
2020
).
46.
D. R.
Yarkony
,
Rev. Mod. Phys.
68
,
985
(
1996
).
47.
F.
Bernardi
,
M.
Olivucci
, and
M. A.
Robb
,
Chem. Soc. Rev.
25
,
321
(
1996
).
48.
B. G.
Levine
and
T. J.
Martínez
,
Annu. Rev. Phys. Chem.
58
,
613
(
2007
).
49.
S.
Matsika
and
P.
Krause
,
Annu. Rev. Phys. Chem.
62
,
621
(
2011
).
50.
A. W.
Kohn
,
Z.
Lin
, and
T.
Van Voorhis
,
J. Phys. Chem. C
123
,
15394
(
2019
).
51.
D. J.
Griffiths
and
D. F.
Schroeter
,
Introduction to Quantum Mechanics
(
Cambridge University Press
,
2018
).
52.
P. F.
Barbara
,
T. J.
Meyer
, and
M. A.
Ratner
,
J. Phys. Chem.
100
,
13148
(
1996
).
53.
J. R.
Reimers
,
J. Chem. Phys.
115
,
9103
(
2001
).
54.
R. P.
Fornari
,
J.
Aragó
, and
A.
Troisi
,
J. Chem. Phys.
142
,
184105
(
2015
).
55.
S.
Refaely-Abramson
,
R.
Baer
, and
L.
Kronik
,
Phys. Rev. B
84
,
075144
(
2011
).
56.
T.
Körzdörfer
and
J.-L.
Bredas
,
Acc. Chem. Res.
47
,
3284
(
2014
).
57.
E.
Cho
,
V.
Coropceanu
, and
J.-L.
Brédas
,
J. Chem. Theory Comput.
16
,
3712
(
2020
).
58.
G. D.
Gillispie
,
Chem. Phys. Lett.
63
,
193
(
1979
).
59.
B. D.
Wagner
,
M.
Szymanski
, and
R. P.
Steer
,
J. Chem. Phys.
98
,
301
(
1993
).
60.
D.
Tittelbach-Helmrich
and
R. P.
Steer
,
Chem. Phys.
197
,
99
(
1995
).
61.
M.
Montalti
,
A.
Credi
,
L.
Prodi
, and
M. T.
Gandolfi
,
Handbook of Photochemistry
(
CRC Press
,
2006
).
62.
M. J.
Frisch
,
G. W.
Trucks
,
H. B.
Schlegel
,
G. E.
Scuseria
,
M. A.
Robb
,
J. R.
Cheeseman
,
G.
Scalmani
,
V.
Barone
,
G. A.
Petersson
,
H.
Nakatsuji
,
X.
Li
,
M.
Caricato
,
A. V.
Marenich
,
J.
Bloino
,
B. G.
Janesko
,
R.
Gomperts
,
B.
Mennucci
,
H. P.
Hratchian
,
J. V.
Ortiz
,
A. F.
Izmaylov
,
J. L.
Sonnenberg
,
Williams
,
F.
Ding
,
F.
Lipparini
,
F.
Egidi
,
J.
Goings
,
B.
Peng
,
A.
Petrone
,
T.
Henderson
,
D.
Ranasinghe
,
V. G.
Zakrzewski
,
J.
Gao
,
N.
Rega
,
G.
Zheng
,
W.
Liang
,
M.
Hada
,
M.
Ehara
,
K.
Toyota
,
R.
Fukuda
,
J.
Hasegawa
,
M.
Ishida
,
T.
Nakajima
,
Y.
Honda
,
O.
Kitao
,
H.
Nakai
,
T.
Vreven
,
K.
Throssell
,
J. A.
Montgomery
, Jr.
,
J. E.
Peralta
,
F.
Ogliaro
,
M. J.
Bearpark
,
J. J.
Heyd
,
E. N.
Brothers
,
K. N.
Kudin
,
V. N.
Staroverov
,
T. A.
Keith
,
R.
Kobayashi
,
J.
Normand
,
K.
Raghavachari
,
A. P.
Rendell
,
J. C.
Burant
,
S. S.
Iyengar
,
J.
Tomasi
,
M.
Cossi
,
J. M.
Millam
,
M.
Klene
,
C.
Adamo
,
R.
Cammi
,
J. W.
Ochterski
,
R. L.
Martin
,
K.
Morokuma
,
O.
Farkas
,
J. B.
Foresman
, and
D. J.
Fox
, Gaussian 16, Revision C.01 (Gaussian Inc., Wallingford, CT,
2016
).
63.
S.
Miertuš
,
E.
Scrocco
, and
J.
Tomasi
,
Chem. Phys.
55
,
117
(
1981
).
64.
J.
Tomasi
,
B.
Mennucci
, and
R.
Cammi
,
Chem. Rev.
105
,
2999
(
2005
).
65.
N. S.
Hush
,
Prog. Inorg. Chem.
8
,
391
(
1967
).
66.
R.
Send
and
F.
Furche
,
J. Chem. Phys.
132
,
044107
(
2010
).
67.
W. P.
Bricker
,
P. M.
Shenai
,
A.
Ghosh
,
Z.
Liu
,
M. G. M.
Enriquez
,
P. H.
Lambrev
,
H.-S.
Tan
,
C. S.
Lo
,
S.
Tretiak
,
S.
Fernandez-Alberti
, and
Y.
Zhao
,
Sci. Rep.
5
,
13625
(
2015
).
68.
P. M.
Shenai
,
S.
Fernandez-Alberti
,
W. P.
Bricker
,
S.
Tretiak
, and
Y.
Zhao
,
J. Phys. Chem. B
120
,
49
(
2016
).
69.
D. B.
Lingerfelt
,
D. B.
Williams-Young
,
A.
Petrone
, and
X.
Li
,
J. Chem. Theory Comput.
12
,
935
(
2016
).
70.
Y.
Sugihara
,
N.
Inai
,
M.
Taki
,
T.
Baumgartner
,
R.
Kawakami
,
T.
Saitou
,
T.
Imamura
,
T.
Yanai
, and
S.
Yamaguchi
,
Chem. Sci.
12
,
6333
(
2021
).
71.
A.
Manian
,
R. A.
Shaw
,
I.
Lyskov
, and
S. P.
Russo
,
J. Chem. Theory Comput.
18
,
1838
(
2022
).
72.
F.
Neese
,
Wiley Interdiscip. Rev.: Comput. Mol. Sci.
2
,
73
(
2012
).
73.
F.
Neese
,
Wiley Interdiscip. Rev.: Comput. Mol. Sci.
8
,
e1327
(
2018
).
74.
C.
Lee
,
W.
Yang
, and
R. G.
Parr
,
Phys. Rev. B
37
,
785
(
1988
).
75.
A. D.
Becke
,
J. Chem. Phys.
98
,
1372
(
1993
).
76.
P. J.
Stephens
,
F. J.
Devlin
,
C. F.
Chabalowski
, and
M. J.
Frisch
,
J. Phys. Chem.
98
,
11623
(
1994
).
77.
R.
Grotjahn
and
M.
Kaupp
,
J. Chem. Theory Comput.
16
,
5821
(
2020
).
78.
J.
Wang
and
B.
Durbeej
,
J. Comput. Chem.
41
,
1718
(
2020
).
79.
W. J.
Hehre
,
R.
Ditchfield
, and
J. A.
Pople
,
J. Chem. Phys.
56
,
2257
(
1972
).
80.
F.
Weigend
and
R.
Ahlrichs
,
Phys. Chem. Chem. Phys.
7
,
3297
(
2005
).
81.
J. S.
Binkley
,
J. A.
Pople
, and
W. J.
Hehre
,
J. Am. Chem. Soc.
102
,
939
(
1980
).
82.
Y.
Zhao
and
D. G.
Truhlar
,
Theor. Chem. Acc.
120
,
215
(
2008
).
83.
C.
Wang
,
J.
Wang
,
N.
Wu
,
M.
Xu
,
X.
Yang
,
Y.
Lu
, and
L.
Zang
,
RSC Adv.
7
,
2382
(
2017
).
84.
P.
Habitz
and
C.
Votava
,
J. Chem. Phys.
72
,
5532
(
1980
).
85.
B. R.
Henry
and
W.
Siebrand
,
J. Chem. Phys.
54
,
1072
(
1971
).
86.
F.
Metz
,
S.
Friedrich
, and
G.
Hohlneicher
,
Chem. Phys. Lett.
16
,
353
(
1972
).
87.
V.
Lawetz
,
G.
Orlandi
, and
W.
Siebrand
,
J. Chem. Phys.
56
,
4058
(
1972
).

Supplementary Material