The role of X–H bonds (X = C, N and O) in internal conversion processes: dibenzoterrylene as an example

Abstract

We have developed a theoretical framework for calculating rate constants of internal conversion (k_IC) in the Franck–Condon (FC) and Herzberg–Teller (HT) approximations. The method accounts for anharmonic vibrational contributions and the Duschinsky effect. Our approach employs recursive dynamic programming to sum over multiple vibrational quantum number combinations and uses a Lagrange-multiplier technique with dispersion broadening to improve the accuracy of the calculated rate constants. We validate the methods by performing calculations on dibenzoterrylene (DBT), which is a molecule emitting in the near-infrared spectral range. The calculations confirm that anharmonic vibrational effects are the main contribution to k_IC, while the Duschinsky effect is significant only for molecules whose lowest excitation energy exceeds 22 000 cm⁻¹. The contributions of the individual X–H bonds are quantified by using the X–H mode approximation (k_IC-XH) and the XH bond approximation (k_IC-proton). The calculations show that the CH bonds of the tetracene moiety of DBT have the largest contribution to k_IC. Deuteration of these bonds leads to a significant decrease in k_IC with complete deuteration resulting in the largest overall effect. The calculated rate constants highlight the important role of the X–H bonds as acceptors of electronic excitation energy, offering strategies for modulating the k_IC through selective substitution of the hydrogen atoms with heavier atoms such as D, F or Cl.

---

1 Introduction

The fluorescence quantum yield (Φ_fl) is the primary indicator of the luminescence efficiency of fluorophores.^1,2 It shows how a large fraction of the electronic excitation energy is emitted as photons in comparison to other nonradiative deactivation pathways. When intermolecular interactions are weak, intramolecular nonradiative channels such as internal conversion (IC) and intersystem crossing (ISC) are the main processes responsible for quenching the fluorescence.¹ For organic molecules for which the energy gap between the first excited electronic state (S₁) and the ground electronic state (S₀) is in the near-infrared (NIR) region, IC is the main nonradiative deactivation pathway.³ More complex cluster systems, comprising organic and inorganic elements—such as aptamer complexes with luminescent metal–ion clusters—may have intersystem crossing as the dominant nonradiative decay channel rather than internal conversion.⁴

For organic NIR molecules, Φ_fl is determined by the competition between the radiative channel, i.e., the fluorescence rate constant (k_r) and the IC rate constant (k_IC). Molecules emitting in the NIR range usually exhibit a low fluorescence quantum yield of <0.2–0.3.^5–7

Dibenzoterrylene (DBT) is an NIR emitter^8–10 with a large Φ_fl of 0.58.¹¹ DBT can be used as a single-photon source¹² in nanophotonics¹¹ and in single-molecule spectroscopy.¹⁰ DBT also has a large Stark effect.¹³ Its rate constant of intersystem crossing (k_ISC) is small making IC the main quenching mechanism.^9,11 Recently, all deactivation rate constants from the S₁ state including the Φ_fl were measured with high accuracy.¹¹ The photophysical properties of deuterated DBT analogs were also investigated, demonstrating that varying the degree of X–H (where X = C, N, or O) deuteration significantly decreases k_IC. However, the specific positions of the deuterium atoms were not determined.¹¹

Measurements of the deuteration effect on k_IC have shown that the main IC acceptors of excited electronic energy are vibrational modes associated with X–H bond vibrations.¹⁴

In this work, we employ our previously developed methods and the methods developed here for calculating k_IC of DBT to elucidate the role of the X–H bonds in the IC process. We use a general theoretical framework that considers the IC contributions from all vibrational modes as well as more specialized approaches that focus only on the contributions from the X–H bonds.^15–18

In the next section, we give a brief overview of the developed techniques for calculating rate constants for internal conversion. The methodology developed in this work is presented in Section 3. The computational details are given in Section 4. The results of the calculations on DBT are reported in Section 5. The main results of the study are summarized in Section 6.

2 Methods for calculating IC rate constants

Internal conversion is a nonradiative process in which the electronic excitation energy is converted into molecular vibrations of the final electronic state without changing the electronic spin state.¹ To the best of our knowledge, it was first shown by Adirovich that IC is described by the nonadiabatic coupling operator.¹⁹

The theory of calculating the IC rate constants (k_IC) was initially formulated by Bixon and Jortner, where they presented the computational formalism.²⁰ A more general expression for evaluating k_IC was reported by Plotnikov and Konoplev.^14,21–23 Their expression for the IC rate constant is

where i is the initial electronic state, f is the final electronic state, and m and n are vibrational levels of i and f, respectively. At low temperatures, one can assume that m = 0. Γ_fn is the relaxation width of the vibrational level |fn〉, and Δ_if = |E_i0 − E_fn| is the energy difference between the vibrational level of the initial electronic state and vibrational level of the final electronic state, and V_i0,fn is the matrix element of a perturbation operator. Eqn (1) can be applied not only to IC but also to ISC processes. In the latter case, the perturbation operator is the spin–orbit coupling operator.¹⁴ Analysis of eqn (1) shows that the electronic excitation energy is converted into vibrational states of f, which does not involve all vibrational modes but only specific accepting and promoting vibrational modes are excited.^1,14 The bands of the involved vibrational modes are broadened by low-frequency vibrational modes. The bandwidth of the transition can be described by a Lorentzian function according to eqn (1). The transition is a two-step process: the energy of the initial state (i0) excites the accepting and promoting modes of the final state (fn) and then the energy irreversibly dissipates into excitations of the remaining vibrational modes that have a lower energy.¹ The second irreversible step leads to eqn (1).¹ The promoting modes are the ones contributing most to the nonadiabatic matrix elements, whereas the accepting modes have the largest overlap of the vibrational wave functions of the initial and final electronic states.^1,14

Plotnikovs analysis of eqn (1) showed that the most important promoting modes are the X–H (X = C, N, O) vibrational modes with frequencies of about 3300 cm⁻¹, while the accepting modes are the vibration modes of the C–C bonds.¹⁴ Typically, one or two C–C modes with energies of (∼1400 cm⁻¹ and ∼400 cm⁻¹) play the main role.^14,21–23 Based on these findings, Plotnikov derived a practical expression consisting of an expansion of the nonadiabatic coupling matrix elements (NACME) in contributions from vibrational modes of the X–H bonds.^14,21–23 The anharmonicity of the vibrational modes was considered using Morse functions. However, they did not calculate the NACME and spin–orbit coupling matrix elements (SOCME) but estimated them.^14,21–23

Artyukhov and Mayer calculated the electronic NACME and SOCME at the semi-empirical INDO level,²⁴ while the integrals involving vibrational functions were estimated from the fitting curve by Plotnikov et al.¹⁴ Artyukhov performed calculations in the Franck–Condon (FC) approximation. However, they did not calculate integrals of vibrational functions nor performed any summation over vibrational levels.

Valiev et al. recently adapted eqn (1) for modern quantum-chemical methods, where all parameters except Γ_fn were calculated from first principles.¹⁵ Γ_fn is in the femtosecond range and does not depend much on vibrational quantum numbers when the energy difference is in the visible range of the spectrum.^1,3,15 They chose a Γ_fn value of 10¹⁴ s⁻¹.¹⁵ Γ_fn can also be calculated using the Lax–Pekar model,^25–28 which also yields Γ values in the femtosecond range.²⁹ The summation in eqn (1) was carried out over normal modes of the X–H vibrations as well as for 1 to 10 accepting modes depending on the molecule.¹⁵

Valiev et al. showed that the Herzberg–Teller (HT) contribution to k_IC can be significant. It may be of the same order of magnitude as the FC contribution.¹⁶ Kovarskii et al. were the first ones to point out that the HT contribution can be substantial.^30–32 Valiev et al. reported very recently that vibrational anharmonicity can significantly affect k_IC values.^17,18,33 They showed that the anharmonic effects originating from X–H vibrations are important. By considering the anharmonicity, they reproduced experimental IC rate constants with a one order of magnitude accuracy for molecules with large energy gaps.³³ The importance of the anharmonicity was also verified by repeating the calculations on deuterated molecules and comparing the calculated rate constants with experimental data.³⁴ Makshantsev also showed that the anharmonicity plays an important role for IC when the energy gap is ∼3 eV.^35,36

Valiev et al. adapted eqn (1) to modern quantum-chemical methods and used an approximate expression for the nonadiabatic coupling matrix elements (NACME), which involves only X–H bonds.¹⁸ The methodology has been successfully applied to calculations of IC rate constant for a variety of molecules.^4,37–44 The methodology was recently extended for studying the influence of external magnetic and electric fields on k_IC values.^45,46

Calculating nonradiative rate constants using eqn (1) has computational advantages and disadvantages. Although it enables explicit computations of combinations of excited vibrational modes of the final electronic state, it is not feasible to calculate all combinations. For example, considering that benzene with 30 vibrational modes has a maximum vibrational quantum number of 10, the total number of combinations is 10³⁰, which is impossible to handle. However, all these combinations are not needed. Since the energy conservation condition must be met, the number of possible combinations is much smaller. Furthermore, considering only the effective vibrational levels contributing to the IC rate constant results in a small subset of all possible combinations.^1,14 One can use the method of undetermined Lagrange multipliers to identify the contributing combinations,^1,17,33 which has been reported by Valiev et al.¹⁷ and Manian et al.^29,47 The method of Lagrangian multipliers is an efficient approach to finding the vibrational levels whose combinations satisfy the energy conservation condition.

The rate-constant expression in eqn (1) belongs to the time-independent (TI) formalism.^1,2 An alternative form of eqn (1) can be derived using the correlation-function formalism in the time-dependent TD domain. Using a Fourier transform, the summation over the quantum numbers can be replaced by integration over time.² This approach was originally developed by Kubo and Toezuawa,⁴⁸ Lax,²⁵ and Pekar.^49,50 The methodology has more recently been refined and adapted to modern quantum-chemical methods by Marian et al.,² Shuai et al.,^51–53 and Santoro et al.⁵⁴ Valiev et al.⁴⁶ used a slightly different TD approach based on temperature-dependent quantum Green's functions, which is a methodology developed by Tyablikov and Moskalenko.⁵⁵

The TD methods avoid summation over a large number of vibrational states. However, the integrand oscillates rapidly requiring an accurate numerical integration approach.^1,14,56 The integral diverges in most cases requiring regularization by using an appropriate damping factor such as Lorentzian, Gaussian, or Voigt broadening.^1,14,56 A fitting parameter is then introduced in the exponent of the broadening function. The parameter must be chosen to be large enough to ensure the convergence of the integral but small enough to avoid that it significantly affects the calculated k_IC. The fitting parameter is associated with vibrational relaxation, but in many cases it is simply a fitting parameter.¹ The TI method using eqn (1) incorporates intramolecular vibrational relaxation, which depends only weakly on the environment and can in principle be calculated using the Lax-Pekar model.^25,49,50

In this work, we develop a new TI approach based on eqn (1). The summation is performed with a dynamic programming approach that employs recursive calculations using a two-array method. Two arrays of dimension E are created, where E is the integer part of the energy in cm⁻¹ of the electronic transition. The dynamic programming approach enables calculations of all combinations that yield the transition energy. For a molecule with E = 40 000 cm⁻¹, 30 vibrational modes and considering 20 vibrational levels of each mode, there are only 24 × 10⁶ possible combinations, which is manageable even on a laptop. A similar approach in the harmonic approximation was recently reported by Manian et al.²⁹ We also show how various methods to calculate the IC rate constant of DBT using eqn (1) lead to the same conclusion about the significance of X–H bonds.

3 Methods based on Plotnikov's expression

The nonradiative rate constant (k_nr) can be calculated using eqn (1) when the transition is much slower than the vibrational relaxation (Γ_fn).¹⁴ The expression can be simplified when the temperature is low, since the initial vibrational state is in the lowest level (m = 0) at ambient temperature.^14–18 The sum over n satisfies the energy conservation condition. Γ_fn is typically ∼10¹⁴ s⁻¹ and depends weakly on n, and the energy difference Δ_if. The expression for the nonradiative rate constant can then be written as .^14–18 Γ_fn is calculated using the Lax–Pekar model using^{25–28,49,50}where y_j is the Huang–Rhys factor of the j-th accepting mode, ω_j is its frequency, and σ_j are the Bose–Einstein occupation numbersT is temperature and k_B is the Boltzmann constant. The Huang–Rhys factors y_j were calculated aswhere K_j is the displacement of the harmonic oscillator along the j-th normal coordinate, given byIn eqn (5), f_νq is the projection of the force acting on the ν-th atom along the q direction in the final electronic state, evaluated at the geometry of the initial state, m_ν is the mass of the ν-th atom, and L_νqj are matrix elements of the linear transformation between the Cartesian (R) and the normal coordinates (Q),^25,28 which is R_νq − R_0νq = M_ν^−1/2L_νqjQ_j.

In the adiabatic Herzberg–Teller (HT) approximation and neglecting the Duschinsky transformation, k_IC can be written as³³

D, P_j and P_jj′ considering the NACME between the i-th and f-th electronic states are given bywhere Ψ is the wave function of the i-th and f-th electronic states, and are the NACMEs of the first and second order, respectively. M_ν is the mass of the ν-th atom.

g_j, t_j and b_j depend on the choice of the type of the vibrational wave function ψ. In the harmonic approximation, they are given by

where n_j is the vibrational quantum number and ω_j is the frequency of vibrational mode j. The anharmonicity of the j-th mode can be considered by introducing the Morse vibrational wave functions³³where z = 2β exp(−α(R − R_e)), and b_n = 2β − 2n − 1, . D_e is the dissociation energy at the equilibrium distance (R_e) and α is the anharmonicity constant. The normalization factor N_n is is the n-th Laguerre polynomial, and Γ(b_n + n + 1) is the Gamma function.

D_e and α can be derived from the anharmonicity parameter (χ_e) and the vibrational energy (ω_e) as D_e = ω_e/4χ_e.^33,54 The dissociation energy of X–H bonds is ∼4.5 eV and ω_e of X–H bonds is ∼3300 cm⁻¹, resulting in a χ_e value of 0.02,³³ which yields good estimates for anharmonic effects.

In the anharmonic approximation at the Morse level, the expressions for g_j, t_j and b_j are

where Δ = exp(−αK), K is the displacement of the equilibrium position of the oscillator, and I_n(A,B,C) isThe derivative of I_n(A,B,C) in eqn (16) can be calculated numerically and the analytical expression of I_n(A,B,C) is

We apply the Morse model to vibrational modes whose energy is larger than 2000 cm⁻¹. Their vibrational wave functions associated to X–H(D) vibrational modes are replaced by a Morse function.

3.1 Duschinsky transformation

The Duschinsky effect was previously introduced using perturbation theory with the perturbation operator whereas the nuclear wave function was unperturbed.³³ The Duschinsky transformation can then be written as Q⃑′ = JQ⃑ + K⃑, where Q⃑ and Q⃑′ are the normal coordinates of the final and initial electronic states, respectively, J is the Duschinsky rotation matrix, and K is the displacement. In the anharmonic approximation, the Duschinsky rotation has a small effect on k_IC,^51–53 whereas in the harmonic approximation the Duschinsky effect on k_IC is significant.³³ In the previous implementation at the FC level by Valiev et al., eqn (6) was modified by replacing P_i with . At the HT level, P_i was replaced by and P_ij was replaced by .

We have here developed an alternative approach to consider the Duschinsky effect in k_IC calculations. Since the Duschinsky matrix is orthogonal i.e., J^TJ = I, it can be written as J = I + A, where I is the unit matrix and the a_ij matrix elements of A are generally small. The Duschinsky rotation can then be written as yielding

Expanding the vibrational wave functions (ψ_n) of the final state in a Taylor series expansion and considering the two first terms, one obtainsSince the a_ij elements are much smaller than 1, higher-order contributions can be neglected. Inserting eqn (19) and (20) into eqn (6) and taking only the second term into account yields the working expression in FC approximationThe new matrix elements (T_i) are in the harmonic approximation given byand in the anharmonic approximation, they areThe quadratic terms in eqn (21) are the main contributions as in eqn (6). The largest contribution originates from the second term with T_i. The new matrix elements appearing in the HT approximation areandwhich are given byandwhere . In the anharmonic approximation:whereThe expression for k_IC in the HT and Duschinsky approximations is

3.2 The dynamic programming algorithm

The summation in eqn (6), (21) and (37) is performed over all combinations of the accepting modes fulfilling the energy conservation condition Δ_if = n⃑·[small omega, Greek, vector]. We consider the combinations using a recursive scheme similar to the one used in ref. 29. Such recursive schemes are well known in the field of dynamic programming and are applied to combinatorial problems.⁵⁷ We use a two-array algorithm and consider anharmonic effects. The dimension of the arrays is obtained by first discretizing Δ_if + δ into steps of 1 cm⁻¹, where δ is 200 cm⁻¹. The final dimension of the arrays is obtained by multiplying the integer of the discretized energy (E + δ) range by the number of vibrational modes times the number of considered vibrational levels of each mode.

We use the following dynamic programming algorithm

1. Initialize an array of probabilities A(E + δ).

2. Set A(0) = 1 and all other elements are set to zero.

3. Iterate over vibrational modes:

(a) Loop over all vibrational modes i = 1,…,N.

(b) Initialize a temporary array B(E + δ) with all elements set to 0.

(c) Loop over the discrete energy values E_j = 0,…,E + δE.

(d) For each vibrational quantum number n_j = 0,…,n_max, do:

i. Compute the vibrational energy ε_nj.

ii. Calculate the energy:

E_new,j = E_j + ε_nj

iii. Update probability:

B(E_new,j) = B(E_new,j) + g_j(n_j)²·A(E_j).

4. Final step: update the probability array:

A(E) = B(E).

5. Add the probabilities within the energy range δ.

The dynamic programming approach enables computation of all combinations in the configuration space specified by the discretized energy range E + δ, the number of vibrational modes (N), and the number of considered vibrational levels of each mode (n_max). The number of such combinations is E × N × n_max, which is the order of millions of combinations, when N is ∼100–200 and n_max is ∼10 and E is ∼40 000 cm⁻¹. This is a very small subset of all possible combinations.

3.3 The Lagrange multipliers technique

The maximum probability of energy-conserving combinations can also be estimated using Lagrange's method of undetermined multipliers,^17,33 which works only when the probabilities are calculated for accepting modes having n_j ≥ 1 but not for modes with n_j < 1. We used this scheme in our previous studies, where only accepting modes were considered, and the X–H vibrations were anharmonic modes.

We extend the method of Lagrange multipliers to energy dispersion contributions considering effective modes. Vibrational quantum numbers of the accepting modes (n_k ≥ 1) fulfilling the energy condition with n_k = y_k exp ω_kλ and the corresponding Lagrange multipliers (λ) are determined. The maximum FC factor (g_λ) can then be calculated for the accepting modes using¹

However, lots of n_k combinations can give FC factors values close to g_λ. We correct g_λ for dispersion contributions and estimate a corrected FC factor using

The FC factors in eqn (38) and (39) were obtained in the harmonic approximation using the calculated vibrational levels of the acceptive modes.

3.4 The X–H approximation

We assume that only the X–H vibrations are promoting modes and approximate the combinations of C–C modes by excitations of two modes at 1400 and 400 cm⁻¹ by one. The k_IC rate constant can then be calculated usingwhere P_ν is calculated only for the ν-th proton and the X_ν atom of the X_ν–H_ν moiety usingL_νq is calculated for a single X_ν–H_ν bond using , where q_Xν, q_Hν are the Cartesian coordinates of the X_ν and H_ν atoms, R_XH is the bond length of the X–H bond, L is the relationship between the Cartesian coordinates and the normal coordinates of the atoms of the X–H moiety. L is almost equal to one because X is much heavier than H. t_ν is calculated using eqn (15) at the anharmonic level for X–H bonds with the typical vibrational energy (ω) of 3300 cm⁻¹ and the χ_e value is 0.02. The term is calculated using a fitting curve that depends only on the energy difference (E_pq) between the initial (Ψ_p) and the final (Ψ_q) electronic states.³³

Calculations of the NACME can be avoided by estimating it from the expansion coefficients of the molecular orbitals in the basis functions (χ) at the X–H moieties. An approximate IC rate constant is then given by

k_IC-XH = 1.6 × 10⁹ × 〈Ψ_p|[capital Lambda, Greek, circumflex]|Ψ_q〉², (42)

where the NACME is obtained usingand the electronic factor of each X–H bond isa and b denote virtual molecular orbitals (MO) and i, j are the indices of the occupied MOs. c_ki, c_k'j, c_ka and c_k'b are MO coefficients of the occupied and virtual MOs, respectively. We assume that is equal to 0.1 for each X–H bond. A^p_ia and A^q_jb are configurational interaction coefficients of states p and q. The X–H approximations enable calculations of k_IC for large molecules. It also identifies which specific X–H bonds provide the largest contributions to the IC rate constant.

We employ the X–H approximation, the dynamic programming algorithm, and the method of Lagrange multipliers technique to calculate the rate constant for internal conversion (k_IC) of DBT and deuterated DBT.

4 Computational details

DBT has three conformers of which the symmetrically twisted one belonging to the C_2h point group has the lowest energy,⁵⁸ which was confirmed by calculations at the density functional theory (DFT) level using the B3LYP functional, def2-TZVP basis sets and the D4 dispersion correction term.^59–62 The energy of the studied conformer calculated with Turbomole version 7.8^63,64 at the B3LYP/def2-TZVP/D4 level is 3.8 kcal mol⁻¹ and 5.7 kcal mol⁻¹ below the asymmetrically twisted and the saddle-shaped conformer, respectively. Calculations at the complete active space self-consistent field (CASSCF) level with 10 electrons in 10 orbitals using the 6-31G(p,d) basis set⁶⁵ were performed with the FIREFLY software package.⁶⁶ The CASSCF calculations showed that the ground electronic state of DBT is dominated by one Slater determinant. Calculations were also performed with Turbomole at the singles and doubles coupled-cluster level with a perturbed treatment of the triple excitations (CCSD(T)) using def2-TZVP basis sets and the reduced-virtual-space (RVS) approach.^67–70 The occupation numbers of the frontier orbitals, the D₂ diagnostics, the T₂(CCSD) norm, and the C₀² show that DBT has a small multireference character, which is also expected because the open-shell singlet character of polyacenes increases with their length.⁷¹ We have used the MN15 functional because it has been successfully used in calculations of rate constants for molecules whose ground electronic state (S₀) has a weak multireference character.^42,72 Molecular structure optimizations of the S₀ state and of the lowest excited singlet state (S₁) as well as calculations of the harmonic vibrational frequencies (ω), the nonadiabatic coupling matrix element (NACME), and the Huang–Rhys factors y were performed with the MN15 functional and the 6-31G(d,p) basis set using the Gaussian 16 software package.⁶⁵

All vibrational modes were considered in the calculations of the rate constants in the harmonic (k_IC-HARM-FC) and anharmonic (k_IC-ANHARM-FC) FC approximation as well as in the calculation of the rate constants in the harmonic (k_IC-HARM-HT) and anharmonic (k_IC-ANHARM-HT) HT approximation. The Duschinsky effect was considered in the calculations of the rate constants using unperturbed (k_{IC-(AN)HARM-FC/HT-D1}) and perturbed vibrational wave functions (k_{IC-(AN)HARM-FC/HT-D2}). The IC rate constant was also calculated at approximate levels considering only NACME contributions of the atoms of the X–H moiety (k_IC-proton) using eqn (40) and (42).

The contributions of the promoting modes are analyzed using the t_j values obtained with eqn (6). The contribution to k_IC-proton and k_IC-XH from each X–H bond is identified. We consider two deuteration patterns of the X–H bonds. We deuterate the X–H bonds with the largest contribution to k_IC and the X–H bonds with the smallest contribution to k_IC, which are here called the rapid and slow IC transition, respectively.

5 Results and discussion

5.1 Internal conversion for nondeutrated DBT_h20

The vibrational relaxation width (Γ) calculated using eqn (2) is 7.8 × 10¹³ s⁻¹, which is very close to the previously used Γ value of 10¹⁴ s⁻¹. The calculated k_IC rate constants for the nondeuterated DBT_h20 are summarized in Table 1. The calculated k_IC shows that the Herzberg–Teller approximation is important because it increases the rate constant leading to a good agreement with the experimental data.¹¹ The anharmonic and Duschinsky contributions are very small in the HT approximation, whereas in the Franck–Condon approximation, the anharmonic effect increases the rate constant by a factor of 4. In the harmonic FC approximation, the Duschinsky effect using perturbed vibrational wave functions increases k_IC by an order of magnitude. In the anharmonic FC approximation, the Duschinsky effect is small. The k_{IC-ANHARM-FC-D1} rate constant is almost equal to k_IC-HARM-HT, even though k_IC-HARM-FC is a factor of ∼40 smaller than k_IC-HARM-HT, whereas the k_IC-HARM-HT rate constant is equal to k_{IC-ANHARM-HT-D1}.

Table 1 Calculated vibrational relaxation width (Γ in s⁻¹), de-excitation energy E (in cm⁻¹), and k_IC rate constants (in s⁻¹) for DBT_h20 at different levels of approximation using the TDDFT/MN15/6-31G(d,p) level of the electronic structure calculations. The rate constants are calculated using the dynamic programming and Lagrange multipliers algorithms

	Calculated		Experiment^11
E (cm^−1)	11 700		12 050–13 450
Γ (s^−1)	7.8 × 10^13		—
	Dynamic	Lagrangian
	programming	multipliers
kIC-HARM-FC	4.6 × 10^6	3.0 × 10^6	1.1–7.3 × 10^8
kIC-HARM-FC-D1	5.1 × 10^6	3.1 × 10^6
kIC-HARM-FC-D2	2.9 × 10^7	8.8 × 10^6
kIC-ANHARM-FC	1.9 × 10^7	3.1 × 10^6
kIC-ANHARM-FC-D1	2.2 × 10^7	3.1 × 10^6
kIC-ANHARM-FC-D2	2.2 × 10^7	9.0 × 10^6
kIC-HARM-HT	1.7 × 10^8	9.4 × 10^6
kIC-HARM-HT-D1	1.8 × 10^8	9.0 × 10^6
kIC-HARM-HT-D2	2.1 × 10^8	1.1 × 10^7
kIC-ANHARM-HT	1.7 × 10^8	8.4 × 10^6
kIC-ANHARM-HT-D1	1.7 × 10^8	7.9 × 10^6
kIC-ANHARM-HT-D2	1.7 × 10^8	7.7 × 10^6
kIC-proton	5.0 × 10^8
kIC-XH	1.5 × 10^9

---

The probably best calculated IC rate constant (k_{IC-ANHARM-HT-D2}) of 1.7 × 10⁸ s⁻¹ agrees very well with the IC rate constant of 1.1 × 10⁸ s⁻¹ measured for DBT in cyclohexane and with the k_IC of 2.7 × 10⁸ s⁻¹ measured for DBT in toluene.¹¹ The dichloromethane molecules of the solvent seem to enhance the IC relaxation leading to a larger experimental k_IC of 7.3 × 10⁸ s⁻¹.¹¹ A recent study showed that solvent molecules can affect rate constants of nonradiative transitions and the quantum yield of luminescence.⁴

The rate constant calculated using the X–H approximation (k_IC-XH) of 1.5 × 10⁹ s⁻¹ is somewhat larger than the experimental value, whereas the k_IC-proton rate constant of 5.0 × 10⁸ s⁻¹ obtained by calculating the NACME only for the X–H bonds agrees very well with the experimental one. The rate constants calculated using the Lagrange multipliers method and the dynamic programming algorithm agree well. However, they are more than an order of magnitude smaller than the rate constants calculated at the two other levels of approximations probably because the vibrational energies calculated in the FC and harmonic approximation were used in the other levels of approximation (Table 2).

Table 2 The deuteration effect on the rapid and slow IC relaxation pathways and the corresponding IC rate constants (k_IC in s⁻¹)

Rapid	DBTh20	DBTd4-rapid	DBTd8-rapid	DBTd12-rapid	DBTd20
kIC-ANHARM-FC	1.9 × 10^7	1.5 × 10^7	4.8 × 10^6	5.0 × 10^6	3.0 × 10^6
kIC-ANHARM-FC-D1	2.2 × 10^7	1.9 × 10^7	5.3 × 10^6	5.8 × 10^6	3.0 × 10^6
kIC-ANHARM-FC-D2	2.2 × 10^7	1.9 × 10^7	5.3 × 10^6	5.8 × 10^6	3.0 × 10^6
kIC-ANHARM-HT	1.7 × 10^8	1.6 × 10^8	1.2 × 10^8	6.9 × 10^7	1.8 × 10^7
kIC-ANHARM-HT-D1	1.7 × 10^8	1.6 × 10^8	1.2 × 10^8	7.0 × 10^7	1.8 × 10^7
kIC-ANHARM-HT-D2	1.7 × 10^8	1.5 × 10^8	1.0 × 10^8	6.8 × 10^7	1.8 × 10^7

---

Slow	DBTh20	—	DBTd8-slow	DBTd12-slow	DBTd20
kIC-ANHARM-FC	1.9 × 10^7		1.8 × 10^7	1.8 × 10^7	3.0 × 10^6
kIC-ANHARM-FC-D1	2.2 × 10^7		2.0 × 10^7	2.0 × 10^7	3.0 × 10^6
kIC-ANHARM-FC-D2	2.2 × 10^7		2.0 × 10^7	2.0 × 10^7	3.0 × 10^6
kIC-ANHARM-HT	1.7 × 10^8		1.3 × 10^8	7.5 × 10^7	1.8 × 10^7
kIC-ANHARM-HT-D1	1.7 × 10^8		1.2 × 10^8	7.3 × 10^7	1.8 × 10^7
kIC-ANHARM-HT-D2	1.7 × 10^8		1.1 × 10^8	7.2 × 10^7	1.8 × 10^7
kIC-Exp.^11 (toluene)	2.7 × 10^8				1.4 × 10^8
kIC-Exp.^11 (cyclohexane)	1.1 × 10^8				0.6 × 10^7
kIC-Exp.^11 (dichloromethane)	7.3 × 10^8				4.0 × 10^7

---

We studied the energy-dependence of the rate constants by varying the de-excitation energy from 100 cm⁻¹ to 30 000 cm⁻¹. The calculated rate constants as a function of the de-excitation energy in Fig. 1 show that IC rate constants calculated at various approximation levels qualitatively agree in the whole energy range. More figures showing the energy dependence of the rate constants are given in the ESI.† The Duschinsky effect on the IC rate constant is very large when the de-excitation energy exceeds 23 000 cm⁻¹, where k_{IC-ANHARM-FC-D2} is several orders of magnitude smaller than k_IC-ANHARM-FC. The Duschinsky effect becomes more significant when the de-excitation energy is large because the electronic energy is transferred to highly excited vibrational modes leading to a strong mixing of the vibrational wavefunctions.⁵⁴ The Duschinsky approximation assumes that the matrix elements of the A matrix are small, which might not hold when the de-excitation energy exceeds 22 000 cm⁻¹. Rate constants calculated at the Duschinsky level using de-excitation energies larger than 22 000 cm⁻¹ may therefore be underestimated, whereas for de-excitation energies between 100 and 21 000 cm⁻¹ the calculated Duschinsky contribution to the IC rate constant of DBT is reliable.

Fig. 1 The computed IC rate constant of DBT calculated in the FC and HT approximation. The rate constants were calculated in the harmonic and the anharmonic approximation. The Duschinsky contributions (D2) were considered. The rate constants in a logarithmic scale are shown as a function of the de-excitation energy from S₁ to S₀. The energy dependence of the k_IC-proton rate constant calculated in the X–H approximation is also shown. The energy-dependence of rate constants calculated at other approximation levels are reported in the ESI.†

The HT approximation is important for the IC rate constant of DBT. Considering HT terms increases the rate constant by a factor of almost 40 yielding a rate constant that agrees well with the experimental one. The anharmonic and Duschinsky effects are very small in the HT approximation, whereas at the harmonic FC level, the Duschinsky rotation increases the rate constant by a factor of 6 and anharmonicity increases it by a factor of 4.

It should also be noted that all methods explicitly computing the contribution of accepting modes to the FC factors show a slight decrease in k_IC at ∼2000 cm⁻¹, the IC rate constant as a function of the de-excitation energy passes a maximum, which happens when Gibbs free energy is equal to the reorganization energy according to Marcus theory.⁷³

Fig. 2 shows the vibrational modes with the largest contribution to t_j. The main promoting mode is the 165 A_u mode with a vibrational energy of 3260 cm⁻¹, which involves the two X–H bonds belonging to the tetracene moiety. The calculations of the k_IC-proton rate constant suggest that the main contribution originates from the X–H bonds localized on the tetracene fragment, whereas calculations of the k_IC-XH rate constant suggest that also the X–H bonds of the naphthalene moieties contribute to the IC rate constant. All the employed methods show that the main contribution to the IC rate constant of DBT involves the vibrational modes of the X–H bonds of the tetracene moiety. The largest contributions involve combinations of vibrational modes at wavenumbers of ∼3000 cm⁻¹. They also provide the largest contribution to the NACME.

Fig. 2 The promoting vibrational modes of the X–H bonds are encircled: (a) the A_u vibrational mode at 3260 cm⁻¹, (b) the X–H vibrational modes with the largest contributions to k_IC-proton and (c) the X–H vibrational modes with the largest contributions to k_IC-XH.

5.2 Internal conversion for deuterated DBT

Mishra et al.¹¹ could not determine which hydrogen atoms were replaced by deuterium because the synthesis yielded a mixture of compounds with different degrees of deuteration. The IC rate constants were measured for the partially deuterated DBT_d12, which was a mixture of molecules containing 10–12 deuterium atoms, and for DBT_d20, which was a mixture consisting of molecules with 18–20 deuterium atoms.¹¹ The measured k_IC in toluene solution is 2.7 × 10⁸ s⁻¹, 1.7 × 10⁸ s⁻¹ and 1.4 × 10⁸ s⁻¹ for DBT_h20, DBT_d12 and DBT_d20, respectively.¹¹ Perdeuteration decreases k_IC by a factor of ∼2.0. The same ratio was obtained for rate constants measured in cyclohexane solution. The experimental k_IC for DBT_h20 and DBT_d20 is 1.1 × 10⁸ s⁻¹ and 6.0 × 10⁷ s⁻¹, respectively.

We are able to computationally determine the effect of replacing hydrogen atoms with deuterium. We studied the rapid IC transition case, where the C–H moieties associated with the promoting modes of the tetracene fragment were deuterated. First, we deuterated the two C–H moieties (DBT_d4-rapid) shown in Fig. 2a and then the eight C–H moieties (DBT_d8-rapid) on the same moiety as shown in Fig. 2b. Finally, we deuterated all promoting C–H moieties (DBT_d12-rapid) as shown in Fig. 2c. In the slow IC transition case, we deuterated the C–H moieties (DBT_d8-slow) that are not marked in Fig. 2c and all C–H moieties (DBT_d12-slow) that do not belong to tetracene. The perdeuterated molecule is denoted DBT_d20.

The k_IC rate constant decreases systematically when increasing the number of deuterium atoms. The largest effect is obtained for perdeuterated DBT showing that the X–H moieties play an important role for the IC process. In the FC approximation, deuteration of the rapid IC transition channel leads to a smaller k_IC as compared to the slow IC transition case. In the HT approximation, the difference is smaller. Calculations on the deuterated molecules indicate that various X–H moieties have different contributions to the IC rate constant because the probability of accepting the electronic excitation energy differs.

The calculated rate constant for DBT_h20 agrees well with the experimental one, whereas the calculated rate constant for DBT_d20 is smaller than the measured one because it consists of a mixture of molecules with 18–20 deuterium atoms. The nondeuterated C–H moieties can still be a fast IC transition channel.

6 Conclusions

We have developed and implemented new methods to calculate internal conversion rate constants (k_IC) in the Franck–Condon (FC) and Herzberg–Teller (HT) approximations. The computational methods can also account for anharmonic and Duschinsky effects. The summation over energy-conserving combinations of the vibrational modes is performed using a recursive dynamic programming approach or alternatively using a method based on Lagrange multipliers that also includes dispersion-broadening contributions. The methods were tested on the energetically lowest conformer of dibenzoterrylene (DBT), which emits light in the near-infrared region (NIR) of the electromagnetic spectrum.

The anharmonic and Duschinsky effects are important when calculating k_IC in the Franck–Condon (FC) approximation. The anharmonic and Duschinsky effects increase k_IC by a factor of ∼5. However, when considering the two effects k_IC calculated in the FC approximation is still an order of magnitude smaller than the experimental IC rate constant. The IC rate constant calculated in the harmonic Herzberg–Teller (HT) approximation agrees well with the experimental k_IC and the anharmonic and Duschinsky contributions are very small at the HT level. Calculations of the energy dependence of the IC rate constant show that k_IC decreases with increasing de-excitation energy and that the Duschinsky effect is significant when the de-excitation energy exceeds 22 000 cm⁻¹.

The computational methods were used to assess the accuracy of the X–H and proton approximations, which are computationally efficient methods that consider only contributions from X–H moieties to k_IC. The IC rate constants calculated in the X–H and proton approximations are somewhat larger than those obtained at more accurate levels of theory and experimentally. Calculations of rate constants considering only X–H contributions showed that k_IC mainly originates from vibrations of the C–H bonds of the tetracene moiety. Deuteration of these C–H bonds leads to a smaller k_IC than when other parts of the molecule are deuterated. The largest isotope effect is obtained for perdeuterated DBT, which confirms that the X–H moieties are the main IC acceptors of the electronic excitation energy. Calculations of the k_IC-XH and k_IC-proton rate constants enable determination of the individual contributions from various C–H moieties, which offers an opportunity to adjust the internal conversion rate constant by replacing H atoms with a substituent or a heavier atom such as D, F or Cl.

When the final electronic wave function is expanded in nuclear coordinates (eqn (6)), we assume that the series converges and that higher-order terms are much smaller than the first-order term. However, it is possible that for molecules with large energy gaps higher-order Herzberg–Teller effects may become significant, which is an interesting direction for future research.

Author contributions

RRV suggested the theory and the algorithms. BSM derived the expressions. RRV and BSM implemented the algorithms and carried out calculations. RRV and RTN wrote the first version of the manuscript. DS, RRV and RTN discussed and validated the results. All authors contributed by editing and reviewing the manuscript.

Conflicts of interest

There are no conflicts to declare.

Data availability

The data supporting this article have been included as part of the ESI.†

Acknowledgements

This work was supported by the Academy of Finland through projects 340582 (RRV), 340583 (DS) and 346369 (TK).

Notes and references

1. E. S. Medvedev and V. I. Osherov, Radiationless Transitions in Polyatomic Molecules, Springer-Verlag, Berlin, 1995.
2. T. J. Penfold, E. Gindensperger, C. Daniel and C. M. Marian, Chem. Rev., 2018, 118, 6975–7025.
3. R. Englman and J. Jortner, Mol. Phys., 1970, 18, 145–164.
4. R. R. Ramazanov, R. T. Nasibullin, D. Sundholm, T. Kurtén and R. R. Valiev, J. Phys. Chem. Lett., 2024, 15, 10710–10717.
5. M. S. K. Rurack, Anal. Chem., 2011, 83, 1232–1242.
6. A. Alessi, M. Salvalaggio and G. Ruzzon, J. Lumin., 2013, 134, 385–389.
7. E. J. Peterson, W. Qi, I. N. Stanton, P. Zhang and M. J. Therien, Chem. Sci., 2020, 11, 8095–8104.
8. A. Pscherer, M. Meierhofer, D. Wang, H. Kelkar, D. Martn-Cano, T. Utikal, S. Götzinger and V. Sandoghdar, Phys. Rev. Lett., 2021, 127, 133603.
9. A. A. Nicolet, C. Hofmann, M. A. Kol'chenko, B. Kozankiewicz and M. Orrit, ChemPhysChem, 2007, 8, 1215–1220.
10. F. Jelezko, P. Tamarat, B. Lounis and M. Orrit, J. Phys. Chem., 1996, 100, 13892–13894.
11. K. Mishra, Z. Wu, C. Erker, K. Müllen and T. Basché, Chem. Sci., 2025, 16, 90–97.
12. C. Toninelli, I. Gerhardt, A. Clark, A. Reserbat-Plantey, S. Götzinger, Z. Ristanović, M. Colautti, P. Lombardi, K. Major and I. Deperasińska, et al., Nat. Mater., 2021, 20, 1615–1628.
13. A. Moradi, Z. Ristanović, M. Orrit, I. Deperasińska and B. Kozankiewicz, ChemPhysChem, 2019, 20, 55–61.
14. V. G. Plotnikov, Int. J. Quantum Chem., 1979, 16, 527–541.
15. R. R. Valiev, V. N. Cherepanov, G. V. Baryshnikov and D. Sundholm, Phys. Chem. Chem. Phys., 2018, 20, 6121–6133.
16. R. R. Valiev, V. N. Cherepanov, R. T. Nasibullin, D. Sundholm and T. Kurtén, Phys. Chem. Chem. Phys., 2019, 21, 18495–18500.
17. 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., 2020, 22, 22314–22323.
18. R. R. Valiev, R. T. Nasibullin, V. N. Cherepanov, A. Kurtsevich, D. Sundholm and T. Kurtén, Phys. Chem. Chem. Phys., 2021, 23, 6344–6348.
19. E. I. Adirovich, Nekotorye voprosy teorii lyuminestsentsii kristallov, Gostekhizdat, Moscow, 1951.
20. M. Bixon and J. Jortner, J. Chem. Phys., 1968, 48, 715–726.
21. V. G. Plotnikov and G. G. Konoplev, Fiz. Tverd. Tela, 1973, 15, 680–688.
22. V. G. Plotnikov and G. G. Konoplev, Zh. Eksp. Teor. Fiz., 1973, 65, 960–972.
23. V. G. Plotnikov and G. G. Konoplev, Institute for Theoretical Physics, Acad. Sci. USSR, 1974.
24. V. Y. Artyukhov, A. I. Galeeva, G. V. Maier and V. V. Ponomarev, Opt. Spectrosc., 1997, 83, 685–690.
25. M. Lax, J. Chem. Phys., 1952, 20, 1752–1760.
26. P. V. Yurenev, M. K. Kretov, A. V. Scherbinin and N. F. Stepanov, J. Phys. Chem. A, 2010, 114, 12804–12816.
27. M. K. Kretov, I. M. Iskandarova, B. V. Potapkin, A. V. Scherbinin, A. M. Srivastava and N. F. Stepanov, J. Lumin., 2012, 132, 2143–2150.
28. M. Kretov, A. Scherbinin and N. Stepanov, Russ. J. Phys. Chem. A, 2013, 87, 245–251.
29. A. Manian, R. A. Shaw, I. Lyskov, W. Wong and S. P. Russo, J. Chem. Phys., 2021, 155, 054108.
30. V. A. Kovarskii, Kinetika bezizluchatel'nykh protsessov, AN MSSR, Chisinau, 1968.
31. V. A. Kovarskii, Mnogokvantovye perekhody, Shtiintsa, Chisinau, 1974.
32. V. A. Kovarskii, N. F. Perel'man, I. S. Averbukh, S. A. Baranov and S. S. Todirashku, Neadiabaticheskie perekhody v sil'nom elektromagnitnom pole, Shtiintsa, Chisinau, 1980.
33. R. Valiev, B. Merzlikin, R. Nasibullin, A. Kurtzevitch, V. Cherepanov, R. Ramazanov, D. Sundholm and T. Kurtén, Phys. Chem. Chem. Phys., 2023, 25, 6406–6415.
34. V. L. Ermolaev, Opt. Spectrosc., 2016, 121, 567–584.
35. B. I. Makshantsev, Spectrosc. Lett., 1972, 5, 1–6.
36. B. I. Makshantsev, Spectrosc. Lett., 1972, 5, 7–12.
37. R. R. Valiev, Y. He, T. Weltzin, A. Zhu, D. Lee, E. Moore, A. Gee, G. Drozd and T. Kurtén, Phys. Chem. Chem. Phys., 2025, 27, 998–1007.
38. G. T. Drozd, T. Weltzin, S. Skiffington, D. Lee, R. Valiev, T. Kurtén, L. R. Madison, Y. He and L. Gargano, Environ. Sci.: Atmos., 2024, 4, 509–518.
39. B. Szymański, S. R. Sahoo, R. Valiev, O. Vakuliuk, P. Aski, K. N. Jarzembska, R. Kamiński, G. Baryshnikov, M. B. Teimouri and D. T. Gryko, New J. Chem., 2024, 48, 2416–2420.
40. B. Szymański, S. R. Sahoo, O. Vakuliuk, R. Valiev, R. Ramazanov, P. Aski, K. N. Jarzembska, R. Kamiński, M. B. Teimouri, G. Baryshnikov and D. T. Gryko, Chem. Sci., 2025, 16, 2170–2179.
41. K. Górski, D. Kusy, S. Ozaki, M. Banasiewicz, R. Valiev, S. R. Sahoo, K. Kamada, G. Baryshnikov and D. T. Gryko, J. Mater. Chem. C, 2024, 12, 1980–1987.
42. N. Ibrayev, E. Seliverstova, R. Valiev, A. Aymagambetova and D. Sundholm, Phys. Chem. Chem. Phys., 2024, 26, 25986–25993.
43. L. I. Valiulina, V. N. Cherepanov and K. Khoroshkin, Phys. Chem. Chem. Phys., 2024, 26, 22337–22345.
44. I. Sahalianov, R. R. Valiev, R. R. Ramazanov and G. Baryshnikov, J. Phys. Chem. A, 2024, 128, 5138–5145.
45. B. S. Merzlikin, V. N. Cherepanov, K. Khoroshkin and R. R. Valiev, Chem. Phys. Lett., 2024, 836, 141032.
46. R. Valiev, B. Merzlikin, R. Nasibullin, V. Cherepanov, D. Sundholm and T. Kurtén, Phys. Chem. Chem. Phys., 2024, 26, 4151–4158.
47. A. Manian, Z. Chen, H. Sullivan and S. Russo, ChemRxiv, 2025, preprint, chemrxiv/assets/orp/resource/item/67b7c567fa469535b9379533 DOI:10.26434/chemrxiv-2025-kp7jf.
48. R. Kubo and Y. Toyozawa, Prog. Theor. Phys., 1955, 13, 160–182.
49. S. I. Pekar, Research in Electron Theory of Crystals, Moscow, 1951.
50. M. A. Krivoglaz and S. I. Pekar, Proc. Inst. Phys., Acad. Sci. UkrSSR, 1953, 4, 37–70.
51. Y. Niu, Q. Peng and Z. Shuai, Sci. China, Ser. B: Chem., 2008, 51, 1153–1158.
52. Y. Niu, Q. Peng, C. Deng, X. Gao and Z. Shuai, J. Phys. Chem. A, 2010, 114, 7817–7831.
53. Q. Peng, Y. Niu, Q. Shi, X. Gao and Z. Shuai, J. Chem. Theory Comput., 2013, 9, 1132–1143.
54. A. Humeniuk, M. Bužančić, J. Hoche, J. Cerezo, R. Mitrić, F. Santoro and V. Bonačić-Koutecký, J. Chem. Phys., 2020, 152, 054107.
55. S. V. Tyablikov and V. A. Moskalenko, Proc. Steklov Inst. Math., 1961, 64, 267–283.
56. M. Bracker, C. M. Marian and M. Kleinschmidt, J. Chem. Phys., 2021, 155, 014102.
57. M. Mankowski and M. Moshkov, Dynamic Programming Multi-Objective Combinatorial Optimization, Springer, 2021, vol. 331, pp. 1–210.
58. PubChem Compound Summary for CID 135991, Available at https://pubchem.ncbi.nlm.nih.gov/compound/135991#section=3D-Conformer, Accessed 2024.
59. A. D. Becke, J. Chem. Phys., 1993, 98, 5648–5652.
60. C. Lee, W. Yang and R. G. Parr, Phys. Rev. B: Condens. Matter Mater. Phys., 1988, 37, 785–789.
61. F. Weigend and R. Ahlrichs, Phys. Chem. Chem. Phys., 2005, 7, 3297–3305.
62. E. Caldeweyher, C. Bannwarth and S. Grimme, J. Chem. Phys., 2017, 147, 034112.
63. Y. J. Franzke, C. Holzer, J. H. Andersen, T. Begušić, F. Bruder, S. Coriani, F. Della Sala, E. Fabiano, D. A. Fedotov, S. Fürst, S. Gillhuber, R. Grotjahn, M. Kaupp, M. Kehry, M. Krstić, F. Mack, S. Majumdar, B. D. Nguyen, S. M. Parker, F. Pauly, A. Pausch, E. Perlt, G. S. Phun, A. Rajabi, D. Rappoport, B. Samal, T. Schrader, M. Sharma, E. Tapavicza, R. S. Treß, V. Voora, A. Wodyński, J. M. Yu, B. Zerulla, F. Furche, C. Hättig, M. Sierka, D. P. Tew and F. Weigend, J. Chem. Theory Comput., 2023, 19, 6859–6890.
64. TURBOMOLE V7.8 2023, a development of University of Karlsruhe and Forschungszentrum Karlsruhe GmbH, 1989–2007, TURBOMOLE GmbH, since 2007, available from https://www.turbomole.org.
65. 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, D. Williams-Young, 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. J. A. Montgomery, 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.
66. A. Granovsky, Firefly version 8, Available at https://classic.chem.msu.su/gran/firefly/index.htm, Accessed 2024.
67. I. Purvis, D. George and R. J. Bartlett, J. Chem. Phys., 1982, 76, 1910–1918.
68. K. Raghavachari, G. W. Trucks, J. A. Pople and M. Head-Gordon, Chem. Phys. Lett., 1989, 157, 479–483.
69. C. Hättig, D. P. Tew and A. Köhn, J. Chem. Phys., 2010, 132, 231102.
70. R. Send, V. R. I. Kaila and D. Sundholm, J. Chem. Phys., 2011, 134, 214114.
71. M. Bendikov, H. M. Duong, K. Starkey, K. N. Houk, E. A. Carter and F. Wudl, J. Am. Chem. Soc., 2004, 126, 7416–7417.
72. H. S. Yu, X. He, S. L. Li and D. D. G. Truhlar, Chem. Sci., 2016, 7, 5032–5051.
73. R. A. Marcus, J. Chem. Phys., 1956, 24, 966–978.

---

Footnote

† Electronic supplementary information (ESI) available: The Cartesian coordinates of the molecular structure of the studied molecule. Energy-dependent rate constants of internal conversion calculated at different levels of approximation. See DOI: https://doi.org/10.1039/d5cp01603g

---