Ultrafast and long-time excited state kinetics of an NIR-emissive vanadium(III) complex II. Elucidating triplet-to-singlet excited-state dynamics

Abstract

We report the non-adiabatic dynamics of V^IIICl₃(ddpd), a complex based on the Earth-abundant first-row transition metal vanadium with a d² electronic configuration which is able to emit phosphorescence in solution in the near-infrared spectral region. Trajectory surface-hopping dynamics based on linear vibronic coupling potentials obtained with CASSCF provide molecular-level insights into the intersystem crossing from triplet to singlet metal-centered states. While the majority of the singlet population undergoes back-intersystem crossing to the triplet manifold, 1–2% remains stable during the 10 ps simulation time, enabling the phosphorescence described in Dorn et al. Chem. Sci., 2021, DOI: 10.1039/D1SC02137K. Competing with intersystem crossing, two different relaxation channels via internal conversion through the triplet manifold occur. The nuclear motion that drives the dynamics through the different electronic states corresponds mainly to the increase of all metal–ligand bond distances as well as the decrease of the angles of trans-coordinated ligand atoms. Both motions lead to a decrease in the ligand-field splitting, which stabilizes the interconfigurational excited states populated during the dynamics. Analysis of the electronic character of the states reveals that increasing and stabilizing the singlet population, which in turn can result in enhanced phosphorescence, could be accomplished by further increasing the ligand-field strength.

---

1 Introduction

The development of photoactive coordination compounds largely focuses on 4d or 5d metals because they often feature long-lived excited states that can be exploited in multiple applications. Such long-lived states are due to the large ligand-field splitting that increases the energy of metal-centered (MC) states, which are then inaccessible for relaxation.¹ In contrast, complexes based on 3d metals have weaker ligand-field splitting and typically show rapid non-radiative relaxation through MC states.² For practical applications, however, 4d/5d metals face the severe disadvantage that their rare natural occurrence makes them costly.³ Thus, significant effort is currently directed to the development of photoactive Earth-abundant 3d transition metal compounds, as low-cost alternatives, e.g. of photosensitizers.^4–6

So far, 3d-metal complexes of chromium(0),⁷ iron(III),^8,9 cobalt(III),¹⁰ and copper(I)¹¹ exhibit luminescence in the visible region of the light spectrum. At longer wavelengths, luminescence has been observed for homoleptic chromium(III) complexes using the ddpd,^12,13 H₂tpda¹⁴ and dqp¹⁵ ligands (ddpd = N,N′-dimethyl-N,N′-dipyridine-2-ylpyridine-2,6-diamine; H₂tpda = 2,6-bis(2-pyridylamino)pyridine; dqp = 2,6-di(quinolin-8-yl)pyridine), which emit at the beginning of in the near-infrared region of the light spectrum (NIR) at 775, 782, and 747 nm, respectively. Using the ddpd ligand and vanadium(III) in place of chromium(III), the emission wavelength can be shifted to considerably longer wavelengths. Indeed, for the homoleptic [V^III(ddpd)₂]³⁺ complex, luminescence was found in the NIR-II range around 1100 nm,¹⁶ paving the way for promising applications in different fields, from photovoltaics to bio-imaging.^17,18

Here we present the heteroleptic chlorido-derivative V^IIICl₃(ddpd)¹⁹ (see Fig. 1), another exciting candidate for photoluminescence in the NIR-II range, which is able to emit at wavelengths of 1111 and 1219 nm at room temperature. While in paper I (see ref. 20) the experimental spectroscopic and electrochemical characterization of V^IIICl₃(ddpd) aided by static CASSCF-NEVPT2 calculations is presented, this work focuses on the electronic excited states and their time-dependent non-adiabatic evolution.

Fig. 1 Structure of V^IIICl₃(ddpd) investigated in this work.

A mechanistic understanding of photoluminescence can help identifying the competing non-adiabatic relaxation pathways in a molecule. This can enable a better modulation of quantum yields, which ultimately leads to design principles for improving luminescence in complexes based on Earth-abundant transition metals. In this respect, non-adiabatic molecular dynamics simulations are indispensable because spectroscopic experiments cannot access details on the motion of the atoms of the system. However, simulations of coupled electronic and nuclear degrees of freedom in transition metal complexes are extremely demanding for several reasons. One is the size of the molecules, which prevents dynamical studies in full dimensionality including all nuclear degrees of freedom.^21–23 Another is the large number of closely lying (near-degenerate) electronic states that can participate in the photodynamics and need to be taken into account. In order to overcome these challenges, we recently combined the non-adiabatic trajectory surface hopping method SHARC^24,25 with a linear vibronic coupling (LVC) approach²⁶ to compute the electronic potentials.²⁷ In this way, the photodynamics of a number of transition metal complexes have been investigated,^16,28–30 where time-dependent density functional theory (TDDFT) was used to parametrize the LVC potentials. Given its positive trade-off between accuracy and computational cost, TDDFT is likely the most popular method to study transition metal complexes,^31,32 and it has also been applied in on-the-fly surface hopping studies of other transition-metal complexes.^33–37 However, TDDFT can become problematic for open-shell systems with degenerate ground states due to its single-reference nature. This is the case for pseudo-octahedral vanadium(III) complexes, such as V^IIICl₃(ddpd), which possesses a d² ground-state electron configuration.

Here we present full-dimensional non-adiabatic photodynamics simulations of V^IIICl₃(ddpd) using trajectory surface hopping molecular dynamics in solution based on complete-active space self-consistent field (CASSCF)-derived LVC potentials. We show that CASSCF is able to describe all low-lying MC states expected from ligand-field theory. The calculation of the entire set of MC states is only possible with a multi-configurational method, as demonstrated by a transition-density matrix analysis. This analysis also reveals the presence of an important higher-order excitation character in low-lying excited states with ligand contributions, which are missing in methods relying on single excitations such as TDDFT. This work emphasizes the need for going beyond single-reference methods to obtain even a qualitatively correct electronic structure in V^IIICl₃(ddpd). The demanding CASSCF-based nonadiabatic dynamics simulations uncover the population pathways from the triplet to the singlet excited MC states via intersystem crossing (ISC), providing new insight into the luminescent states of V^IIICl₃(ddpd) and strategies to enhance the luminescence properties of related 3d metal complexes.

2 Methods

2.1 Trajectory surface hopping simulations

The nonadiabatic dynamics of V^IIICl₃(ddpd) including all 123 normal modes were simulated using an LVC model within the SHARC surface hopping approach.^27,38 Trajectories were set up using 2000 initial coordinates and momenta from a Wigner distribution of the ground-state PES.^39,40 The initial excited states were selected stochastically within an energy range of 0.5 eV around the maximum of the second absorption band (3.0–3.5 eV) based on their oscillator strength.⁴¹ Accordingly, 2000 trajectories were started in the T₆ (664), T₇ (891), T₈ (346), T₉ (55), T₁₀ (24), T₁₁ (19), and T₁₂ (1). The trajectories were then propagated for 10 ps using a nuclear time step of 0.5 fs and an electronic time step of 0.02 fs within the local diabatization method.⁴² In total, this amounts to 20 ns of non-adiabatic dynamics carried out with the PySHARC driver,²⁷ which allows for efficient in-memory communication between a Python interface and the routines of SHARC. An energy-based decoherence correction with a constant of C = 0.1 a.u. was used.⁴³ During the surface hops, the kinetic energy was adjusted by re-scaling the velocity vectors. Surface hopping probabilities were approximated using the wave function overlaps described in ref. 44.

2.2 Linear vibronic coupling model

In an LVC model, the potential energy surfaces (PESs) are approximated in terms of the ground-state PES V₀ and the first-order (linear) vibronic coupling term W on the basis of the mass-frequency scaled normal mode coordinate Q

V(Q) = V₀(Q)1 + W(Q) (1)

The ground-state PESs are approximated as harmonic oscillators with frequencies ω_i, i.e.,

The coupling term reads

where ε_n are the vertical excitation energies at the Franck–Condon (FC) geometry, while κ_i⁽ⁿ⁾ and λ_i^(n,m) are the intrastate and interstate coupling elements for the normal mode coordinate Q_i. The coupling elements were obtained numerically from CASSCF calculations (see below) on geometries displaced by ±0.05 units from the ground-state equilibrium geometry for each normal mode. The intrastate coupling κ_i⁽ⁿ⁾ values were obtained as numerical gradients and the interstate couplings λ_i^(n,m) were obtained from the change in the wavefunction overlaps. The spin–orbit couplings (SOCs) in the LVC model potential were approximated by constant SOCs obtained at the ground-state equilibrium geometry.

The absorption spectrum and the density-of-states were calculated from a Wigner distribution of 20 000 geometries using the LVC model potential.

2.3 Electronic structure calculations

The triplet ground state of V^IIICl₃(ddpd) in acetonitrile (MeCN) was optimized using unrestricted B3LYP,^45,46 the def2-TZVP basis set,⁴⁷ and the D3 dispersion correction using Becke–Johnson damping (D3BJ)⁴⁸ as implemented in the ORCA4.1 program package.⁴⁹ A numerical frequency calculation at the same level of theory confirms that no imaginary frequencies are obtained. Solvent effects were modeled using the conductor polarizable continuum model (C-PCM). For the self-consistent-field (SCF) calculations, the resolution-of-identity approximation (RIJCOSX), tight SCF convergence criteria (TightSCF), and the Grid4 integration grid were used. For the geometry optimization, tight convergence criteria (TightOpt) were used. Scalar relativistic effects are taken into account with the zeroth-order regular approximation⁵⁰ (ZORA). Cartesian coordinates are given in Table S1 in the ESI.†

The electronic excited states of V^IIICl₃(ddpd) in MeCN were calculated with CASSCF⁵¹ and the ANO-RCC-VDZP basis set^52,53 using the v18.0.o180122-0800 version of the OpenMolcas program package.⁵⁴ Thereby, 16 triplet and 15 singlet states were calculated at the state-averaged CASSCF level of theory with equal weights. An active space comprising 10 electrons in 13 orbitals was used (10,13). These orbitals correspond mainly to the five metal d orbitals, four ligand π orbitals, and four ligand π* orbitals; they are shown in Fig. S1 in the ESI.† A closer analysis also shows contributions from the p orbitals of the chlorido ligands in the active orbitals, which result in a small participation of the chlorido ligands in the charge flow during the excitation. Solvent effects of MeCN were described using the C-PCM model. SOCs between all states were calculated using the atomic mean field (AMFI) approximation⁵⁵ in the restricted active space state interaction (RASSI) program in OpenMolcas.

In addition, to assess the performance of the CASSCF results, the energies of the 16 triplet and 15 singlet states were also calculated with the multi-state complete active space second-perturbation theory (MS-CASPT2)⁵⁶ using an imaginary level shift of 0.1 a.u. (ref. 57) and IPEA shift values^58,59 of 0.0 a.u. and 0.25 a.u. The corresponding results are discussed in Section S3 in the ESI.†

2.4 Quantitative wavefunction analysis

Electronically excited states are often described in terms of the canonical orbitals that characterize the most important configurations in the excited-state wave functions.⁶⁰ In the presence of many configurations, this description becomes tedious and often misleading. A more comprehensible and quantitative analysis of the wavefunctions can be obtained by analyzing the one-particle transition-density matrix.⁶¹ For electronic structure methods including only single excitations, the one-particle transition-density matrix contains the full information of the excited-state wave functions. For electronic structure methods including higher excitation levels, such as CASSCF, this description is only approximate. A singular-value decomposition of the transition-density matrix provides a set of hole and electron natural transition orbitals that give a compact representation of the excited-state wave function – particularly in cases where a single pair of hole and electron natural-transition orbitals represents ∼99% of the character of the excitation. This orbital representation is convenient when analyzing the electronic structure at a single geometry. For multiple geometries, it is not possible to compute meaningful averaged orbitals in a straightforward manner. In these cases, it is possible to partition the molecule into several fragments and compute the contributions of the hole and electron parts of the fragments to the full transition density matrix. The fragments can comprise several atoms, e.g., a ligand, or correspond to individual atoms; importantly, their representation is invariant to the displacement of the nuclei. This procedure provides an easily understandable illustration of the charge flow that occurs during the dynamics in the different electronic states. The analysis of the transition-density matrix has been performed using the TheoDORE routine⁶² integrated in the WFA module in OpenMolcas.⁵⁴

3 Results and discussion

3.1 Electronic structure calculations

Fig. 2 shows the triplet (a) and singlet (b) electronic states of V^IIICl₃(ddpd) calculated at the CASSCF(10,13) level of theory. As the vanadium metal center possesses a d² ground-state electron configuration, the ground state of V^IIICl₃(ddpd) is a triplet state, T₀. Only slightly above T₀, we find two further triplet states, T₁ at 0.04 eV (332 cm⁻¹) and T₂ 0.10 eV (799 cm⁻¹). Since in V^IIICl₃(ddpd) the vanadium center has an almost octahedral coordination sphere, the three states T₀–T₂ correspond to the three components of the ³T₁ term – the ground state term predicted by ligand-field theory (see the Tanabe–Sugano diagram in Fig. S2 in the ESI†). The CASSCF(10,13) natural orbital occupation numbers (Table S4 in the ESI†) show that in all three states the d electrons occupy two of the d_xy, d_yz, and d_xz orbitals, and so they can be assigned the (t_2g)² configuration. To avoid confusion between the T_n nomenclature of triplet states ordered by energy and the ³T_n term symbol, here we also label these three states as ³MC₁ states. Thereby, MC refers to the metal-centered character of these states, and we implicitly retain the threefold spatial degeneracy of the ³T₁ term.

Fig. 2 Characters and energies of the (a) triplet and (b) singlet states calculated with CASSCF(10,13). Note that the energy axis is not on scale. ^3,1DL denotes d-assisted ligand excited states, while other state labels follow ligand-field terms of octahedral complexes.

The higher-lying triplet states are described likewise. Between 1.54 and 1.81 eV (12 460–14 577 cm⁻¹), we find the three components of the ³T₂ term of the (t_2g)¹(e_g)¹ configuration, that will be referred to as (threefold nearly-degenerate) ³MC₂ states. At energies of 3.16–3.53 eV (25 505–28 507 cm⁻¹), there are four nearly degenerate states. Three of these states (T₆–T₈) also possess the (t_2g)¹(e_g)¹ configuration: these states can be assigned to the term (threefold nearly-degenerate ³MC₃ states). The fourth state, T₉, possesses a (e_g)² configuration; it is assigned to the ³A₂ term (non-degenerate ³MC₄ state). Finally, we find two sets of three states each, T₁₀–T₁₂ between 4.38 and 4.47 eV (35 303–36 046 cm⁻¹)) and T₁₃–T₁₅ between 4.53 and 4.62 eV (36 546–37 285 cm⁻¹), respectively. Based on the natural orbital occupation numbers (Table S4 in the ESI†), these states correspond to the transition of an electron from a ddpd π orbital to a ddpd π* orbital. However, these states are actually a linear combination of equal parts of π → π* single excitations and (π, d) → (π*, d) double excitations, as becomes apparent when analyzing the one-particle transition density (see discussion below). Thus, we label these sets of states as the d-assisted ligand-centered states ³DL₁ and ³DL₂, respectively.

The S₁–S₁₅ singlet states are shown in Fig. 2(b). As can be seen, there is a group of five states (S₁–S₅) at energies between 1.53 and 1.66 eV (12 329–13 362 cm⁻¹). These states all share a (t_2g)² electron configuration. They correspond to the three components of the ¹T₂ term and the two components of the ¹E term, which are predicted as the lowest-energy singlet states in the Tanabe–Sugano diagram (Fig. S2 in the ESI†), where they are almost degenerate. In this case, it is harder to distinguish which of the states S₁–S₅ corresponds to each of the two terms due to the small energy gap (0.13 eV, 1033 cm⁻¹) as well as their shared electron configuration. Therefore, and for the sake of simplicity, we denote these five states only as ¹MC₁ states.

At higher energy, a group of four states appears between 3.25 and 3.49 eV (26 221–28 167 cm⁻¹). Three of them (S₆, S₈, and S₉) share the (t_2g)¹(e_g)¹ electron configuration. The remaining S₇ state is described by a (t_2g)^1.65(e_g)^0.35 configuration. For this state, the non-integer occupation numbers of the natural orbitals originate from the multi-configurational nature of the CASSCF wave function. We note that the natural orbital occupation numbers for all other electronic states lie much closer to 0, 1, or 2, i.e., usually differing by less than 0.1 from these integer values. Following their different electronic configurations, we ascribe the S₆, S₈, and S₉ states to the ¹T₂ term (threefold nearly-degenerate ¹MC₂ states), while S₇ is assigned to the ¹A₁ term (non-degenerate ¹MC₃ state).

Between 3.71 and 3.91 eV (29 913–31 509 cm⁻¹), we find a further group of three singlet states that share the (t_2g)¹(e_g)¹ configuration, which we assign to the ¹T₁ term (threefold nearly-degenerate ¹MC₄ states). Finally, we reach a group of three singlet states between 4.39 and 4.48 eV (35 429–36 149 cm⁻¹) with similar features as the ³DL_1,2 triplet states. Based on the natural orbital occupation numbers, these states correspond to the π → π* excitation; however, a closer analysis reveals additional (π, d) → (π*, d) double excitation character, which will be discussed in more detail below.

The analysis of the CASSCF natural orbitals provides a detailed, yet only qualitative description of the electronic states of V^IIICl₃(ddpd). One strategy to obtain a quantitative description is to analyze the transition density matrix (TDM) γ_IJ, which describes the change in electron density associated with a transition from the initial state Ψ_I—here the ground state Ψ₀—to the final state Ψ_J. This can be easily done for electronic structure methods including only single excitations, such as TDDFT, as in this case one needs to compute only the one-particle TDM γ_IJ⁽¹⁾. For methods like CASSCF that include higher-order excitations, one should also analyze higher-order TDMs γ_IJ⁽ⁿ⁾. Unfortunately, there is no computer program available yet that offers this possibility. Thus, here we are content with scrutinizing the one-particle TDM γ_IJ⁽¹⁾ of the CASSCF electronic states to infer the limitations of the single-excitation character in the different electronic states. For simplicity, henceforth, we drop the superscript and refer to the one-particle TDM just by γ_IJ.

In order to verify to which extent the analysis of the one-particle TDM is appropriate for the CASSCF electronic states, we first determined the sum of the single excitations ∑Ω of the CASSCF wave functions with respect to the ground state of V^IIICl₃(ddpd), i.e., the T₀ state, which is the lowest-energy component of the nearly-degenerate ³MC₁ states. This quantity, ∑Ω for γ_0J, for all the triplets is shown in Fig. 3(a) in blue. Remarkably, only for the states T₁, T₄, and T₇, the sum of single excitations ∑Ω is close to 1. For two additional states (T₂ and T₃), ∑Ω still lies between 80 and 90%; however, for the remaining triplet states, ∑Ω is significantly smaller than one. Clearly, only a small subset of states can be accurately described as single excitations from the T₀ reference state.

Fig. 3 (a) Sum of single excitation character ∑Ω for the triplet states calculated from the transition density matrices γ_0J, γ_1J, and γ_2J of the T₀, T₁, and T₂ states, respectively. Note that ∑Ω is set to 1 for γ_II to indicate that a state can be described quantitatively using itself as a reference when taking into account only single excitations. (b) Fragment partition and classification of excited-state characters in V^IIICl₃(ddpd). (c) Transition density matrix analysis of the triplet states. Analysis used the transition-density matrix with the largest ∑Ω out of γ_0J, γ_1J, and γ_2J. (d) Characterization of DL excited states.

We recall that based on the natural orbital occupation numbers, the ground state of V^IIICl₃(ddpd) adopts the (t_2g)² configuration, while a number of excited states (³MC_2,3) possess the (t_2g)¹(e_g)¹ configuration. At first glance, the latter states seem to be the result of promoting a single electron from a t_2g orbital to an e_g orbital. However, as revealed by the one-particle TDM analysis, their wave functions must also include higher-order excitations, leading to ∑Ω values significantly smaller than one. Yet, it is possible to describe the (t_2g)¹(e_g)¹ excited states mainly by single excitations if a different reference state other than the T₀ state is chosen. In particular, this can be achieved by selecting the T₁ and T₂ states, which are the complementary components of the triply nearly-degenerate ³MC₁ term. This is shown by the orange and red bars in Fig. 3(a) that denote the TDM γ_1J and γ_2J, which use the T₁ and T₂ states as a reference, respectively. Taking together all three one-particle TDMs γ_IJ (I = 0, 1, and 2), we now see that each of the states T₀–T₈ can be described by single excitations with ∑Ω ≈ 1 for at least one of the reference states (at least one of the three bars is close to 1). Thus, each of the states T₀–T₈ is reasonably described by the one-particle TDM, if the appropriate reference state is selected. For the T₉ state, we see in Fig. 3(a) that ∑Ω for all three γ_IJ is close to zero. This is because the T₉ state (³MC₄) corresponds to the (e_g)² configuration, which naturally cannot be described by single excitations from (t_2g)².

For the states T₁₀–T₁₅ (³DL_1,2), ∑Ω amounts only to ca. 0.5 when using the appropriate reference state. Based on the natural orbital occupation numbers, these states feature the (t_2g)² configuration at the vanadium center and the π → π* excitation at the ddpd ligand. However, these states do not correspond to the simple π → π* excitation as, otherwise, ∑Ω would amount to 1 for one of the TDMs. In order to clarify the character of the states T₁₀–T₁₅, we use charge transfer numbers⁶³ and partition V^IIICl₃(ddpd) into three fragments (Fig. 3(b)): the metal center M, the ddpd ligand (ligand L1), and the three chlorido ligands (ligand L2). The CASSCF configurations can then be assigned to one of the nine different types of excitations sketched in Fig. 3(b): local excitations – either centered at the metal atom (MC) or at one of the two ligands (L1C, L2C) – metal-to-ligand charge transfer excitations (ML1CT, ML2CT), ligand-to-metal charge-transfer excitations (L1MCT, L2MCT), and ligand-to-ligand excitations (L1L2CT, L2L1CT). Fig. 3(c) collects the resulting charge-transfer numbers for each electronic state, see also Tables S6–S8 in the ESI.† As can be seen, the majority of the single excitation character of the T₁₀–T₁₅ states corresponds to L1C excitations, i.e., excitations centered at the ddpd ligand. Combining this observation with the natural orbital occupation numbers leads to the conclusion that these states correspond to linear combinations of 50% π → π* single excitations and 50% (π, d) → (π*, d) double excitations, see Fig. 3(d). Accordingly, we label these states d-assisted ligand-centered (DL) states.

From the remaining states in Fig. 3(c) (the T₉ state is excluded as it corresponds to a double excitation), one can see that the states T₀–T₈ are MC excitations, with small admixtures of ML1CT and ML2CT characters. In passing, we note that the admixtures of chloride orbitals with the vanadium d orbitals can only be observed in the orbitals using an iso value lower than 0.05 (cf. Fig. S1(a and b) in the ESI†). This further highlights the benefits of using the quantitative TDM analysis compared to the qualitative orbital-based description.

In the accompanying paper I,²⁰ the higher-energy part of the absorption spectrum of V^IIICl₃(ddpd) has been calculated using LR-TDDFT. The TDDFT calculations predicted excited states with L1MCT and L1C characters following the manifold of MC states. Based on the double excitation character found for the ³DL states, the excitation energies of the L1C states are likely to be overestimated by LR-TDDFT – as LR-TDDFT includes only single excitations – while standard hybrid functionals such as the employed B3LYP functional also tend to underestimate the energies of CT states. It is also interesting to discuss whether other DFT-based methods such as spin-flip TDDFT⁶⁴ (SF-DFT) or a ΔSCF approach⁶⁵ should, in principle, be able to provide a better description of the electronic states of V^IIICl₃(ddpd) than “standard” LR-TDDFT. While SF-DFT can be used to describe multi-configurational ground states of open-shell systems⁶⁴ by applying a LR-TDDFT approach from a spin-flipped reference state, this is only reasonable if the reference state itself can be described as a single configuration. As shown in Fig. 2, this is problematic due to the near-degeneracy of the excited states, both among the triplet and the singlet manifold. Even if it was possible to choose one of the two excited states with rather unique electronic configurations, i.e., the S₇ state or the T₉ state, the subsequent success in SF-DFT calculations would be questionable. Due to the (e_g)² configuration of the T₉ state, SF-DFT in a linear-response framework would not be capable of describing any of the ground-state components T₀, T₁, and T₂. Furthermore, the (t_2g)^1.65(e_g)^0.35 configuration of the S₇ state highlights its multi-reference character, making it inaccessible to single-reference methods. In ΔSCF methods, it is possible to compute single components of spatially degenerate excited states of different symmetry.⁶⁵ However, as soon as this degeneracy is lifted, e.g., by moving away from a highly-symmetric reference geometry, these states become inaccessible as well. Thus, we conclude that the correct description of the electronic states of V^IIICl₃(ddpd) requires multi-configurational methods including higher-order excitations.

3.2 Excited-state dynamics simulations

The time evolution of the diabatic electronic states, i.e., states that are characterized by the same electronic configuration, is shown in Fig. 4(a). Since in the LVC model no intruder states can enter the dynamics simulation at any geometry, all electronic states can be described on the basis of the diabatic states of the equilibrium geometry at all simulation times. For the following analysis, we combine the population of the different components of the nearly degenerate ³MC₁, ³MC₂, and ³MC₃ states into one population each. In addition, the populations of all singlet states were combined into one (¹MC); individual diabatic singlet state populations can be found in Fig. S5 in the ESI.† For completeness, we also show the time evolution of the adiabatic states in Fig. S4 in the ESI.†

Fig. 4 (a) Time-evolution of the diabatic electronic state populations. Thin lines represent state populations and thick lines represent fitted curves. (b) Time-resolved diabatic state character of 50 sample trajectories. Colors indicate the dominant character, summing up to at least 60% of the total population for each trajectory (this threshold allowed an average assignment of 95%, i.e. 19 000 of the 20 000 time steps for each trajectory; see Fig. S7† for a more strict threshold of 70%, which accounts for 87% of the simulation time). Blank spots indicate that no dominant state could be assigned at a particular time step. (c) Kinetic model used to fit the time evolution of the populations. (d) Fitted time constants including error estimates obtained using the bootstrap method with 100 copies.⁶⁶

At t = 0 fs, the ³MC₃ state is populated by 95%, while the remaining population is in the ³MC₄ and ³DL₁ states (2.5% each). Within the first 1 ps, the population of the ³MC₃ state decreases to ca. 35%. Concomitantly, the population in the ³MC₄ state increases up to ca. 45% until 300 fs and then decreases again. The population of ³MC₂ also increases – albeit at a slower rate – reaching ca. 40% after 1 ps. Additionally, the ³MC₁ and ³DL₁ states are also populated in small amounts. Most importantly, the population is also transferred to the singlet states (¹MC), reaching ca. 10% after 100 fs and remaining constant until 1 ps. At later times, all states depopulate strongly except ³MC₁, illustrating how important it is to simulate dynamics beyond 1 ps. After 10 ps, the population of ³MC₁ has reached ca. 80%, and the remaining ca. 20% is mainly in the ³MC₂ state. The population of the singlet ¹MC states decreases from its peak of 10% to values of 1–2% after 4 ps. Noticeably, the singlet population is almost entirely in the lowest singlet state ¹MC₁, as shown in Fig. S5 in the ESI.†

The time-evolution of the diabatic populations of 50 random trajectories is shown in Fig. 4(b) up to 1 ps (see also Fig. S6 in the ESI† for 10 ps). This allows us to easily follow the individual pathways of the trajectories through the different electronic states and propose a kinetic model to explain their time evolution. We see that trajectories in the ³MC₃ state (dark blue) tend to change to either the ³MC₄ state (light blue), the ³DL₁ state (orange), or the singlet ¹MC state (purple). From either of these states, trajectories then usually transfer to the ³MC₂ state (green), from where they undergo internal conversion (IC) to the ³MC₁ ground state (red). The ³DL₂ states (brown) are rarely visited and thus henceforth disregarded. Based on these observations, we propose the following simplified kinetic model,

for which, fits for the populations are shown in Fig. S8 in the ESI.† As some of the fitted curves did not follow closely the state populations, we described the ³MC₃ → ³MC₄ and ³MC₂ → ³MC₁ reactions by a biexponential function, addingwhich improves the fit considerably: now all fit curves (thick lines in Fig. 4(a)) closely follow the population curves (thin lines). The complete kinetic model is shown in Fig. 4(c) and the fitted time constants are listed in Fig. 4(d).

Accordingly, IC via³MC₃ → ³MC₄ occurs with a fast time constant of τ₁ = 210 ± 20 fs and a slow time constant of . Competing with that is both, ISC to the ¹MC singlet states with a single time constant of τ_ISC = 1682 ± 318 fs and IC via the ³DL₁ states with a time constant of τ_a = 1393 ± 638 fs. The time constant of ISC is thereby similar to the time constant τ = 1.5 ps that describes the rise of a long-lived signal in the transient absorption experiments in the accompanying paper I,²⁰i.e., which describes the population of the singlet states via ISC. Based on the time constants fitted in the simulations, we can calculate the yields of the three deactivation pathways that start from the ³MC₃ state. In total, 79.2% of the excited-state population decays via³MC₃ → ³MC₄, 11.4% of the excited population decays via³MC₃ → ³DL₁, and 9.4% of the excited state population undergoes ISC to the ¹MC states.

From these results, it follows how one can increase ISC to the ¹MC states and stabilize the population in the singlet manifold, thus enhancing phosphorescence in compounds related to V^IIICl₃(ddpd). As discussed above, the main competing reaction to ³MC₃ → ¹MC ISC, that populates the singlet states, is IC via³MC₃ → ³MC₄. The ³MC₃ and ³MC₄ states are characterized by (t_2g)¹(e_g)¹ and (e_g)² configurations, respectively, while the majority of the singlet population is in the ¹MC₁ states described by the (t_2g)² configuration. The t_2g orbitals are non-bonding in ligand-field theory while the e_g orbitals possess an antibonding character. If the ligand-field strength was increased, it would lead to a destabilization of the e_g orbitals, which, in turn, would increase the energy of all electronic states where the e_g orbital is occupied. Thus, increasing the ligand-field strength would increase the energy of the ³MC₃ and ³MC₄ states while leaving the ¹MC₁ states unaffected. Furthermore, the increase would be more pronounced for the ³MC₄ state than for the ³MC₃ state, since in the ³MC₄ state, the e_g orbitals are occupied twice. This could lead to a situation where ³MC₄ is inaccessible during the dynamics, thus quenching the main competing pathway to ISC. Additionally, increasing the ligand-field strength can also stabilize the population in the singlet manifold, since the depopulation of the singlet states occurs via¹MC → ³MC₂. ³MC₂ is characterized by the (t_2g)¹(e_g)¹ configuration. Thus, it will be shifted to higher energies compared to ¹MC₁ and may become inaccessible for depopulation of the singlet states as well.

We now turn to discuss the internal coordinates that show most distinct changes during the course of the simulations. These coordinates were identified by analysing the average motion of all trajectories on the basis of the normal modes of the reference equilibrium geometry,^67,68 see also Section S4.6 in the ESI.† The identified coordinates include the distances V–N and V–Cl between the central vanadium atom and the ligating N and Cl atoms, as well as the angles X–V–Y between the trans-coordinated ligand atoms (X, Y = N, Cl). The average values of these internal coordinates from the 2000 trajectories are shown in Fig. 5(a–i) in blue. For comparison, the respective value at the equilibrium geometry is shown as a red line. Atom numbering is shown in Fig. 5(j). Note that the initial average values at t = 0 fs do resemble their reference values at the equilibrium geometry, as demonstrated by 0–200 fs zoom in Fig. S12 in the ESI.†

Fig. 5 Time-evolution of important internal coordinates: (a–c) V–N distances, (d–f) V–Cl distances and (g–i) trans angles. Average values of the 10 000 trajectories in blue, the average of trajectories in singlet states in green, the value at the equilibrium geometry in red, and the average value at the ISC hopping geometries in orange. (j) Atom numbering. (k) Superposition of structures at the ISC hopping geometries.

As seen in Fig. 5(a–c), the distances between the vanadium atom and the ligating N atoms in the ddpd moiety increase rapidly at the beginning of the dynamics from ca. 2.11–2.12 Å to 2.3–2.4 Å. The effect is most pronounced for the V–N₂ distance – the vanadium and the nitrogen atom trans to a chlorido ligand – which reaches values close to 2.5 Å. At later times, all V–N distances decrease again, however, ending up at 2.31 Å (V–N₂) and 2.25–2.27 Å (V–N_3,4), i.e., distances that are still considerably larger than in the equilibrium geometry. The evolution of the V–Cl distances, Fig. 5(d–f), also shows a rapid initial increase from values of 2.36–2.38 Å to the range of 2.4–2.6 Å. Thereby, the increase is smaller for the chlorido ligand that is coordinated trans to a ligating nitrogen atom (V–Cl₄₁) than for the other two chlorido ligands. However, in contrast to the V–N distances, the V–Cl distances increase in average with longer simulation times, reaching 2.49–2.52 Å after 10 ps. The X–V–Y angles between the trans-coordinating atoms in V^IIICl₃(ddpd) are shown in Fig. 5(g–i). The N₃–V–N₄ angle decreases in the beginning of the dynamics from 169° to 150–160°, while at longer simulation times it increases again to values above 160° – however, still below the ground-state reference value. For the Cl₄₁–V–N₂ and Cl₄₂–V–Cl₄₃ angles, the initial decrease with respect to the starting geometries is smaller. However, both angles continue to decrease slightly at later simulation times and reach values of 168 and 171° after 10 ps, i.e., also smaller than their reference values at the ground-state minimum-energy geometry which are close to 180°.

The observations of the increase of all metal–ligand bond distances and the decrease of the trans angles can be nicely rationalized by the nature of the electronic states accessed during the dynamics. Compared to the ³MC₁ ground state, that possesses the (t_2g)² configuration, the excited states ³MC₂, ³MC₃, and ³MC₄ (which account for the majority of the population during the excited-state dynamics) are described by (t_2g)¹(e_g)¹ and (e_g)² configurations. These configurations feature electrons in the antibonding e_g orbitals, which are stabilized at increased metal–ligand bond distances. At the same time, the energy of these excited states decreases when the system moves further away from the ideal octahedral coordination geometry, as the reduced overlap between the metal and ligand orbitals results in a smaller ligand-field strength that decreases the energy of these excited states. Accordingly, the system moves to longer metal–ligand bond distances and smaller trans angles during the dynamics.

Fig. 5(a–i) also display the average value of the respective internal coordinates at the ISC hopping geometries as an orange line. Such geometries include both, hops from triplet to singlet states as well as hops back from singlet to triplet states. A superposition of the ISC hopping geometries is shown in Fig. 5(k). As can be seen in Fig. 5(a–c), the average V–N distances at the ISC hopping geometries are 2.34 Å for V–N₂ and 2.28–2.29 Å for V–N₃ and V–N₄. These values are only reached by the ensemble of trajectories (blue curve) in the first few ps after the initial V–N bond stretching. Thus, during this time window, both ISC from the triplet to the singlet states as well as back-ISC from the singlet to the triplet states is possible. At longer simulation times, both types of ISC are quenched, which leads to the singlet population being stable after ca. 4 ps. The same behavior is observed for the N₃–V–N₄ angle in Fig. 5(g): ISC hops occur at an average angle of 157° – a value which is only realized in the ensemble of trajectories during the first few ps. At later simulation times, the N₃–V–N₄ angle increases, moving the system away from the ISC hopping zone and stabilizing the remaining singlet population. For the other trans angles shown in Fig. 5(h and i), the ensemble average stays close to the average of the ISC hopping geometries for longer simulation times, before finally moving away by ca. 4°. The V–Cl distances of the ensemble stay most of the time close to the ISC hopping averages. Only the V–Cl₄₁ ensemble distance increases towards the end of the simulation time to ca. 2.50 Å compared to the ISC hopping average of 2.46 Å.

Finally, it is also useful to analyze the evolution of the coordinates associated with the singlet states only, green curves in Fig. 5(a–i). Note that because the majority of the population is in the triplet states at all simulation times, the averages of the exclusive triplet trajectories (not shown) closely resemble the evolution of the full ensemble. Compared to the full ensemble of trajectories (blue curve), the main difference of the V–N distances and the N₃–V–N₄ angle in the singlet trajectories is noticeable during the first few ps: the V–N bonds are smaller and the N₃–V–N₄ angle is larger than in the triplet trajectories. This is due to the fact that the majority of the singlet population is in the ¹MC₁ states, which are characterized by the (t_2g)² configuration with no electrons in the e_g orbitals. Thus, in the singlet states, the V–N distances do not tend to increase much, while keeping the N₃–V–N₄ angle closer to the octahedral angle of 180°. At later times, the full ensemble curve – corresponding to ca. 98% of triplet trajectories – resembles the average of the singlet curves, as now more and more electronic population reaches the ³MC₁ triplet ground state that also possesses the (t_2g)² configuration.

4 Conclusion

Here we investigated the electronic excited states of the phosphorescent complex V^IIICl₃(ddpd) and its photodynamics using trajectory surface hopping and CASSCF potentials in 123 degrees of freedom. By analyzing the electronic wave functions and transition-density matrices, we evidence the complicated electronic structure of the MC excited states as well as the presence of excitations involving the organic ddpd ligand that requires multi-configurational methods including higher-order excitations.

Our dynamical study reveals that after excitation from the triplet ground state possessing a (t_2g)² configuration to a higher-lying MC state of (t_2g)¹(e_g)¹ character, the system evolves through different decay pathways. The main deactivation pathway involves IC to a different MC state with the (e_g)² configuration. Minor pathways include ISC to the singlet MC manifold and an additional IC pathway through triplet states involving excitations to the ddpd ligand. The time constant τ_ISC = 1.7 ± 0.3 ps obtained for the ISC in the present simulation is in good agreement with that of τ = 1.5 ps obtained for the rise of a long-lived signal from transient absorption experiments in the accompanying paper I.²⁰ All three mechanisms end in a further triplet MC state with the (t_2g)¹(e_g)¹ configuration, which eventually relaxes back to the triplet ground state.

The dynamics is driven by elongation of all metal–ligand bond distances and a decrease of the angle between all pairs of trans-coordinated ligands further away from the octahedral 180° ideal. Both features result in a decrease of the ligand-field strength, which stabilizes the excited triplet states visited during the dynamics. The ISC leads to a temporary 10% population in the singlet states, that subsequently decreases to 1–2% at later simulation times up to 10 ps. Based on the analysis of the electronic character of all states involved in the dynamics, it is suggested that the overall population of the singlet states can be increased by increasing the ligand-field strength through suitable modification of the ligands of V^IIICl₃(ddpd). Such modifications are expected to result both in a more favorable population of the singlet states compared to the competing pathway to the triplet MC states, as well as a stabilization of the population in the singlet states, thus leading to larger phosphorescence yields for derivatives of V^IIICl₃(ddpd).

Data availability

The data that support the findings of this study are available from the corresponding author upon reasonable request.

Author contributions

Conceptualization, J. P. Z. and L. G.; calculations, J. P. Z. and T. K.; analysis, J. P. Z.; writing—original draft preparation, J. P. Z.; writing—review and editing, L. G.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

The authors thank S. Mai and M. Heindl for fruitful discussion during the course of this work. This work was supported by Deutsche Forschungsgemeinschaft [DFG, Priority Program SPP 2102 “Light-controlled reactivity of metal complexes” (GO 1059/8-1)]. The authors thank the Vienna Scientific Cluster for allocation of computer time.

Notes and references

1. K. Kalyanasundaram and M. Grätzel, Coord. Chem. Rev., 1998, 77, 347–414.
2. G. Auböck and M. Chergui, Nat. Chem., 2015, 7, 629–633.
3. CRC Handbook of Chemistry and Physics, ed. D. R. Lide, CRC Press, 2005.
4. O. S. Wenger, J. Am. Chem. Soc., 2018, 140, 13522–13533.
5. B. Bozic-Weber, E. C. Constable and C. E. Housecraft, Coord. Chem. Rev., 2013, 257, 3089–3106.
6. C. Förster and K. Heinze, Chem. Soc. Rev., 2020, 49, 1057–1070.
7. L. A. Büldt, X. Guo, R. Vogel, A. Prescimone and O. S. Wenger, J. Am. Chem. Soc., 2017, 139, 985–992.
8. P. Chábera, Y. Liu, O. Prakash, E. Thyrhaug, A. El Nahhas, A. Honarfar, S. Essén, L. A. Fredin, T. C. B. Harlang, K. S. Kjær, K. Handrup, F. Ericson, H. Tatsuno, K. Morgan, J. Schnadt, L. Häggström, T. Ericsson, A. Sobkowiak, S. Lidin, P. Huang, S. Styring, J. Uhlig, J. Bendix, R. Lomoth, V. Sundström, P. Persson and K. Wärnmark, Nature, 2017, 543, 695–699.
9. K. S. Kjær, N. Kaul, O. Prakash, P. Chábera, N. W. Rosemann, A. Homarfar, O. Gordivska, L. A. Fredin, K.-E. Bergquist, L. Häggström, T. Ericsson, L. Lindh, A. Yartsev, S. Styring, P. Huang, J. Uhlig, J. Bendix, D. Strand, V. Sundström, P. Persson, R. Lomoth and K. Wärnmark, Science, 2019, 363, 249–253.
10. A. K. Pal, C. Li, G. S. Hanan and E. Zysman-Colman, Angew. Chem., Int. Ed., 2018, 57, 8027–8031.
11. R. Hamze, J. L. Peltier, D. Sylvinson, M. Jung, J. Cardenas, R. Haiges, M. Soleilhavoup, R. Jazzar, P. I. Djurovich, G. Bertrand and M. E. Thompson, Science, 2019, 363, 601–606.
12. S. Otto, M. Grabolle, C. Förster, C. Kreitner, U. Resch-Genger and K. Heinze, Angew. Chem., Int. Ed., 2015, 54, 11572–11576.
13. C. Wang, S. Otto, M. Dorn, E. Kreidt, J. Lebon, L. Sršan, P. Di Martino-Fumo, M. Gerhards, U. Resch-Genger, M. Seitz and K. Heinze, Angew. Chem., Int. Ed., 2018, 57, 1112–1116.
14. S. Otto, C. Förster, C. Wang, U. Resch-Genger and K. Heinze, Chem.–Eur. J., 2018, 24, 12555–12563.
15. J.-R. Jiménez, B. Doistau, C. M. Cruz, C. Besnard, J. M. Cuerva, A. G. Campaña and C. Piguet, J. Am. Chem. Soc., 2019, 141, 13244–13252.
16. M. Dorn, J. Kalmbach, P. Boden, A. Päpcke, S. Gómez, C. Förster, F. Kuczelinis, L. M. Carrella, L. A. Büldt, N. H. Bings, E. Rentschler, S. Lochbrunner, L. González, M. Gerhards, M. Seitz and K. Heinze, J. Am. Chem. Soc., 2020, 142, 7947–7955.
17. A. Zampetti, A. Minotto and F. Cacialli, Adv. Funct. Mater., 2019, 29, 1807623.
18. Z. Guo, S. Park, J. Yoon and I. Shin, Chem. Soc. Rev., 2014, 43, 16–29.
19. C. Förster, M. Dorn, T. Reuter, S. Otto, G. Davarci, T. Reich, L. Carrella, E. Rentschler and K. Heinze, Inorganics, 2018, 6, 86.
20. M. Dorn, J. Kalmbach, P. Boden, A. Kruse, C. Dab, C. Reber, G. Niedner-Schatteburg, S. Lochbrunner, M. Gerhards, M. Seitz and K. Heinze, Chem. Sci., 2021 10.1039/D1SC02137K.
21. K. Falahati, H. Tamura, I. Burghardt and M. Huix-Rotllant, Nat. Commun., 2018, 9, 4502.
22. M. Fumanal, E. Gindensperger and C. Daniel, J. Phys. Chem. Lett., 2018, 9, 5189–5195.
23. M. Pápai, G. Vankó, T. Rozgonyi and T. J. Penfold, J. Phys. Chem. Lett., 2016, 7, 2009–2014.
24. M. Richter, P. Marquetand, J. González-Vázquez, I. Sola and L. González, J. Chem. Theory Comput., 2011, 7, 1253–1258.
25. S. Mai, P. Marquetand and L. González, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2018, 8, e1370.
26. H. Köppel, W. Domcke and L. S. Cederbaum, Advances in Chemical Physics Vol. 57, John Wiley & Sons, 1984, pp. 59–246.
27. F. Plasser, S. Gómez, M. F. S. J. Menger, S. Mai and L. González, Phys. Chem. Chem. Phys., 2019, 21, 57–69.
28. S. Mai, M. F. S. J. Menger, M. Marazzi, D. L. Stolba, A. Monari and L. González, Theor. Chem. Acc., 2020, 139, 65.
29. J. P. Zobel, O. S. Bokareva, P. Zimmer, C. Wölper, M. Bauer and L. González, Inorg. Chem., 2020, 59, 14666–14678.
30. M. Heindl, J. Hongyan, S.-A. Hua, M. Oelschlegel, F. Meyer, D. Schwarzer and L. González, Inorg. Chem., 2021, 1672–1682.
31. C. J. Cramer and D. G. Truhlar, Phys. Chem. Chem. Phys., 2009, 11, 10757–10816.
32. S. Mai and L. González, Angew. Chem., Int. Ed., 2020, 59, 16832–16846.
33. A. J. Atkins and L. González, J. Phys. Chem. Lett., 2017, 8, 3840–3845.
34. X.-Y. Liu, Y.-H. Zhang, W.-H. Fang and G. Cui, J. Phys. Chem. A, 2018, 122, 5518–5532.
35. S. Mai and L. González, Chem. Sci., 2019, 10, 10405–10411.
36. Y.-G. Fang, L.-Y. Peng, X. Y. Liu, W.-H. Fang and G. Cui, Comput. Theor. Chem., 2019, 1155, 90–100.
37. F. Talotta, M. Boggio-Pasqua and L. González, Chem.–Eur. J., 2020, 26, 11522–11528.
38. S. Mai, M. Richter, M. Heindl, M. F. S. J. Menger, A. Atkins, M. Ruckenbauer, F. Plasser, L. M. Ibele, S. Kropf, M. Oppel, P. Marquetand and L. González, SHARC2.1: Surface Hopping Including Arbitrary Couplings –Program Package for Non-Adiabatic Dynamics, http://sharc-md.org, 2019.
39. M. Barbatti and R. Crespo-Otero, Top. Curr. Chem., 2016, 368, 415–444.
40. L. Sun and W. L. Hase, J. Chem. Phys., 2010, 133, 044313.
41. M. Barbatti, G. Grannuci, M. Persico, M. Ruckenbauer, M. Vazdar, M. Eckert-Maksić and H. Lischka, J. Photochem. Photobiol., A, 2007, 190, 228–240.
42. G. Granucci, M. Persico and A. Toniolo, J. Chem. Phys., 2001, 114, 10608–10615.
43. G. Granucci, M. Persico and A. Zoccante, J. Chem. Phys., 2010, 133, 134111.
44. F. Plasser, M. Ruckenbauer, S. Mai, M. Oppel, P. Marquetand and L. González, J. Chem. Theory Comput., 2016, 12, 1207–1219.
45. A. D. Becke, J. Chem. Phys., 1993, 98, 5648–5652.
46. C. Lee, W. Yang and R. G. Parr, Phys. Rev. B, 1988, 37, 785–789.
47. F. Weigend and R. Ahlrichs, Phys. Chem. Chem. Phys., 2005, 7, 3297–3305.
48. S. Grimme, S. Ehrlich and L. Goerigk, J. Comput. Chem., 2011, 32, 1456–1465.
49. F. Neese, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2012, 2, 73–78.
50. C. van Wüllen, J. Chem. Phys., 1998, 109, 392–399.
51. B. O. Roos, P. R. Taylor and P. E. M. Siegbahn, Chem. Phys., 1980, 48, 157–173.
52. B. O. Roos, R. Lindh, P.-Å. Malmqvist, V. Veryazov and P.-O. Widmark, J. Phys. Chem. A, 2004, 108, 2851–2858.
53. B. O. Roos, R. Lindh, P.-Å. Malmqvist, V. Veryazov and P.-O. Widmark, J. Phys. Chem. A, 2005, 109, 6575–6579.
54. I. F. Galván, M. Vacher, A. Alavi, C. Angeli, F. Aquilante, J. Autschbach, J. J. Bao, S. I. Bokarev, N. A. Bogdanov, R. K. Carlson, L. F. Chibotaru, J. Creutzberg, N. Dattani, M. G. Delcey, S. S. Dong, A. Dreuw, L. Freitag, L. M. Frutos, L. Gagliardi, F. Gendron, A. Giussani, L. González, G. Grell, M. Guo, C. E. Hoyer, M. Johansson, S. Keller, S. Knecht, G. Kovačević, E. Källman, G. L. Manni, M. Lundberg, Y. Ma, S. Mai, J. P. Malhado, P. Å. Malmqvist, P. Marquetand, S. A. Mewes, J. Norell, M. Olivucci, M. Oppel, Q. M. Phung, K. Pierloot, F. Plasser, M. Reiher, A. M. Sand, I. Schapiro, P. Sharma, C. J. Stein, L. K. Sørensen, D. G. Truhlar, M. Ugandi, L. Ungur, A. Valentini, S. Vancoillie, V. Veryazov, O. Weser, T. A. Wesołowski, P.-O. Widmark, S. Wouters, A. Zech, J. P. Zobel and R. Lindh, J. Chem. Theory Comput., 2019, 15, 5925–5964.
55. B. Schimmelpfennig, AMFI, An Atomic Mean-Field Spin-Orbit Integral Program, University of Stockholm, 1996.
56. J. Finley, P.-Å. Malmqvist, B. O. Roos and L. Serrano-Andrés, Chem. Phys. Lett., 1998, 288, 299–306.
57. N. Forsberg and P.-Å. Malmqvist, Chem. Phys. Lett., 1997, 274, 196–204.
58. G. Ghigo, B. O. Roos and P.-Å. Malmqvist, Chem. Phys. Lett., 2004, 396, 142–149.
59. J. P. Zobel, J. J. Nogueira and L. González, Chem. Sci., 2017, 8, 1482–1499.
60. M. Jäger, L. Freitag and L. González, Coord. Chem. Rev., 2015, 304–305, 146–165.
61. F. Plasser and H. Lischka, J. Chem. Theory Comput., 2012, 8, 2777–2789.
62. F. Plasser, J. Chem. Phys., 2020, 152, 084108.
63. S. Mai, F. Plasser, J. Dorn, M. Fumanal, C. Daniel and L. González, Coord. Chem. Rev., 2018, 361, 74–97.
64. Y. Shao, M. Head-Gordon and A. Krylov, J. Chem. Phys., 2003, 118, 4807–4818.
65. T. Ziegler, A. Rauk and E. J. Baerends, Theor. Chim. Acta, 1977, 43, 261–271.
66. S. Nangia, A. W. Jasper, T. F. Miller III and D. G. Truhlar, J. Chem. Phys., 2004, 120, 3586–3597.
67. K. Kurtz, A. Hofmann and R. de Vivie-Riedle, J. Chem. Phys., 2001, 114, 6151–6159.
68. F. Plasser, M. Barbatti, A. J. A. Aquino and H. Lischka, J. Phys. Chem. A, 2009, 113, 8490–8499.

---

Footnote

† Electronic supplementary information (ESI) available: Optimized ground-state geometry (Section S1), additional details of the CASSCF calculations (Section S2), comparison to benchmark CASPT2 calculations and experimental absorption spectra (Section S3), additional analysis of the surface hopping dynamics (Section S4), and an additional file containing the parameters of the LVC model as well as the ground-state normal modes (lvc.zip). See DOI: 10.1039/d1sc02149d

---