Saikat Mukherjee
email: saikat.mukherjee@univ-amu.fr.

Sergey Dmitry A Fedorov
email: dfedorov@nevada.unr.edu

Sergey A Varganov
email: svarganov@unr.edu.

Modeling Spin-Crossover Dynamics

Keywords:

INTRODUCTION

Nonradiative transitions between electronic states with different spin multiplicity play important roles in many areas of molecular science, as evident from the different names used for these processes, including spin-crossovers [START_REF] Halcrow | Spin-Crossover Materials: Properties and Applications[END_REF][START_REF] Swart | Spin States in Biochemistry and Inorganic Chemistry[END_REF], intersystem crossings (3), spin-forbidden reactions [START_REF] Harvey | Spin-forbidden reactions: computational insight into mechanisms and kinetics[END_REF], and two/multi-state reactivity [START_REF] Shaik | Two-state Reactivity: Personal Recounting of Its Conception and Future Prospects[END_REF]. These transitions, we will call them spincrossovers for consistency, can be initiated by various stimuli, such as light, temperature, pressure, and an external magnetic field. The concepts of spin-crossover and internal conversion (IC), transitions between electronic states with the same spin, are rooted in two approximations central to quantum chemistry. The Born-Oppenheimer approximation (BOA) separates electronic and nuclear motions, while the spin-free Hamiltonian in the non-relativistic Schrödinger equation decouples electronic motion from electronic spin. These approximations lead to nuclei propagating adiabatically on a single electronic state with a well-defined spin. Nonadiabatic processes characterized by the energy transfer between the nuclear, electronic and spin degrees of freedom are forbidden. However, if an energy gap between electronic states becomes small enough to be comparable with the nuclear kinetic energy or the relativistic spin-dependent energy, these approximations break down, and transitions between different electronic states become allowed. The probability of transition between the electronic states with the same spin depends on the nonadiabatic coupling (NAC) defined in the next section. The NAC between the electronic states with different spin vanishes due to the orthogonality of the spin eigenfunctions.

Spin crossovers are mediated by the spin-dependent couplings (SDC), which mix the electronic states of different spin and can arise from the spin-orbit, spin-spin, hyperfine and an external magnetic field interactions. Here, we will mostly focus on the spin-orbit coupling (SOC), which is usually the largest contribution to SDC.

Modeling spin-crossovers can provide insight into photochemical reaction mechanisms, help to interpret complex vibronic spectra, calculate the rates of formally spin-forbidden reactions, and predict the lifetimes of excited electronic states. While time-independent statistical theories can be used to study spin-crossover kinetics (see the sidebar Nonadiabatic Statistical Theories), the non-equilibrium nature of spin-crossovers often requires the use of nonadiabatic molecular dynamics (NAMD), which propagates nuclear and electronic degrees of freedom (DOF) in time [START_REF] Agostini | Different flavors of nonadiabatic molecular dynamics[END_REF]. The NAMD methods can be separated into two types. The grid-based methods propagate nuclear DOF on precomputed electronic potential energy surfaces (PESs), while in the direct dynamics methods, the electronic properties are calculated "on-the-fly" as nuclei propagate. Both types rely on electronic structure methods to obtain energies, energy gradients, and couplings for multiple electronic states at different nuclear geometries. In general, because spin-crossovers and ICs can happen on a similar time scale, both types of nonadiabatic transitions have to be considered in NAMD simulations.

In NAMD, nuclear DOF can be propagated in time on electronic states obtained in several different ways (Figure 1). Most electronic structure calculations solve the spin-free timeindependent Schrödinger equation (TISE) within BOA producing a separate manifold of adiabatic electronic states for each spin multiplicity. As their name suggests, these states correspond to electrons instantaneously adjusting to nuclear motion and can have mixed character (for example, / character). However, because the Hamiltonian has no spindependent terms, the states from different spin manifolds do not mix. We will call this representation spin-diabatic [START_REF] Granucci | Surface hopping trajectory simulations with spinorbit and dynamical couplings[END_REF]; it also has been called the molecular Coulomb Hamiltonian (MCH) representation [START_REF] Mai | A general method to describe intersystem crossing dynamics in trajectory surface hopping[END_REF]. The breakdown of BOA introduces NAC between the spin-diabatic states with the same spin, and the addition of the spin-dependent terms into the Hamiltonian leads to SDC between the states of different spin (Figure 1b). While the spin-diabatic states are relatively easy to obtain, the crossings of -dimensional PESs of the same-spin states produce the -dimensional conical intersections with singular NAC (avoided crossings for ; is the number of nuclear DOF). These features present challenges for propagating nuclear DOF, especially if the PESs and coupling surfaces have to be fitted for grid-based NAMD. Smooth PESs and couplings can be obtained by transformation into the diabatic representation, where electronic states preserve their character, and the singular NAC is replaced by the smooth interstate coupling V (Figure 1a). However, for polyatomic molecules, only an approximate diabatization can be achieved [START_REF] Mead | Conditions for the definition of a strictly diabatic electronic basis for molecular systems[END_REF], and at least for some NAMD methods, the diabatic representation leads to less accurate results [START_REF] Mai | A general method to describe intersystem crossing dynamics in trajectory surface hopping[END_REF][START_REF] Tully | Mixed quantum-classical dynamics[END_REF]. Spin-diabatic states can be also transformed into the fully adiabatic representation by eliminating SDC at the expense of introducing the non-zero NAC between the resulting spin-mixed states (Figure 1c). While adiabatic PESs and NAC surfaces are not suitable for fitting, the NAC localized at conical intersections and the ability to describe each component of spin multiplets separately (see states 3, 4 and 5 in Figure 1c) make this representation desirable for direct dynamics [START_REF] Granucci | Surface hopping trajectory simulations with spinorbit and dynamical couplings[END_REF][START_REF] Mai | A general method to describe intersystem crossing dynamics in trajectory surface hopping[END_REF]. In general, performing NAMD simulations in fully adiabatic representation requires the implementations of analytical energy gradient and NAC for spin-mixed states, which is still a challenge; however, an approximate "on-the-fly" adiabatization has been implemented [START_REF] Mai | A general method to describe intersystem crossing dynamics in trajectory surface hopping[END_REF].

In section 2, we review the NAMD and some aspects of electronic structure methods for modeling spin-crossover dynamics. While many different flavors of NAMD have been developed to describe IC dynamics at conical intersections, and some of them potentially could be extended to model spin-dependent processes, we limit our discussion to the methods that have been applied to spin-crossovers. In section 3, we describe a few representative applications of these methods to molecular systems of different complexity. We conclude by posing several fundamental questions related to the spin-dependent processes. These questions should be possible to address with the future methodological developments in NAMD.

METHODS FOR MODELING SPIN-CROSSOVER DYNAMICS

In this section, we describe three types of NAMD used to model spin-crossover dynamics. The stochastic trajectory surface hopping (TSH) method solves the time-dependent Schrödinger equation (TDSE) for electrons while propagating classical nuclei. In the multiple spawning method, the nuclear wave function is expanded in the adaptive basis of frozen Gaussians propagating in time. In the multiconfiguration time-dependent Hartree (MCTDH) method, the nuclear wave function is represented as a superposition of multidimensional products of basis functions, where both the superposition coefficients and the basis functions are time-dependent. In principle, all three types of NAMD can be performed on the precomputed (usually diabatic) PESs or as direct dynamics (usually using the spin-diabatic or fully adiabatic representations). However, the strength of TSH and multiple spawning methods is in direct dynamics, while MCTDH was originally developed as a grid-based method. Finally, since any NAMD simulation requires accurate electronic energies, gradients and interstate couplings, we discuss several aspects of electronic structure methods used to model spin-crossover dynamics.

Trajectory surface hopping

In the TSH method, the nuclear wave packet on a particular electronic state is represented as a swarm of independent trajectories obeying the classical Newton equations of motion (EOM).

The force acting on the nuclei is equal to the negative energy gradient of the corresponding electronic state [START_REF] Tully | Trajectory Surface Hopping Approach to Nonadiabatic Molecular Collisions: The Reaction of H + with D 2[END_REF][START_REF] Barbatti | Nonadiabatic dynamics with trajectory surface hopping method[END_REF][START_REF] Wang | Recent Progress in Surface Hopping: 2011-2015[END_REF]. At each timestep, a stochastic process determines whether the system will propagate on the current electronic state or hop to another state leading to population transfer between the states. The classical propagation and purely local nature of nuclei make the TSH method easy to implement and widespread [START_REF] Marquetand | Nonadiabatic ab initio molecular dynamics including spin-orbit coupling and laser fields[END_REF][START_REF] Persico | An overview of nonadiabatic dynamics simulations methods, with focus on the direct approach versus the fitting of potential energy surfaces[END_REF][START_REF] Crespo-Otero | Recent Advances and Perspectives on Nonadiabatic Mixed Quantum-Classical Dynamics[END_REF]. Although TSH was originally formulated to describe ICs, it has been extended to model spin-crossovers [START_REF] Granucci | Surface hopping trajectory simulations with spinorbit and dynamical couplings[END_REF][START_REF] Richter | SHARC: ab Initio Molecular Dynamics with Surface Hopping in the Adiabatic Representation Including Arbitrary Couplings[END_REF][START_REF] Cui | Generalized trajectory surface-hopping method for internal conversion and intersystem crossing[END_REF][START_REF] Martínez-Fernández | Competing ultrafast intersystem crossing and internal conversion: a time resolved picture for the deactivation of 6-thioguanine[END_REF][START_REF] Zaari | Nonadiabatic Transition State Theory and Trajectory Surface Hopping Dynamics: Intersystem Crossing Between 3 B 1 and 1 A 1 States of SiH 2[END_REF]. Gonzalez and coworkers introduced the surface hopping including arbitrary couplings (SHARC) method with "on-the-fly" adiabatization of spin-diabatic states [START_REF] Mai | A general method to describe intersystem crossing dynamics in trajectory surface hopping[END_REF][START_REF] Richter | SHARC: ab Initio Molecular Dynamics with Surface Hopping in the Adiabatic Representation Including Arbitrary Couplings[END_REF]. Persico and coworkers included SOC in the TSH method to study spin-crossovers in both spin-diabatic and fully adiabatic representations [START_REF] Granucci | Surface hopping trajectory simulations with spinorbit and dynamical couplings[END_REF][START_REF] Favero | Dynamics of acetone photodissociation: a surface hopping study[END_REF].

We define the full electron-nuclear Hamiltonian as a sum of the nuclear kinetic energy , the spin-free electronic Hamiltonian and the SOC operator :

where r, R and s are vectors of the electronic, nuclear and spin DOF, respectively. The nuclearnuclear repulsion is included in . The time-dependent electronic wave function is expanded in the basis of time-independent electronic states that depend parametrically on nuclear positions:

Inserting this wave function into the electronic TDSE with the Hamiltonian yields the EOM for the coefficients (we will use the atomic units with ):

where a dot above a variable indicates a time derivative, and are the matrix elements of the spin-free and SOC Hamiltonians. The NAC is defined as , where is the nuclear velocity vector, and is the NAC vector. In the spin-diabatic representation, are the eigenfunctions of , and = , where are the electronic state energies. If states I and J have different spin quantum numbers and , vanishes due to the orthogonality of the spin eigenfunctions. If I and J are both singlet states or , vanishes. A detailed discussion of the SOC selection rules can be found in [START_REF] Marian | Spin-Orbit Coupling in Molecules[END_REF][START_REF] Fedorov | Spin-orbit coupling in molecules: Chemistry beyond the adiabatic approximation[END_REF]. In contrast to NAC, which is mostly localized around conical intersections, SOC only weakly depends on molecular geometry and can lead to interstate transitions far away from the crossing seam of two PESs [START_REF] Zaari | Nonadiabatic Transition State Theory and Trajectory Surface Hopping Dynamics: Intersystem Crossing Between 3 B 1 and 1 A 1 States of SiH 2[END_REF]. In the diabatic representation, the off-diagonal elements of are interstate couplings, while should vanish (in quasi-diabatic representations, usually some residual NAC is still present), and the SOC selection rules are the same as in the spin-diabatic representation. Finally, in the fully adiabatic representation, the spin-mixed states are the eigenfunctions of and are coupled by non-zero NAC. The results of TSH simulations are averaged over a swarm of trajectories with initial conditions sampled from some positionmomentum distribution [START_REF] Sun | Comparisons of classical and Wigner sampling of transition state energy levels for quasiclassical trajectory chemical dynamics simulations[END_REF][START_REF] Barbatti | Effects of different initial condition samplings on photodynamics and spectrum of pyrrole[END_REF].

Among various strategies to compute the probability of transitions between electronic states, Tully's fewest switches surface hopping (FSSH) algorithm [START_REF] Tully | Molecular dynamics with electronic transitions[END_REF] is the most common. In the generalized FSSH algorithm, the hopping probability from state I to state J in the time interval is where are the elements of the electronic density matrix. For a successful hop, two conditions must be simultaneously fulfilled. First, for a random number r selected from the interval [0,1], the following must be true:

Second, because the electronic energy changes as a result of the hop, to conserve the total energy, the nuclear kinetic energy must be adjusted by rescaling the nuclear velocity vector. If velocity rescaling along the NAC vector cannot compensate for the electronic energy change, the hop is rejected, and the velocity component along the NAC vector is reversed [START_REF] Hammes-Schiffer | Proton transfer in solution: Molecular dynamics with quantum transitions[END_REF]. This is known as a frustrated hop. Truhlar and coworkers introduced the fewest switches with time uncertainty (FSTU) method, where a classically frustrated hop can occur if an allowed hopping geometry can be reached within the Heisenberg interval of time uncertainty [START_REF] Jasper | Fewest-switches with time uncertainty: A modified trajectory surface-hopping algorithm with better accuracy for classically forbidden electronic transitions[END_REF]. The TSH dynamics can also be performed using the Landau-Zener [START_REF] Hu | Trajectory Surface Hopping Study of the O( 3 P) + Ethylene Reaction Dynamics[END_REF][START_REF] Rajak | Trajectory surface hopping study of the O( 3 P) + C 2 H 2 reaction dynamics: Effect of collision energy on the extent of intersystem crossing[END_REF] or more sophisticated Zhu-Nakamura formulas [START_REF] Xu | Intersystem crossing-branched excited-state intramolecular proton transfer for o-nitrophenol: An ab initio on-the-fly nonadiabatic molecular dynamic simulation[END_REF][START_REF] Yue | Benchmark Performance of Global Switching versus Local Switching for Trajectory Surface Hopping Molecular Dynamics Simulation: Cis ↔ Trans Azobenzene Photoisomerization[END_REF] to calculate the hopping probability at the crossings between electronic states. Standard FSSH formulation suffers from overcoherence problem [START_REF] Schwartz | Quantum decoherence and the isotope effect in condensed phase nonadiabatic molecular dynamics simulations[END_REF][START_REF] Subotnik | Understanding the Surface Hopping View of Electronic Transitions and Decoherence[END_REF], which implies that the off-diagonal density matrix elements do not necessarily decay during the dynamics (see the sidebar Spin Decoherence in Molecular Magnets and Qubits). Several, decoherence correction schemes were proposed [START_REF] Schwartz | Quantum decoherence and the isotope effect in condensed phase nonadiabatic molecular dynamics simulations[END_REF][START_REF] Jasper | Electronic decoherence time for non-Born-Oppenheimer trajectories[END_REF][START_REF] Granucci | Including quantum decoherence in surface hopping[END_REF][START_REF] Subotnik | A new approach to decoherence and momentum rescaling in the surface hopping algorithm[END_REF][START_REF] Shu | Implementation of Coherent Switching with Decay of Mixing into the SHARC Program[END_REF] to correct for this behavior.

Although the FSSH algorithm is often criticized as an ad hoc theory, there have been attempts to show that the decoherence-corrected surface hopping can be obtained from the mixed quantumclassical Liouville equation [START_REF] Subotnik | Can we derive Tully's surface-hopping algorithm from the semiclassical quantum Liouville equation? Almost, but only with decoherence[END_REF][START_REF] Wang | Fewest Switches Surface Hopping in Liouville Space[END_REF][START_REF] Kapral | Surface hopping from the perspective of quantum-classical Liouville dynamics[END_REF][START_REF] Martens | Surface Hopping by Consensus[END_REF].

Multiple spawning

The multiple spawning method belongs to a family of techniques where the nuclear wave function is represented as a linear combination of the Gaussian basis functions, which are local in nuclear configuration space [START_REF] Heller | Frozen Gaussians: A very simple semiclassical approximation[END_REF]. In the full multiple spawning (FMS) method [START_REF] Martinez | Multi-Electronic-State Molecular Dynamics: A Wave Function Approach with Applications[END_REF][START_REF] Ben-Nun | Ab Initio Quantum Molecular Dynamics[END_REF], these trajectory basis functions (TBFs) propagate in time on multiple electronic states. If TBFs are properly distributed in configuration space, FMS approaches the exact solution of TDSE.

However, for systems with more than a few DOF, propagating a large number of TBFs is too computationally expensive. To minimize the basis set size without significant loss in accuracy, FMS uses the adaptive basis set that is expanded by spawning new TBFs in the regions of strong coupling between electronic states (Figure 2). In principle, FMS can be employed as a gridbased method requiring the global knowledge of PESs and couplings. However, the strength of localized TBFs is fully utilized in the direct dynamics ab initio multiple spawning (AIMS) method [START_REF] Ben-Nun | Ab Initio Quantum Molecular Dynamics[END_REF][START_REF] Levine | Implementation of ab initio multiple spawning in the Molpro quantum chemistry package[END_REF][START_REF] Yang | Ab Initio Multiple Spawning: First Principles Dynamics Around Conical Intersections[END_REF], which uses the saddle-point approximation (SPA) and the independent first generation approximation (IFGA) [START_REF] Mignolet | A walk through the approximations of ab initio multiple spawning[END_REF].

The FMS/AIMS formalism is based on the expansion of the total wave function in the basis of the electronic wave functions :

The nuclear wave functions are superpositions of the multidimensional Gaussian TBFs , where are complex amplitudes and N TBF is a total number of TBFs. Each TBF follows a classical trajectory with the position and momentum centers and propagating on the electronic state I according to the classical Hamilton EOM, while the phase propagates semiclassically. The time-independent Gaussian widths collected in the vector are parametrized for each type of nucleus. The following EOM for amplitudes is obtained by inserting Equations 7 and 8 into the TDSE:

Equation 9 is written in matrix form, with the elements of the overlap and Hamiltonian matrices defined as

The overlap matrix elements and their time-derivatives arise due to the nonorthogonality of TBFs propagating on the same electronic state. The diagonal Hamiltonian matrix elements are electronic state energies, and the off-diagonal elements are interstate couplings. It is important to point out that, in contrast to the independent trajectories in the TSH method, the TBFs in FMS/AIMS are coupled, which prevents overcoherence [START_REF] Ibele | Excited-state dynamics of molecules with classically driven trajectories and Gaussians[END_REF]. The original AIMS method has been developed to study ICs at conical intersections with NAC arising from the kinetic energy matrix elements in Equation 12. To model spin-crossovers between spin-diabatic states, AIMS was generalized by adding SOC to the Hamiltonian (61-63).

The matrix elements in Equation 12 are obtained by integration over all electronic, nuclear and spin DOF, which requires the knowledge of the entire PESs and coupling surfaces.

The SPA makes the direct AIMS dynamics possible by using the Taylor expansion to approximate these matrix elements. In the zeroth-order SPA, the matrix elements between two TBFs centered at and are calculated at the position of their centroid , which removes the problem of nonlocality [START_REF] Martinez | Multi-Electronic-State Molecular Dynamics: A Wave Function Approach with Applications[END_REF][START_REF] Martínez | Non-adiabatic molecular dynamics: Split-operator multiple spawning with applications to photodissociation[END_REF]. For example, the SOC matrix elements are approximated as

It is important to point out that for SOC, the SPA works very well because the magnitude of coupling does not change drastically around the state crossing.

The key feature of AIMS is an adaptive basis set that is expanded by spawning new TBFs in the regions of strong coupling between electronic states. For two TBFs propagating on the states I and J with different spin, the spawning is triggered if a preset threshold becomes smaller than the effective coupling parameter

The denominator, defined as the energy gap between two states, is introduced to avoid an excessive spawning far from the state crossing regions. A new TBF is created with zero amplitude, and therefore zero population. This is an important distinction from the stochastic methods, such as TSH, where a mere fact of a hop means that the population has been transferred between electronic states. In AIMS, the population is transferred through the solution of the TDSE accounting for quantum effects in a natural way. Although a choice of the coupling parameter ensures that spawning occurs only where significant population can be transferred, the event of spawning itself does not automatically translate into population transfer. Therefore, spawning is just an efficient and physically motivated way of reducing the basis set size. During a spawning event, classical energy conservation can be maintained by rescaling both the momentum and the position of the new TBF. The finite size of TBFs allows for significant adjustments of their momenta and positions, potentially reducing the number of frustrated spawning events compared to the number of frustrated hops in the TSH method. A detailed description of the spawning algorithm can be found in [START_REF] Ben-Nun | Ab Initio Quantum Molecular Dynamics[END_REF][START_REF] Levine | Implementation of ab initio multiple spawning in the Molpro quantum chemistry package[END_REF][START_REF] Yang | An "optimal" spawning algorithm for adaptive basis set expansion in nonadiabatic dynamics[END_REF].

The IFGA reduces the cost of AIMS simulations by treating the initial TBFs independently. Every initial TBF stays coupled to all TBFs that it spawned. This is contrasted with FMS, in which every TBF is coupled to every other TBF, even to TBFs spawned by other initial basis functions, resulting in the number of matrix elements to be calculated growing quadratically. While rigorous, the FMS approach is inefficient because initial TBFs often quickly move away from each other, and their overlaps become negligible. Similar to TSH, the results of AIMS dynamics are averaged over multiple simulations, each starting from a single initial TBF, with positions and momenta sampled from some distribution [START_REF] Ben-Nun | Ab Initio Quantum Molecular Dynamics[END_REF][START_REF] Hack | Comparison of full multiple spawning, trajectory surface hopping, and converged quantum mechanics for electronically nonadiabatic dynamics[END_REF]. It is important to point out that even if IFGA is applied and multiple AIMS simulations are run independently, the number of spawned basis functions in individual AIMS simulations can become too large, drastically increasing simulation time or even rendering it intractable. In such cases, a careful tuning of the spawning threshold parameter is required. In addition, this problem can be alleviated by using the recently introduced stochastic-selection ab initio multiple spawning (SSAIMS) method, which essentially separates groups of TBFs into independent simulations if these groups stay uncoupled during dynamics [START_REF] Curchod | SSAIMS-Stochastic-Selection Ab Initio Multiple Spawning for Efficient Nonadiabatic Molecular Dynamics[END_REF]. Such situation is very common in high-dimensional problems.

Multiconfiguration time-dependent Hartree

The TDSE can be solved by expanding nuclear wave function into a time-independent orthogonal basis set, such as the eigenfunctions of a harmonic oscillator, with time-dependent coefficients. While such grid-based approach rigorously describes the motion of a nuclear wave packet, it scales exponentially with the number of nuclear DOF, and therefore, applicable to only small systems. In MCTDH [START_REF] Meyer | The multi-configurational time-dependent Hartree approach[END_REF][START_REF] Manthe | Wave-packet dynamics within the multiconfiguration Hartree framework: General aspects and application to NOCl[END_REF][START_REF] Beck | The multiconfiguration time-dependent Hartree (MCTDH) method: a highly efficient algorithm for propagating wavepackets[END_REF], the computational cost is lowered significantly by nonlinear differential equations that can be solved by a predictor-corrector integration scheme.

However, a suitable integration scheme, called constant mean-field (CMF) integration has been designed to solve these EOM accurately and efficiently [START_REF] Beck | The multiconfiguration time-dependent Hartree (MCTDH) method: a highly efficient algorithm for propagating wavepackets[END_REF][START_REF] Meyer | Studying molecular quantum dynamics with the multiconfiguration time-dependent Hartree method[END_REF].

To describe nonadiabatic transitions, an extra electronic DOF with the number of SPFs equal to the number of electronic states is added. If only one set of SPFs is used to describe all electronic states, this is called a single-set formulation. If the PESs are very different from each other, a multi-set formulation with different SPF sets for each electronic state is preferable [START_REF] Meyer | Studying molecular quantum dynamics with the multiconfiguration time-dependent Hartree method[END_REF].

A vibrational mode combination technique is introduced to further reduce the computation cost.

The idea is to combine the physical coordinates into particles or logical coordinates, and thereby, shorten the size of the expansion coefficient vector. A good practice is to combine strongly correlated modes in one particle and keep all the particle grids similar in size. Alternatively, the DOF with similar vibrational frequencies can be combined. It is advisable not to combine too many modes and not to construct particles with large size. The idea of mode combination is the foundation for multilayer (ML) MCTDH [START_REF] Vendrell | Multilayer multiconfiguration time-dependent Hartree method: Implementation and applications to a Henon-Heiles Hamiltonian and to pyrazine[END_REF], which was originally developed by Wang and Thoss [START_REF] Wang | Multilayer formulation of the multiconfiguration timedependent Hartree theory[END_REF] and later reformulated by Manthe [START_REF] Manthe | A multilayer multiconfigurational time-dependent Hartree approach for quantum dynamics on general potential energy surfaces[END_REF] for an arbitrary number of layers. ML-MCTDH is very computationally efficient and can handle more than 1000 nuclear DOF ( 75), but strict convergence can be difficult to achieve.

Another important MCTDH development is the variational multiconfigurational Gaussian wave packet (vMCG) method [START_REF] Worth | A novel algorithm for non-adiabatic direct dynamics using variational Gaussian wavepackets[END_REF][START_REF] Richings | Quantum dynamics simulations using Gaussian wavepackets: the vMCG method[END_REF][START_REF] Lasorne | Direct quantum dynamics using variational multi-configuration Gaussian wavepackets. Implementation details and test case[END_REF]. However, to the best of our knowledge, it has not yet been applied to spin-crossovers. The idea was conceived from the G-MCTDH method [START_REF] Burghardt | Approaches to the approximate treatment of complex molecular systems by the multiconfiguration time-dependent Hartree method[END_REF][START_REF] Worth | Using the MCTDH wavepacket propagation method to describe multimode non-adiabatic dynamics[END_REF],

where some of SPFs are replaced by Gaussian functions. Replacing all SPFs by multidimensional frozen Gaussians, as done in vMCG, leads to the same wave function as in the FMS method. However, the Dirac-Frenkel variational principle used in vMCG yields the EOM for the expansion coefficients and the Gaussian parameters that are different from the ones in FMS. Similar to FMS/AIMS, the direct dynamics version of vMCG has been also implemented [START_REF] Lasorne | Direct quantum dynamics using variational multi-configuration Gaussian wavepackets. Implementation details and test case[END_REF].

Because the grid-based MCTDH works with non-local PESs and couplings, it is necessary to remove the singular NAC by transforming the adiabatic or spin-diabatic electronic states to the diabatic representation [START_REF] Mukherjee | Beyond Born-Oppenheimer theory for spectroscopic and scattering processes[END_REF][START_REF] Naskar | ADT: A Generalized Algorithm and Program for Beyond Born-Oppenheimer Equations of[END_REF][START_REF] Baer | Beyond Born-Oppenheimer: Conical Intersections and Electronic Nonadiabatic Coupling Terms[END_REF]. Vibronic coupling models are often used to construct the diabatic Hamiltonians [START_REF] Penfold | Spin-Vibronic Mechanism for Intersystem Crossing[END_REF][START_REF] Köppel | Multimode Molecular Dynamics Beyond the Born-Oppenheimer Approximation[END_REF]. For the model with electronic states, the -dimensional Hamiltonian matrix is expanded around some, usually the Frank-Condon (FC), geometry :

where is the ground state potential, and is the diagonal matrix of the vertical excitation energies. Truncating the series at and produces the linear vibronic coupling (LVC)

and quadratic vibronic coupling (QVC) models, respectively. The LVC parameters, the diagonal forces and the off-diagonal interstate couplings in , for each normal mode are where is the adiabatic energy of the electronic state . The spin-vibronic Hamiltonian can be obtained as a sum of the spin-free vibronic Hamiltonian and the spin-orbit Hamiltonian (85-89). However, it is possible to construct the spin-vibronic Hamiltonian directly from the fully adiabatic spin-mixed states (90-94).

Electronic structure methods

Here, we discuss some practical aspects of employing different electronic structure methods to model spin-crossover dynamics. For comprehensive review of electronic structure methods used in NAMD the reader is referred to [START_REF] Crespo-Otero | Recent Advances and Perspectives on Nonadiabatic Mixed Quantum-Classical Dynamics[END_REF][START_REF] Yang | Ab Initio Multiple Spawning: First Principles Dynamics Around Conical Intersections[END_REF][START_REF] Ibele | Excited-state dynamics of molecules with classically driven trajectories and Gaussians[END_REF][START_REF] Lischka | Multireference Approaches for Excited States of Molecules[END_REF]. The success of NAMD simulations depend on the accuracy of the electronic structure methods used to obtain the energies, energy gradients, and couplings for multiple electronic states at different molecular geometries, as has been demonstrated in the recent TSH study of the spin-crossover dynamics in thioformaldehyde [START_REF] Mai | The Influence of the Electronic Structure Method on Intersystem Crossing Dynamics. The Case of Thioformaldehyde[END_REF]. On the other hand, because the number of electronic structure calculations limits the timescale of direct NAMD simulations and the number of nuclear DOF in grid-based dynamics, the electronic structure methods must be computationally efficient. For example, a 2 ps direct TSH simulation with the average 0.2 fs time step (often a variable time step is used) and a small sampling over 100 trajectories will require 10 6 electronic structure calculations. Assuming the trajectories can be run in parallel, to complete them within one week will require calculating energies, gradients and couplings for multiple electronic states approximately every minute. A similar AIMS simulation, which have to calculate matrix elements between TBFs, will require even more electronic structure calculations. For grid-based NAMD, the same 10 6 electronic structure calculations will be required to build the global PESs and coupling surfaces for a system with 3 nuclear DOF, assuming the calculations generate a dense 3-dimensional grid with 100 points per each DOF. Modern surface interpolation techniques aimed to overcome the "curse of dimensionality" by using the sparse grids and focusing on the relevant parts of PESs can push the number of DOF to 24-39 (8-13 atoms), for roughly the same number of electronic structure calculations [START_REF] Conte | Efficient Generation of Permutationally Invariant Potential Energy Surfaces for Large Molecules[END_REF][START_REF] Zhu | Fitting coupled potential energy surfaces for large systems: Method and construction of a 3-state representation for phenol photodissociation in the full 33 internal degrees of freedom using multireference configuration interaction determined data[END_REF][START_REF] Shen | Construction of Quasi-diabatic Hamiltonians That Accurately Represent ab Initio Determined Adiabatic Electronic States Coupled by Conical Intersections for Systems on the Order of 15 Atoms. Application to Cyclopentoxide Photoelectron Detachment in the Full 39 Degrees of Freedom[END_REF]. Constructing the vibronic Hamiltonian models commonly used with the MCTDH method is less computationally expensive. Assuming the same 100 electronic structure calculations per vibrational mode, the 15mode LVC and QVC models will require roughly 15 100=1500 and 15 2 100/2=11250 calculations, respectively. However, these models typically describe only small parts of PESs in the vicinity of the FC geometry. The electronic structure methods used in the direct NAMD must be also very reliable. A failure to produce smooth electronic energies, gradients and couplings during time propagation can lead to the failure of the entire simulation.

Because of their computational efficiency, the density functional theory (DFT) and multiconfigurational self-consistent field (MCSCF) based methods are widely used in direct NAMD. The conventional linear-response time-dependent DFT (LR-TDDFT) method is a popular choice. However, the LR-TDDFT reliance on a single-reference closed-shell wave function for the ground state and singly excited configurations for excited states limits its applicability and accuracy. The situations with strongly multi-configurational electronic states, such as diradical states in organic molecules, electronic states of transition metal complexes, and crossing seams between the ground and excited electronic states cannot be described accurately.

The excited states dominated by double excitations also cannot be modeled by conventional LR-TDDFT. The spin-flip TDDFT capable of describing the multi-configurational and doubly excited electronic states (100) could become an important method for modeling spin-crossover dynamics.

The MCSCF methods, including the complete active space self-consistent field (CASSCF) version, are well-suited for the direct NAMD simulations of the processes where multiple electronic configurations play an important role, assuming relatively small active spaces can describe these processes. Efficient implementations of analytical energy gradient, NAC and SOC for state-averaged (SA)-CASSCF are widely available. However, neglecting the dynamic electron correlation in CASSCF can lead to incorrect energy gaps between electronic states of different character. While covalent states often have a significant contribution from multiple electron configurations, and therefore, are stabilized by the static correlation recovered by CASSCF, for ionic states dominated by a single configuration, little correlation energy is recovered. This leads to the overstabilization of covalent states with respect to ionic states, which in turn can affect the location of state crossings and overall outcome of NAMD simulations.

Another important issue is related to the selection of the active space (101), which, ideally, should describe all possible outcomes of an NAMD simulation. However, such active space is usually too large to be practical. Therefore, a smaller active space must be selected to describe the electronic transitions of interest and the parts of a molecule where chemical bonds are expected to break and/or form. This approach could lead to NAMD oversampling the reaction pathways best described by the chosen active space at the expense of other pathways. In addition, if dynamics leads to molecular geometries that are not correctly described by the selected active space, the active space orbitals could rotate out of the active space. The results of such sudden active space change include discontinuities in electronic energies, gradients and couplings, as well as failed CASSCF convergence often leading to the failure of the entire NAMD simulation.

While using the post-MCSCF methods, such as multi-reference perturbation theory and configuration interaction, to account for dynamic electron correlation in direct NAMD is possible, it is often impractical due to the high computational cost. This led to the development of parametrized approaches, such as scaled-CASSCF (102), -CASSCF (103), and floating occupation molecular orbitals CASCI (104-106). Despite the fact that these methods were developed to model IC dynamics, they are also applicable to spin-crossovers. In the future, other emerging multireference methods with analytical energy gradients, NAC and SOC (95) are expected to find applications in modeling spin-crossover dynamics.

For constructing the vibronic coupling Hamiltonian models commonly used in the gridbased dynamics, the requirements for electronic structure methods are less stringent than for direct dynamics. A much smaller number of electronic structure calculations is required, and the failure of any particular calculation does not lead to the failure of the NAMD simulation. Also, implementations of analytical energy gradient and couplings, while helpful, are not required. As a result, the post-MCSCF methods, such as MS-CASPT2 and MR-CISD, can be used to obtain the vibronic Hamiltonian parameters. However, these reduced requirements come at the price of restricting dynamics to the relatively small regions of PESs where vibronic coupling models are valid.

The most common approach to calculate SOC between spin-diabatic states is the firstorder perturbation theory in combination with the Breit-Pauli or Douglas-Kroll-Hess spin-orbit

Hamiltonians [START_REF] Fedorov | Spin-orbit coupling in molecules: Chemistry beyond the adiabatic approximation[END_REF]107,108). Such calculations produce the SOC matrix elements between the individual components of the spin states. The root-mean-square of these matrix elements for the two states with spin S and S' is defined as the SOC constant representing the effective coupling between two spin states [START_REF] Fedorov | Spin-orbit coupling in molecules: Chemistry beyond the adiabatic approximation[END_REF]. The SOC matrix elements have been implemented for the multireference [START_REF] Fedorov | Spin-orbit coupling in molecules: Chemistry beyond the adiabatic approximation[END_REF]109,110), DFT (111-114) and single-reference (115) methods. For the systems with heavy elements where S is not a good quantum number, including complexes of third-row transition metals, lanthanides and actinides, it is desirable to carry out NAMD simulations in the spin-adiabatic basis with the NAC between the spin-mixed states driving the interstate population transfer. While such spin-adiabatic simulations are starting to emerge [START_REF] Martínez-Fernández | Competing ultrafast intersystem crossing and internal conversion: a time resolved picture for the deactivation of 6-thioguanine[END_REF](116)(117)(118), the progress requires interfacing NAMD with the electronic structure methods capable of calculating the analytical NAC between the spin-mixed states (119, 120).

APPLICATIONS

In this section, we describe several selected applications of NAMD to model spincrossovers in systems of different complexity. First, we focus on the SO 2 molecule, which due to its small size and complex photophysics has been studied with the grid-based methods and several flavors of TSH dynamics. Second, we describe two test applications of the AIMS method to spin-crossovers in the GeH 2 and H 2 CS molecules. Finally, we highlight the capabilities of NAMD to model spin-crossovers in large systems by discussing the full-dimensional direct TSH simulations of the [Ru(bpy) 3 ] 2+ complex and the MCTDH modeling of the active center of myoglobin protein using a reduced dimensionality LVC model.

SO 2 molecule

The SO 2 molecule has been studied with different NAMD approaches, including fulldimensional grid-based quantum dynamics (85, 121) and direct TSH dynamics (122, 123). Xie et al. constructed three-dimensional PESs for the lowest singlet ( ) and triplet ( ) states using the MRCI+Q/aug-cc-pVTZ level of theory with the full-valence active space (121). Only two singlet states and the lowest triplet state were included in the NAMD performed in the quasi-diabatic representation with constant SOC. The Chebyshev propagation to more than 1200 fs was carried out in the Jacobi coordinates using the discrete variable representation (DVR) for the radial DOF and the finite basis representation (FBR) for the angular DOF. This study showed a reasonable agreement between the calculated and experimental absorption spectra and demonstrated the effect of the lowest triplet state on the excited state dynamics. The calculated electronic state population showed that the spin-crossover between and states occurs within 100 fs. Köppel and coworkers also investigated the SO 2 photoexcitation dynamics with triplet states using the DVR based quantum dynamics [START_REF] Lévêque | Communication: Theoretical prediction of the importance of the 3 B 2 state in the dynamics of sulfur dioxide[END_REF]. They obtained the PESs of the lowest six states, which are symmetry-isolated from the manifold, with MRCI+Q/cc-PVTZ and the full-valence active space. The triplet 1 components were combined into the symmetric (+) and antisymmetric (-) states (Figure 3). Threedimensional SOC surfaces were calculated using the Breit-Pauli Hamiltonian. After diabatization of the six-state model Hamiltonian, the NAMD simulations were initiated from the vertically excited 1 B 1 state. After 1 ps, the populations of the 1 B 1 and 1 A 2 states were 24% and 41% states were included in the 50-trajectory simulation. Again, a fast S 2 S 1 IC, followed by the S 1 T 2 and S 1 T 3 spin-crossovers were observed. These were followed by population transfer to T 1 via successive ICs. The liquid phase QM/MM simulations predicted the effective spincrossover rate that is about twice faster than that in the gas phase. Recently, Gonzalez and coworkers parameterized an LVC model for SO 2 using only a single excited state electronic structure calculation and a ground state vibrational frequency calculation. The model was used to carry out a 200-trajectory TSH simulation on four singlet and three triplet diabatic states (124).

This work demonstrated that, at least for this small system, the main timescales predicted by direct dynamics simulations can be reproduced with the LVC model dynamics while drastically reducing computational expenses.

Applications of AIMS to GeH 2 and H 2 CS

The AIMS method has been tested on the spin-crossover between the lowest excited 3 B 1 and ground 1 A 1 states of GeH 2 [START_REF] Fedorov | Ab Initio Multiple Spawning Method for Intersystem Crossing Dynamics: Spin-Forbidden Transitions between 3 B 1 and 1A1 States of GeH 2[END_REF][START_REF] Fedorov | Predicting Intersystem Crossing Rates with AIMS-DFT Molecular Dynamics[END_REF]. The spin-diabatic states were obtained using the full valence active space CASSCF and unrestricted DFT (B3LYP) methods with 6-31G* basis set, as implemented in the GAMESS suite of programs (125). The TBFs were propagated using the state-specific energies and gradients, while SOC was calculated perturbatively using the SA(2)-CASSCF orbitals and the Breit-Pauli Hamiltonian (126). In the DFT-based dynamics, SOC was also evaluated using the high-spin unrestricted DFT orbitals. The 48 initial TBFs, which started form the FC region with the initial conditions sampled from the Wigner distribution, generated about 1200 TBFs during the 150 fs simulations. Due to the presence of a heavy atom in GeH 2 , the SOC between the two states is relatively large (350 cm -1 ) leading to a fast decay of the excited state population. After 150 fs, around 60% of population was transferred in both CASSCF and DFT simulations. Surprisingly, the CASSCF and DFT-based dynamics predicted essentially the same lifetime for the 3 B 1 state (186 fs and 182 fs). However, this agreement is believed to be due to a fortunate cancelation of two effects: the higher spin-crossover energy barrier and stronger SOC predicted by DFT compared to the corresponding CASSCF values.

Another implementation of the generalized AIMS method has been used to study the spin-crossover dynamics in thioformaldehyde [START_REF] Curchod | Communication: GAIMS-Generalized Ab Initio Multiple Spawning for both internal conversion and intersystem crossing processes[END_REF]. The simulations were performed at the SA(4)-CASSCF(4,3)/6-31G* level of theory and included four electronic states (S 0 , S 1 , T 1 , T 2 ).

The 20 initial TBFs sampled from the Wigner distribution were started from the FC region on the S 1 state. During the 200 fs dynamics, 326 TBFs were spawned with 306 of them on the triplet states (Figure 4). The population transferred to T 1 was negligible, whereas T 2 had a significant population of 8% after 200 fs. These results can be explained by El-Sayed's rule, which predicts a stronger SOC between the states of different character. Therefore, the S 1 state with nπ character is coupled stronger to the ππ T 2 state than to nπ T 1 . As the C=S bond time evolution indicates, the dynamics of TBFs on the nπ states S 1 and T 1 are similar, whereas the T 2 TBFs have a longer average C=S bond consistent with the ππ character of the T 2 state. While the 200 fs simulations were too short to calculate the accurate spin-crossover rates, they provided valuable qualitative insight into the H 2 CS nonadiabatic dynamics.

Limitations of the direct dynamics when applied to spin-crossovers arise from the generally longer timescale of these processes compared to ICs. In the AIMS method, longer trajectories required for the spin-crossover simulations are challenging because of the growing size of the adaptive basis set. The need to calculate the matrix elements between each TBF pair results in a superlinear growth of computational cost. Techniques to reduce the number of TBFs, such as removal of TBFs with small contributions to the total wave function and careful selection of spawning criteria can significantly accelerate AIMS simulations. For high-dimensional systems, SSAIMS method (67) can drastically reduce the computational cost by uncoupling groups of TBFs in the simulation without sacrificing the accuracy. Another way to overcome the quickly growing basis set problem is to accelerate the electronic structure calculations through parallelization and the use of graphical processing units (127).

Large systems

Due to its simplicity and computational efficiency, the direct TSH method is the most widely used type of NAMD to model spin-crossovers in various systems [START_REF] Cui | Generalized trajectory surface-hopping method for internal conversion and intersystem crossing[END_REF][START_REF] Martínez-Fernández | Competing ultrafast intersystem crossing and internal conversion: a time resolved picture for the deactivation of 6-thioguanine[END_REF][START_REF] Favero | Dynamics of acetone photodissociation: a surface hopping study[END_REF](128)(129)(130)(131)(132)(133)(134).

The grid-based TSH dynamics on the precomputed PESs is also used to study spin-crossovers in small molecules (135-140). A recent study of the [Ru(bpy) 3 ] 2+ complex is an excellent example demonstrating the capability of direct TSH dynamics to model spin-crossovers in large systems (132) . This complex serves as a prototype for a class of compounds where ultrafast spincrossovers occur between the singlet and triplet metal-to-ligand-charge transfer (MLCT) states.

The electronic energies, gradients, NAC and SOC were calculated with LR-TDDFT using the PBE functional. Since the density of excited states is very high, 101 trajectories were initialized from different singlet excited states and propagated for 30 fs on the manifold of 15 singlet and 15 triplet states (Figure 5). Almost 65% of the total population was transferred to the triplet states corresponding to the time constant of 26 3 fs, which is in good agreement with the experimental value of 15 10 fs. Two interesting observations came out of this study. First, in contrast to Kasha's rule, the spin-crossovers occur between the high-lying singlet and triplet states around the FC geometry. Second, in addition to the high density of states and strong SOC (up to 350 cm - 1 ), the vibrational motions of the N and Ru atoms play a major role in promoting the spincrossovers.

There

CONCLUSIONS AND FUTURE DIRECTIONS

In the last decade, spin-crossover dynamics evolved from a somewhat exotic niche within the field of nonadiabatic dynamics primarily concerned with ICs at conical intersections to a very active area of research. This rapid development has been driven by the realization that, despite being mediated by relativistic effects, spin-crossovers can occur on a time scale similar to that of ICs, even in molecules without heavy atoms where the relativistic effects have been traditionally

neglected. Expanding the capability of NAMD to model spin-crossovers lead to a strong interest in nonadiabatic dynamics in the transitional metal and even lanthanide complexes, where spindependent processes play a central role. These new methodological developments also raised multiple fundamental questions: How strong the SOC can be before we need to completely abandon the spin-diabatic representation, and therefore, the concept of spin-crossover? When the spin-dependent relativistic effects can be accounted for using the perturbation theory and when they must be treated variationally at the molecular orbital level (interestingly, the perturbation approach seems to work reasonably well for lanthanide complexes)? Can the higher-order effects, such as the spin-spin and hyperfine couplings, drive spin-crossovers, especially if SOC is small as between states with >1? Can an external magnetic field be used to control the spincrossover rates in photochemistry, thermally activated spin-forbidden reactions, and molecular magnets? The future development of the NAMD and electronic structure methods, some of which we attempted to describe in this review, should help to answer these questions.

SUMMARY POINTS

1. Spin-crossovers, transitions between electronic states with different spin multiplicity, are common and can be initiated by different stimuli, including light, temperature, pressure, and magnetic field.

2. While spin-crossovers mediated by the relatively weak spin-orbit and higher-order couplings often occur on a slower time scale than ICs, these two types of nonadiabatic processes can compete with each other.

3. The NAMD methods originally developed to simulate ICs have been extended to model spin-crossover dynamics, with the goals of predicting the spin-crossover rates and gaining insights into the reaction mechanisms involving electronic states with different spin.

4. Because building the full-dimensional PESs and interstate coupling surfaces is not practical for systems with more than a few nuclear DOF, the direct and reduceddimensionality NAMD methods are the most common approaches to model spincrossover dynamics.

5. The NAMD simulations require accurate and computationally efficient electronic structure methods. These methods must produce accurate energies, energy gradients and couplings for multiple electronic states at different molecular geometries.

FUTURE ISSUES

1. At what point relativistic effects become too strong to justify working in the spin-diabatic basis?

2. How the slow spin-crossovers can be modeled using NAMD? Perhaps some combination of the short NAMD simulations and the nonadiabatic statistical theories can be developed.

3. Can the higher-order effects, such as the spin-spin and hyperfine couplings, influence the spin-crossover rates? This is especially important for the situations when the SOC is small.

4. Can the effect of an external magnetic field on the spin-crossover rates be modeled with NAMD? This should help to understand how the spin-crossover dynamics can be controlled with a magnetic field.

5. Can computationally efficient multireference electronic structure methods be developed for direct NAMD? Such methods must account for dynamic electron correlation and have analytical energy gradient, NAC and SOC.

  representing a nuclear wave function as a linear combination of f-dimensional products of the single particle functions (SPFs) : where are the nuclear coordinates, f is the number of nuclear DOF, is the number of SPFs for the k th DOF, and are the time-dependent expansion coefficients. Each SPF is expressed as a linear combination of the time-independent primitive basis functions with the time-dependent coefficients: Therefore, the MCTDH ansatz employs standard wave packet expansion, but with timedependent SPFs. This reduces the number of basis functions required for a converged calculation by providing the variationally determined basis for the optimal description of the evolving wave packet. There are two constraints invoked into the formal MCTDH derivation: (a) initially orthonormal SPFs remain orthonormal at all times; (b) the constraint operator ensuring the uniqueness of the nuclear wave function is Hermitian (70). The wave function defined by Equations 15 and 16 is inserted into the Dirac-Frenkel variational principle to obtain the EOM for the time-dependent expansion coefficients and SPFs: In Equation 17, composite indices and the configuration function are defined as and . In Equation 18, is the vector of SPFs, and is the -dimensional identity matrix. The Hamiltonian is split into the separable terms acting only on a single DOF, and that includes all correlations between DOF. and are the mean-field and density matrices. The projector operator ensures that the SPF time-derivative is orthogonal to the space spanned by SPFs. The eigenfunctions of the density matrix are called natural orbitals, and the eigenvalues correspond to populations of these orbitals. As the space spanned by the natural orbitals is equivalent to that of the original SPFs, the natural orbital populations provide a measure for the quality of the MCTDH wave function. If the population of the highest natural orbital is negligibly low, this orbital (SPF) is redundant, and hence, the MCTDH wave function is of good quality. The time-dependent Hartree method is a limiting case of MCTDH with all On the other hand, the standard wave packet expansion into a primitive basis corresponds to When SPFs do not form the complete basis set, the variational method ensures that the available SPFs provide the best possible basis set to describe the wave function at each timestep. The EOM (Equations 17 and 18) are coupled

6 .

 6 photoexcitation to the bright singlet Q band, the complex experiences the Jahn-Teller symmetry breaking and the population transfer to the singlet 1 MLCT states within 26 fs. Since the 1 MLCT states are strongly spin-orbit coupled to the triplet3 MLCT manifold, a rapid singlet-triplet spincrossover takes place. Finally, the population is transferred from 3 MLCT to the quintet 5 MLCT

Figure 1 .

 1 Figure 1. Electronic states in the diabatic (a), spin-diabatic (b), and fully adiabatic (c)

Figure 3 .

 3 Figure 3. The lowest energy states of SO 2 . One-dimensional cut through the PESs along the

Figure

  Figure reproduced with permission from Reference (85); copyright 2014 AIP Publishing.

Figure 4 .

 4 Figure 4. Generalized AIMS dynamics of thioformaldehyde after photoexcitation to the S 1 state.

ACKNOWLEDGEMENTS S.A.V. thanks Todd Martinez, Basile Curchod, Edward Hohenstein, Benjamin Levine, and Dmitry G. Fedorov for fruitful discussions. This material is based on work supported by the National Science Foundation (CAREER award CHE-1654547). Acknowledgment is made to the Donors of the American Chemical Society Petroleum Research Fund for partial support of this research (award 60481-ND6).

respectively, while 25% of the population was distributed over the two 3 B 2 components, and remaining 10% decayed to 3 B 1 (-) and 3 A 2 (-) states. The overall dynamics was described as a three-step process with the characteristic time intervals 0-16 fs, 55-65 fs and 150-165 fs. In the first step, the 1 B 1 3 A 2 spin-crossover occurs in the FC region around the crossing point labeled as ISC3 in Figure 3, followed by the rapid 3 A 2 3 B 1 IC through the conical intersection CI2.

During this time, most of the singlet population is transferred to 1 A 2 via CI1. In the second step, the population on 1 A 2 propagates into the region of smaller bending angle, and the 1 A 2 3 B 1 spin-crossover takes place at ISC5. In the third step, the wave packet returns to the FC region repeating the first step population transfer. This study confirmed the important role of the triplet states in the deactivation process of SO 2 and demonstrated the complexity of nonadiabatic dynamics, even in small molecules.

The deactivation pathways in SO 2 were also studied by Gonzalez and coworkers using the direct TSH dynamics (122). They have run 111 trajectories on the four lowest singlet and three lowest triplet states calculated at the MR-CIS/ANO-RCC-VDZP level of theory with SOC calculated using the effective Fock-type spin-orbit operator. The simulations showed a rapid 1 B 1 1 A 2 IC, followed by a significant population transfer from 1 A 2 to 3 B 2 with small amount of population going to 3

Sidebars NONADIABATIC STATISTICAL THEORIES

Nonadiabatic statistical theories (NASTs) are based on the assumption that the rate of intramolecular energy distribution is much faster than the spin-crossover rate [START_REF] Lykhin | Nonadiabatic transition state theory: Application to intersystem crossings in the active sites of metalsulfur proteins[END_REF][START_REF] Jasper | Multidimensional Effects in Nonadiabatic Statistical Theories of Spin-Forbidden Kinetics: A Case Study of 3 O + CO → CO 2[END_REF][START_REF] Harvey | Understanding the kinetics of spin-forbidden chemical reactions[END_REF]. These theories can account for quantum effects, such as tunneling and zero-point vibrational energy [START_REF] Lykhin | Intersystem crossing in tunneling regime: T 1 → S 0 relaxation in thiophosgene[END_REF], and are ideally suited to study kinetics of slow spin-crossovers in large complex systems for which long molecular dynamics simulations are not feasible. NASTs require electronic structure information only at very few nuclear geometries, making these theories compatible with the high-level electronic structure methods [START_REF] Lykhin | Intersystem crossing in tunneling regime: T 1 → S 0 relaxation in thiophosgene[END_REF][START_REF] Pokhilko | Spin-Forbidden Channels in Reactions of Unsaturated Hydrocarbons with O( 3 P)[END_REF] and molecular fragmentation techniques [START_REF] Kaliakin | Locating Minimum Energy Crossings of Different Spin States Using the Fragment Molecular Orbital Method[END_REF].

SPIN DECOHERENCE IN MOLECULAR MAGNETS AND QUBITS

Electron spin and orbital angular momentum are responsible for the magnetic properties of the atomic and molecular systems that have been proposed for high-density memory and quantum information applications [START_REF] Coronado | Molecular magnetism: from chemical design to spin control in molecules, materials and devices[END_REF][START_REF] Giménez-Santamarina | Exploiting clock transitions for the chemical design of resilient molecular spin qubits[END_REF]. In molecular magnets, spin relaxation (population transfer between spin states) is responsible for the loss of magnetization, while spin decoherence is related to the loss of quantum information in spin qubits. The ability to predict the spin relaxation (T 1 ) and decoherence (T 2 ) times is critical for the development of practical spinbased molecular magnets and qubits with long lifetimes. Depending on the system and operating temperature, the spin-vibronic and hyperfine interactions can be responsible for electron spin relaxation and decoherence [START_REF] Ullah | In Silico Molecular Engineering of Dysprosocenium-Based Complexes to Decouple Spin Energy Levels from Molecular Vibrations[END_REF][START_REF] Upadhyay | Spin coherence and optical properties of alkali-metal atoms in solid parahydrogen[END_REF][START_REF] Upadhyay | Enhanced spin coherence of rubidium atoms in solid parahydrogen[END_REF].

SPIN EFFECTS IN ATOM-DIATOM COLLISIONS

It is important to mention the full-dimensional nonadiabatic quantum dynamics calculations on collision reactions between an atom and a diatomic molecule that include the SOC effects (142-145). For example, considerable efforts have been devoted to construct the spin-vibronic PESs for the F+H 2 reaction (146-148), where the excited spin-orbit state F * ( 2 P 1/2 ) plays a major role in low energy collisions. The quantum dynamics calculations carried out on these PESs successfully reproduce the experimental collision cross-sections, reaction rate constants and branching ratios for the F+H 2 /HD/D 2 reactions (145, 148-151).