Investigating XC-functionals towards describing experimentally relevant excited-state properties of NIR-BODIPY derivatives

Abstract

The predictive and analytical power of time-dependent density functional theory (TD-DFT) has been instrumental in the design and mechanistic understanding of numerous organic chromophores. Yet, the widely popular boron-dipyrromethene (BODIPY) dye class suffers from notorious TD-DFT accuracy issues, undermining the serviceability of the technique. Highly correlated wave function approaches are much better at reproducing photophysical properties but become computationally unviable when making the push towards larger near-infrared (NIR) active structures. In an effort to find the protocol most capable of helping experimentalists design and analyze novel NIR BODIPYs, we have benchmarked 11 global or range-separated hybrid exchange–correlation functionals (XCFs) with different amounts of Hartree–Fock exchange. By relating both transition energies and oscillator strengths, first through a set of resolution-of-the-identity second-order coupled cluster (riCC2) calculations and then directly to experimental data, it is revealed that M06-2X and M06-HF behave most consistently for singlet and triplet excitations. To optimize accuracy across states, we recommend a hybrid approach where singlets are obtained through full TD-DFT and triplets are treated using the Tamm–Dancoff approximation.

---

1 Introduction

Due to their strong absorption, high quantum yields, narrow absorption/emission bands, and rich synthetic library, 4,4-difluoro-4-bora-3a,4a-diaza-s-indacene (BODIPY) chromophores and their aza-analogues have been able to penetrate and proliferate across a wide range of fields.^1–3 Several applications, including photovoltaics,⁴ fluorescence bio-imaging,⁵ chemosensing,^6,7 and photodynamic therapy,^8–10 have invoked the development of BODIPY dyes active in the near-infrared (NIR), more particularly within their respective biological (650–1350 nm) and phototherapeutic (600–800 nm) windows. The most common modifications to push optical properties to higher wavelengths include substitution of the meso-carbon with nitrogen, placing aryl groups at the 1,3,5,7-positions (Fig. 1), fusing aryl groups to the pyrrole moiety, or fusing multiple BODIPY units to each other.^1,11–13 The latter two strategies often work best in retaining fluorescence quantum yields, as little to no extra rotational/vibrational freedom is introduced. Besides tuning absorption wavelengths, a significant portion of BODIPY research is dedicated to implementing additional functionalities, many relying on mechanisms involving singlet and triplet excited states having localized or charge transfer (CT) character. Typical examples encompass the spin–orbit charge transfer intersystem crossing (SOCT-ISC) used in photodynamic therapy, where a locally excited singlet (¹LE) is converted to its triplet counterpart (³LE) via a charge transfer state in an orthogonally arranged donor–acceptor design,¹⁴ as well as triplet upconversion phenomena like thermally activated delayed fluorescence (TADF)¹⁵ and triplet–triplet annihilation¹⁶ that are leveraged to increase the quantum yield of fluorescence. These fields and others have taken advantage of the important contribution of Christel Marian, who is celebrated by this special issue.^17,18

Fig. 1 Structure of the BODIPY core.

The challenge of understanding processes like these lies in the full characterization of each intermediate state and the mechanisms by which they interact. Here, computational chemistry has become an indispensable component of material research by providing a means to determine a plethora of useful molecular properties. Various levels of theory are available ranging from highly correlated wavefunction methods to the widely used time-dependent density functional theory (TD-DFT), each weighing computational cost against chemical accuracy. The choice of method is, however, not straightforward for NIR-BODIPY researchers as the larger size of NIR-absorbing compounds limits them to cheaper methods like adiabatic TD-DFT with the vertical approximation, while charge transfer excitations and the cyanine-like class, to which BODIPYs belong, have well-known issues in TD-DFT. Range-separated exchange–correlation functionals (XCF) have dealt somewhat with the former problem, but an accurate description of the lowest BODIPY singlet (S₁) and triplet (T₁) excited states relies heavily on the inclusion of double excitations which are not accounted for in the adiabatic approximation of TD-DFT, resulting in typical overestimations of 0.3 eV or more for S₁ and underestimations up to 0.9 eV for T₁.^19–21 Applying the Tamm–Dancoff approximation (TDA) somewhat attenuates T₁-errors, especially for functionals with high amounts of Hartree–Fock (HF) exchange. Unfortunately, singlet–triplet gaps are not improved as S₁ excitation energies experience a systematic upshift as well.²¹ Therefore, Boulanger et al.²² and Chibani et al.²³ proposed joint TD-DFT/ab initio approaches where vertical excitations are calculated using the Bethe–Salpeter (BSE) formalism or the scaled opposite spin configuration interaction singles perturbative doubles (SOS-CIS(D)). Related works by Momeni and Brown^20,24 and Feldt and Brown²⁵ also demonstrated good accuracy of other high-level methods such as CASPT2, SAC-CI, LCC2*, and DLPNO-STEOM-CCSD. Using the compound set of Momeni, Helal et al.²⁶ benchmarked the increasingly popular double-hybrid (DH) XC-functionals alongside a wide range of TD-DFT XCFs and found two empirical dispersion-corrected, spin-component-scaled, DH XC-functionals with mean absolute errors (MAEs) comparable to those of the mentioned wavefunction methods. Another work by Toffoli et al.²⁷ presented time-independent Δ-self-consistent-field DFT as an inexpensive and accurate way to screen BODIPYs excited-state properties after finding out that it aided in the description of other boron-containing compounds with double-excitation-related TD-DFT issues.^28,29

Although these and past^30–33 studies have uncovered the challenges associated with smaller BODIPYs and suggested several remedies, the experimentalist remains unfulfilled in his pursuit towards a good all-round approach to efficiently screen NIR BODIPY structures for the properties that are relevant to their research. Firstly, there are three promising fused NIR BODIPY derivatives, i.e., (aza-)BOIMPYs,^34–36 pyrrolopyrrole (aza-)BODIPYs,^37,38 and α,α-linked BisBODIPYs,³⁹ that have yet to be investigated in any benchmark, and secondly, for a method to be sound in its description of photophysical phenomena, it should go beyond the lowest vertical excitation energies, determining the characteristics of several local and charge transfer states with accuracy. Towards the latter end, Sbai and Guthmuller found MN15 and SOS-PBE-QIDH to be balanced in their treatment of excitations belonging to three BODIPY photocatalysts.⁴⁰ Additionally, the inclusion of oscillator strengths (f) in a benchmark against experimental data provides a broader basis on which to evaluate method performance.⁴¹

This study aims to discover which TD-DFT method provides the most reliable results on NIR BODIPY-type chromophores whilst being cheap enough for utilization in the screening of a typical candidate batch encountered in BODIPY research. 7 Literature compounds^{34,35,39,42–45} (Fig. 2) were selected to represent a set of highly fluorescent BODIPY structures spanning the green to NIR range, one of which (compound 2) is known to be a performant SOCT-ISC photosensitizer. As only the experimental S₀ → S₁ excitation energies and oscillator strengths are available in the literature, our investigation starts off with a S₁-focused assessment of the resolution-of-the-identity second-order approximated coupled cluster (riCC2) method with neither, either, or both the conductor-like screening model⁴⁶ (COSMO) and spin-component scaling^47,48 (SCS) enabled. The best-performing combination then serves as the reference to test a total of 11 XCFs covering pure, global hybrid, and long-range corrected (LC) hybrid XCFs, thereby considering the first two singlet (S) and first three triplet (T) states alongside the S₁–T₁ and S₁–T₂ energy gaps. Both full TD-DFT and its Tamm–Dancoff approximated version are taken into consideration as well as the polarizable continuum model⁴⁹ (PCM) when comparing with the COSMO results. The analysis is additionally supported by a suite of excitation-induced CT characterization tools, i.e., electron density difference (EDD) plots, CT amplitudes, and two different CT distance metrics.^50–53

Fig. 2 BODIPY derivative structures considered in the benchmark.

In summary, this article starts off with a technical overview of the employed methodologies, followed by an evaluation of the riCC2 methods in their ability to estimate the experimental results. The ensuing main discussion covers the benchmarking of various XC-functionals in relation to riCC2 and experiments, focusing on the lowest singlet and triplet states, and the effects of the Tamm–Dancoff approximation.

2 Computational methods

2.1 Vertical excitation energies and oscillator strengths

2.1.1 Experimental results. Experimental oscillator strengths have been extracted from the literature via a procedure inspired by the work of Tarleton et al.⁵⁴ As suggested by the authors and others,⁵⁵ a solvent refractive index correction factor (n) was applied for f^exp to improve correlation coefficients (R²) and regression slopes. UV-vis absorption spectra were reconstructed from images using WebPlotDigitizer⁵⁶ and the total areas of the presumed S₁-peaks were converted to oscillator strengths using eqn (1).

2.1.2 Quantum chemistry for excitation energies and oscillator strengths. All DFT and TD-DFT(/TDA) calculations were carried out on the Gaussian 16 platform.⁵⁷ Compounds 5 and 6 had their octoxyl moieties reduced to methoxyls to save on computational costs. Ground-state equilibrium geometries were optimized under a tight convergence criterion using the ultrafine pruned (99, 590) grid, the M06-2X functional, and the 6-311G(d) basis set. Their true minimal position is confirmed by the absence of imaginary frequencies in the subsequent vibrational frequency calculations, save for a small methyl group rotation on compound 7. Given the reluctance of this compound to stabilize any further and the negligible impact of the vibration, the geometry was taken as is. Then, the 30 lowest singlet and triplet vertical excitation properties were computed for each of the 11 XC-functionals (Table 1) by means of full TD-DFT and with the Tamm–Dancoff approximation. The latter was employed owing to its recognized ability to improve the accuracy on the triplet excitation energies, alleviating triplet instability problems. Indeed, Tozer and co-workers^58,59 have shown that, when combined with range-separated hybrids, TDA yields accurate results for local (in contrast to charge transfer) triplet excitations.

Table 1 Categorized overview of the investigated XC-functionals, detailing the amounts of HF exchange. For the long-range corrected functionals, the percentages at interelectronic distances 0 and ∞ are reported, alongside values for the range-separating parameter ω (Bohr⁻1)

	XC-functionals	% HF exchange
GGA	BLYP^61,62
	PBE^63
Global hybrid GGA	B3LYP^62,64–66	20
	PBE0^67	25
Global hybrid meta-GGA	M06^68	27
	M06-2X^68	54
	M06-HF^69	100
LC hybrid GGA	LC-ωPBE17^70,71	0–100 (ω = 0.17)
	LC-BLYP17^72,73	0–100 (ω = 0.17)
	LC-BLYP20^72,73	0–100 (ω = 0.20)
	LC-BLYP33^72,73	0–100 (ω = 0.33)

---

The BLYP/PBE/B3LYP/PBE0 XCFs were included as popular references for the lower rungs of Jacob's ladder, including GGA's and global hybrids. M06-2X was included as a popular hybrid meta-GGA XCF, which has demonstrated good accuracy for BODIPYs and CT excitations.^32,60 The M06 and M06-HF variants are there to observe the effects of the HF-exchange amount. Then, the long-range corrected (LC) XCFs, or range-separated hybrids, were added as these should be able to deal with CT excited states. Again, varying amounts of HF-exchange were tested by tuning the range-separating parameter ω. The 6-311G* basis set was maintained throughout despite other works calling out for the use of a larger basis set in TD-DFT.^30,31 Our additional calculations at the M06-2X/6-311+G(2d,p) level gave a mean absolute deviation of only 0.02 eV with isolated errors never exceeding 0.08 eV. All steps were performed both in gas-phase and under a non-equilibrium linear-response PCM regime using the corresponding experimental solvents.

A set of four riCC2 methods were compared with the experimental results on the basis of the first singlet excited state with the intent of serving as an intermediary basis for the benchmarking of higher TD-DFT singlet and triplet excitations. Preceding works by Brown et al. have employed similar strategies and demonstrated a good linear correlation between the riCC2 model and experiments.^25,74 Vertical excitation properties were computed in the TURBOMOLE v7.5.1 package,⁷⁵ using both the regular and spin-component-scaled version of riCC2, with and without the inclusion of COSMO solvent effects. To avoid geometry-related dependencies, the DFT-optimized structures were assigned to the gas-phase and solvated riCC2 calculations accordingly. A subset of calculations was carried out with both the double and triple-ζ augmented Dunning basis sets: aug-cc-pVDZ and aug-cc-pVTZ. The observed differences in vertical excitation energies were never more than 0.03 eV; hence, we opted to use aug-cc-pVDZ as the regular and auxiliary basis set. Similar observations were made in previous investigations.^24,25,76

2.2 Charge transfer characteristics from quantum chemistry calculations

Charge transfer properties were computed from TD-DFT(/TDA) densities generated by the cubegen utility of Gaussian16. CT-amplitudes (q_CT) and CT-distances (D_CT) were determined according to a method developed by Jacquemin et al.^50,51

Δρ(r⃑) = ρ^ES(r⃑) − ρ^GS(r⃑) (2)

Subtracting the ground- and excited-state densities (eqn (2)) yields a map of the total density change upon excitation that can be partioned into the electron (eqn (3)) and hole (eqn (4)) distributions.The q_CT is obtained by integrating either distribution (eqn (5)), whilst D_CT is acquired by measuring the distance between the barycenters of both (eqn (6)).

A well-known limitation of the D_CT index is that it cannot properly describe the degree of locality for excitations in symmetrical molecules such as compounds 5 and 6. To remove the influence of system shape on the quantity, three very similar indices were recently developed by Wang et al.,⁵² ourselves,⁵³ and Lieberherr et al.,⁷⁷ all making use of an optimal transport scheme to determine the average distance across which electronic density is displaced during an excitation.⁷⁸ The so-called earth mover's CT distance (^EMD_CT, eqn (7)) index, as described in our work, focuses on facilitating the demanding computation for larger systems by using the CHelpG charge model,⁷⁹ making it the best option for use in this work.

3 Results and discussion

3.1 Coupled-cluster methods versus experiment

Extracting the first singlet transition energies at maximum absorption (ΔE^exp) and oscillator strengths (f^exp) from experimental data was straightforward as all compounds exhibited a distinct S₁-peak with negligible interference from neighbouring excitations (Table 2). The riCC2 method, which already shows good linear correlation with the experimental data and, in accordance with the literature,⁸⁰ the implementation of spin-component scaling creates even better R² values for both ΔE and f at the cost of consistent ΔE overestimation (Fig. 3). Likewise, the introduction of solvent effects via COSMO was accompanied by a consistent oscillator strength increase and a reduced overestimation of the excitation energies (Table 2). SCS-riCC2 is thus naturally chosen as the basis for the ensuing TD-DFT benchmark.

Table 2 SCS-riCC2 first singlet vertical excitation energies and oscillator strengths with COSMO solvation alongside the experimental data

	ΔE^exp (eV)	ΔE^calc (eV)	f ^exp	f ^calc	Solvent	n ^81	n × f^exp
1	2.470	2.595	0.432	0.605	CH2Cl2	1.4242	0.615
2	2.450	2.590	0.471	0.601	CH2Cl2	1.4242	0.670
3	2.098	2.311	0.574	0.848	CH2Cl2	1.4242	0.818
4	2.073	2.312	0.475	0.763	CH2Cl2	1.4242	0.677
5	1.943	2.134	0.590	0.992	CHCl3	1.4459	0.853
6	1.813	1.929	0.683	1.013	CHCl3	1.4459	0.988
7	1.651	1.832	0.579	1.006	Toluene	1.4961	0.867

---

Fig. 3 S₁ vertical excitation energies (top) and oscillator strengths (bottom) of different riCC2 methods plotted against their respective λ_max and peak surface derived experimental counterparts. The dashed lines represent perfect coincidence of the theoretically and experimentally obtained values.

3.2 TD-DFT and TD-DFT/TDA versus SCS-riCC2

As a preliminary step, COSMO SCS-riCC2 transitions were correlated with PCM TD-DFT(/TDA) results based on the similarity in the overall character, prioritizing resemblance in nature over excitation energies. Indeed, DFT and HF molecular orbitals (MO) of the given character could often rank differently in terms of energy; therefore, their natures were compared based on the MO topology. For some cases, this cross-method assignment was complicated by the presence of two TD-DFT(/TDA) excited states whose characters appeared to be a mixture of two SCS-riCC2 excitations. At every one of these occurrences, the second of the two SCS-riCC2 states was either S₃ or T₄ (states not considered in the benchmark) with the former always having a negligible f compared to S₂. To prevent any oscillator strengths of the second SCS-riCC2 singlets from “diluting” across two TD-DFT(/TDA) states, we opted to compare them with the sum of the TD-DFT(/TDA) oscillator strengths and the thereon weighted average excitation energies. Whenever two TD-DFT(/TDA) excitation energies differed significantly (>0.2 eV), only one state was selected based on the trends in the overall character across XC-functionals. This process has been repeated for all other riCC2 methods, the results of which are summarized in the ESI.‡ Regression data for the subsequent correlations can likewise be found there. The mean absolute errors and determination coefficients of the discussed excitations and singlet–triplet gaps are displayed in Table 3.

Table 3 Mean absolute errors (MAE) and determination coefficients (R²) of PCM TD-DFT(/TDA) versus COSMO SCS-riCC2 regressions. TD-DFT results are displayed in the top half, with the Tamm–Dancoff approximated analogs rendered in the bottom half

XCF	ΔES0–S1 (eV)		f S0–S1		ΔES0–S2 (eV)		f S0–S2		ΔES0–T1 (eV)		ΔES0–T2 (eV)		ΔES0–T3 (eV)		ΔES1–T1 (eV)		ΔES1–T2 (eV)
	MAE	R ^2	MAE	R ^2	MAE	R ^2	MAE	R ^2	MAE	R ^2	MAE	R ^2	MAE	R ^2	MAE	R ^2	MAE	R ^2
BLYP	0.18	0.89	0.27	0.51	0.89	0.56	0.15	0.02	0.38	0.97	0.79	0.46	0.93	0.43	0.38	0.37	0.79	0.59
PBE	0.18	0.88	0.30	0.52	0.88	0.56	0.15	0.02	0.38	0.97	0.79	0.46	0.94	0.42	0.39	0.38	0.80	0.61
B3LYP	0.17	0.95	0.11	0.72	0.43	0.55	0.12	0.12	0.48	0.98	0.60	0.93	0.61	0.61	0.61	0.50	0.73	0.88
PBE0	0.18	0.95	0.09	0.77	0.31	0.55	0.12	0.22	0.55	0.97	0.60	0.94	0.57	0.67	0.72	0.54	0.76	0.91
M06	0.14	0.96	0.09	0.80	0.32	0.57	0.11	0.29	0.52	0.94	0.58	0.93	0.52	0.67	0.64	0.57	0.70	0.94
M06-2X	0.19	0.98	0.03	0.96	0.17	0.66	0.02	1.00	0.51	0.92	0.28	0.98	0.19	0.96	0.70	0.70	0.47	0.99
M06-HF	0.15	0.98	0.02	0.98	0.44	0.90	0.04	0.99	0.68	0.64	0.16	0.96	0.17	0.78	0.83	0.79	0.23	0.77
LC-BLYP17	0.14	0.94	0.09	0.88	0.25	0.48	0.04	1.00	0.51	0.99	0.50	1.00	0.53	0.98	0.62	0.59	0.62	0.96
LC-ωPBE17	0.14	0.94	0.09	0.89	0.24	0.49	0.04	1.00	0.57	0.99	0.55	1.00	0.56	0.98	0.70	0.59	0.67	0.96
LC-BLYP20	0.14	0.95	0.07	0.93	0.19	0.58	0.02	0.99	0.57	0.98	0.47	1.00	0.44	0.98	0.70	0.65	0.59	0.98
LC-BLYP33	0.17	0.98	0.02	1.00	0.26	0.90	0.04	0.99	0.97	0.71	0.42	0.98	0.22	0.71	1.14	0.65	0.60	0.96
BLYP	0.47	0.87	0.26	0.07	0.88	0.50	0.19	0.00	0.34	0.97	0.71	0.39	0.90	0.41	0.74	0.27	1.11	0.53
PBE	0.42	0.83	0.36	0.06	0.90	0.42	0.17	0.00	0.35	0.97	0.71	0.39	0.91	0.41	0.70	0.16	1.07	0.59
B3LYP	0.42	0.92	0.24	0.62	0.38	0.53	0.15	0.03	0.35	0.99	0.51	0.90	0.55	0.55	0.76	0.49	0.93	0.76
PBE0	0.44	0.94	0.30	0.70	0.28	0.53	0.14	0.06	0.36	0.98	0.47	0.94	0.48	0.59	0.81	0.51	0.91	0.81
M06	0.38	0.95	0.28	0.71	0.28	0.55	0.14	0.09	0.34	0.96	0.45	0.94	0.43	0.60	0.72	0.54	0.84	0.86
M06-2X	0.43	0.98	0.32	0.91	0.27	0.48	0.03	1.00	0.29	0.96	0.18	0.96	0.12	0.95	0.72	0.66	0.61	0.97
M06-HF	0.39	1.00	0.23	0.84	0.61	0.81	0.12	0.98	0.24	0.82	0.19	0.93	0.29	0.82	0.64	0.80	0.28	0.80
LC-BLYP17	0.36	0.94	0.22	0.91	0.25	0.31	0.04	0.99	0.35	1.00	0.44	0.98	0.48	0.95	0.71	0.59	0.80	0.87
LC-ωPBE17	0.38	0.94	0.24	0.90	0.24	0.32	0.04	0.99	0.38	1.00	0.45	0.99	0.50	0.97	0.75	0.59	0.83	0.87
LC-BLYP20	0.36	0.96	0.23	0.93	0.23	0.42	0.01	1.00	0.36	0.99	0.38	0.98	0.38	0.99	0.72	0.63	0.74	0.93
LC-BLYP33	0.39	0.99	0.21	0.88	0.39	0.79	0.09	0.99	0.37	0.91	0.20	0.98	0.10	0.76	0.76	0.75	0.59	0.97

---

For each compound, COSMO SCS-riCC2 mainly describes the S₀ → S₁ transition as a local single electron excitation from the highest occupied to the lowest unoccupied molecular orbital (HOMO → LUMO) with doubles contributions ranging from 11.6 to 13.1 percent. As seen in Fig. 4, XCF trends for the first singlet excitation energy differ between smaller (1–4) and larger (5–7) BODIPY derivatives. Whilst most XC-functionals stay within ∼0.1 eV intervals for the smaller four, the larger three clearly show increased ΔE_S0–S1 deviations whenever more HF exchange is introduced. Following this, the global and range-separated hybrids possessing higher amounts of HF exchange (M06-2X, M06-HF & LC-BLYP33) prove themselves as most consistent since the S₁ excitation energies are generally overestimated for compounds 1–4. In terms of oscillator strengths, M06-2X, M06-HF, and LC-BLYP33 also outperform the rest, as for all compounds besides 7, proportional underestimations are observed when less HF exchange is present. Applying the Tamm–Dancoff approximation results in a rather uneven inflation of the f values, compromising the regression qualities of the GGA functionals and strongly increasing MAEs for the others. Changes in excitation energies under the TD-DFT/TDA regime come down to systematic upshifts of about 0.23 eV for the higher-level XC-functionals, leading to regressions with roughly identical R² values but considerably worse overestimations. The lower-level XCFs (BLYP, B3LYP, PBE & PBE0) experience the same effects but to a less consistent degree.

Fig. 4 Correlation plots of PCM TD-DFT(/TDA) against COSMO SCS-riCC2 for the vertical S₀ → S₁ excitation energies (top-left) and oscillator strengths (top-right), the vertical S₀ → S₂ excitation energies (middle-left) and oscillator strengths (middle-right), the vertical S₀ → T₁ excitation energies (bottom-left), and the vertical S₀ → T₂ excitation energies (bottom-right). Each graph displays full TD-DFT on the left and TD-DFT/TDA on the right. The dashed lines represent perfect coincidence of PCM TD-DFT(/TDA) with COSMO SCS-riCC2.

Despite exhibiting higher 13.4–15.5% doubles contributions in SCS-riCC2, the TD-DFT(/TDA) second singlet vertical excitation energies express more uniform trends across the compound spectrum (Fig. 4). Starting from severe underestimations at the BLYP (MAE: 0.89 eV) level, the errors diminish as the amount of HF exchange rises, eventually transitioning over to an overestimation when passing beyond the mid-HF range around M06-2X & LC-BLYP20. The S₂ excitations of compounds 2 at 3.37 eV and 6 at 3.04 eV are notably more underestimated by the pure XC-functionals (1.95 eV & 1.85 eV for PBE and 1.93 & 1.86 for BLYP). Likely caused by the stronger CT-character (^EMD_CT: 3.14 Å and 2.86 Å), their negative deviations are noticably improved with the addition of HF exchange. Meanwhile, the implementation of TDA does not help much in this regard, only shifting some results and without positively influencing the MAE or R² values. On the oscillator strength side, most of the f_S0–S2 values fall below 0.1, providing too small of a basis to make any statements here regarding method accuracy. Of the two excitations with noteworthy f values, compound 4 demonstrates good consistency, whilst compound 7 is strongly dependent on the presence of HF exchange. Although in the latter case, each XC-functional exhibits some degree of mixing among S₂ and S₃, no real relationship is found between the character and oscillator strength.

The SCS-riCC2 lowest vertical triplet excitations (S₀ → T₁) all exhibit local character related to those of the first singlet transitions, also with a dominant HOMO → LUMO single-electron promotion. The weight of double excitations is on average 3% lower for triplets, raising expectations of improved linear response TD-DFT performance. Regardless, markedly large underestimations are observed with MAEs ranging between 0.38 eV (BLYP) and 0.97 eV (LC-BLYP33). The rampant decreases of first triplet excitation energies following the increasing share of HF exchange is a clear sign of the triplet instability problem, which is known to be stronger in cyanine-type molecules.^58,59,82 This issue is frequently remedied by the TDA approach and accordingly we observe a nullification of most HF exchange related errors. Even so, almost every transition energy stands about 0.3 eV below their reference value. Other than that, the GGA functionals and the hybrids with low to median amounts of HF exchange hold good linear regressions across pure TD-DFT and TD-DFT/TDA with slopes and R² values close to 1 (Fig. 4).

The second (Fig. 4) and third (Fig. 5) triplet states behave similar to S₂ with the same upward shift following rising levels of HF exchange. Two T₂-excitations deviate from this general trend: the first at 2.46 eV, almost completely localized on the anthracene moiety of compound 2, and the second at 2.19 eV belonging to compound 7, Table 5 (COSMO SCS-riCC2 values). Both have a relatively low lying position and show signs of triplet instability. On the T₃ side, two molecules (5 at 3.31 eV & 6 at 3.17 eV) are found to display enhanced sensitivity to HF exchange on account of their partial charge transfer nature (^EMD_CT: 1.76 Å and 2.77 Å, as determined with M06-2X). Here, an optimal balance is achieved around M06-2X and LC-BLYP20 with determination coefficients falling as soon as too much or too little HF exchange is present. The mean absolute errors are minimal at the M06-HF level given it is the sole XC-functional that does not completely underestimate the excitation energies. When the TDA approximation is active, it has to give way to LC-BLYP33 and M06-2X because the slight upshift lessens their underestimation. Nonetheless, the XC-functionals with a high degree of HF exchange lack reliability as their slopes are amongst the worst of the hybrid XCFs. Considering everything, M06-2X, LC-ωPBE17, LC-BLYP17, and LC-BLYP20 possess the most amenable T₂/T₃ linear regression in both regimes (TD-DFT(/TDA)) even with the irregularities that arise for the T₂ excitation energies of compounds 2 and 7 (Fig. 4).

Fig. 5 Correlation plots of PCM TD-DFT(/TDA) against COSMO SCS-riCC2 for the vertical S₀ → T₃ excitation energies (top-left), the energy gaps from S₁ to T₁ (top-right) and to T₂ (bottom-left), and the TD-DFT first singlet to TD-DFT/TDA (first and second) triplet energy gaps (bottom-right). The first three graphs display full TD-DFT on the left and TD-DFT/TDA on the right. The dashed lines represent perfect coincidence of PCM TD-DFT(/TDA) with COSMO SCS-riCC2.

Moving on to the TD-DFT S₁ → T₁ plots (Fig. 5), a coincidental dependency on HF exchange is observed stemming from the opposing trends of the smaller and larger compounds in the first singlet and triplet results. Where vertical S₁ excitation energies ascend in relation to HF exchange for 5–7 and remain more condensed for 1–4, the T₁ excitations descend in the latter group with smaller changes for the former. The outcome hereof is MAEs starting from 0.38 eV at the pure GGA BLYP level and jumping up to 1.14 eV for the LC-BLYP33 case. Furthermore, determination coefficients ΔE_S1–T1 are poor across the board despite the individual ΔE_S0–S1 and ΔE_S0–T1 results showing high R² values for most XC-functionals. This is likewise a consequence of an unfortunate meeting of good singlet estimates with bad triplet estimates and vice versa. Under the TD-DFT/TDA formalism, MAEs are evened out at about 0.7 eV following the attenuation of HF-dependent T₁-errors and amplification of pure XCF S₁-errors. Regression qualities stay poor with only M06-HF and LC-BLYP33 experiencing a slight improvement. Better results are obtained for ΔE_S1–T2 where higher level XC-functionals (M06, M06-2X, LC-ω PBE17, LC-BLYP17, LC-BLYP20, and LC-BLYP33) yield good linear correlations, thanks to their consistency in both excited states. The errors are again mainly affected by the shift in HF exchange of the triplet state that now gives M06-HF the lowest MAE (0.23 eV). The Tamm–Dancoff approximation generally inflates the excitation energies in S₁ to a larger extent than for T₂, causing all XCFs except LC-BLYP33 to lose accuracy.

ISC and reversed ISC rates are often estimated as a function of spin–orbit coupling (SOC) constants and singlet–triplet energy gaps. The second item is potentially problematic in predictive or mechanistic BODIPY studies as considerable errors on ΔE_S1–T1 are brought about by the compounded over- and underestimations of the first singlet and triplet states. However, the general depth of the first BODIPY triplet makes the (local) S₁ → T₁ transition an unlikely participant in ISC. More probable pathways are those via the less problematic second and third triplet states, rendering the TD-DFT method still viable for computational investigations in an experimental context.

Also, following the observed tendency of TDA to alleviate triplet related issues whilst being detrimental in the singlet domain, it becomes sensible to compare pure TD-DFT singlets and TD-DFT/TDA triplets in the hope of obtaining more reliable energy gaps (Table 4). As expected, this approach outperforms the non-hybrid methods with lower MAEs across the board and slope values closer to one. The smallest improvements belong to the pure GGA functionals (BLYP & PBE), which is primarily a consequence of their triplet results being quite insenstive to TDA. In contrast, the high HF exchange segment (M06-HF and LC-BLYP33) benefits most effectively from the combination of stable TD-DFT singlet energies and TDA-corrected triplet energies. Other than that, no big changes are observed on the side of determination coefficients, keeping M06-HF & LC-BLYP33 as the best options for ΔE_S1–T1 and M06-2X & LC-BLYP20 for ΔE_S1–T2.

Table 4 Regression data of PCM TD-DFT singlet to TD-DFT/TDA triplet gaps versus COSMO SCS-riCC2. MAE: mean average error, R²: determination coefficient, M: slope, and B: y-intercept

XCF	ΔES1–T1 (eV)				ΔES1–T2 (eV)
	MAE	R ^2	M	B	MAE	R ^2	M	B
BLYP	0.35	0.36	0.94	−0.38	0.72	0.59	0.69	−0.52
PBE	0.36	0.37	0.96	−0.38	0.72	0.59	0.70	−0.52
B3LYP	0.48	0.46	1.05	−0.45	0.65	0.86	0.94	−0.60
PBE0	0.53	0.49	1.10	−0.47	0.64	0.90	1.01	−0.64
M06	0.46	0.52	1.06	−0.43	0.57	0.93	1.05	−0.60
M06-2X	0.48	0.66	1.10	−0.42	0.37	0.98	1.03	−0.39
M06-HF	0.40	0.83	1.10	−0.34	0.19	0.74	1.02	−0.05
LC-BLYP17	0.47	0.55	1.00	−0.46	0.55	0.94	0.99	−0.54
LC-ωPBE17	0.50	0.55	1.04	−0.48	0.58	0.94	1.02	−0.59
LC-BLYP20	0.49	0.60	1.03	−0.47	0.50	0.97	1.03	−0.52
LC-BLYP33	0.55	0.76	1.12	−0.48	0.37	0.94	1.15	−0.47

---

3.3 TD-DFT and TD-DFT/TDA versus experimental data

When replicating the SCS-riCC2 results under a solvated environment, the mid- or high-HF exchange XCFs are continuously among the top performers, yet, in order to verify their effectiveness in aiding experimental research, a final comparison with the true S₁ is necessary. Upon plotting the TD-DFT vertical excitation energies against their experimental counterparts (Fig. 6), it is revealed that the better XC-functionals practically correlate as well as the riCC2 methods. The determination coefficients are very favorable and LC-BLYP33, in particular has a really good regression slope. Overestimations of the vertical excitation energies with respect to the maximum absorption energies generally fall within a rather high 0.3–0.4 eV range, but this is to be expected in a vertically approximated scheme. For the rigid pigment violet 19 molecule, TD-DFT calculations predict that the vertical excitation energy of the dominant excitation is 0.18 eV larger than the excitation energy corresponding to the maximum absorption when accounting for vibronic coupling.⁸³ For a more flexible cyanine, this difference amounts to 0.26 eV,⁸⁴ while for the chromophores of the green fluorescent protein, these differences range between 0.12 and 0.17 eV as a function of its protonation state.⁸⁵ Finding a rationale for these differences as a function of the molecular structure, without explicitly assessing the excited-state geometry, therefore does not appear straightforward. Lastly, the oscillator strength regressions are found to be better than those on the riCC2 side with lower MAEs and slopes closer to one.

Fig. 6 S₁ vertical excitation energies (top) and oscillator strengths (bottom) of different TD-DFT methods plotted against their respective λ_max and peak surface derived experimental counterparts. The dashed lines represent perfect coincidence of the theoretically and experimentally obtained values.

Table 5 TD-DFT (M06-2X/6-311G*) vertical excitation energies, oscillator strengths, density-based CT-amplitudes and CT-dipoles, and CHelpG-based earth mover's charge transfer distances. Transitions are in the original, non-riCC2-correlated order

Compound	State	ΔE (eV)	f	q CT	^EM D CT	μ CT
1	S1	2.869	0.605	0.36	1.40	0.97
	S2	3.812	0.000	1.16	2.86	14.50
	T1	1.551	—	0.54	1.41	1.76
	T2	3.188	—	0.57	1.51	0.55
	T3	3.389	—	0.48	1.48	0.33
2	S1	2.868	0.614	0.36	1.44	1.01
	S2	2.970	0.000	1.21	3.14	15.80
	T1	1.529	—	0.54	1.44	1.84
	T2	2.188	—	0.36	1.40	0.14
	T3	2.963	—	1.21	3.14	15.78
3	S1	2.500	0.846	0.35	1.69	1.23
	S2	3.731	0.026	0.60	1.65	2.75
	T1	0.945	—	0.68	1.59	2.94
	T2	3.048	—	0.55	1.63	2.15
	T3	3.058	—	0.52	1.59	2.16
4	S1	2.502	0.751	0.37	1.90	1.93
	S2	3.341	0.371	0.69	3.21	9.37
	T1	1.077	—	0.68	1.80	4.53
	T2	2.737	—	0.67	2.80	8.12
	T3	3.023	—	0.50	1.55	1.49
5	S1	2.214	0.928	0.36	1.47	0.08
	S2	3.341	0.007	0.55	1.95	0.33
	T1	1.132	—	0.43	1.70	0.14
	T2	2.768	—	0.48	1.49	0.10
	T3	3.067	—	0.57	1.76	0.32
6	S1	2.078	0.939	0.37	1.49	0.14
	S2	3.097	0.007	0.82	2.86	1.60
	T1	1.010	—	0.51	1.56	0.36
	T2	2.613	—	0.49	1.63	0.12
	T3	2.979	—	0.83	2.77	1.55
7	S1	2.000	1.032	0.43	1.44	1.43
	S2	3.076	0.923	0.34	1.44	0.63
	T1	1.009	—	0.66	1.41	1.44
	T2	1.779	—	0.40	1.39	0.69
	T3	2.872	—	0.60	1.94	5.23

---

3.4 Further comments on charge transfer

Upon consulting the EDDs, D_CT and ^EMD_CT, we came across some additional oddities that were excluded from the main discussion of the benchmark and are consequently considered here. The first case is compound 1, whose HOMO−2 and HOMO−3 orbitals are located both on the mesityl and BODIPY moieties. The presence of either orbital across those perpendicular moieties shifts to a perceptible degree between XC-functionals (Table S8, ESI‡). The S₂ and T₂ states, which are mixtures of HOMO−3 → LUMO and HOMO−2 → LUMO, are accordingly of more local or CT character depending on the method (Table 5). Nonetheless, we opted to preserve the basis of the correlation on the nature of the excitation and the topology of the orbitals involved rather than the EDD topology. A second example is compound 4 whose first few frontier orbitals (…HOMO−1, HOMO, LUMO, LUMO+1…) are completely or mostly localized on either the pyrrole or benzene side (Table S11, ESI‡). As a result, several excitations exhibit reasonably high charge transfer amplitudes and distances, e.g., S₂ with a q_CT of 0.69 e and D_CT of 2.81 Å. Even its first singlet returns an ^EMD_CT of 1.90 Å despite having a q_CT (0.37e) that indicates local character and is usually accompanied by earth mover's CT-distances around the 1.4–1.5 Å range. Its aza-analog, compound 3, has a similar off-balance frontier MO distribution, also leading to a noticeable but less pronounced inflation of the charge transfer character. Regardless, the (aza-)BOIMPY excitation energies and oscillator strengths behaved in line with the local transitions belonging to the two regular BODIPY structures (1–2) and showed no signs of CT-related issues.

4 Conclusions

The present work has covered the benchmarking of several adiabatic TD-DFT methods towards describing the first few singlet and triplet transitions with the objective of providing a robust basis upon which complex NIR-BODIPY dyes can be designed and analyzed. Using 7 BODIPY structures spanning the green to NIR range, we tested a collection of XC-functionals under solvated conditions, first by comparing with the spin-component-scaled resolution-of-the-identity second-order coupled cluster (SCS-riCC2) results as a means to investigate a broader set of excited states, hererin taking great care to ensure that the transition characters match across methods. Then, the XCFs were evaluated in direct correspondence with experimentally derived S₁ excitation energies and oscillator strenghts.

In light of the entire study, the M06-2X and M06-HF XC-functionals have shown themselves to be the most reliable in terms of accuracy and correlation quality. The local and singlet excitations are treated better by the moderate amount of Hartree–Fock exchange in M06-2X, whilst the higher HF-degree in M06-HF helps to counteract the tendency towards the underestimation in charge transfer and triplet excitations. The LC-BLYP20 and LC-BLYP33 are also respectable choices, but have often ended up lagging slightly behind one of their Minnesota counterparts. Observations pertaining to the Tamm–Dancoff approximation are found to be positive only in the first triplet case, where it succesfully offsets the HF-related instabilities. In contrast, the already overestimated results of the lowest singlet state are pushed further off by the TDA regime, while all other transitions typically experience a rather small upshift, never really doing any big favors in terms of regression performance. The recommended approach therefore is a hybrid procedure in which singlets are treated under the full TD-DFT scheme and the Tamm–Dancoff approximation is employed to handle the triplet domain. The choice between M06-2X and M06-HF depends on the focus of the investigation at hand. M06-2X is a better option when on overall consistency and local singlets are of the essence, whilst M06-HF is best resorted to if the focus lies on energy gaps and charge transfer states. Lastly, the obtained excitation energies will remain rather over and underestimated, but in our observations they do not particularly exceed the usual vertical approximation errors.

Conflicts of interest

The authors have no conflicts to declare.

Data availability

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

Acknowledgements

This work has been realized through the support of the University of Namur and the Special Research Fund of Hasselt University (BOF20DOCNA01). W. Maes thanks the Research Foundation Flanders (FWO Vlaanderen) for financial support (projects G0D1521N and W000620N). B. Champagne thanks the Research Foundation Flanders (FWO Vlaanderen) for financial support (projects G0D1521N). The calculations were performed on the computers of the Consortium des Équipements de Calcul Intensif (CÉCI, https://www.ceci-hpc.be) and particularly those of the Technological Platform of High-Performance Computing, for which the authors gratefully acknowledge the financial support of the FNRS-FRFC, of the Walloon Region, and of the University of Namur (Conventions No. U.G006.15, U.G018.19, U.G011.22, RW1610468, RW/GEQ2016, RW1117545, and RW2110213).

References

1. A. Loudet and K. Burgess, Chem. Rev., 2007, 107, 4891–4932.
2. H.-B. Cheng, X. Cao, S. Zhang, K. Zhang, Y. Cheng, J. Wang, J. Zhao, L. Zhou, X.-J. Liang and J. Yoon, Adv. Mater., 2023, 35, 2207546.
3. D. Wang, X. Wang, S. Zhou, P. Gu, X. Zhu, C. Wang and Q. Zhang, Coord. Chem. Rev., 2023, 482, 215074.
4. J. Huang, X. Wang, Y. Xiang, L. Guo and G. Chen, Adv. Energy Sustainability Res., 2021, 2, 2100016.
5. B. Zhao, L. Liao, Y. Zhu, Z. Hu and F. Wu, J. Lumin., 2023, 263, 120099.
6. P. Kaur and K. Singh, J. Mater. Chem. C, 2019, 7, 11361–11405.
7. L. Gai, Y. Liu, Z. Zhou, H. Lu and Z. Guo, Coord. Chem. Rev., 2023, 481, 215041.
8. V.-N. Nguyen, Y. Yan, J. Zhao and J. Yoon, Acc. Chem. Res., 2021, 54, 207–220.
9. X. Zhang, Z. Wang, Y. Hou, Y. Yan, J. Zhao and B. Dick, J. Mater. Chem. C, 2021, 9, 11944–11973.
10. J. Deckers, T. Cardeynaels, S. Doria, N. Tumanov, A. Lapini, A. Ethirajan, M. Ameloot, J. Wouters, M. Di Donato, B. Champagne and W. Maes, J. Mater. Chem. C, 2022, 10, 9344–9355.
11. J. Wang, N. Boens, L. Jiao and E. Hao, Org. Biomol. Chem., 2020, 18, 4135–4156.
12. V. K. Shukla, G. Chakraborty, A. K. Ray and S. Nagaiyan, Dyes Pigm., 2023, 215, 111245.
13. M. Kaur, A. Janaagal, N. Balsukuri and I. Gupta, Coord. Chem. Rev., 2024, 498, 215428.
14. M. A. Filatov, Org. Biomol. Chem., 2019, 18, 10–27.
15. A. Endo, M. Ogasawara, A. Takahashi, D. Yokoyama, Y. Kato and C. Adachi, Adv. Mater., 2009, 21, 4802–4806.
16. T. N. Singh-Rachford, A. Haefele, R. Ziessel and F. N. Castellano, J. Am. Chem. Soc., 2008, 130, 16164–16165.
17. C. M. Marian, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2012, 2, 187–203.
18. T. J. Penfold, E. Gindensperger, C. Daniel and C. M. Marian, Chem. Rev., 2018, 118, 6975–7025.
19. B. Le Guennic and D. Jacquemin, Acc. Chem. Res., 2015, 48, 530–537.
20. M. R. Momeni and A. Brown, J. Chem. Theory Comput., 2015, 11, 2619–2632.
21. V. Postils, F. Ruipérez and D. Casanova, J. Chem. Theory Comput., 2021, 17, 5825–5838.
22. P. Boulanger, S. Chibani, B. Le Guennic, I. Duchemin, X. Blase and D. Jacquemin, J. Chem. Theory Comput., 2014, 10, 4548–4556.
23. S. Chibani, A. D. Laurent, B. Le Guennic and D. Jacquemin, J. Chem. Theory Comput., 2014, 10, 4574–4582.
24. M. R. Momeni and A. Brown, J. Phys. Chem. A, 2016, 120, 2550–2560.
25. M. Feldt and A. Brown, J. Comput. Chem., 2021, 42, 144–155.
26. W. Helal, Q. Alkhatib and M. Gharaibeh, Comput. Theor. Chem., 2022, 1207, 113531.
27. D. Toffoli, M. Quarin, G. Fronzoni and M. Stener, J. Phys. Chem. A, 2022, 126, 7137–7146.
28. D. Toffoli, A. Ponzi, E. Bernes, M. d Simone, C. Grazioli, M. Coreno, M. Stredansky, A. Cossaro and G. Fronzoni, J. Chem. Phys., 2019, 151, 134306.
29. D. Toffoli, C. Grazioli, M. Monti, M. Stener, R. Totani, R. Richter, L. Schio, G. Fronzoni and A. Cossaro, Phys. Chem. Chem. Phys., 2021, 23, 23517–23525.
30. B. Le Guennic, O. Maury and D. Jacquemin, Phys. Chem. Chem. Phys., 2012, 14, 157–164.
31. S. Chibani, B. Le Guennic, A. Charaf-Eddin, O. Maury, C. Andraud and D. Jacquemin, J. Chem. Theory Comput., 2012, 8, 3303–3313.
32. S. Chibani, B. Le Guennic, A. Charaf-Eddin, A. D. Laurent and D. Jacquemin, Chem. Sci., 2013, 4, 1950–1963.
33. A. Charaf-Eddin, B. Le Guennic and D. Jacquemin, RSC Adv., 2014, 4, 49449–49456.
34. L. J. Patalag, P. G. Jones and D. B. Werz, Angew. Chem., Int. Ed., 2016, 55, 13340–13344.
35. L. J. Patalag, P. G. Jones and D. B. Werz, Chem. – Eur. J., 2017, 23, 15903–15907.
36. T. Freese, L. J. Patalag, J. L. Merz, P. G. Jones and D. B. Werz, J. Org. Chem., 2021, 86, 3089–3095.
37. L. Wang, Z. Xiong, X. Ran, H. Tang and D. Cao, Dyes Pigm., 2022, 198, 110040.
38. S. Shimizu, Chem. Commun., 2019, 55, 8722–8743.
39. Q. Wu, G. Jia, B. Tang, X. Guo, H. Wu, C. Yu, E. Hao and L. Jiao, Org. Lett., 2020, 22, 9239–9243.
40. A. Sbai and J. Guthmuller, Phys. Chem. Chem. Phys., 2024, 26, 25925–25935.
41. X. Wu, X. Xie and A. Troisi, J. Mater. Chem. C, 2024, 12, 18886–18892.
42. J. Tao, D. Sun, L. Sun, Z. Li, B. Fu, J. Liu, L. Zhang, S. Wang, Y. Fang and H. Xu, Dyes Pigm., 2019, 168, 166–174.
43. Z. Wang and J. Zhao, Org. Lett., 2017, 19, 4492–4495.
44. G. M. Fischer, M. Isomaki-Krondahl, I. Gottker-Schnetmann, E. Daltrozzo and A. Zumbusch, Chem. – Eur. J., 2009, 15, 4857–4864.
45. S. Shimizu, T. Iino, Y. Araki and N. Kobayashi, Chem. Commun., 2013, 49, 1621–1623.
46. A. Klamt and G. Schüürmann, J. Chem. Soc., Perkin Trans. 2, 1993, 799–805.
47. S. Grimme, J. Chem. Phys., 2003, 118, 9095–9102.
48. A. Hellweg, S. A. Grün and C. Hättig, Phys. Chem. Chem. Phys., 2008, 10, 4119–4127.
49. J. Tomasi, B. Mennucci and R. Cammi, Chem. Rev., 2005, 105, 2999–3094.
50. T. Le Bahers, C. Adamo and I. Ciofini, J. Chem. Theory Comput., 2011, 7, 2498–2506.
51. D. Jacquemin, T. L. Bahers, C. Adamo and I. Ciofini, Phys. Chem. Chem. Phys., 2012, 14, 5383–5388.
52. Z. Wang, J. Liang and M. Head-Gordon, J. Chem. Theory Comput., 2023, 19, 7704–7714.
53. M. Fraiponts, W. Maes and B. Champagne, J. Chem. Theory Comput., 2024, 20, 2751–2760.
54. A. S. Tarleton, J. C. Garcia-Alvarez, A. Wynn, C. M. Awbrey, T. P. Roberts and S. Gozem, J. Phys. Chem. A, 2022, 126, 435–443.
55. J. C. Garcia-Alvarez and S. Gozem, J. Chem. Theory Comput., 2024, 20, 7227–7243.
56. F. Marin, A. Rohatgi and S. Charlot, WebPlotDigitizer, a polyvalent and free software to extract spectra from old astronomical publications: application to ultraviolet spectropolarimetry, arXiv, 2017, preprint, arXiv:1708.02025 DOI:10.48550/arXiv.1708.02025.
57. 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. 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, Gaussian16 Revision C.01, 2016, Gaussian Inc., Wallingford CT.
58. M. J. G. Peach, M. J. Williamson and D. J. Tozer, J. Chem. Theory Comput., 2011, 7, 3578–3585.
59. M. J. G. Peach and D. J. Tozer, J. Phys. Chem. A, 2012, 116, 9783–9789.
60. T. Cardeynaels, S. Paredis, J. Deckers, S. Brebels, D. Vanderzande, W. Maes and B. Champagne, Phys. Chem. Chem. Phys., 2020, 22, 16387–16399.
61. A. D. Becke, Phys. Rev. A: At., Mol., Opt. Phys., 1988, 38, 3098–3100.
62. C. Lee, W. Yang and R. G. Parr, Phys. Rev. B: Condens. Matter Mater. Phys., 1988, 37, 785–789.
63. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865–3868.
64. A. Becke, J. Chem. Phys., 1993, 98, 5648–5652.
65. S. H. Vosko, L. Wilk and M. Nusair, Can. J. Phys., 1980, 58, 1200–1211.
66. P. J. Stephens, F. J. Devlin, C. F. Chabalowski and M. J. Frisch, J. Phys. Chem., 1994, 98, 11623–11627.
67. C. Adamo and V. Barone, J. Chem. Phys., 1999, 110, 6158–6170.
68. Y. Zhao and D. G. Truhlar, Theor. Chem. Acc., 2008, 120, 215–241.
69. Y. Zhao and D. G. Truhlar, J. Phys. Chem. A, 2006, 110, 13126–13130.
70. O. A. Vydrov and G. E. Scuseria, J. Chem. Phys., 2006, 125, 234109.
71. H. Sun, C. Zhong and J.-L. Brédas, J. Chem. Theory Comput., 2015, 11, 3851–3858.
72. H. Iikura, T. Tsuneda, T. Yanai and K. Hirao, J. Chem. Phys., 2001, 115, 3540–3544.
73. Y. Tawada, T. Tsuneda, S. Yanagisawa, T. Yanai and K. Hirao, J. Chem. Phys., 2004, 120, 8425–8433.
74. I. A. Elayan, L. Rib, R. A. Mendes and A. Brown, J. Chem. Theory Comput., 2024, 20, 3879–3893.
75. S. G. Balasubramani, G. P. Chen, S. Coriani, M. Diedenhofen, M. S. Frank, Y. J. Franzke, F. Furche, R. Grotjahn, M. E. Harding, C. Hättig, A. Hellweg, B. Helmich-Paris, C. Holzer, U. Huniar, M. Kaupp, A. Marefat Khah, S. Karbalaei Khani, T. Müller, F. Mack, B. D. Nguyen, S. M. Parker, E. Perlt, D. Rappoport, K. Reiter, S. Roy, M. Rückert, G. Schmitz, M. Sierka, E. Tapavicza, D. P. Tew, C. van Wüllen, V. K. Voora, F. Weigend, A. Wodyński and J. M. Yu, J. Chem. Phys., 2020, 152, 184107.
76. B. Le Guennic, G. Scalmani, M. J. Frisch, A. D. Laurent and D. Jacquemin, Phys. Chem. Chem. Phys., 2017, 19, 10554–10561.
77. A. Z. Lieberherr, P. Gori-Giorgi and K. J. Giesbertz, J. Chem. Theory Comput., 2024, 20, 5635–5642.
78. N. Bogo and C. J. Stein, Phys. Chem. Chem. Phys., 2024, 26, 21575–21588.
79. C. M. Breneman and K. B. Wiberg, J. Comput. Chem., 1990, 11, 361–373.
80. I. Knysh, F. Lipparini, A. Blondel, I. Duchemin, X. Blase, P.-F. Loos and D. Jacquemin, J. Chem. Theory Comput., 2024, 20, 8152–8174.
81. P. Winget, D. M. Dolney, D. J. Giesen, C. J. Cramer and D. G. Truhlar, Minnesota Solvent Descriptor Database, Department of Chemistry and Supercomputer Institute, Minneapolis, MN, 1999.
82. B. Moore II, H. Sun, N. Govind, K. Kowalski and J. Autschbach, J. Chem. Theory Comput., 2015, 11, 3305–3320.
83. B. Champagne, V. Liégeois and F. Zutterman, Photochem. Photobiol. Sci., 2015, 14, 444–456.
84. B. Champagne, M. Guillaume and F. Zutterman, Chem. Phys. Lett., 2006, 425, 105–109.
85. F. Zutterman, V. Liégeois and B. Champagne, ChemPhotoChem, 2017, 1, 281–296.

---

Footnotes

† This work is dedicated to Christel Marian.

‡ Electronic supplementary information (ESI) available. See DOI: https://doi.org/10.1039/d5cp01830g

---