Deactivation of excited states in transition metal complexes: Insight from computational chemistry

Investigation of the excited state decay dynamics of transition metal systems is a crucial step for the development of photoswitchable molecular based materials with applications in growing fields as energy conversion, data storage or molecular devices. The photophysics of these systems is an entangled problem arising from the interplay of electronic and geometrical rearrangements that take place on a short time scale. Several factors play a role in the process: various electronic states of di↵erent spin and chemical character are involved, the system undergoes important structural variations and several nonradiative processes can occur. Computational chemistry is a useful tool to get insight into the microscopic description of the photophysics of these materials since it provides unique information about the character of the electronic spin states involved, the energetics and time evolution of the system. In this review article, we present an overview of the state of the art methodologies available to address the several aspects that have to be incorporated to properly describe the deactivation of excited states in transition metal complexes. The most recent developments in theoretical methods are discussed and illustrated with examples. 2 10.1002/chem.201801990 Ac ce pt ed M an us cr ip t Chemistry A European Journal This article is protected by copyright. All rights reserved.


Introduction
Understanding the dynamics of excited electronic states of molecular transition metal (TM) complexes is a paramount objective of current research in the fields of spectroscopy, photochemistry and photophysics. [1] The interest in these molecular-based systems is prompted by their potential applications as sensitizers of solar cells for energy conversion, phosphorescent dyes for organic light emitting diodes (OLED), photoswitchable devices for high-density data storage, luminiscence-based sensors, molecular electronic devices, photocatalysts or activators of electron transfer processes in biological systems. [2][3][4][5][6][7] After photoexcitation of the system into an electronic excited state having significant radiative coupling with the ground state, a sequence of interconnected processes takes place during the excited state evolution that leads to a relatively stable occupation of the relevant final excited state. [8] Several factors play a role in the dynamics of the excited state relaxation and control the deactivation to the proper electronic state. Among these factors are the spin, spatial symmetry and character of the electronic states involved in the process. Indeed, various electronic states of very di↵erent nature, occasionally lying in a narrow energy range, are involved. Those include states in which the excited electron is basically localized on the transition metal (metal centered, MC), centered on the ligand (LC), states involving charge separation, either charge transfer from the metal to the ligand (MLCT) or conversely, ligand to metal charge transfer states (LMCT) and ligand to ligand charge transfer (LLCT) states.
Charge transfer states are usually more accessible from the ground state because of their larger oscillator strength, but metal centered states can also become involved in the photocycle and they actually play an important role in the photochemistry of numerous complexes. [9] In many cases, the character of the excited states is readily accessed from a visual inspection of the orbitals and their occupation numbers. However, it is not uncommon that the molecular orbitals are strongly delocalized over ligand and metal or that the wave function adopts a high multiconfigurational character and that the character of the di↵erent electronic states is less obvious. Among the di↵erent (semi-)automatic alternatives for the visual inspection, we mention here the orthogonal valence bond interpretation of the multiconfigurational wave function, [10,11] the analysis in terms of so-called charge transfer numbers extracted from the one-electron transition density matrix recently presented by Mai et al., [12] and the density (di↵erence) based indexes as descriptors for the character of the excited state character by Ciofini and co-workers. [13] The photoinduced intramolecular electron transfer can be accompanied by important structural changes both in the molecular system and its nearest environment, which include cooperative e↵ects, solvation e↵ects or thermal distortions. The extent of the geometrical relaxation within the molecular complex as a result of the deactivation process relies on the particular system. To quote one typical example, Figure 1 shows how a quasi-octahedral complexes with a Fe II N 6 core, like the Fe II (bpy) 3 complex, undergoing a spin crossover transition from a singlet to a quintet state, experience an enlargement of the Fe-N distance of around 0.2Å, [14,15] because of the occupation with two electrons of the anti-bonding orbitals in the quintet state. Conversely, the isoelectronic complex Ru II (bpy) 3 only shows very small variations in the metal-ligand distance along the deactivation process, [16] where the triplet and quintet MC states lie higher in energy and the system remains in the 3 MLCT state.
A host of intramolecular nonradiative processes can potentially occur, such as intersystem crossing (ISC) between electronic states of di↵erent spin multiplicities, internal conversion (IC) between states of the same spin quantum number, and vibrational relaxation and/or intramolecular vibrational redistribution. Moreover, in many cases all these processes can be extremely fast, and even take place in a shorter timescale than molecular vibrations. [17,18] Approaching the study of these systems by experimental techniques has been possible thanks to the development in the last 15 years of ultrafast optical and X-ray spectroscopies. [19][20][21][22] These techniques give insight into both the electronic and geometrical structure of the photophysics of transition metal complexes, achieving a time resolution in the femtosecond scale.
In view of all these di↵erent characteristics that are at play in the excited state dynamics, it is clear that a theoretical description of photocatalysis, light-induced magnetism, phosphorescence, electron transfer, etc. puts very high requirements on the computational method to be used. In fact, there is not a unique state-of-theart methodology able to capture all the aspects in one shot. Since the pioneering  [23][24][25] important progress has been made in the development of new approaches to allow for accessing excited electronic states, [13,[26][27][28] spin-orbit coupling interactions and intersystem crossing processes, [29][30][31] and dynamic contributions. Several reviews concerning the ability of these methodologies together with several applications have been published in the last years. [32][33][34][35][36][37][38] In the present review, the di↵erent components needed for a proper description of the dynamics of excited states are described and the pertinent theoretical approaches discussed. A first crucial factor is the computation of the energy di↵erences between the various spin-states involved in the process. These include the absorption spectra from the initial ground state, the adiabatic energy di↵erence between initial and final states, and the variation of the energy of the relevant electronic states with the geometrical changes, that is, the potential energy surfaces (PES). Additionally, core-level excitations are the basis of many spectroscopic techniques and the accurate computational description can be of help for a correct interpretation of the measurements.
These questions will be addressed in Section 2. Since intersystem crossings between electronic states of di↵erent spin are expected along the potential energy surfaces, inclusion of spin-orbit (SO) coupling is mandatory. In Section 3 the importance of spin-orbit e↵ects and how to treat them will be discussed. Both direct spin-orbit interaction and higher-order contributions will be analyzed.
The aforementioned factors give a static description of the system, disregarding the e↵ect of time and temperature in the dynamic evolution of the excited state. Inclusion of temperature and time by molecular dynamics simulations allows to incorporate thermal geometrical distortions, while surface hopping techniques open the possibility to simulate conical intersections and hence calculate intersystem crossing rate constants. The interplay of nuclear movement and electron density can only be achieved from a full quantum mechanical treatment of both the nuclei and the electrons, which is available in quantum dynamics, such as the multi-configuration time-dependent Hartree (MCTDH) method. These approaches will be illustrated in Section 4.

The energetics
One of the most important ingredients for an accurate account of the photochemistry of transition metal complexes is the relative energy of the di↵erent electronic states involved in the deactivation cascade. There are several aspects that deserve a close inspection as graphically explained in Figure 2. In the first place, it is of fundamental interest to have a good estimate for the adiabatic energy di↵erence of the di↵erent electronic states relevant to the photophysical phenomenon under study. This often implies states with di↵erent spin moment and di↵erent equilibrium geometry and is shortly discussed in Sec. 2.1. In spin-crossover complexes this is known as the highspin low-spin energy di↵erence, E 0 HL . Photo-induced processes are triggered by the absorption of photons, and hence, vertical excitation energies ( E F C , FC=Franck-Condon) are relevant to get information about the initial population of excited states (Sec. 2.2). The deactivation of these excited states is largely determined by the potential energy surfaces on which the nuclei can move to adapt the nuclear configuration to the change in the electron distribution. In addition to minima on the excited state energy surface, other points are of special interest such as those that locate conical

Adiabatic energy di↵erence
Ever since the first computational studies of the spin crossover (SCO) phenomenon, [39,40] there has been a constant flow of publications in which di↵erent computational schemes are tested for their ability to reproduce the adiabatic energy di↵erence between HS and LS state, E 0 HL . Usually, the E 0 HL is a small quantity with typical values of less than 2000 cm 1 for SCO systems. Combined with the fact that the geometry of the two spin states is usually quite di↵erent, the accurate computation of this parameter is a hard task and the inclusion of the zero-point energy correction is mandatory.
In most cases, attention is focused on the singlet-quintet energy di↵erence in the quasi-octahedral complexes with a Fe II N 6 core, but it goes without saying that other transition metals, other oxidation states and other coordination modes have also been looked at. Furthermore, adiabatic energy di↵erences are also essential in the study of luminescent materials. Analogous to the HS-LS energy di↵erence, one also has to locate the optimal geometry of the lowest emissive state and precisely determine the energy di↵erence with the ground state.
From the earlier density functional theory (DFT) studies in which nearly all existing density functionals have been tested, [41][42][43][44] it can be concluded that E 0 HL strongly depends on the functional of choice but that there are a few that perform remarkably well and reproduce the experimental data in most cases. One of the key factors for correct relative energies of HS and LS states is the amount of exact Fock exchange in the density functional. This was recognized by Reiher [45] and based on the observation that the standard B3LYP functional has a clear tendency to overstabilize the HS state, the weight of the Fock exchange was reduced to 10%. The resulting functional, known as B3LYP*, has been successfully applied in studies of many transition metal complexes. [46][47][48] A second widely used functional is the TPPSh functional. This hybrid variant of the meta-GGA TPSS functional was triggered as the most accurate functional in a careful analysis of the di↵erent factors that play a role in the relative stability of di↵erent spin states. [49] Apart from the electronic energy di↵erence, the author also included zero-point energy corrections, entropy e↵ects, dispersion corrections and (scalar) relativistic e↵ects in the final stability comparison of high-spin and low-spin states. Not considering these e↵ects may seriously a↵ect the outcome of the calculations and lead to wrong conclusions about the performance of the functional for predicting the HS-LS stability. The disadvantage of these two functionals is that they belong to the so-called class of hybrid functionals and hence require the evaluation of the exact Fock-exchange. Since this is a relatively costly operation, the application of these functionals to large systems could become cumbersome. Moreover, analytical gradients for the meta-GGA functionals are not available in all standard computational packages. As an alternative, Swart explored di↵erent combinations of standard GGA exchange and correlation functionals to come up with the OPBE (combination of the OPTX exchange and the PBE correlation functionals) as one of the best performing pure functionals for the calculation of the relative energies of di↵erent spin states in transition metal complexes. [50,51] The reparametrization of the OPBE functional to handle weak interactions (dispersion) in a more accurate manner gave rise to the S12g functional, which can be considered to be one of the most complete GGA functionals for dealing with spin states in transition metal complexes. [52,53] A somewhat di↵erent route towards an accurate yet computationally e cient method was taken by Vela and co-workers, [54] who applied the GGA+U method after having carefully benchmarked the U-value against experimental HS-LS energy di↵erences.
Apart from the calculations based on density functional theory, the relative stability of spin states can also be addressed with computational schemes that use the N -electron wave function as central entity. Among the many possible wave function based approaches the CASSCF/CASPT2 method has emerged as one of the most accurate schemes. The complete active space self-consistent field (CASSCF) wave function is constructed by performing a full configuration interaction in a small set of (valence) orbitals to capture the main static electron correlation e↵ects. The simultaneous optimization of the orbital expansion coe cients and the configuration interaction coe cients to minimize the energy leads to a reference wave function for the subsequent treatment of the dynamic electron correlation by complete active space second-order perturbation theory (CASPT2) to obtain a consistent description of the electronic structure. The standard active space first proposed by Pierloot and Vancoillie [55] includes the TM-3d orbitals and two occupied ligand orbitals directed along the TM-ligand bonds plus a second shell of TM-d orbitals (the exact number of orbitals can vary depending on the coordination mode and number of d-electrons). Combined with a reasonably large basis set, very accurate estimates of the energy di↵erence can be obtained. A full review of the ins-and-outs of the CASSCF/CASPT2 approach in its application to spin state energetics can be found in Ref. [56] It is important to note that the lack of analytical CASPT2 gradients (not to mention the Hessian) makes it impossible to optimize the geometries or calculate the zero-point energy correction. Hence, the CASSCF/CASPT2 should always be combined with DFT calculations when one aims at an accurate estimate of E 0 HL . [57] In this aspect it is very interesting to mention the last developments in combining a multiconfigurational wave function directly with the speed of DFT in multiconfigurational pair-density functional theory (MC-pDFT). [58] The main advantage of this method is that it produces results with similar accuracy as the CASSCF/CASPT2 approach at the cost of a CASSCF calculation.
After the initial testing on small systems, there are now some recent studies in which the ability of the method has been tested to correctly describe the spin state energetics of transition metal complexes. [59] A final remark in this section concerns the role of the environment on the energy di↵erence. In most applications the calculations have been done in a gas phase setting, that is, an isolated complex with no environment. There are however some studies in which doubt is casted on the validity of this approximation. The above-mentioned GGA+U method is ideal for implementation in a periodic approach of the electronic structure, and application on systems with translational symmetry. The spin state energetics of isolated models compared to those of a complex in a lattice can be rather di↵erent, as shown by Vela et al. [54] and very recently also by Phung and co-workers in a study of the spin state energetics in bi-iron complexes. [60] The periodic GGA+U approach was also applied to explain the (non-)occurrence of SCO in a series of closely related systems [61] based on the e↵ect of the intermolecular interactions on E 0 HL . Obviously, such a rationalization is not possible when one only considers the SCO complex itself. Staying within the molecular approach, Radoń et al. found an important e↵ect on the HS-LS energy di↵erence by explicitly including the first shell of solvent molecules for the Fe III (H 2 O) 6 complex [62] and crystal packing e↵ects were put forward to explain the di↵erent thermal SCO behaviour of the two Fe sites in Fe II (methyl-tetrazole) 6 (BF 4 ) 2 in a cluster model study by Rudavskyi et al. [63]

Vertical excitation energies
From the very first beginning of CASPT2 in the early 1990s, [64] the method has been intensively applied to calculate vertical excitation energies in transition metal complexes. In general, accurate results can be obtained provided that one uses large be considered (see Figure 3). The generalization of CASPT2 to restricted active space second-order perturbation theory (RASPT2) [66] made that virtually any complex with one transition metal atom became within reach of multiconfigurational perturbation theory. [67,68] The rapid increase of the size of the active space severely hinders the extension to complexes with more than one metal center, but recent developments in quantum Monte Carlo [69] and density matrix renormalization group (DMRG) theory linked to CASPT2 have opened the door to treat polynuclear complexes with multiconfigurational approaches. [70,71] A computationally cheaper alternative is provided by time-dependent DFT (TD-DFT) [72] in the linear response version developed by Casida. [73,74] There exist by now a rather large amount of experience on how to simulate and accurately reproduce absorption spectra for organic molecules, [75,76] but the amount of data that has been published for transition metal complexes is less abundant and much more scattered.
Early work in the group of Baerends [77][78][79][80] showed that reasonable agreement with experimental data can be obtained using GGA functionals for a series of transition metal carbonyls and a collection of Zr and Ni tetrapyrrole (sandwich) complexes. The lower energy excitations were reproduced with good accuracy whereas the deviation for the higher excited states grows larger. In a follow-up study by Hummel et al. it was found that the excitation energies are rather strongly dependent on the applied functional and B3LYP was tagged as the most reliable one. [81]  bipyridine complexes. [82] On the contrary, PW91 gave a better account of the Ir(III) complex studied by Brahim and Daniel, [27,83] although it was argued that this may be caused by a cancellation of errors due to the lack of solvent e↵ects, whose inclusion may bring B3LYP back in closer agreement with the experimental data. Atkins and González constructed the absorption spectrum of [Ru(bpy) 3 ] 2+ by performing 1500 single point calculations on di↵erent conformations taken from a molecular dynamics simulation [84] and found that TD-DFT with the PBE functional gives a 0.3 eV underestimation of the first absorption band (higher excitations were not considered), which was argued to be in line with previous findings on the performance of the PBE functional for excited states in organic systems. [85] The systematic study of Latouche and co-workers [86] on the excitation energies in Pt II and Ir III complexes show that standard hybrid functionals behave better than the GGA and range-separated functionals, which tend to under-and overestimate the relative energies, respectively. The overestimation of the range-separated functionals can be remedied by a first-principles tuning of the ! parameter as extensively reviewed by Bokarev et al. [34] To end this small (and necessarily incomplete) overview of the performance of TD-DFT using di↵erent functionals, we mention the study of Pápai et al. addressing the excitation energies and potential energy surfaces of three Fe II complexes relevant for spin crossover. It was established that the B3LYP* functional gives relative energies for the metal-centered states that are in close agreement with CASPT2 results. [46] Furthermore, it should be noticed that most studies use the Tamm-Danco↵ approximation to full TD-DFT. This gives in general better results and avoids problems with spurious low-lying triplet states.
The third method for calculating vertical excitation energies combines the speed of DFT with the intrinsic accuracy of multiconfigurational methods, the DFT/MRCI approach of Grimme and Waletzke. [87] The DFT part of the calculation accounts for the dynamic (short-range) electron correlation, while multiconfigurational reference configuration interaction (MRCI) ensures a correct treatment of the non-dynamic electron correlation. The method counts with some (fixed) scaling parameters that largely eliminate the double counting of the correlation e↵ects inherent to mixed approaches. The method was shown to yield excitation energies in organic systems with an accuracy of ±0.2 eV, [88] and Escudero and Thiel addressed the accuracy of the DFT/MRCI excitation energies for TM complexes. [89,90] The comparison to TD-DFT and CASPT2 estimates shows that the method is also reliable in these inorganic systems. It systematically reproduces the same energetic ordering of the excited states as found in CASPT2 and the di↵erences in relative energies are significantly smaller than those found with TD-DFT, which were reported to be larger than 0.7 eV in some cases. Similar accuracy of the DFT/MRCI method was reported in the applications of Marian and co-workers to various Ir III -pyridyl complexes. [91,92]

Potential energy surfaces
Another aspect where computations can be of help for understanding the photochemical processes in transition metal complexes is in the determination of geometrical rearrangements that drive the changes of the electronic structure. While ground state geometries can be routinely determined with experimental techniques like x-ray di↵raction, the precise details of excited electronic states are more elusive. Concerning the calculation of ground state geometries, standard DFT can be used in most cases and virtually any density functional reproduces with good accuracy the experimental structures and vibrational frequencies. In addition, time-dependent DFT o↵ers a unique possibility to obtain detailed information on the geometries of excited states. In the first place, one has the possibility to find the optimal structure of the low-lying excited states of di↵erent spin multiplicities, and subsequently vibrational frequencies and IR/Raman intensities of these states can be addressed. Complementary information can be obtained by studying how the energies of ground and excited states evolve along certain predefined geometry distortions of the complex. Often a single coordinate defined as the interpolation between the optimal geometries of initial and final spin states gives a simple yet reasonable description of the process and allows one to generate an ab initio version [46,93,94] of the qualitative picture introduced by Hauser to explain the light-induced spin crossover. [95] Energies can be computed with one of the previously commented schemes, DFT-or wave function-based.
Moreover, potential energy surfaces can be mapped out along (a selection of) the vibrational normal modes [82,[96][97][98][99] to be used as input for explicit excited state dynamics simulations as further described in Section 4.3. Among the critical points on the PES of the di↵erent spin states, the minimal energy crossing points (MECPs) are of special relevance not only for luminescent properties [12,100] and photochemical reactions in transition metal complexes, [101] but also to spin-forbidden reactions. [102] In the latter case, MECPs determine the lowest energy where two electronic states of different spin are degenerate, and hence, take us to the lowest energy barrier in two-state reactivity. [103] In the former cases, the precise location of the MECP can be decisive whether a radiative or non-radiative deactivation mechanism is taken. In addition to the traditional algorithm introduced by Bearpark et al., [104] we also mention the automated global mapping procedure developed by Maeda and co-workers, [105] recently applied to the photochemistry of Re-complexes. [106] Vibrational frequencies and intensities of the IR/Raman transitions can also be extracted from the autocorrelation function of the time derivative of the dipole moment. [107,108] This requires long (ab initio) molecular dynamics runs, but has the advantage that the calculated vibrational spectrum explicitly includes the e↵ects of the environment, which is not possible in the time-independent approach of the vibrations calculated from the second derivative of the energy.

Core level spectroscopy
An important part of the experimental information on the excited state dynamics is based on x-ray spectroscopy in its many di↵erent variants: [109] x-ray absorption (near-edge) spectroscopy (XAS, XANES), x-ray photoelectron spectroscopy (XPS, ionization), [110] resonant inelastic x-ray scattering (RIXS), [111,112] Auger electron spectroscopy (AES), and x-ray emission spectroscopy (XES) [113] to name a few. The semi-empirical description of core level spectroscopy has a long standing tradition and goes back to the landmark papers by Thole and co-workers. [114,115] The method described therein has been applied in uncountable occasions, mostly for solid-state compounds. A satisfactory ab initio modelling of the di↵erent x-ray techniques has only emerged recently. [112,[116][117][118][119] The method is based on the previously mentioned RASSCF/RASPT2 approach in which the core orbital is also included in the active space and a large number of roots is calculated to address the core excitations. In the case of XAS, the intensities of the di↵erent transitions are accessible via the standard calculation of the dipole transition moments between initial and final states. The calculation of the intensities is somewhat more involved when initial and final states have a di↵erent number of electrons, for example in XPS. For these cases, Grell and co-workers outlined a strategy based on the Dyson orbital formalism, [120] which was shown to give very reliable predictions of the experimentally measured spectra, not only with respect to the peak positions, but also for the relative intensities of the di↵erent transitions. [118] 3 Spin-orbit coupling

Higher-order SO coupling
In addition to the direct spin-orbit coupling described above, there is also the possibility of higher-order coupling between electronic states through the spin-orbit operator. This higher order coupling is intermediated by other (excited) states and is conceptually nicely described with the following equation derived from second-order quasi-degenerate perturbation theory for the interaction between states 'a' and 'b' by spin-orbit coupling: where the first term on the right hand side is the direct coupling between the states and the second term represent the higher-order coupling intermediated by other states.
This second-order coupling has been invoked by Iuchi and Koga to investigate the interaction between quintet and singlet states in [Fe(bpy) 3 ] 2+ by SO coupling. [121] There is no direct interaction between these two spin states within the one-electron SO operator approximation. However, the perturbative treatment of the interaction leads to sizeable couplings through triplet states, which are likely to play an important role in the relaxation of the HS state in thermal spin crossover processes.
This perturbative approach can only be applied to address the higher-order coupling between the lowest states of each spin multiplicity. Its application to higher lying electronic states (as has to be done in the study of deactivation processes of excited states) leads to problems when the energy di↵erence between state 'b' and 'µ' becomes close to zero. Instead, the variational approach based on e↵ective Hamiltonian theory [122] does not show this limitation and can be applied to derive the e↵ective coupling (direct plus all higher-order couplings) between all electronic states that play a role in the phenomenon under study. The method has recently been applied to calculate the e↵ective coupling between the quintet and the triplet MLCT states in [Fe(bpy) 3 ] 2+ ; [123,124] two high-lying excited states that play a fundamental role in the light-induced magnetism of this compound.

Introducing time and temperature
So far, we have only considered those properties related to the electronic structure of transition metal complexes in a static framework, that is, time and temperature are not defined. Such a static description of the deactivation is inherently limited as nuclear motion may influence both the photoexcitation from the ground state and may also occur during the deactivation process. A relatively simple, yet interesting strategy to account for the nuclear motion caused by thermal disorder in the structure of transition metal complexes consists of a large number of static electronic calculations performed on a series of geometries (snapshots) extracted along the trajectory of a molecular dynamics simulation. [84,125,126] In this way, one can get a first impression of the influence of vibrational motion on the electronic properties of the ground state of the complex and its photoexcitation. Nevertheless, the description of the decay process from the photoexcited state requires more sophisticated methods that will be briefly covered in the forthcoming sections. The method has been applied to Fe II,III -polypyridyl complexes [47,123,124,128] , a Cu I complex [129] and several OLED-related Ir III complexes, [92] among many other applications including organic molecules. A similar approach, based on second-order perturbation theory, was introduced by Peng et al. [130] and has recently been applied to study the non-radiative decay processes in Ir III OLEDs. [131] The results of the aforementioned studies indicate that the golden rule, despite its approximate nature, leads to intersystem crossing rates that are at least of the correct order of magnitude in those cases where experimental data are available. Hence, it o↵ers an interesting alternative to the approaches discussed in the next sections, which are more expensive but introduce real time resolution, and therefore, a much more detailed description of the excited state dynamics.

Non-adiabatic on-the-fly dynamics
Ultrafast non-radiative decay processes going through conical intersections might involve large geometrical distortions pushing the system outside of the Frank-Condon portion of the phase space to achieve statistical significance, but they must also be capable of describing electronic states of rather di↵erent character. In the following, we overview recent developments using molecular dynamics techniques to describe ultrafast phenomena on TM complexes. In-depth reviews are also available covering the use of molecular dynamics and DFT for ultrafast ISC processes on TM complexes [27] and the reactivity that might derive from such ultrafast processes. [37] The review article of Penfold et al. gives an excellent in-depth description of the computational schemes that are currently being used and further developed for the study of excited state dynamics. [38] Ab initio molecular dynamics (AIMD) propagates the nuclei classically, while the electronic potential is calculated on-the-fly at each time step by quantum mechanics.
Trajectory surface hopping [132,133] (TSH) is a variant of AIMD, where the nuclei are classically propagated over an adiabatic PES generated on-the-fly for a particular electronic state, which incorporates non-adiabatic processes by introducing a stochastic method to trigger the transition between di↵erent electronic states. [134] TSH simulations can incorporate all degrees of freedom of the system in a e cient manner that allows treating medium to relatively large systems with hundreds of atoms. The semi-classical propagation method developed by Nakamura and co-workers [135] can be used to introduce quantum e↵ects like tunnelling in the time evolution of the nuclear movement.
Nonetheless, these methods rely on the use of (semi-)classical nuclei and hence exclude quantum e↵ects related to nuclear vibration. The limitations on the electronic part will depend on the method of choice for the electronic structure calculations and usually require some previous assessment. [136] The original hopping algorithm uses various criteria such as populations and the SO coupling to push the simulation through a conical intersection on-the-fly when the appropriate conditions are met, but over the years several other techniques have been developed. [137] Capano et al. [138] describe with TSH the photodynamics of [Cu(dmp) 2 ] + (dmp=2,9-dimethyl-1,10-phenantroline), where photoexcitation leads to a bright singlet excited state that quickly decays to a singlet MLCT with the electron localized on a single ligand, a phenomenon that requires a time-dependent method to be captured in the simulation. In this case, the TSH dynamics used a hopping algorithm allowing for non-adiabatic transitions between singlet states on-the-fly, whereas ISC rates were calculated a posteriori on the resulting TSH trajectories. The results show a fast deactivation process of 100 fs that, surprisingly, not only passes through several IC, but also several ISC involving a manifold of triplet excited states, even though the initial and metastable states have both singlet spin multiplicity. This convoluted decay path highlights the critical role that ISC can have even in those cases where it is not anticipated. A consistent treatment of the IC and ISC along the whole simulation can be achieved with the SHARC method. [139] This method calculates the hopping probabilities through an equation of motion based on a diagonalized electronic Hamiltonian including SO coupling and hence, in combination with appropriate quantum chemistry methods, it o↵ers TSH dynamics with non-adiabatic couplings and SO coupling computed on-the-fly along the full trajectories. Atkins and González [84] studied the ultrafast deactivation of the [Ru II (bpy) 3 ] 2+ complex with the SHARC method, following the decay from a bright potentials. [140][141][142] AIMS with CASSCF has been applied to a dimethylnitramine Fe complex to explore the di↵erent decomposition paths after electronic excitation. [143] Nonetheless, application to TM complexes subject to spin crossover requires including the spin-orbit coupling associated with conical intersections. In this regard, the generalized AIMS method developed by Curchod et al. is promising as it is capable of describing IC and ISC events. [144]

Quantum dynamics with MCTDH
Spin crossover systems are characterized by strong vibronic coupling between the states involved in the deactivation process, especially in proximity of conical intersections.
Even at low temperatures, a spin crossover model system quickly breaks the Born-Oppenheimer (BO) approximation due to the strong vibronic motion triggered by the photoexcitation. [145] The following decay may reach or avoid the minimum of the HS state depending on the strength of the tunnelling constant arising from spin-orbit coupling. Furthermore, time dependent calculations on the ultrafast spin crossover (20 fs) of a Fe-Co Prussian blue showed that spin-orbit coupling can rapidly vary during the deactivation process. [146] Therefore, propagating the nuclei in classical trajectories does not su ce to treat systems with strong vibronic coupling and/or strong spinorbit coupling and it becomes necessary to break the separation between nuclear and electronic motion imposed by the BO approximation by means of quantum dynamics.
Recent experimental and theoretical developments on the spin-vibronic mechanism related to ISC are reviewed in detail by Penfold et al. [38] Quantum dynamics can be executed with the Multi-Configuration Time-Dependent Hartree [147][148][149]  Another example is the description by Huix-Rotllant et al. [152] of the light induced mechanism of carbon monoxide release on the heme complex, which not only includes the photolysis step but also a spin crossover process. The Multi-Layer MCTDH simulation included 179 electronic states and 10 vibrational modes reaching 1 ps of length, a short time span for common biologic processes, but long enough to describe ultrafast processes. The heme-CO photolysis is dominated by symmetry breaking JT distortions, occurring within the first 15-20 fs after photo excitation to the 1 Q band and transfer to a singlet 1 MLCT. At this point the SCO triggers, driving the system towards the HS quintet state. First, a singlet-triplet ISC populates the 3 MLCT at ⇠70 fs and later, at 300 fs, a triplet-quintet ISC populates the 5 MLCT state.

Summary and Outlook
The theoretical description of the electronic properties of transition metal complexes related with the deactivation of excited states is a rapidly expanding field. Especially the recent advances in the description of the excited state dynamics are providing the scientific community with new tools for further understanding the details of the complex mechanisms that rule the photophysics of these complexes. All the properties without explicit time dependence, such as the relative energies of the di↵erent electronic states, optimal geometries, critical points on potential energy surfaces, spin-orbit coupling, etc., are readily computed with any of the well-established quantum chemical strategies as discussed in the first part of this overview. In most cases, there is enough experience in the literature to decide on the most appropriate computational scheme, either based on a multiconfigurational wave function or alternatively with density functional theory. Among the few points that still have not received too much attention and would certainly be worthwhile to explore are the excited state absorption (ESA) and the possibility to optimize geometries within the MC-pDFT scheme. The theoretical study of ESA can be very useful to interpret the results of the widely used pump-probe experiments in time-resolved spectroscopy. Recently, a promising approach has been published based on the RASSCF approach considering a large number of roots. [153] Good results are reported for the organic chromophore benzophenone and it would be very interesting to investigate the performance of this approach for transition metal complexes. The second new aspect could o↵er a multiconfigurational geometry optimizer that is competitive with the standard DFT schemes, both in accuracy and speed. This can possibly be advantageous in situations where the assumption of a single Kohn-Sham determinant becomes arguable as in systems with strong biradicalar character. MC-pDFT analytical gradients can be used for single state calculations [154] and when the state averaged variant is available, the method can in principle be tested for excited state geometry optimizations.
The calculation of intersystem crossing rates through the application of Fermi's golden rule or the perturbative expression derived by Peng and co-workers [130] has shown its usefulness in several applications. However, it would be desirable to be able to treat the internal conversions in a similar fashion so that (rough) lifetime estimates can be obtained for all stages in the photocycle without having to rely on an explicit time propagation of the system. An interesting study in this aspect is the one by Valiev et al., [155] in which a procedure is presented to calculate ISC and IC rates based on the same starting equation for non-radiative decay.
Given the success of the ab initio calculations of the di↵erent x-ray spectra, a straightforward theoretical simulation of the time resolved core-level spectroscopic measurements can possibly be obtained by performing the calculations on a regular time interval along the trajectory of a molecular dynamics simulation of the system.
A similar procedure could be used for excited state absorptions.
The direct simulation of the dynamics of excited states in transition metal complexes just set o↵ seriously. TSH methods and MCTDH dynamics are very powerful techniques, capable of accurately describing all electronic and nuclear changes involved in an ultrafast deactivation. Since these calculations are very CPU intensive at the moment (important improvements can still be expected from more e cient algorithms and developments of the computer hardware), it is not possible to just go for the method with all the bells and whistles though. Therefore, their application is not straightforward, the methods have to be tailored to the system at hand and hence, some previous knowledge of the system is necessary. Are large molecular distortions to be expected during the deactivation? Does the decay path go through ISC? How many states are involved? Unfortunately, that information might not be available and the nature of the system might hide the underlying complexity of its deactivation process.
Nonetheless, the implementations of these methods is constantly improved as more experience is gained with these simulations.