Photoinduced electron transfer in a bilirubin–oxygen complex: the triplet-reactant computational challenge

Abstract

Although bilirubin photochemistry is central to neonatal jaundice phototherapy, the mechanism of bilirubin photooxidation remains unclear. Here, we use a comprehensive computational approach to investigate whether this mechanism might be initiated by photoinduced electron transfer (PET) to molecular oxygen (O₂), generating superoxide (O₂˙⁻). For this purpose, we employed a simplified bilirubin model compound—tetramethyldipyrrinone (TMD). We use a combination of multireference, coupled cluster methods, and density functional theory techniques to assess the feasibility of the TMD–O₂ complex formation and PET from TMD to O₂. Our results support the feasibility of PET in the TMD–O₂ complex and suggest that PET could potentially initiate bilirubin photooxidation. Beyond the bilirubin case, this work underscores the need for efficient and accurate protocols to compute binding free energies of weak O₂ encounter complexes with organic chromophores in solution, and it highlights the broader challenge of modeling triplet photochemistry with dense manifolds of near-degenerate states, crossings, and strongly state-specific solvent response.

---

1 Introduction

The bilirubin molecule (Fig. 1) is a product of heme catabolism in mammals.¹ Its accumulation in newborns causes neonatal jaundice, which can be treated with phototherapy, serendipitously discovered already in 1956.² However, the understanding of the underlying molecular mechanism remains incomplete.³ The primary degradation mechanism of bilirubin is its photoisomerization, yielding products with solubilities that differ from the native pigment. This process occurs on a sub-picosecond timescale and has been investigated multiple times using classical photochemical techniques,⁴ femtosecond spectroscopies,^5–10 as well as ab initio methods.^10,11 In addition, photocyclization has been examined as an alternative light-induced degradation pathway of bilirubin.^3,7,12

Fig. 1 Stereochemical formula of 4Z,15Z-bilirubin.

There is, however, another possible reaction that has recently come to the forefront: photooxidation, specifically the light-induced reaction between bilirubin and molecular oxygen. Oxygen is ubiquitous in mammalian organisms and is known to bind strongly to some biomolecules, notably porphyrin-based proteins.¹³ Provided that oxygen is present in close proximity to bilirubin, new reaction pathways may open. Two mechanisms have been proposed (Fig. 2): (i) reaction with singlet oxygen generated via bilirubin–oxygen triplet–triplet annihilation.¹⁴ This pathway is most likely inefficient because it requires formation of triplet bilirubin, which is improbable given bilirubin's very short excited-state lifetime.³ (ii) Photoinduced electron transfer (PET) from excited bilirubin to triplet oxygen, yielding the superoxide radical O₂˙⁻.¹⁴ The charge-transfer route might be plausible as superoxide formation upon photoexcitation is well documented for other biomolecular chromophores,¹⁵ and superoxide generally plays important roles in biological redox networks.¹⁶ PET to oxygen was also recently proposed as a relaxation mechanism of excited cyanine dyes.¹⁷

Fig. 2 Proposed two types of bilirubin (BR) photooxidation mechanism.

To investigate the bilirubin photooxidation reaction, a bilirubin model compound, tetramethyldipyrrinone (TMD), was recently synthesized,¹⁴ and its photooxidation was studied in solution. Upon irradiation of TMD in methanol in the presence of triplet oxygen, a hydroperoxide product forms (Fig. 3). The same product was observed in reactions with singlet oxygen, yet no singlet oxygen was detected during direct photooxidation.¹⁴

Fig. 3 Proposed mechanism of the hydroperoxide formation by photooxidation. The investigated electron-transfer step from the S₁ state of TMD to the contact radical ion pair (CRIP) is depicted in red.

These experimental observations raise a key mechanistic question: Could the reaction proceed via electron transfer from the S₁ state of TMD to triplet O₂? The S₁ state is expected to have a lifetime of only a few picoseconds (a similar bilirubin model compound, (Z)-vinylneoxanthobilirubic acid methyl ester, exhibits an S₁ lifetime of approximately 6 ps⁷). Therefore, the electron transfer from the S₁ state of TMD has to be very fast if it plays any role, and the oxygen molecule must already be in close proximity to TMD at the moment of excitation since the diffusion is much slower compared to the S₁ lifetime.

Because the PET mechanism has not yet been verified, neither experimentally nor theoretically, the following questions must be addressed: (i) How probable is it that oxygen appears in close proximity to the bilirubin chromophore so that the photochemical reaction can take place? (ii) Is the photoinduced electron-transfer reaction energetically allowed? (iii) Is the photochemical reaction kinetically feasible on the relevant excited-state timescales?

In this work, we aim to provide new insight using techniques of computational photochemistry and photodynamics, focusing on the plausibility and efficiency of the PET mechanism. A central underlying issue is the suitable choice of electronic-structure methods¹⁸ and the efficient inclusion of environmental effects,¹⁹ which are crucial for charge-transfer reactivity. The present system is relatively large and requires treatment of multiple electronic states, which constrains the choice of feasible methods. We therefore employ efficient local coupled-cluster theory (DLPNO-CCSD(T)²⁰), correlated multireference theory (NEVPT2²¹), constrained density functional theory (CDFT²²) and time-dependent density functional theory (TDDFT²³). It is the combination of these approaches that allows us to report results with increased confidence despite individual methodological limitations.

Beyond the bilirubin case, understanding the reactivity of excited organic molecules with O₂ is a general problem: such reactions can yield multiple pathways and products, with outcomes shaped by subtle electronic and environmental factors.²⁴ Understanding PET is also crucial for photovoltaics and photocatalysis.²⁵ However, modeling photochemical processes in large triplet systems, where multiple excited states of different character (such as charge-transfer and locally excited) are involved and state-specific solvation effects are important, restricts the range of applicable methods and remains a challenge for computational chemistry.

2 Methods

2.1 Free energy calculation of the TMD–O₂ complex formation

We used a similar procedure to that successfully applied to association free energies of noncovalently bound host–guest complexes:²⁶ (i) we calculated electronic energies with the accurate DLPNO-CCSD(T) method extrapolated to the complete basis set and pair natural orbital space. (ii) We obtained thermostatistical corrections from energy to Gibbs free energy in vacuum at the ωB97X-D3/6-31+G* level. (iii) We included a solvent correction using B3LYP/6-31G* and the SMD solvent model as implemented in Q-Chem v6.1.²⁷ (iv) In addition, we considered multiple possible oxygen binding sites (optimized at the ωB97X-D3/6-31+G* level), which contribute to the configurational entropy.²⁸ We performed all calculations in ORCA v6.0²⁹ unless stated otherwise. Further details of the computational protocol are provided in the SI.

2.2 Modeling electronic absorption spectra

The ground state of the TMD was sampled with adiabatic molecular dynamics in the ABIN v1.1 code.³⁰ The electronic structure was described at the OM3 level in vacuum. Temperature was maintained at 298.15 K with the Nosé–Hoover thermostat. The dynamics were propagated for 100 000 steps with the velocity Verlet integrator and a time step of 10 a.u. (≈0.24 fs). The first half of the trajectory was discarded as a thermal equilibration period. The ground state density was mimicked with 100 configurations selected from the second half of the trajectory at constant intervals of 500 steps. To model the spectrum, we used the nuclear ensemble approach.³¹ We employed coupled-cluster, multireference, and TDDFT electronic structure methods. Further computational details are provided in the SI.

2.3 Ionization energies and electron affinities

We calculated ionization energies (IEs) as the difference between the energies of ionized TMD and neutral TMD, either at the neutral ground-state geometry (vertical IE) or allowing relaxation to the optimized ionic geometry (adiabatic IE). Similarly, the vertical and adiabatic electron affinities (EAs) were calculated as the differences in the electronic energies of neutral triplet oxygen and the superoxide radical anion with or without the geometry relaxation. The water environment influencing IEs and EAs was modeled with the integral equation formalism polarizable continuum model (IEFPCM³²). For the vertical ionization energy, only the fast component of polarization was included (non-equilibrium IEFPCM), while for the adiabatic ionization energy both the fast and slow components of the polarization were considered (equilibrium IEFPCM).^33,34 Calculations were done at ωB97X-D3/6-31+G* level in the Q-Chem v6.2 code.²⁷

2.4 Mapping the potential energy surface

We optimized the ground-state geometries of TMD (S₀) and the TMD–O₂ complex (T₁) with DFT using ωB97X-D3/6-31+G*. The locally excited state was optimized with TDDFT at the ωB97X-D3/6-31+G* level (for TMD as S₁ and for the complex as T₆), and the electron-transfer state with constrained density functional theory (CDFT) using B3LYP/6-31G*. We performed the DFT and TDDFT calculations in ORCA v6.0,²⁹ and the CDFT optimizations in NWChem v6.8.³⁵

To estimate the reaction coordinate between the critical points on the PES, we used the geodesic interpolation, as an efficient alternative to the linear interpolation in internal coordinates (LIIC), implemented in the Python package geodesic-interpolate.³⁶

The choice of electronic structure method was dictated by the molecular size and the presence of multiple excited states of different character which require computationally efficient approaches, particularly for dynamical simulations.

As a reference method, we used the multireference quasi-degenerate NEVPT2 (QD-NEVPT2)/aug-cc-pVDZ with an active space of 6 electrons in 5 orbitals and 14 electronic states, as implemented in ORCA v6.0.²⁹ The active space comprised the π orbitals (HOMO−1 and HOMO) of TMD, two π* orbitals of O₂, and the π* orbital (LUMO) of TMD (orbital plots are provided in Fig. S2). This setup provided an accurate description of both the bright LE state and the CT state. Along the geodesic coordinates, we additionally evaluated the ground-state curve in vacuum with DLPNO-CCSD(T)/aug-cc-pVDZ using TightSCF and TightPNO in ORCA v6.0²⁹ to provide an accurate single-reference reference for the ground-state energetics. For comparison, we also employed XMS-CASPT2(6,5)/6-31G* in OpenMolcas v24.06,³⁷ CDFT/B3LYP/def2-TZVPD and CDFT-CI/B3LYP/def2-TZVPD in Q-Chem v6.2,²⁷ and TDDFT with ωB97X-D3/def2-TZVPD in ORCA v6.0.²⁹ as a faster method suitable for subsequent simulations.

Because state-specific solvation is not available for QD-NEVPT2 in ORCA, we incorporated solvent effects by adding state-specific solvent shifts computed at lower levels to the QD-NEVPT2 vacuum energy profiles. Nonequilibrium shifts were taken from TDDFT ωB97X-D3/def2-TZVPD calculations in ORCA v5.03 using CPCM water with ε_opt = 1.78. Equilibrium shifts for the ground and charge-transfer states were obtained from CDFT/B3LYP/def2-TZVPD with PCM water using ε_stat = 78.39, while the equilibrium shift for the LE state was computed with TDDFT ωB97X-D3/def2-TZVPD in ORCA v5.03 using CPCM with ε_stat = 78.39. This additive scheme provides an approximate description of solvation effects on the multireference reference curves used for the PET energetics.

2.5 Nonadiabatic simulations

To model the nonadiabatic transitions following the photoexcitation, we used the Landau–Zener surface hopping algorithm³⁸ in the ABIN v1.1 code.³⁰ The simulations were performed on time-dependent density functional theory (TDDFT) potential energy surfaces within the Tamm–Dancoff approximation (TDA) at the ωB97X-D3/def2-TZVPD level. The use of TDDFT for excited state simulations is questionable for multiple reasons: (i) poor description of bond breaking³⁹ (ii) possible appearance of spurious charge transfer states⁴⁰ (iii) issues with topology and topography of S₀/S₁ intersections.^41,42 We have partially mitigated the problems by using range-separated functionals and benchmarking them against the QD-NEVPT2 method. Still, the results have to be interpreted with caution.

We generated the initial conditions for the nonadiabatic dynamics from a ground-state adiabatic molecular dynamics simulation of the complex in vacuum using the OM3 method, as implemented in the ABIN v1.1 code.³⁰ The simulation was propagated for 100 000 steps with a time step of 10 a.u. (≈0.24 fs) using the velocity Verlet integrator and the quantum thermostat based on generalized Langevin equation (GLE)⁴³ at 298.15 K. We selected 400 geometries and velocities along the trajectory at constant intervals of 250 steps and discarded the first 100 to ensure thermal equilibration, which resulted in a set of 300 initial conditions. Finally, we randomly selected 35 initial conditions for launching the simulations.

The nonadiabatic dynamics were started in the state with highest oscillator strength (local HOMO–LUMO excitation on TMD) mimicking the experimental conditions. The excited-state index corresponding to the local state was highly sensitive to the complex geometry. We had to consider 9 excited states. A time step of 5 atomic units (approximately 0.12 fs) was used throughout the simulations. The simulations were complicated by the presence of spin-contaminated states. We classified a state as spin-contaminated when 〈Ŝ²〉 > 2.2.

3 Results and discussion

The results are organized according to the sequence of elementary steps in the proposed PET mechanism: (i) we examine the feasibility of TMD–O₂ complex formation, and (ii) we assess the feasibility of electron transfer from the bright excited state of TMD to O₂.

3.1 TMD–O₂ complex formation

We estimate the feasibility of TMD–O₂ complex formation by calculating the Gibbs free energy of complex formation. To do so, we first need to perform a conformational analysis of TMD and the TMD–O₂ complex.

3.1.1 Conformational analysis of TMD. The presence of a methine bridge in the TMD molecule leads to two types of isomers: Z and E configurations relative to the double bond and s-cis, s-trans conformers arising from the rotation around the single bond (Fig. 4).

Fig. 4 Configurations and conformations of the TMD with their relative energies (s-cis Z serving as a reference) and estimated Boltzmann populations at 298 K. Groups relevant for different types of isomers are highlighted: red color for E/Z isomerization, blue color for s-cis/s-trans interconversion. All calculations were performed at the ωB97X-D3/6-31+G* level.

Calculations at the ωB97X-D3/6-31+G* level of theory indicate that the s-cis Z has the lowest energy (Fig. 4) and therefore we study the PET only on this isomer. Thermal populations of E configurations are negligible at 298 K.

To confirm the stability of the s-cis Z conformer, we ran a 24-ps-long molecular dynamics simulation in vacuum at 298.15 K, employing a quantum GLE thermostat. The molecule remains in its s-cis Z conformation throughout the whole simulation, validating our choice of conformer.

Note that under irradiation, a photostationary state containing an approximately 1 : 1 mixture of the Z and E isomers is formed in experiments.¹⁴ However, for simplicity we decided to focus solely on the s-cis Z conformer.

3.1.2 Conformational analysis of TMD–O₂ complex. We have identified 29 unique local minima on the TMD–O₂ ground-state potential energy surface at the ωB97X-D3/6-31G* level, corresponding to different O₂ binding positions and small variations in the TMD geometry (details of the procedure in the SI). The lowest-energy conformation is shown in Fig. S7 in the SI, with O₂ located above the pyrrole ring. Several similar minima share this binding motif. In these structures, the O⋯H distance to the NH group in the lactam ring is 2.6–2.7 Å, and together they account for about 39% of the population at 298 K (Boltzmann distribution). In other conformations, O₂ is located further above the pyrrole ring or on the opposite side of the pyrrole or lactam ring relative to the lowest-energy structure. The energies of the local minima relative to the global minimum and corresponding populations calculated from the Boltzmann distribution are presented in Table S4 in the SI. Overall, the conformational analysis shows that O₂ can bind at several positions, but it preferentially interacts with the aromatic ring. In particular, conformations with O₂ located above the pyrrole ring are the most favorable.

The calculations with implicit solvent models imply that the solvent effects are not significant. The calculated solvation Gibbs energies with implicit solvent models are small, see Table S3 in the SI. We also located minimum geometries in water and methanol at the PCM level, and the differences relative to the vacuum geometry were negligible. However, implicit models may not capture all relevant interactions, and the energetics and population could still differ, for example in water, where oxygen might hinder hydrogen bond formation.

3.1.3 Gibbs free energy of complex formation. The calculated Gibbs free energy of TMD–O₂ complex formation is 9.24 kJ mol⁻¹ in methanol and 8.22 kJ mol⁻¹ in water (methods used for the calculation are in Section 2.1). These values can be decomposed into the following contributions: electronic interaction energy in vacuum of −10.11 kJ mol⁻¹, thermostatistical correction in vacuum of 27.90 kJ mol⁻¹ (including the standard-state correction), a configurational entropic term of −8.77 kJ mol⁻¹, and solvation reaction Gibbs energies of −0.80 kJ mol⁻¹ in water and 0.22 kJ mol⁻¹ in methanol. The sum of these contributions yields the final Gibbs energies. A full benchmark of the electronic interaction energies and solvation Gibbs energies is provided in the SI (Tables S2 and S3).

To estimate the equilibrium constant and the fraction of the complex, we assume concentration of TMD c_TMD = 10 µmol dm⁻³ (concentration used in experiments with a similar compound⁹), concentration of oxygen c_oxygen = 1 mmol dm⁻³ corresponding to oxygen solubility in water⁴⁴ and 10 mmol dm⁻³ for oxygen solubility in methanol.⁴⁵

Solving the equilibrium equations yields the complex fractions 0.004% in water and 0.024% in methanol (Table 1). The photooxidation quantum yield of a similar compound—the bilirubin subunit (Z)-isovinylneoxanthobilirubic acid methyl ester—is reported to be between 0.0011 and 0.0014, depending on the wavelength.⁹ This suggests that at least approximately 0.1% of TMD molecules must form the complex for the photooxidation to be explainable only with the PET, considering that the ET efficiency would be 100% assuming the most favorable case. However, since the calculation is approximate—due to the sensitivity of the exponential function and several other limitations (details in the SI)—the reaction channel remains feasible if the efficiency of the electron transfer is high.

Table 1 Calculated binding Gibbs energies of TMD to oxygen in water and methanol at 298 K with the corresponding equilibrium constant and the fraction of the complex

Solvent	ΔGbinding (kJ mol^−1)	K	fcomplex (%)
Water	8.22	0.036	0.004
Methanol	9.24	0.024	0.024

---

3.2 Photoinduced electron transfer from TMD to oxygen

We assess the feasibility of electron transfer from the S₁ state of TMD to O₂ by comparing the energies of bright and charge-transfer states, first for the isolated molecules and then for the complex. Finally, we perform exploratory nonadiabatic dynamical simulations to estimate the quantum yield of PET.

3.2.1 Calculations on the isolated molecules. We first study PET based on the properties of isolated TMD and oxygen molecules: the energy of the bright S₁ state and the asymptotic charge-transfer state.
3.2.1.1 Position of bright excited state. Excitation to the bright state of TMD is the first step in the PET process. The excitation corresponds to the S₀ → S₁ transition. It is mainly described by the HOMO–LUMO transition with a π → π* character, as shown in Fig. 5. At the DLPNO-STEOM-CCSD/def2-TZVP level, the vertical excitation energy at the minimum S₀ geometry is 3.63 eV in vacuum and 3.44 eV in water, with high oscillator strengths around 0.6. The minimum S₁ state geometry is structurally very similar to the S₀ minimum (Fig. S1). The adiabatic excitation energy is 3.20 eV in vacuum and 3.08 eV in water. The vertical excitation energy at the S₁ minimum is 2.75 eV and 2.60 eV, respectively.

Fig. 5 Comparison of experimental absorption spectrum of Z-TMD in methanol (digitized from ref. 14) and calculated spectra with DLPNO-STEOM-CCSD, DLPNO-NEVPT2(2,2), MRCISD/OM3(16,16), LR-TDDFT/TDA/B3LYP, LR-TDDFT/TDA/ωB97X-3c (ω = 0.1a₀⁻¹) and LR-TDDFT/TDA/ωB97X-D3 (ω = 0.25a₀⁻¹) methods. See Table S1 for the employed basis sets. The Gaussian broadening parameter was set to 0.11 eV for DLPNO-STEOM-CCSD, 0.27 eV for TDDFT, 0.15 eV for MRCISD/OM3, and 0.049 eV for NEVPT2. Methanol solvent was included with CPCM apart from MRCISD/OM3 where the calculation was done in vacuum. Orbitals contributing to the S₀ → S₁ transition responsible for the main absorption band were calculated at the DLPNO-STEOM-CCSD/def2-TZVP level. Isovalue was set to 0.2.

The experimental absorption spectrum in methanol has a maximum around 3.1 eV. We calculated linear electronic absorption spectrum of s-cis Z-TMD with different electronic structure methods and compared them with the experimental spectrum, see Fig. 5. The highly correlated DLPNO-STEOM-CCSD method provided the best agreement with the experimental spectrum in terms of the position, intensity of the absorption maximum, and the peak width. Unfortunately, the DLPNO-STEOM-CCSD method in its present implementation does not allow the calculation of triplet excited states and is therefore not applicable to the complex.

The DLPNO-NEVPT2 method reproduces the experimental band position well but clearly overestimates the peak intensity of the main peak on the absolute scale. This discrepancy is most likely related to the limited active space and to the fact that only a single excited state was computed. Because of its multireference character and its ability to describe local and charge-transfer excitations on the same footing, NEVPT2 is suitable for the further description of the states in the TMD–O₂ complex.

Among the density functional methods, the B3LYP functional reproduces the experimental spectrum well. However, global hybrids such as B3LYP generally perform poorly for charge-transfer excitations, predicting their positions to be artificially low.⁴⁰ Within range-separated functionals, ωB97X-3c (ω = 0.1a₀⁻¹) performs best, with a blue-shift around 0.1 eV compared to the experiment. Its computational cost is low, but excited-state gradients are not available. We therefore also tested the ωB97X-D3 functional with the default range-separation parameter ω = 0.25a₀⁻¹ which gives a larger blue shift compared to the experiment. We still use the ωB97X-D3 functional with the ω = 0.25a₀⁻¹ later in the exploratory molecular dynamics simulations because the vacuum calculations yield LE and CT energetics similar to the NEVPT2 results in water.

The MRCISD/OM3(16,16) method provides us with a peak that is blue-shifted by 0.4 eV relative to the experiment even with this large active space. A limitation of this benchmark is that excited-state calculations were performed near the Franck–Condon region, so the comparison is restricted to a small part of configuration space. Further details of the benchmark are provided in the SI.

3.2.1.2 Position of charge-transfer state. We determine the asymptotic position of CT state from the ionization energy of TMD and the electron affinity of O₂.

The experimental adiabatic electron affinity (AEA) of O₂ in the gas phase is 0.45 eV.⁴⁶ The DLPNO-CCSD(T) method with extrapolation to the complete PNO space and complete basis set (DLPNO-CCSD(T)) yields an AEA of 0.41 eV, in good agreement with experiment. We therefore adopt this method as a reference. At the DLPNO-CCSD(T) level, the adiabatic CT energy calculated as AIE–AEA (Table 2) is 6.59 eV. The vertical charge-transfer (CT) energy at infinite separation, calculated as VIE–VEA, is 7.48 eV.

Table 2 Vertical ionization energies (VIE) and adiabatic ionization energies (AIE) of TMD, and vertical (VEA) and adiabatic (AEA) electron affinities of oxygen. Geometries were optimized at ωB97X-D3/6-31+G* in vacuum, and single-point energies were evaluated with DLPNO-CCSD(T) extrapolated to CBS and complete PNO space. The ωB97X-D3 calculations in water (PCM) are also shown, using nonequilibrium solvation for vertical and equilibrium solvation for adiabatic quantities

Method	Solvent	VIE (eV)	VEA (eV)	AIE (eV)	AEA (eV)
DLPNO-CCSD(T)/CBS	Vacuum	7.27	−0.21	7.00	0.41
ωB97X-D3/def2-TZVPD	Vacuum	7.02	−0.19	6.75	0.28
ωB97X-D3/6-31+G*	Vacuum	7.03	−0.15	6.76	0.42
ωB97X-D3/6-31+G*	Water	6.35	1.21	5.00	3.64

---

The IE and EA values in water were calculated with the ωB97X-D3/6-31+G* functional, which performed well for the vacuum IE/EA (Table 2) and allows nonequilibrium solvation for vertical quantities, using a PCM description of bulk water. The corresponding vertical CT energy at infinite separation in water is 5.14 eV, and the adiabatic CT energy is 1.36 eV. Solvation stabilizes both the oxidized donor and the reduced acceptor.

Finally, we estimated the vertical CT energies in the complex by adding a Coulomb stabilization term to the asymptotic value in vacuum obtained at the DLPNO-CCSD(T) level and in water obtained at the ωB97X-D3 level. The coulombic term was calculated in a point-charge approximation using an effective separation of 3.4 Å (the distance between the geometry centers, defined as the arithmetic mean of atomic coordinates, of TMD and O₂ in the minimum-energy T₁ complex geometry). The Coulomb term is −4.24 eV in vacuum and −2.38 eV in water. This yields vertical CT energies of 3.24 eV (vacuum) and 2.76 eV (water) at the T₁ minimum geometry.

3.2.1.3 Comparing bright state and charge-transfer state energetics. The comparison of vertical and adiabatic quantities provides a first energetic estimate of the PET process. At the T₁ minimum of the complex, the vertical CT state lies below the bright S₁ state of TMD in both vacuum and water, indicating that PET is energetically feasible (Table 3). Formation of fully separated ions, however, is only favorable in solution: in water, the energy of the separated, relaxed ions (1.36 eV) is lower than the estimated vertical CT energy at the complex geometry (2.76 eV), whereas in vacuum the opposite holds (6.59 eV vs. 3.24 eV, Table 3). The quantitative values in water must be interpreted with caution, since the EA of O₂ is described here by a continuum solvent model; more accurate solvation energies would likely require a cluster–continuum treatment or explicit solvent.⁴⁷ Finally, this analysis assumes that the S₁ state of TMD is not strongly perturbed by the presence of oxygen, which is confirmed by calculations of the full complex in the following section. A concise summary of the PET energetics is provided in Table 3.

Table 3 Summary of PET energetics (eV) in vacuum and water. VCT(∞) = VIE–VEA and ACT(∞) = AIE–AEA, obtained at the DLPNO-CCSD(T)/CBS level in vacuum and at ωB97X-D3/6-31+G*/PCM in water using the IE/EA values from Table 2. VCT(T₁ min) is the vertical CT energy at the T₁ minimum geometry, estimated from these IE/EA values plus a point-charge Coulomb term. VEE(S₁) is the vertical excitation energy of the bright S₁ state at the DLPNO-STEOM-CCSD/def2-TZVP level

Environment	VCT(∞)	ACT(∞)	VCT(T1 min)	VEE(S1)
Vacuum	7.48	6.59	3.24	3.63
Water	5.14	1.36	2.76	3.44

---

3.2.2 Calculations on the complex. We then study PET by calculating the excited states of the complex at critical geometries, interpolating between them to obtain an approximate reaction coordinate, and finally performing exploratory nonadiabatic dynamical simulations to estimate the quantum yield of PET.
3.2.2.1 Potential energy surface mapping and interpolations. We first benchmark the methods capable of describing the CT state in the complex by comparing the energies with the coupled cluster asymptotic value 7.48 eV. The quasi-degenerate (QD)-NEVPT2(6,5)/aug-cc-pVDZ (14 states) provided the best agreement with the DLPNO-CCSD(T) method, giving a vertical CT energy of 7.58 eV at a separation of 700 Å. The XMS-CASPT2(6,5)/6-31G* method overestimates the CT energy (8.69 eV), independent of the choice of IPEA and imaginary shift (Table S5 in the SI); adding diffuse functions lowers the energy by about 0.5 eV which still leads to a significant overestimation. The CDFT/B3LYP/def2-TZVPD calculation converged only at a separation of 7.5 Å; the CT state energy at this distance is 5.23 eV, subtracting the Coulomb interaction between molecules estimated from Mulliken charges gives an asymptotic vertical CT energy of 7.03 eV. The TDDFT/TDA approach with ωB97X-D3/def2-TZVPD functional places the CT state at 6.64 eV at 700 Å. Summary of the results is provided in Table 4.

Table 4 Asymptotic vertical charge-transfer (CT) energies. DLPNO-CCSD(T) values are obtained as IE(TMD)–EA(O₂). Other methods use the ground-state-to-CT energy gap at a TMD–O₂ separation of 700 Å (7.5 Å for CDFT with an electrostatic correction)

Method	Vertical CT energy (eV)
DLPNO-CCSD(T)/CBS	7.48
QD-NEVPT2(6,5)/aug-cc-pVDZ	7.58
XMS-CASPT2(6,5)/6-31G*	8.69
CDFT/B3LYP/def2-TZVPD	7.03
ωB97X-D3/def2-TZVPD	6.64

---

Then, we calculate the triplet excited states in the ground-state minimum geometry of the TMD–O₂ complex in vacuum. The nature of the excited states was analyzed to identify the locally-excited (LE) donor state and the intermolecular CT states. Note that there are two nearly-degenerate CT states corresponding to a charge transfer excitation from TMD to the two nearly-degenerate π* orbitals of O₂.

The position of the LE state appears between 3.6 and 4.0 eV, depending on the method used. QD-NEVPT2 gives an energy of about 3.6 eV, while other methods overestimate its position slightly, up to ∼4.0 eV. The LE state always has a large oscillator strength. Its energy is close to the S₀ → S₁ transition of isolated TMD, so the presence of O₂ predictably has only a minor effect on the donor excitation. In contrast, the CT states have very small oscillator strengths (on the order of 0.002). Their excitations populate the π* orbitals of O₂, consistent with charge transfer.

We used the same electronic-structure methods as for the asymptotic calculations. Their energies are strongly method dependent, in line with the asymptotic calculations. QD-NEVPT2 places the lowest CT state at about 4.0 eV at the ground-state minimum. CDFT predicts a similar value, but this agreement is probably due to error cancellation, because at infinite separation CDFT underestimates the CT energy by ∼0.5 eV relative to QD-NEVPT2. CDFT-CI lowers the CT energy to about 3.5 eV, consistent with this expected shift. TDDFT methods place the CT state around 3.0 eV. XMS-CASPT2 gives the highest CT energy, close to 5.0 eV, as expected from its behavior in the CT benchmark. The CASPT2 and TDDFT curves are provided in Fig. S3 in the SI.

The minima of the LE and the CT states are key points on the PES for the PET mechanism. The minimum of the LE state is geometrically very similar to the ground-state complex minimum (Fig. 6). The geometry of TMD in the complex closely resembles the S₁ minimum of isolated TMD, indicating that the O₂ molecule has only a minor influence on the LE-state structure.

Fig. 6 QD-NEVPT2(6,5)/aug-cc-pVDZ energies of the ground state (GS), locally excited (LE) state, and charge-transfer (CT) state of the TMD–O₂ complex along geodesic interpolation coordinates connecting the GS and LE minima and the LE and CT minima in vacuum (left figure) and in water (right figure). We also show the DLPNO-CCSD(T) ground-state curve in vacuum. Note that there are always two nearly-degenerate CT states shown, corresponding to a charge transfer to the two nearly-degenerate π* orbitals of O₂.

The CT minimum was optimized with CDFT, which provides a convenient way to optimize a state with explicit charge separation (Fig. 6). In this geometry, the O₂ molecule is located between the two NH groups of TMD, and the rings of TMD are essentially coplanar. A similar CT-state geometry was also found at the TDDFT level.

These three stationary points (ground state minimum, LE minimum, and CT minimum) define the endpoints of the geodesic interpolations³⁶ used below to construct approximate relaxation pathways for PET. The interpolation between these points allowed us to understand the relaxation of the TMD–oxygen complex after photoexcitation, representing a first estimate of a reaction coordinate for PET. We assume that, after excitation to the LE state, the system initially relaxes toward the minimum of this state. Along such a path, close approach of the LE and CT states would be a first indication that PET is feasible and potentially fast.

We therefore interpolated energies along the two coordinates in a subsequent fashion: (i) between the ground-state (GS) minimum and the LE minimum, and (ii) between the LE minimum and the CT minimum (Fig. 6). As a reference method, we used QD-NEVPT2(6,5) with 14 states, which reproduces both the bright LE state and the CT state well compared to experiment. The resulting energy profiles in vacuum are shown in Fig. 6. The CT state remains above the LE state along the entire GS–LE path, indicating that PET is unlikely to be feasible in the gas phase.

The GS curve in vacuum can be computed accurately with DLPNO-CCSD(T) (Fig. 6). The curve agrees with QD-NEVPT2 at the GS–LE part of the coordinate, however it diverges at the LE–CT part of the coordinate. In DLPNO-CCSD(T), the GS energy rises at the end of the LE–CT coordinate, yet this behavior might be related to increasing spin contamination of the DLPNO-CCSD(T) calculations observed in this region as highlighted by T₁ diagnostics (Table S7, SI). Nevertheless, CDFT and CDFT-CI yield GS curves comparable to DLPNO-CCSD(T), see Fig. S4 in the SI.

Solvent effects were considered in two limiting regimes. In the equilibrium regime, the solvent is assumed to fully relax to the electronic state of the solute; in the nonequilibrium regime, only the fast electronic polarization of the solvent responds, while the nuclear configuration of the solvent is frozen. Because state-specific solvation is not available for the NEVPT2 method in the ORCA code, we added state-specific solvent shifts obtained from TDDFT and CDFT calculations to the QD-NEVPT2 vacuum curves. Nonequilibrium solvation was modeled with the CPCM approach using the optical dielectric constant of water (ε_opt = 1.78), and equilibrium solvation with a PCM model using the static dielectric constant (ε_stat = 78.39).

The resulting profiles in water are shown in Fig. 6. One limiting pathway corresponds to excitation into the nonequilibrium-solvated LE state, rapid relaxation on this LE surface, crossing to the nonequilibrium-solvated CT state, and subsequent relaxation to the CT minimum on the equilibrium-solvated surface as the solvent reorganizes. The CT state eventually crosses the ground state and becomes the lowest state. Note that the GS/CT crossing likely occurs even closer to the LE minimum geometry because the reference DLPNO-CCSD(T) GS curve rises more than the QD-NEVPT2 curve.

A second limiting mechanism would involve solvent relaxation on the LE state (transition from nonequilibrium- to equilibrium-solvated LE), followed by a crossing to the CT state and subsequent relaxation as above. In reality, the dynamics will likely lie between these two extremes, sampling partially equilibrated solvation environments. In all cases, the presence of an LE/CT crossing near the GS minimum in water provides an indication that PET could be feasible.

Note that we have also calculated singlet charge-transfer state energies along the geodesic interpolation coordinate in addition to the triplet excited states manifold (Fig. S5 in the SI). The energies of the singlet and triplet CT states are almost degenerate due to the weak interaction between TMD and O₂. However, we do not expect the singlet CT state to be involved in the mechanism since bilirubin has a very short excited-state lifetime.

3.2.2.2 Nonadiabatic simulations. To estimate the efficiency of the PET process, we performed nonadiabatic Landau–Zener surface hopping simulations with TDDFT in vacuum. We chose the ωB97X-D3/def2-TZVPD functional (ω = 0.25a₀⁻¹) because it provides energies comparable to the NEVPT2 reference results in water. In particular, the LE state lies at about 4 eV in both descriptions (Fig. 6 and Fig. S3). The CT state positioned about 3 eV is also in reasonable agreement with the NEVPT2 value in equilibrated water (about 2.2 eV). Note that such an approach is necessarily very approximate and its sole goal is to confirm the feasibility of the charge-transfer reaction.

The initial conditions were generated from the long OM3 simulation in the singlet multiplicity. The reason for using a singlet instead of a triplet multiplicity stems from a benchmark of minimal geometries on both OM3 and DFT levels of theory presented in the SI. The singlet OM3 geometry resembled the triplet DFT minimal geometry better than the triplet OM3 minimal geometry, which placed oxygen too far from TMD. A direct DFT sampling of initial conditions was hindered by a lack of computational resources. Since the subsequent TDDFT-LZSH simulations are intended as exploratory feasibility tests rather than quantitatively converged yield calculations, we believe that the OM3 sampling is acceptable.

We conducted 35 exploratory TDDFT-LZSH trajectories. We quantified CT formation at 110 fs, when 24 trajectories were still running. At this time, 12 of the 24 trajectories were in the CT state. The remaining trajectories evolved into states with a high spin contamination, exhibiting an average total energy drift of 0.43 eV. The drift was linked to the propagation on the spin-contaminated states; therefore, we did not consider these trajectories relevant and excluded them from the analysis.

We therefore also analyzed the trajectories at 50 fs, where the average drift was 0.09 eV. At 50 fs, 20 of 32 trajectories were already in the CT state, while the remaining trajectories were spin-contaminated. These numbers correspond to a CT conversion yield of 0.63 at 50 fs and 0.50 at 110 fs. Again, these numbers are only indicative that the process is feasible, rather than providing an accurate estimate of reaction yields.

Despite the known limitations of TDDFT and the presence of spin-contaminated states in these trajectories, the results support that PET is feasible and can occur on a femtosecond timescale.

4 Conclusions

In this work, we investigated the photoinduced electron transfer (PET) from a bilirubin model compound, TMD, to molecular oxygen as a plausible contribution to bilirubin photooxidation during neonatal jaundice phototherapy. The central objective was to assess, using computational chemistry, whether this pathway is thermodynamically and kinetically viable in an aqueous environment and whether it can compete with other pathways on ultrafast timescales.

Our results support the feasibility of PET in water. The CT state is stabilized relative to the LE state at the NEVPT2 level (evaluated for the minimum geometry of the ground-state encounter complex), consistent with a driving force for electron transfer. Interpolations between key points on the PES at the NEVPT2 level revealed (i) LE–CT crossing in the vicinity of the minimum ground state complex geometry and (ii) a CT-ground state crossing near the minimum of the CT state, both supporting a physically consistent reaction scenario in which population transfer to the CT manifold can occur and subsequently relax toward the ground state.

Complementary to this, the TDDFT/Landau–Zener surface-hopping simulations indicate that the LE → CT conversion can be efficient and ultrafast, yielding a CT population of about 0.6 within 50 fs. While the simulations are necessarily approximate, their purpose here is not to provide an accurate reaction yield but to demonstrate that an efficient electron-transfer channel is compatible with the computed electronic-state topology.

At the same time, our equilibrium estimate for the fraction of the pre-reactive TMD–O₂ encounter complex suggests weak binding in solution (Table 1). Because the predicted fraction depends exponentially on the binding free energy and because ΔG_binding was obtained using approximate statistical-thermodynamic models, a small calculated bound fraction does not by itself exclude PET as a relevant pathway. It rather highlights the sensitivity of this ingredient and the need for improved free-energy methodologies for weak, open-shell encounter complexes in solution.

To further support our results, future experimental validation is crucial. Ultrafast spectroscopic techniques such as time-resolved Raman spectroscopy could directly probe the occurrence and timescale of PET and its competition with the possible alternative pathways involving singlet oxygen generation.¹⁴ Detection of superoxide would serve as a strong indicator of the electron transfer.

Beyond the specific conclusions drawn for the present bilirubin–oxygen system, several more general observations emerge.

First, photochemical interactions of dissolved O₂ with organic chromophores are closer to a rule than an exception. In solution, O₂ efficiently perturbs excited-state dynamics through (i) triplet–triplet energy transfer leading to singlet oxygen^48,49 and (ii) photoinduced electron transfer leading to superoxide⁵⁰ and related reactive oxygen species. At the same time, systematic theoretical descriptions and benchmarks for ground-state O₂ encounter complexes and the ensuing local/charge-transfer manifolds in condensed phase remain scarce, with most studies focusing on small prototype complexes.^51,52

Second, the present system highlights a broader methodological need for robust and efficient protocols to compute binding free energies of weak O₂ encounter complexes in solution. In contrast to classical tight complexes, O₂ association with organic chromophores is often characterized by multiple shallow minima and substantial orientational degeneracy, so that ΔG_binding depends critically on configurational sampling, solvent reorganization, and the precise definition of the bound ensemble and standard state. Related difficulties are well known in the bioinorganic context of hemoglobin.⁵³

Third, the present problem is challenging not because of a single bottleneck, but because several of them coincide: local and CT states are energetically close, the intersections between ground and excited states are important, solvent response is essential and state-dependent and the ground state of the system lies in the triplet electronic manifold. This combination pushes the limits of widely used excited-state protocols. In particular, standard linear-response TDDFT is known to struggle for CT excitations, and the description of near-degeneracies and crossings in open-shell manifolds calls for methods that remain balanced across spin and multiconfigurational regimes. Advancing such approaches is also practically important for predicting transient absorption spectra and other time-resolved observables in oxygenated environments.

Finally, it is encouraging that modern locally correlated wave-function methods now make it feasible to treat oxygen–chromophore complexes beyond functional-dependent DFT at realistic system sizes. In particular, DLPNO-based coupled-cluster approaches provide a systematic route toward high accuracy at a manageable cost for complex molecules, making them a natural backbone for future benchmarks of O₂ binding, CT energetics, and the coupled role of solvation and state crossings in photochemistry.

Conflicts of interest

There are no conflicts to declare.

Data availability

Additional data for this article, including input files, are available in the Zenodo repository at https://doi.org/10.5281/zenodo.19133513. The data supporting this article have been included as part of the supplementary information (SI). See DOI: https://doi.org/10.1039/d6cp00995f.

Acknowledgements

We thank the Czech Science Foundation for the financial support (Grant No. 26-22810S). L. P. and J. J. are students of the International Max Planck Research School Quantum Dynamics and Control. This work was supported by the project “The Energy Conversion and Storage”, funded as project No. CZ.02.01.01/00/22_008/0004617 by Programme Johannes Amos Comenius, call Excellent Research.

References

1. L. Vítek and J. D. Ostrow, Bilirubin chemistry and metabolism; harmful and protective aspects, Curr. Pharm. Des., 2009, 15, 2869–2883.
2. M. Maisels, Sister Jean Ward, phototherapy, and jaundice: a unique human and photochemical interaction, J. Perinatol., 2015, 35, 671–675.
3. A. S. Tatikolov, P. G. Pronkin and I. G. Panova, Bilirubin: Photophysical and photochemical properties, phototherapy, analytical methods of measurement. A short review, Biophys. Chem., 2025, 318, 107378.
4. R. Sloper and T. Truscott, The quantum yield for bilirubin photoisomerisation, Photochem. Photobiol., 1982, 35, 743–745.
5. B. Zietz and T. Gillbro, Initial photochemistry of bilirubin probed by femtosecond spectroscopy, J. Phys. Chem. B, 2007, 111, 11997–12003.
6. C. Carreira-Blanco, P. Singer, R. Diller and J. L. P. Lustres, Ultrafast deactivation of bilirubin: dark intermediates and two-photon isomerization, Phys. Chem. Chem. Phys., 2016, 18, 7148–7155.
7. J. Janoš, D. Madea, S. Mahvidi, T. Mujawar, J. Švenda, J. Suchan, P. Slavíček and P. Klán, Conformational control of the photodynamics of a bilirubin dipyrrinone subunit: Femtosecond spectroscopy combined with nonadiabatic simulations, J. Phys. Chem. A, 2020, 124, 10457–10471.
8. D. Madea, S. Mahvidi, D. Chalupa, T. Mujawar, A. Dvořák, L. Muchová, J. Janoš, P. Slavíček, J. Švenda and L. Vítek, Wavelength-dependent photochemistry and biological relevance of a bilirubin dipyrrinone subunit, J. Org. Chem., 2020, 85, 13015–13028.
9. D. Madea, T. Mujawar, A. Dvořák, K. Pospíšilová, L. Muchová, P. Čubáková, M. Kloz, J. Švenda, L. Vítek and P. Klán, Photochemistry of (Z)-Isovinylneoxanthobilirubic Acid Methyl Ester, a Bilirubin Dipyrrinone Subunit: Femtosecond Transient Absorption and Stimulated Raman Emission Spectroscopy, J. Org. Chem., 2022, 87, 3089–3103.
10. R. Pu, Z. Wang, R. Zhu, J. Jiang, T. C. Weng, Y. Huang and W. Liu, Investigation of Ultrafast Configurational Photoisomerization of Bilirubin Using Femtosecond Stimulated Raman Spectroscopy, J. Phys. Chem. Lett., 2023, 14, 809–816.
11. G. Granucci, M. Mazzoni, M. Persico and A. Toniolo, A computational study of the excited states of bilirubin IX, Phys. Chem. Chem. Phys., 2005, 7, 2594–2598.
12. A. F. McDonagh, G. Agati, F. Fusi and R. Pratesi, Quantum yields for laser photocyclization of bilirubin in the presence of human serum albumin. Dependence of quantum yield on excitation wavelength, Photochem. Photobiol., 1989, 50, 305–319.
13. E. C. Niederhoffer, J. H. Timmons and A. E. Martell, Thermodynamics of oxygen binding in natural and synthetic dioxygen complexes, Chem. Rev., 1984, 84, 137–203.
14. D. Madea, J. Penakova, J. Mehara, R. Akisaka, M. Martinek, J. Roithova and P. Klan, Photooxidation of Dipyrrinones: Reaction with Singlet Oxygen and Characterization of Reaction Intermediates, J. Org. Chem., 2025, 90, 2403–2420.
15. G. Y. Fraikin, Photosensitized Reactions of Oxidative Damage to Biomolecules: Role in Genotoxic and Cytotoxic Processes, Moscow Univ. Biol. Sci. Bull., 2024, 79, 115–129.
16. C. C. Winterbourn, Biological chemistry of superoxide radicals, ChemTexts, 2020, 6, 7.
17. G. Glotz, J. Polena, N. M. Khan, A. Mukherjee, M. Kloz, P. Slavíček and P. Klán, The first microseconds of the life of excited heptamethine cyanine revealed by femtosecond stimulated Raman spectroscopy, Commun. Chem., 2025, 9, 43.
18. J. Janoš and P. Slavíček, What controls the quality of photodynamical simulations? Electronic structure versus nonadiabatic algorithm, J. Chem. Theory Comput., 2023, 19, 8273–8284.
19. J. M. Toldo, M. T. Do Casal, E. Ventura, S. A. Do Monte and M. Barbatti, Surface hopping modeling of charge and energy transfer in active environments, Phys. Chem. Chem. Phys., 2023, 25, 8293–8316.
20. Y. Guo, C. Riplinger, U. Becker, D. G. Liakos, Y. Minenkov, L. Cavallo and F. Neese, Communication: An improved linear scaling perturbative triples correction for the domain based local pair-natural orbital based singles and doubles coupled cluster method [DLPNO-CCSD (T)], J. Chem. Phys., 2018, 148, 011101.
21. C. Angeli, R. Cimiraglia and J. P. Malrieu, n-electron valence state perturbation theory: A spinless formulation and an efficient implementation of the strongly contracted and of the partially contracted variants, J. Chem. Phys., 2002, 117, 9138–9153.
22. B. Kaduk, T. Kowalczyk and T. Van Voorhis, Constrained Density Functional Theory, Chem. Rev., 2012, 112, 321–370.
23. M. A. Marques and E. K. Gross, Time-dependent density functional theory, Annu. Rev. Phys. Chem., 2004, 55, 427–455.
24. R. Schmidt, Photosensitized Generation of Singlet Oxygen, Photochem. Photobiol., 2006, 82, 1161–1177.
25. D. Escudero, Revising intramolecular photoinduced electron transfer (PET) from first-principles, Acc. Chem. Res., 2016, 49, 1816–1824.
26. R. Sure and S. Grimme, Comprehensive benchmark of association (free) energies of realistic host-guest complexes, J. Chem. Theory Comput., 2015, 11, 3785–3801.
27. E. Epifanovsky, A. T. Gilbert, X. Feng, J. Lee, Y. Mao, N. Mardirossian, P. Pokhilko, A. F. White, M. P. Coons and A. L. Dempwolff, et al., Software for the frontiers of quantum chemistry: An overview of developments in the Q-Chem 5 package, J. Chem. Phys., 2021, 155, 084801.
28. H. X. Zhou and M. K. Gilson, Theory of free energy and entropy in noncovalent binding, Chem. Rev., 2009, 109, 4092–4107.
29. F. Neese, F. Wennmohs, U. Becker and C. Riplinger, The ORCA quantum chemistry program package, J. Chem. Phys., 2020, 152, 224108.
30. D. Hollas, J. Suchan, J. Janoš, M. Ončák and P. Slavíček, ABIN. Zenodo, 2026 DOI:10.5281/zenodo.19740845.
31. S. Sršeň, J. Sita, P. Slavíček, V. Ladányi and D. Heger, Limits of the nuclear ensemble method for electronic spectra simulations: Temperature dependence of the (E)-azobenzene spectrum, J. Chem. Theory Comput., 2020, 16, 6428–6438.
32. J. Tomasi, B. Mennucci and R. Cammi, Quantum mechanical continuum solvation models, Chem. Rev., 2005, 105, 2999–3094.
33. B. Jagoda-Cwiklik, P. Slavíček, L. Cwiklik, D. Nolting, B. Winter and P. Jungwirth, Ionization of imidazole in the gas phase, microhydrated environments, and in aqueous solution, J. Phys. Chem. A, 2008, 112, 3499–3505.
34. E. Pluharova, P. Slavicek and P. Jungwirth, Modeling photoionization of aqueous DNA and its components, Acc. Chem. Res., 2015, 48, 1209–1217.
35. M. Valiev, E. J. Bylaska, N. Govind, K. Kowalski, T. P. Straatsma, H. J. J. Van Dam, D. Wang, J. Nieplocha, E. Aprà and T. L. Windus, et al., NWChem: A comprehensive and scalable open-source solution for large scale molecular simulations, Comput. Phys. Commun., 2010, 181, 1477–1489.
36. X. Zhu, K. C. Thompson and T. J. Martínez, Geodesic interpolation for reaction pathways, J. Chem. Phys., 2019, 150, 164103.
37. I. Fdez Galván, M. Vacher, A. Alavi, C. Angeli, F. Aquilante, J. Autschbach, J. J. Bao, S. I. Bokarev, N. A. Bogdanov and R. K. Carlson, et al., OpenMolcas: From source code to insight, J. Chem. Theory Comput., 2019, 15, 5925–5964.
38. J. Suchan, J. Janoš and P. Slavíček, Pragmatic approach to photodynamics: mixed Landau–Zener surface hopping with intersystem crossing, J. Chem. Theory Comput., 2020, 16, 5809–5820.
39. K. Giesbertz, E. Baerends and O. Gritsenko, Charge Transfer, Double and Bond-Breaking Excitations with Time-Dependent Density Matrix Functional Theory, Phys. Rev. Lett., 2008, 101, 033004.
40. A. Dreuw and M. Head-Gordon, Failure of time-dependent density functional theory for long-range charge-transfer excited states: the zincbacteriochlorin- bacteriochlorin and bacteriochlorophyll- spheroidene complexes, J. Am. Chem. Soc., 2004, 126, 4007–4016.
41. J. M. Herbert and A. Mandal, Time-dependent density functional theory, Jenny Stanford Publishing, 2022, pp. 361–404.
42. L. Xu, V. M. Freixas, F. Aleotti, D. G. Truhlar, S. Tretiak, M. Garavelli, S. Mukamel and N. Govind, Conical Intersections Studied by the Configuration-Interaction-Corrected Tamm-Dancoff Method, J. Chem. Theory Comput., 2025, 21, 3600–3611.
43. M. Ceriotti, G. Bussi and M. Parrinello, Nuclear quantum effects in solids using a colored-noise thermostat, Phys. Rev. Lett., 2009, 103, 030603.
44. N. J. Turro, Modern molecular photochemistry, University Science Books, 1991.
45. I. Golovanov and S. Zhenodarova, Quantitative Structure-Property Relationship: XXIII. Solubility of Oxygen in Organic Solvents, Russ. J. Gen. Chem., 2005, 75, 1795–1797.
46. E. S. Chen, W. E. Wentworth and E. C. M. Chen, The electron affinities of NO and O₂, J. Mol. Struct., 2002, 606, 1–7.
47. J. R. Pliego and J. M. Riveros, The cluster- continuum model for the calculation of the solvation free energy of ionic species, J. Phys. Chem. A, 2001, 105, 7241–7247.
48. M. A. Filatov, S. Baluschev and K. Landfester, Protection of densely populated excited triplet state ensembles against deactivation by molecular oxygen, Chem. Soc. Rev., 2016, 45, 4668–4689.
49. M. Stanitska, D. Volyniuk, B. Minaev, H. Agren and J. V. Grazulevicius, Molecular design, synthesis, properties, and applications of organic triplet emitters exhibiting blue, green, red and white room-temperature phosphorescence, J. Mater. Chem. C, 2024, 12, 2662–2698.
50. R. Weinkauf and J. Schiedt, Energetics of Photoinduced Electron Transfer in the Indole-O2 Cluster in Gas Phase: Possible Consequences for Photoexcited Tryptophan in Solution, Photochem. Photobiol., 1997, 66, 569–575.
51. F. Thorning, K. Strunge, F. Jensen and P. R. Ogilby, The complex between molecular oxygen and an organic molecule: modeling optical transitions to the intermolecular charge-transfer state, Phys. Chem. Chem. Phys., 2021, 23, 15038–15048.
52. B. F. Minaev, K. V. Mikkelsen and H. Ågren, Collision-induced electronic transitions in complexes between benzene and molecular oxygen, Chem. Phys., 1997, 220, 79–94.
53. M. Radon and K. Pierloot, Binding of CO, NO, and O₂ to heme by density functional and multireference ab initio calculations, J. Phys. Chem. A, 2008, 112, 11824–11832.

---