The Quest to Simulate Excited-State Dynamics of Transition Metal Complexes

This Perspective describes current computational efforts in the field of simulating photodynamics of transition metal complexes. We present the typical workflows and feature the strengths and limitations of the different contemporary approaches. From electronic structure methods suitable to describe transition metal complexes to approaches able to simulate their nuclear dynamics under the effect of light, we give particular attention to build a bridge between theory and experiment by critically discussing the different models commonly adopted in the interpretation of spectroscopic experiments and the simulation of particular observables. Thereby, we review all the studies of excited-state dynamics on transition metal complexes, both in gas phase and in solution from reduced to full dimensionality.


INTRODUCTION
Transition metal complexes provide a rich photochemistry that can be utilized in applications from solar energy conversion to medicine. 1 This is possible due to the large variety of electronic states of distinct nature that transition metal complexes have to offer. For instance, while long-lived states of metal-to-ligand charge transfer character are key to applications in dye-sensitized solar cells, 2−4 short-lived metal-centered states can mediate dissociation processes in biology. 5,6 In general, the behavior of transition metal complexes after light irradiation is controlled by the presence or absence of radiationless reaction pathways. These can either enable efficient transfer between electronic states or facilitate long-lived excited states that at last may emit. Unveiling nonradiative reaction pathways is therefore key to understanding and ultimately tuning the photochemistry of transition metal compounds.
In general, radiationless transitions between electronic states in molecules are categorized into two types: internal conversionfor transitions between electronic states of the same spin multiplicityand intersystem crossing, which connects electronic states of different spin multiplicity. The nonradiative behavior of photoactivated molecules is driven by the motion of the nuclei. Upon excitation, the molecules in their equilibrium geometry are lifted from the electronic ground-state potential into an excited-state one. This brings the molecule to a non-equilibrium situation that induces nuclear motion toward other conformational regions. On its following journey, the molecule then can pass regions with high probability of nonradiative transitions to other electronic states that ultimately decay back to the electronic ground state. Alternatively, the molecule can end up in regions with no possibility of (further) nonradiative transfers, from where it can only return to the ground state via luminescence.
Based on the nature of the processes, the study of radiationless pathways in experiment has two prerequisites. One is identifying a property that changes during the reaction in order to detect which species are present during the reaction. Since radiationless processes involve a change of the electronic state, these properties may directly or indirectly relate to the electronic potentialsthis is the case of following changes in absorption intensities or shifts in vibrational frequencies. The second is to be able to monitor this property on the same time scale as the reaction occurs. This time scale is dictated by the nuclear motion of the molecule which takes place typically in the femtosecond regime. 7 Thus, it is not surprising that electronic spectroscopy techniques, such as transient-absorption spectroscopy or timeresolved emission spectroscopy as well as vibrational spectroscopic techniques, such as time-resolved infrared spectroscopy with sub-picosecond resolution, have made incredible progress to access the nature of ultrafast processes during nonradiative reactions. 8 While substantial information can be obtained from spectroscopic experiments, often this might not be sufficient to derive a detailed description of the excited-state dynamics. Furthermore, with increasing system size, a larger number of vibrational degrees of freedom and greater density of electronic states may participate in the photodynamics, making it more and more difficult to interpret experimental signals. Most transition metal complexes fall into this category: even small coordination compounds already possess dozens of atoms which add up to hundreds of vibrational degrees of freedom and possess many close-lying electronic states due to the only partially filled d shell. For this reason, experimental studies are almost routinely accompanied by theoretical calculations. 9−13 Most of these efforts involve quantum-chemical calculations of electronic states and relevant points in their potential energy surfaces (PES). 14,15 For example, with the help of calculated absorption spectra, one can infer the electronic states that may play a role in the photodynamics. Naturally, it would be better to go beyond this static description and carry out explicit dynamics simulations that directly monitor the time evolution of the electrons and nuclei in the molecule occurring after light irradiation.
One could go as far as to assert that theory not only helps the interpretation of experiments but is a predictive tool on its own. 16,17 In practice, however, photodynamics simulations of transition metal complexes are still very much limited due to their computational cost. This limitation restricts the size of the molecules that can be studied and compromises its accuracy because it necessitates diverse approximations. In this predicament, one should ponder the benefits of experiment versus theory as follows. The experiment can be considered an exact instrument (within experimental resolution) that probes the full system of the molecule interacting with its real environment. In contrast, theory is an approximate instrument with an inherent methodological error and is often restricted to a truncated system, i.e., reduced molecular model system with or without an environment. The experiment, however, gives only limited information. By monitoring a signal, one observes only the global response of the molecule to light resulting from all (simultaneously occurring) processes. In contrast, theory is able to yield every single deactivation pathway.
Clearly, these contrasting advantages and shortcomings ask for synergy between experiment and theory. To efficiently collaborate, practitioners of each side should understand the mindsets and workflows of the other side to adequately assess the information that can be derived from experiment and theory, respectively. This is particularly important in the developing research field of transition metal photodynamics, where each single discipline is left in the dark when on its own.
The motivation of this Perspective is illuminating what can be done from the theoretical front, explain which are the current strategies of simulating the photodynamics of transition metal complexes and its limitations, help in the interpretation of theoretical simulations, and pose the challenges still present in the field, thereby reviewing the studies done until now. If the gap between theory and experiment can be made smaller, more efficient synergies between the two "approaches" can better meet the current and future challenges in transition metal photochemistry. To this ambitious goal, we present the basic theory and typical ingredients that enter the simulation of excited-state dynamics simulations, with their strengths and weaknesses, in relation to transition metal complexes.

Time-Dependent Schrodinger Equation
In order to understand the limitations of current theories in the simulation of the photodynamics of transition metal complexes, we deem necessary to introduce some underlying working equations. The key equation to solve is the time-dependent Schrodinger equation 18 tot (1) that describes the time evolution of a molecule. In general, the time-dependent molecular wave function Ψ(r, R, t) depends on the time t, the positions r of the electrons, and the positions R of the nuclei. The total Hamiltonian̂ tot can be written as ( , ) ( ) tot mol (2) wherê r R ( , ) mol is the molecular Hamiltonian that contains the kinetic and potential energy terms of the electrons and nuclei in the molecule, as well as a term̂ r R ( , ) SOC that couples electronic states of different multiplicity via relativistic spin− orbit coupling (SOC) and will be discussed later. The interaction between the molecule and the light is given in the semiclassical approximation via the dipole operator μ̂r R ( , ) and a time-dependent electric field ε(t). Given the complexity of eq 1, usually one separates the motion of the nuclei from that of the electrons. To do this, we can express the molecular wave function as a product of electronic wave functions Φ el (r; R) and nuclear wave functions Ψ nuc (R, t): The electronic wave functions Φ el (r; R) are eigenfunctions of the electronic Hamiltonian̂ el that can be obtained through the time-independent electronic Schrodinger equation: 19,20 The time-independent electronic Schrodinger equation can be solved for fixed selections of nuclear coordinates R. Thus, the electronic wave functions (and so its properties) depend only parametrically on the nuclear coordinates R, and their eigenvalues E el (R) yield the potential energy for this specific geometry. As a result, the motion of the nuclear wave functions Ψ i nuc in the electronic potential E i el is described bŷ wherê nuc is the nuclear kinetic energy operator and̂ ii DBOC is the diagonal Born−Oppenheimer correction or diagonal kinetic coupling. 21 The electronic potentials are coupled through offdiagonal terms that include The term̂ ij NAC are the so-called nonadiabatic couplings (NACs), 22 that furthermore couple the motion of the nuclei and electrons. μ̂i j is the transition-dipole moment between electronic states Φ i el and Φ j el , and ∇ R A denotes the gradient with respect to the nuclear coordinates R A .

Nonadiabatic Couplings
In electronic structure theory, the molecular Schrodinger equation is solved (yet approximately) only by neglecting the NACs of eq 7. In this Born−Oppenheimer approximationor adiabatic approximation if the diagonal Born-Opppenheimer correction̂ ii DBOC is includedthe motion of the nuclei and electrons is completely decoupled and the nuclear wave functions are then constrained to a single electronic potential. Naturally, this approximation is only useful when describing processes that take place in a single electronic state. This can be the case, when the electronic state in question is well-separated in energy from all other electronic states, e.g., as is usually found for reactions that occur in the electronic ground state. However, if we are interested in reactions including several electronic states, it is mandatory to include NACs to allow the nuclear wave function to transfer between the different electronic states (internal conversion). In eq 7, the second term in the sum can be written as This shows that the NACs become large when the corresponding electronic states are close in energy (E j − E i ≈ 0). Since it is the NACs that enable the transfer between different electronic states, state transfer is most efficient when the corresponding states are close in energy.

Spin−Orbit Couplings
The SOC couples electronic states of different spin multiplicities and, thus, allows both radiative transitions (phosphorescence 23 ) and nonradiative transitions (intersystem crossing 12,24 ) between them. SOC is a relativistic effect, as it occurs naturally in a formulation of quantum mechanics that includes the principles of the theory of special relativity. Phenomenologically, SOC is explained as the interaction of the magnetic moment of the spin angular momentum with the magnetic field that is induced by the electron orbiting around the nuclei as well as in the field of the other electrons. 25,26 In nonrelativistic quantum theory based on the Schrodinger equation, SOC, has to be introduced ad hoc.
Using the Breit-Pauli operator,̂ SOC can be expressed as The first sum in eq 9 describes the interaction of each electron's orbital angular momentum̂×r p i i A (orbiting around the nucleus A) with its spinŝ i . This interaction depends on the charge of the nucleus Z A . For valence electrons in many-electron systems, SOCs typically scale as Z A 2 , 27 which manifests strongly in heavy atoms. In transition metal complexes, SOCs are often strong enough to allow ultrafast intersystem crossing, 28 and it is thus necessary to include them during nonadiabatic dynamics simulations. The strength of the SOC not only is dependent on the presence of heavy atoms in the molecule but also depends on the character of the electronic states that are coupled. This is rationalized by the generalized El-Sayed rules 24 that state that in order to provide large SOCs for efficient intersystem crossing, (i) the coupled electronic states should differ only by a single excitation, which (ii) involves a change in orbital type, and (iii) the orbitals should be localized at the same site in the molecule.

GENERAL STRATEGY TO SIMULATE
EXCITED-STATE DYNAMICS From the previous section, it is clear that in order to follow the time-evolution of a molecule, we first need to collect several electronic ingredients: the electronic energies E R ( ) i el (eq 4) and implicitly their gradients, the NACs (eq 7), the SOCs (eq 9) and eventually also the transition-dipole moment μ ij (R) if the light− matter interaction is explicitly simulated (eq 6). In the second step, the motion of the molecule can be simulated by Step 1 refers to the ingredients that must be calculated from electronic structure theory.
propagating the nuclear motion in the different (coupled) electronic states. This strategy, which guides the structure of the remainder of this Perspective, is collected in Figure 1.
Both the electronic and nuclear steps face practical challenges. On the one side, sufficient accuracy in the electronic structure part is needed. On the other side, the scaling of the computational cost inherent to the size of the molecule needs to be controlled. In the exact limit, this cost scales exponentially with the molecular size for solving both the electronic and the nuclear problem. To alleviate the steep scaling in the electronic part, a large number of electronic structure methods have been developed whose scaling follows different power laws. Thus, when choosing one electronic structure method for dynamics, the selection is strongly motivated by maximizing accuracy for a given computational effort. This situation is different in the nuclear part, where selecting a method is more a fundamental choice of which processes should be qualitatively well-described, with less focus on quantitative accuracy. A closer insight of the available possibilities to solve the electronic and nuclear problems is given in section 4 and sections 5 and 6, respectively.

Practical Considerations
There exist a large number of both commercial and (academically) freely distributed quantum chemistry computer program packages that offer a wide range of electronic structure methods to calculate the electronic ingredients. In practice, however, the selection of a particular method to underlay nonadiabatic nuclear dynamics of transition metal complexes is limited to those implementations that are able to provide the energies, SOCs and possibly gradients and NACs, depending on the approach chosen for dynamics. When gradients are needed, it is desirable to have implementations that allow for analytical rather than numerical gradients. 29 Analytical implementations obtain the derivatives along all nuclear coordinates in a single calculation instead of requiring the "manual" displacement along each nuclear coordinate which makes numerical implementations computationally more expensive; however, analytical gradients are not available for many methods. When NACs are not available, it is possible to approximate the state-to-state transition probabilities by computing the overlap between the electronic wave functions. 30,31 The SOC matrix elements can also be calculated a posteriori for some quantum chemical methods. 32 The choice of a quantum chemical method is further limited if solvent is included implicitly in the dynamics simulation, as all required properties need to be available including the solvent, or if schemes containing the solvent explicitly are not implemented with the electronic structure method of choice. Fortunately, steady effort in the development of quantum chemistry program packages is continuously adding new implementations that can be employed for nonadiabatic dynamics simulations.
Depending on the number of nuclear degrees of freedom considered, the number of individual electronic structure calculations that need to be performed for the excited-state dynamics simulations can easily reach ∼10 5 . Compared to stationary explorations of PEScomprising the search for potential energy minima, crossing points, and minimum-energy or linear-interpolated paths connecting themwith typically only few dozens of calculations, this is an enormous increase in computational cost. Thus, it is often necessary to balance the trade-off between computational cost and accuracy differently between static and dynamic studies of excited-state processes. We caution, however, that this does not imply that anything goes! A dynamics study with wrong electronic ingredients will only produce a collection of meaningless numbers, regardless of how much computational effort has been invested in the dynamics simulation. Instead, the selection of the electronic structure method has to be guided by carefully testing their performance against experimental reference data or other reliable electronic structure methods. 33 Therefore, it gets clear that besides the question of availability and efficiency, the choice of an electronic structure method should be tailored to describe each transition metal complex and the reaction of interest. Since we are interested in electronically excited states, the first experimental reference of choice is often the static absorption spectrum. 34,35 A calculated spectrum with the chosen method should qualitatively match the experimental one in terms of number of observed absorption bands, their relative intensity, and preferably also band shoulders. For the relative position of the calculated and measured absorption bands, one should strive for energetic differences of 0.1−0.5 eV as this is a reasonable accuracy that can be achieved for excitation energies of medium-sized and larger molecules. Smaller average error may be obtained only with very accurate methods; however, such methods are most likely computationally unfeasible for dynamics except for few-atomic molecules. 36 Sometimes, larger errors in the total excitation energies can be acceptable in excited-state dynamics studies when the error in the relative excitation energies is small and the same along the PESs of interest. For example, dynamics excluding relaxation to the ground state might be well-described despite larger errors in the excitation energies if the error is systematic for all excited states, i.e., when all excitation energies are under-or overestimated by a similar amount and can be accordingly scaled. 37,38 An additional cross-check can be achieved by computing time-resolved spectra that can be directly compared with experimental counterparts, coming, e.g., from time-resolved absorption spectroscopy, 39,40 time-resolved emission spectroscopy, 41,42 time-resolved photoelectron spectroscopy, 43−46 or time-resolved X-ray scattering. 47,48 Such calculations can a posteriori validate the selection of the electronic structure methods.
Particularly, if no experimental references are available, it is useful to compare the calculated energies and geometries of selected points in the PESsuch as ground-and excited-state minima, conical intersections, or interstate crossingsagainst other (higher-level) electronic structure methods. Available benchmark studies can also help in this endeavor. However, while the number of benchmark studies that systematically investigate excitation energies 49−56 or oscillator strengths and excited-state dipole moments 57 of organic molecules continuously grows, corresponding studies for transition metal complexes are much more scarce and limited in size. 58−61 The absence of such data is not surprising as benchmarking transition metal complexes requires considerable effort due to the larger size and the greater variety of characters of electronic states in transition metal complexes compared to organic molecules. For some families of wave-function-based electronic structure methods, there exists a hierarchy in terms of accuracy, 33 which can be consulted in the absence of other reference data to decide on the reliability of the computed PES.
Finally, we note that even if the ground state might not play an apparent role in the excited-state dynamics simulations, the method of choice should be capable to describe it. This is because electronic excited states are described either by a configuration interaction or (coupled-) cluster expansion or from perturbation theory/response theory applied to the electronic ground state. Thus, at least a qualitatively correct ground-state wave function is mandatory to obtain a correct description of the excited states. This requirement should be extended to all geometries that may be visited during the dynamics. Additionally, the quality of the electronic ground state also influences the vertical excitation energies, which in turn can affect the excitation wavelengths and even the course of the dynamics.
In general, electronic structure methods can be grouped into two kinds of approaches, multiconfigurational (sometimes also grouped together with multireference) methods and singlereference methods. 62,63 These approaches will be discussed in the following sections 4.2 and 4.3 with a focus on transition metal complexes. Semiempirical methods, which simplify the solution of eq 4 by replacing expensive integrals with experimentally parametrized corrections, 64 deserve an extra category (section 4.4), as they can be formulated in singlereference or multi-reference variants.

Multiconfigurational and Multireference Electronic
Structure Methods 4.2.1. Applicability. In contrast to single-reference (SR) methods, 65 which describe the ground state by a single configuration expressed by a Slater determinant, multiconfigurational (MC) methods 66 consider more than one configuration/ Slater determinant and multireference (MR) methods 67,68 use such linear combination of multiple configurations/Slater determinants for the ground state to generate further excitations. By construction, MR/MC approaches offer a more flexible description than the SR ansatz, and they are also capable of describing open-shell electronic ground states, which can frequently be encountered in transition metal complexes due to the partially filled d-shell of the metal atom, and they are also capable to describe ligand dissociation. For example, among octahedral complexes, most d n configurations (see Figure 2a) feature spatially degenerate electronic ground-state terms that can only be described correctly by using multiple Slater determinants.
4.2.2. Complete Active Space Self-Consistent Field and Beyond. The most popular MC approach for the study of transition metal complexes is the family of complete-activespace self-consistent field (CASSCF)-based methods. 71−73 The CASSCF method itself 74 is capable of achieving a qualitative correct description where a SR method fails; however, CASSCF excitation energies can often bear sizable errors. 75 Therefore, it is often used as the starting point of a second-order perturbation theory treatment, e.g., in the CASPT2 76 or NEVPT2 77 approaches, or in a subsequent configuration interaction (CI) expansion, as in the multireference configuration interaction (MRCI) 78 ansatz. CASSCF-based methods rely on the selection of an active space, which comprises a subset of occupied and unoccupied orbitals of the molecule. For this subset, the best possible wave function is then calculated in a full CI treatment with simultaneous reoptimization of the orbitals. Subsequent perturbation or CI expansions of a CASSCF wave function add energy corrections to the CASSCF states.
CASSCF-based calculations scale factorially with the size of the active space. This scaling imposes a limit to the maximum size of the active space that can be treated nowadays, perhaps of ∼20 orbitals including 20 electrons, as in the Cr 3 example of ref 79 using 4 × 10 9 Slater determinants in a single electronic structure calculation. Therefore, it is important to carefully restrict the active space to the most important orbitals expected to contribute to the problem of interest 71,80 such as the orbitals that characterize the ground and excited states in the nonadiabatic dynamics simulations. For transition metal complexes, this can include the metal d orbitals, metal− ligand−σ and σ* orbitals, as well as π, π*, and nonbonding n ligand orbitals that typically characterize the low-lying excited states. In addition, for complexes featuring 3d metals, it is advised to include a second set of 3d′ orbitals (a double d shell) when studying processes where the occupation of the d shell changes, such as in metal-to-ligand charge-transfer excited states. 71 One realizes that for transition metal complexes, attempting to include all formally important orbitals in the active space quickly exceeds the computational limits of the standard CASSCF ansatz. Extended approaches to treat larger active spaces have been developed as the restricted (RASSCF) 81 and generalized (GASSCF) 82 variants of CASSCF or the density matrix renormalization group (DMRG) 83,84 approach. The RASSCF/GASSCF methods allow the use of larger active spaces of ∼30 orbitals 82 by defining subspaces in the active space with limited interaction between the orbitals in the different subspaces. The active space limit can be pushed further up to the order of 100 orbitals in DMRG, 85 which uses tensor decomposition methods to approximate the CASSCF wave function. While for typical DMRG calculations, active space sizes of ∼50 orbitals are manageable, 86 when applying a subsequent second-order perturbational treatment to the DMRG wave function as in CASPT2 or NEVPT2 approaches to obtain more accurate electronic energies, the active space size that can be handled computationally decreases again to ∼30 orbitals 87 for a single-point energy calculation. Figure 2c shows an example of a so-called entanglement diagram that can be obtained from a DMRG calculation to estimate the importance of different orbitals in an active space. Other methods to include more than one configuration can be found in ref 63 96,104,105 However, in these cases, the PESs were restricted to one or two dissociative coordinates. Similarly, wave packet dynamics studies of small rhenium-(VII) 108 and chromium(III) 109 complexes employed CASSCF by itself and CASSCF combined with quasi-degenerate perturbation theory, respectively, on one-dimensional PES. In comparison, the number of studies using MR or MC methods in ab initio molecular dynamics simulations (see section 6) is much smaller. We could find only one study: an ab initio multiple spawning (see section 6.2) study of a small iron(II) complex, which included the ground state S 0 and the first-excited singlet state S 1 computed with CASSCF and an active space including 11 orbitals. 110 This scarcity is in strong contrast to ab initio molecular dynamics simulations of organic molecules, where MC/MR methods such as CASPT2 or MRCI are used very often 111−115 with active spaces including up to 11 orbitals. 112 As it will be introduced later (section 5.3), it is also possible to use parametrized PES where to run nonadiabatic simulations at a much lower cost. This strategy has been recently employed in a wave packet dynamics study of a heme−CO complex 116 carried out on 15-dimensional CASSCF/CASPT2 PES, as well as in an ab initio molecular dynamics study of a vanadium(III) complex on its full 123-dimensional CASSCF PES. 69 These studies used active spaces including 9 and 13 orbitals, respectively.

Applicability.
When all molecular geometries visited during the dynamics possess an electronic ground state that can be described by a single configuration, then SR electronic structure methods can be employed to calculate the electronic potentials, gradients, and couplings. This condition might be fulfilled if the molecule possesses a closed-shell electronic ground state, and neither dissociates nor undergoes internal conversion from the first excited singlet state S 1 to the ground state S 0 during the dynamics simulation. Note that the opposite is not true; i.e., a SR method does not necessarily describe a closed-shell ground-state molecule well.
Compared to MC/MR, SR methods feature a considerably lower scaling in the computational cost with regard to the system size, 36 making them more amenable to transition metal complexes. Applications revolve around (closed-shell) lowspin d 6 complexes 42,48,117−136 and d 10 systems 137−141 with one single exception. 142 As a further advantage, SR methods can be considered more "black box" than MC/MR methods as they typically depend on fewer (and less critical) parameters, which makes them more user-friendly.

Time-Dependent Density Functional
Theory. The most popular and computationally lowest-scaling SR method is time-dependent density functional theory (TDDFT). 143 It is based on the standard (time-independent) DFT, in which the ground-state energy of the molecule is expressed as a functional of the electron density. The energy can be calculated in the Kohn−Sham formalism 144 that assumes the electron density of the real system to be identical to that of a fictitious system of noninteracting electrons. The electrons in the noninteracting system can conveniently be described by a set of orbitals to calculate the electron density. This gives the ground-state wave function the form of a single configuration. The differences between the system of interacting and noninteracting electrons are combined in a so-called exchange-correlation (XC) functional that also depends on the density. The drawback of the Kohn−Sham formalism, however, is that the exact form of the XC functional is unknown, which has led to numerous efforts to find approximated forms of XC functionals. 145 For excited-state calculations, TDDFT can be cast in a linear-response (LR) formulation 146 that allows the calculation of excitation energies from matrix eigenvalue problems without the need to propagate the time-dependent density explicitly, and it yields state vectors that are analogous to the CI expansions from wave function methods. TDDFT is thereby typically employed in the adiabatic approximation 144 that uses the XC functionals from ground-state DFT.
In practical applications, TDDFT can suffer from a number of problems, 147,148 the most notorious being the dependence of the calculated energy and properties on the choice of XC functional. For excited states of transition metal complexes, hybrid XC functionals, such as B3LYP or PBE0, seems to be preferable compared to generalized-gradient functionals, as discussed previously. 149,150 This preference is also mirrored in excitedstate dynamics studies, that mostly rely on hybrid functionals 42,48,117−130,136,138−141 with some exceptions. [131][132][133][134]139 In addition to the general functional dependence, TDDFT is well-known to suffer from its inability to describe charge transfer (CT) excitations with standard XC functionals. 151 This problem is alleviated using long-range corrected XC functionals. 152,153 For transition metal complexes, the CT problem is more prominent for interligand CT excitations, while CT transitions between metal and ligands are usually described well. 9 When a large number of excited states is calculated, TDDFT can fail to describe high-energy states with energies between the ionization potential and the negative energy of the highest-occupied molecular orbital correctly. This problem, however, can be corrected using special asymptotically corrected XC functionals. 154 Finally, TDDFT describes excited states only using single excitations. This not only disregards double and higherorder excitations (see an example in Figure 2b) but also can fail to describe single excitations when they involve a spin-flip in open-shell molecules. 155 Some states of double excitation character are, in turn, accessible using a spin-flip TDDFT ansatz. 156 As is summarized in ref 147, TDDFT works best for "low-energy one-electron excitations involving little or no charge transfer and that are not too delocalized". In all other cases, care should be exercised.

Coupled Cluster and Related
Methods. In addition to TDDFT, there exist other SR wave function methods, most notable coupled cluster (CC) 157 as well as the algebraic diagrammatic construction (ADC) scheme of the polarization propagator. 158 Both methods have a clear hierarchy that allows to improve the accuracy of the results systematically with a concomitant increase of computational cost. The most economic variants, approximate second-order CC (CC2) 159 and second-order ADC [ADC (2)], 160 scale slightly larger than TDDFT with the system size (N 5 vs N 4 ) and can deliver similar accurate excitation energies 161,162 as well as excited-state geometries 163 at least for organic molecules, provided the excited states are dominated by single excitations. 158 In contrast to TDDFT, both CC2 and ADC(2) variants include double excitations, however, only at a zeroth-order level. This contributes some admixture of single and double excited states, but it is insufficient to describe states with large doubleexcitation character adequately. Probably for this reason, CC2 or ADC(2) studies on transition metal complexes are rare 164 and in general little reliable.

Semiempirical Methods
Semiempirical methods 64 usually neglect or parametrize molecular integrals occurring in the calculation of the electronic Schrodinger equation (eq 4). This way they can deal with large molecules that might be prohibitive with ab initio wave function theory or even DFT. Semiempirical methods can be based either on molecular orbital (MO) theory 165 or on DFT. 166 Further, they can be combined with CI schemes 167 or MRCI ones 168 to be used in ab initio molecular dynamics simulations.
However, while semiempirical MO methods can reproduce excitation energies of small-and medium-sized organic molecules with deviations of 0.4−0.5 eV, 169 benchmarks for excited states of transition metal complexes are missing. Even for the prediction of ground-state energetics of transition metal complexes, semiempirical MO methods can struggle when tested against DFT results. 170 As a matter of fact, we are not aware of any excited-state dynamic study of transition metal complexes using semiempirical MO methods.
A representative of DFT-based semiempirical methods is density functional tight-binding (DFTB). DFTB relies on a truncated Taylor expansion of the DFT energy with respect to the fluctuation of the electron density around a reference density, which is typically given by a superposition of atomic densities. 64 In analogy to the parent DFT approach, DFTB can be used in a linear-response time-dependent formulation to calculate excited states and has been used in this manner in excited-state simulations. 171 Excitation energies of small and medium-sized organic molecules obtained with linear-response time-dependent DFTB can reproduce TDDFT results using the PBE functional with deviations of ∼0.2 eV. 172 Furthermore, DFTB approaches can reproduce DFT geometries of transition metal complexes reasonably. 173 Even if benchmark studies of DFTB excitation energies are missing for transition metal complexes, excited-state dynamics studies of transition metal complex should emerge in the near future.

Wave Packet Propagation on a Grid
Having selected a suitable method to compute the electronicstructure ingredients, we can now tackle the simulation of the nuclear motion given by the time-dependent nuclear Schrodinger equation: nuc tot nuc (10) where the nuclear and electron degrees of freedom are separated, as explained in section 2.1. Integrating this equation, one can obtain 174 the nuclear wave function wherê t t ( , ) 0 is a propagator that evolves the wave function from the initial time t 0 to the final time t, and̂ is the timeordering operator. Often, the so-called split-operator method 175 is used to numerically evaluate this propagator.
The general solution for eq 10 is a nuclear wave function represented as a linear combination of specific timeindependent basis functions that is also known as a wave packet. This wave packet contains all quantum effects and can split in the presence of couplings. It is usually discretized on a spatial grid along the 3N − 6 degrees of freedom of the PES of the N-atomic molecular system (see Figure 3a). Assuming M grid points for each degree of freedom, a full dimensional wave packet propagation requires the precalculation of M 3N−6 grids points, which becomes quickly prohibitively expensive for all but the smallest molecules. This curse of dimensionality enforces the practical application of grid-based wave packet methods to very small systems orin most of the casesto a selection of only few relevant nuclear degrees of freedom.
In the case of transition metal complexes, this meant performing wave packet dynamics simulations along one 92,100,103−105,107−109,121 or two dimensions. [88][89][90][91][93][94][95][96][97][98]102,106 Most of these studies 88−98,100,102−107 were concerned with ultrafast ligand dissociation in hydride and carbonyl complexes upon photoexcitation, for which the natural coordinate of choice was the bond distance between the metal and one or two of the ligands to detach. How to come up with an optimal lowdimensional coordinate space, for a general molecule and reaction to study, is not necessarily an easy problem and different approaches are employed; see section 5.4. In passing we note that while some of these studies incorporated explicitly laser pulses to initiate the dynamics or even to guide it, 100,104−107 the laser−matter interaction is excluded in many others and excitation is assumed to take place instantaneously (see section 9).

Multiconfigurational Time-Dependent Hartree
The multiconfigurational time-dependent Hartree method (MCTDH) 176 is an efficient form of quantum dynamics, where the nuclear wave function is expanded in a set of singleparticle functions (SPFs) φ as ( ) are the MCTDH expansion coefficients. This approach is analogous to the MC treatment that was introduced in the electronic structure theory (see section 4.2). One subtle difference between both fields is that for the nuclear problem, (symmetric) Hartree products are used in the configurations, while the electronic problem requires (antisymmetrized) Slater determinants that change sign upon exchanging two particles, following the Pauli principle. Similar to the case of propagating wave packets on a grid (Section 5.1), the computational cost of MCTDH also scales exponentially with the degrees of freedom included. The important difference, however, is that in this case the exponential scaling is given by the number of SPFs per coordinate instead of the number of grid points. As the number of SPFs needed is smaller than the number of grid points, MCTDH is computationally more economic and therefore allows to consider a larger number of degrees of freedom than by propagating on a grid. In addition, a number of strategies have been devised in the last years to increase the number of dimensions, e.g., by combining several individual modes through using multimode SPFs 177 or in the multilayer MCTDH variant. 178,179

Vibronic Coupling Models
An extended strategy to harness the efficiency of MCTDH is to employ vibronic coupling models 180 to describe the PESs on which the wave packets can be propagated. In a vibronic coupling model, the PES are expanded in a Taylor series around a reference geometry R = Q 0 (usually the Franck−Condon geometry) using mass-frequency scaled normal coordinates Q (see Figure 4). The Taylor series is often truncated after the first (linear) termwhat is then known as the linear vibronic coupling (LVC) modelso that the PES E el values are approximated in terms of the ground-state PES E 0 el and linear vibronic coupling terms W The ground-state PESs are approximated as harmonic oscillators with frequencies ω i The coupling terms read where E Q ( ) By definition, the use of LVC models is limited due to the harmonic approximation of the potentials. This approximation neglects anharmonic effects that can be essential in different situations, e.g., to describe torsional motion or dissociation. Furthermore, due to the parametrization, the LVC potentials can only describe nuclear motion in the vicinity of the reference geometry. Thus, LVC models work best in rigid molecules.
Despite these limitations, LVC has become the standard approach to calculate PES in wave packet dynamics simulations u s i n g M C T D H f o r t r a n s i t i o n m e t a l c o mplexes, 48,101,[116][117][118][119][120][122][123][124][125][126][136][137][138]140 making it possible to include up to 16 nuclear degrees of freedom. 136 One example of MCTDH using 15 degrees of freedom in a heme−CO complex is shown in Figure 5. Among these studies, it is worth mentioning that only two include explicitly the excitation by a laser pulse. 48,125 The convenience of LVC models is that the computational effort in the electronic structure step is mostly reduced to determining the coupling parameters, and those can be obtained from a small number of calculations for each normal mode individually.

Choice of Degrees of Freedom
Due to the curse of dimensionality, regardless whether one performs wave packet propagation on a grid or using MCTDH as well as combined with LVC models, a difficult decision is always the selection of how many and which degrees of freedom need be considered, i.e., which are the most important coordinates that describe the problem at hand. Going beyond natural dissociation coordinates and involving the degrees of freedom that connect the Franck−Condon region with one (or more) conical intersections and excited-state intermediates is a task that in most cases goes beyond chemical intuition.
One approach thereby is to select normal modes based on the size of their vibronic coupling terms. 120,136,137,181 As large vibronic coupling elements are needed to efficiently transfer population between the electronic states, these coupling modes are necessary to describe the excited-state dynamics. This selection can be extended by adding tuning modes, which are normal modes that are responsible for the largest displacements in the excited-state dynamics by reaching toward the excitedstate minima 48,123−126 and the excited-state crossing points. [117][118][119]122 A different strategy to identify an optimal coordinate subspace is to use lower-cost dynamics methods that allow including all or a very large amount of degrees of freedom and then identify a posteriori the most important ones. These can be then selectively considered in more accurate quantum dynamical approaches. Ab initio molecular dynamics approaches, that will be introduced next, belong to this category of low-cost methods.

AB INITIO MOLECULAR DYNAMICS
One alternative to wave-packet-based dynamics is ab initio molecular dynamics (AIMD). In AIMD, the nuclei are described as classical particles that follow Newton's (classical) laws of motion on electronic potentials E el obtained by quantumchemical methods This is only an approximate description of the nuclear motion.
In reality, the motion naturally follows the laws of quantum mechanics. Accordingly, by definition, AIMD excludes nuclear quantum effects such as tunneling or coherence in the nuclar motion, and canat bestbe corrected a posteriori. Describing the motion of the nuclei classically in AIMD, however, introduces a huge practical advantage for dynamics. As the motion of the nuclei follows a (classical) trajectory R A (t) that is at each time step determined only by the current molecular geometry, the (exponentially scaling) precomputing part of the entire PES is lifted off. Instead, the electronic-structure calculations can be performed "on-the-fly" during the dynamics simulation, whereby the necessary properties to propagate the nuclear dynamics such as electronic potentials and their gradients are only calculated at the current geometry. We note that, in principle, it is also possible to generate a PES within certain approximations where to run classical trajectories. In one case, a semiglobal PES of a copper(I) complex in solution was obtained from molecular dynamics trajectories, albeit on uncoupled S 0 and S 1 states. 182 As it will be describe later, LVC models can also be used to run AIMD trajectories. In all cases though, the classical nature of the nuclei implies that a swarm of trajectories to be propagated is needed, instead of the one single propagation required in wave packet dynamics.
When several coupled electronic states are considered, two problems appear due to the nature of the classical trajectory approximation. One is that, unlike the wave packet that can split in the presence of couplings (recall Figure 3a), in AIMD a recipe is needed to transfer classical particles between different electronic states. The other is that again, unlike a wave packet that spreads over different electronic states and each portion follows the gradient in its corresponding PES, in AIMD every classical particle is confined to a single point of the PES and follows a single gradient that has be to be decided somehow.
In the following, two of the AIMD methods, which have been used up to now for excited-state dynamics of transition metal complexes will be described.

Surface Hopping
Probably the most popular AIMD approach that includes a mechanism to transfer population between different electronic states is surface hopping (SH). 183,184 In SH dynamics, the electronic wave function is allowed to spread over different electronic states as it is expressed as a linear combination of several electronic states: el j j j elec. states (17) Its time evolution is determined by the time-dependent electronic Schrodinger equation (in analogy to eq 10):  (19) Note that in this standard formulation, the Hamiltonian already excludes the light−matter interaction and is restricted to the electronic states of the system. The first term in the parentheses in eq 19 is the coupling between the different electronic states, while the second term can be computed using the NAC between the electronic states and the velocity of the nuclei v R: The trajectory in SH follows the gradient of a single electronic state, the so-called active state, Φ i . After every time step in the simulation, the trajectory is allowed to change the active state, i.e., to "hop" to a different electronic PES (see Figure 3b) with a certain probability. This probability is often calculated using the fewest-switches criterion, 183 which ensures a minimum number of hops along the propagation. This criterion prevents the system from effectively traveling along an averaged gradient in the unfortunate case of a system hopping every time step. The probability P i→j for a hop from initial state Φ i to final state Φ j can be expressed as This equation shows that the probability for a hop can only become large in the presence of large NACsthrough the second term in brackets (cf. eq 20)and when the electronic wave function of the system has already sizable admixture of the final state−through the coefficient c j (t) and the thus necessary * < c t c t ( ) ( ) 1 i i . After a trajectory hops from one electronic state to another, its potential energy changes instantly. To conserve the total energy, its kinetic energy needs to be adjusted. This is done by rescaling the momenta of the nuclei, which, in practice, is best achieved by rescaling along the direction of the NACs. 185,186 A problem appears when a hop should occur according to the probabilities calculated from the electronic wave functions (eq 21), but the trajectory has insufficient nuclear kinetic energy to compensate for the potential energy change during the hop. In a fully quantum description, such a transition can be allowed due to the tunneling effect. However, in the classical description used in SH, such transitionsreferred to as f rustrated hopsare not allowed. Consequently, standard SH is not able to describe processes involving tunneling effects, although using a modified hopping criterion, it is possible to explore tunneling pathways qualitatively (see section 10.4).
In order to mimic a wave packet, AIMD methods employ a swarm of trajectories, which in SH are independent of each other. While following a single path along the nuclear coordinates, the propagation of the electronic wave function still faces the problem that it is completely coherent. 187,188 This means that all parts of the electronic wave function, even if they JACS Au pubs.acs.org/jacsau Perspective are in different electronic states, are all propagated along the same gradient, which is the gradient of the active state. This description is erroneous. Instead, each part of a wave packet should experience the gradient of the electronic state that it occupies and be moved with individual velocities according to corresponding gradient, thus, losing the coherence of motion among them over time. As a remedy, SH simulations employ different types of so-called decoherence corrections. 185,189 For example, in the easiest from all, the energy-based decoherence correction, 190 the electronic populations on the nonactive states are continuously damped at each time step. The decoherence time that determines the rate of this damping thereby depends on the energy difference between the active and nonactive states as well as the kinetic energy of the trajectory. 191 Thus, the larger the energy gap between the states and the faster the system moves, the faster the electronic populations decohere. 6.1.1. Cost of Surface Hopping Simulations. The computational cost of a SH simulation is basically determined by the underlying on-the-fly electronic structure calculations. The total cost depends on the number of trajectories propagated. However, as in SH all trajectories are independent, their calculation can be well-parallelized. A SH trajectory uses a typical nuclear time step of 0.5 fs. This means, for a total simulation time of, say, 500 fs, we need 1000 time steps, and if a swarm of 100 trajectories is considered, this adds up to 10 5 electronic structure calculations that need to be performed. For reference, this number is comparable to the number of calculations necessary to precompute a five-dimensional PES with 10 grid points in each dimension where to propagate a wave packet. The computational advantage of SH methods explains why it has been extensively used in the last decades 192,193 to study the excited states in a broad variety of organic materials and, to a lesser extent, also in transition metal complexes. 42,99,[130][131][132][133][134][135]141 On passing we note that, from these studies, only one 42 considered the effect of a laser excitation explicitly.
The implementation of the LVC model within SH 194 has enabled further efficiency and therefore a cheaper application to transition metal complexes. 69,[127][128][129]142 In this way, it is possible to deal with systems with more than 100 nuclear degrees of freedom, propagate for several picoseconds, and consider thousands of trajectoriesa previously inaccessible venture for on-the-fly SH. Figure 6 exemplifies the capabilities of SH dynamics using LVC potentials defined for 166 normal modes of a ruthenium(II) polypyridyl complex and using almost 9000 trajectories. 129 On the dark side, however, even with LVC, it is sometimes not possible to include all normal modes. Especially low-frequency modes can experience nonphysically large displacements when describing the motion of a molecule in the basis of normal modes, and these modes need to be excluded. [127][128][129]142 6.1.2. Exploiting Surface Hopping to Find Relevant Degrees of Freedom. As discussed in section 5.4, the selection of nuclear degrees of freedom to be included in wave packet/ MCTDH simulations can be challenging. An interesting approach is then to use a combination of wave packet and SH methods. 139,195 For instance, SH simulations were performed for a copper(I) complex in solution including all vibrational degrees of freedom. 139 Using a principal component analysis, the dominant normal modes activated during the SH excitedstate decay were identified and could be used in a subsequent, more accurate wave packet dynamics simulation.
One way to identify the important normal modes that can be obtained from a single simulation run is to follow the activity of each normal mode during the dynamics. However, there exist also more sophisticated approaches that take into account the coupling of the nuclear motion to the evolution of the electronic state population. For example, normal mode coherence or correlation analyses include the comparison of the motion of AIMD trajectories in excited states and in the ground state or monitor the effect of normal modes on excitation energies, energy gaps, and the overlaps between electronic state wave functions. 196 Furthermore, the FrozeNM algorithm can be used to freeze normal modes and observe the effect that their exclusion has on the time evolution of the electronic states. 197 Finally, a machine-learning algorithm has been developed that can identify global reaction coordinates in excited-state reactions from AIMD simulations in an automatic manner, given that the AIMD simulations provide sufficiently large data sets for statistical evaluation. 198 As an alternative avenue, SH simulations performed using LVC potentials are so efficient that they can be gradually repeated reducing its dimensionality until the differences to the full dimensional calculation are acceptably small. In this way, a minimum set of degrees of freedom can be identified, as illustrated in a small platinum(IV) complex that could be reduced from its 15-dimensional space initially considering 200 electronic states to a nine-dimensional problem with 76 electronic states without loss of accuracy. 195

Ab Initio Multiple Spawning
A number of other methods to study excited-state dynamics exist, which in terms of cost and approximations can be placed formally between the MCTDH and SH formalisms. 63 From them, only the ab initio multiple spawning (AIMS) approach 199 has been employed in the excited-state dynamics of transition metal complexes, an iron(II) complex 110 (see Figure 7). AIMS is derived from full multiple spawning (FMS), 200 in which Gaussian functionsalso referred to as trajectory basis functions (TBF)are propagated on classical trajectories. When TBFs enter regions with high probability to transfer population between different electronic states (regions with strong NACs), new TBFs are spawned on the electronic states that are encountered, and both the initial as well as the spawned TBFs are propagated further (see Figure 3c). In contrast to other methods using classical trajectories, FMS requires the precomputation of the complete PESs to mediate the coupling between the TBFs, and it is formally exact.
In AIMS, the couplings between the TBFs are calculated only locally around the regions where the TBFs are getting close to each other. Using this approximation, AIMS simulation can also be performed on-the-fly. Furthermore, AIMS simulations usually employ the independent-first-generation approximation, in which all initial (parent) TBFs are independent; only the spawned (child) TBFs stay coupled to the initial TBFs. This approximation can be justified by assuming that the nuclear wave packet will usually spread rapidly in phase space in the beginning of the dynamics, which then would allow to neglect the coupling between the parent TBFs. This approximation is exaggerated further in SH − substituting spawns by hops−in which there is no coupling at all between the trajectories. It could likewise be justified by assuming that the nuclear wave packet spreads rapidly also after nonadiabatic events.
With the parent and child TBFs coupled in AIMS, the computational demand scales quadratically with the number of TBFs, thus making AIMS more costly than SH when comparable numbers of trajectories and TBFs are used. The higher scaling may be alleviated by introducing a stochasticselection scheme, in which spawned TBFs can be removed. 201 When the coupling between TBFs becomes small, one of the coupled TBFs is selected, the population of the other TBFs is collapsed into the selected TBF, and only the selected TBF is propagated further. In this manner, approximate AIMS simulations could be run at similar cost as SH simulations with results close to that of standard AIMS 202 so far without

JACS Au
pubs.acs.org/jacsau Perspective spin−orbit couplings 203 and not yet applied to any transition metal complex. Furthermore, AIMS simulations can be made to describe tunneling dynamics by allowing to spawn TBFs in the same electronic state, e.g., when the distance between the tunneling particle and its donor particle surpasses a certain threshold. 204

ENVIRONMENTAL EFFECTS
Very often, the phenomena one is interested in occur in an environment, being either a solvent or a biological surrounding structure. Accordingly, dynamics simulations of transition metal complexes should be simulated in the same media. In the following, we discuss two approaches that are readily used for excited-state dynamics: in one the environment is included explicitly and in the other, typically a solvent is only accounted for implicitly (see Figure 8).

Explicit Environments: QM/MM Partitioning
In section 4.4, we discussed the introduction of semiempirical approximations in the electronic structure calculations in order to make calculations of large molecules feasible. When a large amount of, for example, solvent molecules should be included in the calculation, semiempirical methods are rarely enough. In these cases, the system can be further approximated with a partition in two (or more) regions, where one of themat least the transition metal complexis treated quantum mechanically (QM) and the rest is only accounted for with molecular mechanics (MM), i.e., replaced by parametrized force fields. 206 Force fields contain classical energy expressions for bond lengths, bond angles, and dihedral angles as well as long-range interaction terms such as van-der-Waals and electrostatic interactions. This classical approach is computationally very economic, and it allows to simulate dynamics of systems with >10 6 atoms. Therefore, combined with the QM methods to describe the electronic excited states, it is ideal to treat transition metal complexes in solution or in an biological environment. Depending how the interaction between the regions is defined, different QM/MM methods exist. 207,208 In all schemes, the computational effort is mostly due to the expense of the QM part of the calculation. Therefore, the size limits for a typical QM region in hybrid QM/MM calculations is not larger than the molecular size limit for the QM computation of the isolated molecule. Due to the structural flexibility of biological environments or solvents, PES cannot be characterized by unique points such as global energy minima or minimum-energy crossing points in QM/MM calculations. Instead, several thermally accessible minima that can be populated exist and should be properly sampled, 209 increasing the complexity of the calculations.
QM/MM approaches have been applied in nonadiabatic AIMD simulations of several transition metal complexes in solution. 42,130,132,139 In all of these studies, the transition metal complex is described alone in the QM region, while the bulk of the solvent (∼10 3 molecules) are described using MM force fields, as in Figure 8b. Although QM/MM-AIMD simulations of organic systems embedded in biological environments exist, 39,210−214 we are not aware of any example with a transition metal complex.
Including environment effects using QM/MM approaches in wave packet dynamics is also possible, however, it requires more elaborate schemes. These schemes can include the modification of precomputed PES of the isolated molecule with energetic shifts calculated from ground-state QM/MM-MD simula-tions 126,215,216 or an iterative update of the solvent effects that is obtained from the simultaneous simulation of the solvent in an ensemble of classical trajectories, 217 with most applications so far treating only organic solutes. 215−217

Implicit Environments
An alternative approach to incorporate environmental effects is to replace their atomistic description by a dielectric continuum. This approach was created and is particularly useful to account for solvation. 218,219 Popular implementations of this approach include the polarizable continuum model (PCM) 220 and related variants, such as conductor-like screening model (COSMO) 221 or the SMD model. 222 In these models, the solute molecule is placed inside a cavity containing charges on its surface, through which the interaction between the solute molecule and the surrounding solvent continuum is described (see Figure 8c).
These models have been applied in excited-state dynamics of solvated transition metal complexes parametrized with LVC potentials using both MCTDH [117][118][119][120][121][122]126,136,138 and AIMD simulations. [127][128][129]142 Although, in principle, there exist stationary studies where also biological environments are approximately modeled with continuum models with very small dielectric constants, such strategy has not been used for dynamics, as it cannot capture the explicit fluctuations of the complex environments.
Something to keep in mind when describing an excitation process using a continuum model is that the solvent effects can be split into two contributions: a fast, dynamical component and a slow, inertial component. 218 The fast component describes the interaction of the electron density of the solvent with the electron density of the solute in its excited state that can be considered changed instantaneously after excitation. The slow component takes into account the situation that the solvent moleculesthough modeled as a dielectric continuumare still oriented referring to the initial electron density of the ground state of the solute. When describing the dynamical evolution after photoexcitation, the solvent effects will still be well-approximated by the slow component referring to the ground state, where they consist mainly on the modulation of the energy differences between the different electronic states. 223 However, at later times in the dynamics, when the solvent molecules start to adopt their orientation to the changed electron density in the excited state, this approximation becomes worse.

INITIAL CONDITIONS, ZERO-POINT ENERGY, AND TEMPERATURE
Excited-state dynamics simulations of nuclear motion requires the selection of an initial condition that describes the state of the nuclei before the excitation process. In wave packet/MCTDH dynamics, the initial wave function is typically the lowest vibrational state of the electronic ground state. When the PES are expressed in an LVC model, the initial wave function can be described by the ground-state wave function of an harmonic oscillator. Starting the system in the vibrational ground state Ψ 0 corresponds to describing the system at zero temperature, and the total energy of the system is then given by its zero-point energy (ZPE). Such situation is adequate to mimic gas phase experiments performed in a cold beam, as for example those described previously. 99,100,104,107 When the system adopts a finite temperature, vibrationally excited states become populated with temperature-dependent probabilities that can be obtained, e.g., from statistical JACS Au pubs.acs.org/jacsau Perspective sampling. 224 At room temperature, the populations of vibrationally excited states are small for high-frequency modes (>1000 cm −1 ); however their population grows at lower frequencies with contributions exceeding 50% around frequencies of 140 cm −1 . Despite its simplicity, there are no wave packet/MCTDH dynamics simulations of transition metal complexes including temperature. In AIMD simulations, the classical nuclei are described by their position and momenta, which need to be selected as initial conditions. Two common approaches exist for this purpose: Wigner sampling (also referred to as quantum sampling) and molecular dynamics sampling (classical sampling). 225,226 In the Wigner sampling, coordinates and momenta are sampled from a Wigner distribution 227 a simultaneous probability distribution of coordinates and momenta, that is a function of the nuclear wave function. For this, one typically employs the harmonicoscillator approximation for which analytical expressions of the nuclear wave functions are known. 228 Wigner sampling is commonly employed in the zero-temperature formalism, where the total energy of the system is again given by the ZPE. However, the effects of finite temperature can easily be accounted for by allowing the population of excited vibrational states. 229 In the classical sampling, initial position and momenta are taken from snapshots from MD trajectories propagated in the electronic ground state. In the ground-state MD simulation, the system is given a total energy of k B T and the simulations are typically performed at room temperature (T = 300 K). This energy, k B T at room temperature, is much smaller than the ZPE. The inclusion of the different amounts of vibrational energy from both sampling approaches can have an effect, e.g., on reaction rates in excited-state dynamics. 225 Since the total energy at room temperature is better approximated by the ZPE (zero-temperature) Wigner sampling is, thus, usually preferable over MD-based sampling. When studying large systems in a hybrid QM/MM schemewhich cannot be described completely by QMit is also possible to combine Wigner sampling for the QM part and classical sampling for the MM part. 230,231

INCLUDING EXPLICIT LIGHT−MOLECULE INTERACTIONS
One goal of nonadiabatic dynamics simulations is to study the processes occurring after excitation of the molecule into its excited electronic states. In spectroscopic experiments, this excitation is usually realized by irradiation of the molecule with a laser pulse centered around a specific wavelength, that is usually known as the initial pump pulse. Furthermore, time-resolved experiments need a second delayed pulse to probe the excitedstate reaction in time. This pulse sequence is then known as a pump−probe experiment. Probe pulses can ionize the molecule or excite as well as de-excite it to other electronic excited states. Naturally, it would be desirable that the simulations imitate the experiment as close as possible, in particular when a comparison or interpretation of a particular experiment is aimed at. Formally, the laser pulses can be easily added to the Hamiltonian (cf. eq 2) and in the field of quantum dynamics, there exist a number of studies of transition metal complexes using explicit laser pulses. 48,100,[104][105][106]125 One example is shown in Figure 9a, where the dissociation ratio of CO versus CH 3 is manipulated using a sequence of well-timed laser pulses: a pump laser that excites population into an intermediate state and a dump laser that de-excites it to the desired dissociative state. As time-resolved spectroscopic experiments typically use weak laser pulses that only excite a few percent of the total population, the simulation of the excitation process can be cumbersome, particularly in AIMD simulations. In wave packet dynamics propagations, one would choose laser parameters similar to the experimental ones and make sure that the corresponding small populations do not run into numerical problems. In AIMD, although laser interactions can be simulated, 233−235 the use of laser pulses that only excite a tiny fraction of the molecules is impractical 236 because most of the calculated initial conditions are disregarded. For that reason, explicit laser excitations to excite initially the molecule are ignored in AIMD simulations of transition metal complexes− with exception. 42 In general, for simplicity in wave packet dynamics and for economic reasons in AIMD simulations, the excitation is only considered implicitly through the appropriate selection of the initially excited electronic states. We note, however, that the simulation of the time-resolved experimental signal can be computed a posteriori from AIMD simulations, using the available trajectories to excite or ionize the system further. Figure 9b−d shows one example of a simulated timeresolved luminescence spectra signal for a rhenium complex in solution obtained with SH simulations. We also note that in this case, an excellent agreement between theory and experiment only became visible because the global experimental signal was directed simulated. 42 The excitation process in wave packet dynamics is often simulated through a vertical projection of the nuclear wave packet from the electronic ground state into the potential of the absorbing electronic state. This approach is valid in the limit of ultrashort laser pulses (so-called δ-pulses), where the nuclear wave packet in the excited state is simply a copy of the vibrational ground state in the electronic ground state scaled by the transition-dipole moment. 237 At the other extreme, the simulation of long laser pulses (a continuous wave) corresponds to exciting the molecule into a vibrational state of its excited electronic state. In reality, a laser pulse with a finite duration will achieve an intermediate situation, that is best taken into account if explicitly included in the simulation. We note that care has to be taken when including the laser pulse explicitly in the case of exciting into a set of degenerate electronic states, as this can lead to dynamics that are dependent on the polarization of the laser in the simulation. 125 In AIMD simulations, the ultrashort pulse limit for the initial excitation is usually emulated by computing the excited states of a large set of nuclear geometries (sampled as described in section 8) and then selecting the states and associated geometries by the size of the oscillator strength. 236,238 In order to keep a reasonable number of initial conditions, for practical purposes, the initially excited states are selected not at a specific wavelength but within an energy window around the intended excitation energy. As AIMD simulations require a large number of trajectories for statistical averaging, accordingly, a much larger number of trial states need to be generated to start trajectories only in states with (relatively) large oscillator strengths. Typical excitation energy windows range from 0.2 eV 129,142 to 0.5 eV. 69,128,131,134 In some cases, this selection scheme is also approximated by starting trajectories in one specific state based on its large oscillator strength at the Franck−Condon geometry. 110,130,132,133,135 It is also possible to think about an excitation process different from the coherent light coming from the laser pulses typically employed in femtosecond spectroscopy: the interaction of molecules with natural thermal light coming from sources such as the sunas encountered in biological processes. 239 Besides the broader spectral excitation range, natural thermal light is stationary incoherent irradiation and may, thus, trigger other responses than observed when using coherent pulsed laser irradiation. 240 The incoherent nature of natural thermal light during the excitation process can be described through an ensemble of coherent pulses, whereby all effects of incoherence are recast in the postexcitation averaging. 241 The effects of stationary incoherent irradiation can then be implemented in nonadiabatic dynamics simulations through running an ensemble of dynamics simulations, displacing the individual runs in time and using a weighted averaging scheme over all runs. 242 10. BRIDGING TIME SCALES 10.1. State of the Art With the standard techniques described until now, the time scales that can effectively be simulated using nonadiabatic dynamics are basically determined by their computational cost (see an overview in Figure 10). For wave packet propagations (recall section 5), this cost is hidden in the precalculation or parametrization of the PES and scales exponentially with the number of nuclear degrees of freedom. Accordingly, nowadays, simulation times for transition metal complexes range from 1 ps on 15-dimensional PES 117 over 4 ps on five-dimensional PES 48 to 100 ps on two-dimensional PES. 138 In on-the-fly AIMD methods (section 6), the computational cost is dominated by the electronic structure calculation that is performed at every time step. A typical simulation time is 1 ps, e.g., in a ruthenium(II) complex with 138 nuclear degrees of freedom including three singlet and three triplet states. 134 If using parametrized potentials, the time span of AIMD simulations can be substantially extended. For example, using LVC parametrized potentials, SH were run for several picoseconds for iron(II) and vanadium(III) complexes on 244 and 197-dimensional PES including total numbers of 60 and 25 electronic states, respectively. 128,142 In general, the cost of SH simulations on LVC potentials is low enough to perform even longer dynamics simulations. However, since the parametrization of the LVC potentials is usually conducted at the Franck− Condon geometry, the description of the dynamics becomes worse the farther the system moves away from the Franck− Condon geometry. Thus, simulations based on LVC potentials become less reliable for longer simulation times−both using AIMD and wave packet dynamics. In the following, we describe several scenarios that envision how to bridge the gap of large time scales and what can expected to be soon used for the excited-state dynamics of transition metal complexes.

Machine Learning
A promising alternative to extend the time scales is to exploit parametrized PES from machine learning (ML) algorithms, 245−247 although this approach has been only applied to relatively small organic molecules. ML algorithms use parametrized functions that have been fitted to a reference test set− the training set−in order to predict properties of interest for new data points that are not in the training set. For example, ML algorithms were trained on excitation energies and oscillator strengths at a set of different ground-state geometries and used to predict these properties at other ground-state geometries and thus simulate the complete absorption spectrum. 248,249 The first applications of ML algorithms to nonadiabatic dynamics simulations have been reported on on-the-fly computed PES. 243,244,250−254 Studies thereby highlighted difficulties to describe the dynamics in regions of strong NACs or close to conical intersections when using ML with kernel ridge regression algorithms. 251,252 In contrast, dynamics based on ML featuring neural networks seem to handle these problematic regions better. 243,244,253,254 The ability of the ML algorithm to predict properties depending on the training data introduces two problems in the simulation of nonadiabatic dynamics. First, in order to describe molecular properties correctly in all regions of the PES that the molecule visits during the dynamics, all regions must be taken into account in the training of ML algorithm. For large, polyatomic molecules, the knowledge about which regions are important is usually not available a priori. This problem can, however, be solved within the ML algorithm itself through adaptive sampling. 243,244 In this approach, two or more ML algorithms are trained independently and their predictions are compared during the simulation. When the differences between their predictions at a certain molecular geometry surpass a predefined threshold, the geometry is assumed to lie in a region of insufficient training data. The properties at this geometry are then calculated using a reference method, the data are added into the training data, and the ML algorithm can, thus, be improved iteratively. In this way, dynamics on parametrized potentials can be simulated also for long-time scales up to nanosecond in a meaningful manner. 243,244 The second problem concerns the choice of the reference method for the training data. As the ML algorithm is trained to reproduce the results of the reference method, the choice of the reference method defines the ultimate accuracy of the ML algorithm. Since typical training sets need thousands of training points, the compilation of the training data can itself become computationally expensive. However, also this task can be alleviated using ML techniques, i.e., Δ machine learning. 255,256 In Δ-ML, an algorithm can be trained to predict corrections to computationally inexpensive electronic structure methods through comparison to more accurate reference data. This approach, however, has not yet been applied in combination with excited-state dynamics simulations.

Metadynamics
Metadynamics is a technique to accelerate rare events occurring in the electronic ground state. 257 In metasurface hopping, 258 it is possible to speed-up the dynamics and sample faster reaction pathways by accelerating the rate of nonadiabatic transitions. Since this rate depends linearly on the nuclear kinetic energy and quadratically on the NACs, it can be accelerated by scaling up the NACs; this has been first demonstrated in the relaxation rate from the S 1 state of a small organic molecule. 258 A related approach to metasurface hopping is given by multistate metadynamics, 259,260 which extends the original (ground-state) metadynamics approach 257 to multiple electronic states. In metadynamics, the dynamics of a molecule are started in its ground state equilibrium, and successively external potentials are added to drive the molecular dynamics toward a desired region. For multistate metadynamics, this desired region is the conical intersection between the ground and the first excited state. Thus, instead of monitoring the decay to the ground state starting from the excited state, the possible reaction pathways are searched in reverse direction.

Rare Event Sampling
The rate of a process in nonadiabatic dynamics depends on the probability that it occurs. In wave packet dynamics simulations, low-probability processes occurring on a long time scale can be captured by extending the simulation time (within reason), as the wave packet at every time step is allowed to spread over all accessible reaction pathways. In AIMD simulations, such as SH, where trajectories follow a single pathway, describing such rare events can be difficult, as besides the cost of the long simulation, a sufficiently large ensemble of trajectories is required. The margin of error of a random statistical sampling is approximately given by ε = where z γ is the quantile of the standard deviation. This means that for a 95% confidence interval (z 0.95 = 0.98), a process that accounts for 10 or 1% of the reaction yield requires ∼100 or 10000 trajectories, respectively, for statistical meaningful results. In order to sample even rarer events, the army-ants algorithm has been developed for SH, 261 where trajectories are JACS Au pubs.acs.org/jacsau Perspective stochastically driven to explore also low-probability channels. This significantly reduces the number of trajectories needed to converge results for simulating processes with very low probabilities (∼10 −8 ). The army-ants algorithm has also been used to model tunneling paths employing classical trajectories. This is done by pushing the trajectories to go beyond the classical turning points in ground-state molecular dynamics simulations 262 and mean-field Ehrenfest dynamics including excited states. 263,264

Rate Methods
One example of a process that can be slow is intersystem crossing. Fast intersystem crossing events can be monitored by doing excited-state simulations including SOCs (cf. eq 6). However, when intersystem crossing is too slow, time scales can be obtained indirectly from static calculation of intersystem crossing rates, 12 The molecular wave functions Ψ mol are thereby taken as a product of an electronic and a nuclear part, the latter typically approximated by harmonic oscillators, and all coupling elements are calculated at a reference geometry for which usually the equilibrium geometry of the initial state Ψ i mol is selected. This approach has been used to calculate intersystem crossing rates for a number of transition metal complexes. 266−272 Due to the selection of a reference geometry where the intersystem crossing is to occur and the assumption that the system resides in an initial state at the beginning of the intersystem crossing process, these static approaches are limited in the processes that they are able to describe. Thus, most of the applications 266−269 focus on the intersystem crossing connecting the lowest-excited singlet and triplet states S 1 and T 1 , while few others 270,271 also take into account the T 2 state. A different strategy was adopted in ref 272, where intersystem crossing rates were calculated for a network of several singlet, triplet, and quintet states of an iron(II) complex. In this case, the calculations were performed at the ground-state equilibrium geometry, which was deemed reasonable based on the ultrafast nature of the intersystem crossing established from experimental studies, i.e., giving the system only few tens of picoseconds to move away from the ground-state equilibrium geometry after excitation before ISC occurs.

CONCLUSIONS AND CHALLENGES AHEAD
The field of excited-state simulations on transition metal complexes is now about three decades old. It started with exact quantum dynamics simulations in reduced dimensionality (one or two degrees of freedom), but it has been only very recently that it could make the leap to include all the important degrees of freedom or even full dimensionality in selected complexes. The latter was possible either by exploiting parametrized potentials or/and in combination with on-the-fly methods, such as surface hoping or multiple spawning. With these techniques established, one can expect a booming in the field that will help rationalize new experiments and make exciting predictions comparable to the pace that organic photophysics and photochemistry experienced in the last 30 years.
However, the limitations of the methods are still plenty and many are the challenges ahead. Even today, the dynamical methods that can applied to transition metal complexes could be basically grouped at two extremes: one that rely on quantum nuclei but sacrifices degrees of freedom and another that can include all the degrees of freedom but sacrifices the quantum nature of the nuclei. Unfortunately, as of today, the best of both worlds is only possible for the smallest organic systems but not for transition metal complexes. These means that often, depending on the molecule or even its electronic structure (which a priori is difficult to know), one strategy or the other can be more suitable (or the only one viable) and this selection is not a lighthearted decision to take.
One additional aspect is the inclusion of explicit light−matter interactions in the calculations, be it incoherent light sources such as the sun, a coherent excitation with a laser, the simulation of a full pump−probe experiment, or even the application of tailored (optimal) laser pulses aimed to control nuclear and electronic properties of a transition metal complex. While early studies based on wave packets could describe these interactions easily (albeit in reduced dimensionality), the methods based on trajectories that can to date be applied to transition metal complexes still face important difficulties in this endeavor. 273,274 In addition, one should realize that the suitability of the underlying electronic structure method is always key for successregardless of the chosen dynamical strategyand its choice can hamper the efficiency and feasibility of the dynamical method itself, as well as the quality of the final results. For instance, some transition metal complexes might require multireference methods, which, while they might be available for stationary calculations, are not affordable for on-the-fly dynamical methods. Thus, if the applicability of a multireference method in reduced-dimensional models is also too expensive or too simplistic (despite its application also being an intrinsic challenge itself!), one is left with no choice for dynamics to look at.
In this Perspective, we have reviewed every dynamical study of transition metal complexes (to the best of our knowledge), and we could see that the list of "must-haves" for the next years is very long. We are awaiting to see recently developed electronic structure methods, such as DMRG, combined with any method for nuclear dynamics, to see more sophisticated on-the-fly techniques extended to transition metal complexes, to see more systematic studies where spectroscopic observables are included and can be directly compared with experimental signals, and to see the capabilities of machine learning potentials applied to transition metal complexes and studies based on metadynamics or rare event sampling on this field. Of course, simulations should be also be extended to include the description of the environment, beyond solution, as there are a lot of interesting applications of transition metal complexes embedded in DNA, proteins, and other biological environments, as well as molecular machines and other supramolecular assemblies at the nanoscale. A considerable challenge is also posed by polynuclear coordination complexes, where two or more metals are simultaneously light-active and responsible for the photophysics of the system. Besides the growth of the electronic structure problem, such complexes display magnetic exchange interactions, whose description with quantum chemical methods pose a challenge by itself. For such complicated problems, the advent of quantum computation 275 could bring a new perspective by exploiting their unique features of superposition and entanglement. The first steps in the use of quantum algorithms for nonadiabatic dynamics 276 have already been proposed. We thus believe that the next decade will see an increase of fundamental new ideas to harness the excited-state dynamics of transition metal complexes for the benefit of society and humankind in many areas.